Actual source code: ex9.c
2: static char help[] = "Tests MatCreateComposite()\n\n";
4: /*
5: Include "petscmat.h" so that we can use matrices.
6: automatically includes:
7: petscsys.h - base PETSc routines petscvec.h - vectors
8: petscmat.h - matrices
9: petscis.h - index sets petscviewer.h - viewers
10: */
11: #include <petscmat.h>
13: int main(int argc, char **args)
14: {
15: Mat *A, B; /* matrix */
16: Vec x, y, v, v2, z, z2;
17: PetscReal rnorm;
18: PetscInt n = 20; /* size of the matrix */
19: PetscInt nmat = 3; /* number of matrices */
20: PetscInt i;
21: PetscRandom rctx;
22: MatCompositeType type;
23: PetscScalar scalings[5] = {2, 3, 4, 5, 6};
26: PetscInitialize(&argc, &args, (char *)0, help);
27: PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
28: PetscOptionsGetInt(NULL, NULL, "-nmat", &nmat, NULL);
30: /*
31: Create random matrices
32: */
33: PetscMalloc1(nmat + 3, &A);
34: PetscRandomCreate(PETSC_COMM_WORLD, &rctx);
35: MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, n, n / 2, 3, NULL, 3, NULL, &A[0]);
36: for (i = 1; i < nmat + 1; i++) MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, n, n, 3, NULL, 3, NULL, &A[i]);
37: MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, n / 2, n, 3, NULL, 3, NULL, &A[nmat + 1]);
38: for (i = 0; i < nmat + 2; i++) MatSetRandom(A[i], rctx);
40: MatCreateVecs(A[1], &x, &y);
41: VecDuplicate(y, &z);
42: VecDuplicate(z, &z2);
43: MatCreateVecs(A[0], &v, NULL);
44: VecDuplicate(v, &v2);
46: /* Test MatMult of an ADDITIVE MatComposite B made up of A[1],A[2],A[3] with separate scalings */
48: /* Do MatMult with A[1],A[2],A[3] by hand and store the result in z */
49: VecSet(x, 1.0);
50: MatMult(A[1], x, z);
51: VecScale(z, scalings[1]);
52: for (i = 2; i < nmat + 1; i++) {
53: MatMult(A[i], x, z2);
54: VecAXPY(z, scalings[i], z2);
55: }
57: /* Do MatMult using MatComposite and store the result in y */
58: VecSet(y, 0.0);
59: MatCreateComposite(PETSC_COMM_WORLD, nmat, A + 1, &B);
60: MatSetFromOptions(B);
61: MatCompositeSetScalings(B, &scalings[1]);
62: MatMultAdd(B, x, y, y);
64: /* Diff y and z */
65: VecAXPY(y, -1.0, z);
66: VecNorm(y, NORM_2, &rnorm);
67: if (rnorm > 10000.0 * PETSC_MACHINE_EPSILON) PetscPrintf(PETSC_COMM_WORLD, "Error with composite add %g\n", (double)rnorm);
69: /* Test MatCompositeMerge on ADDITIVE MatComposite */
70: MatCompositeSetMatStructure(B, DIFFERENT_NONZERO_PATTERN); /* default */
71: MatCompositeMerge(B);
72: MatMult(B, x, y);
73: MatDestroy(&B);
74: VecAXPY(y, -1.0, z);
75: VecNorm(y, NORM_2, &rnorm);
76: if (rnorm > 10000.0 * PETSC_MACHINE_EPSILON) PetscPrintf(PETSC_COMM_WORLD, "Error with composite add after merge %g\n", (double)rnorm);
78: /*
79: Test n x n/2 multiplicative composite B made up of A[0],A[1],A[2] with separate scalings
80: */
82: /* Do MatMult with A[0],A[1],A[2] by hand and store the result in z */
83: VecSet(v, 1.0);
84: MatMult(A[0], v, z);
85: VecScale(z, scalings[0]);
86: for (i = 1; i < nmat; i++) {
87: MatMult(A[i], z, y);
88: VecScale(y, scalings[i]);
89: VecCopy(y, z);
90: }
92: /* Do MatMult using MatComposite and store the result in y */
93: MatCreateComposite(PETSC_COMM_WORLD, nmat, A, &B);
94: MatCompositeSetType(B, MAT_COMPOSITE_MULTIPLICATIVE);
95: MatCompositeSetMergeType(B, MAT_COMPOSITE_MERGE_LEFT);
96: MatSetFromOptions(B);
97: MatCompositeSetScalings(B, &scalings[0]);
98: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
99: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY); /* do MatCompositeMerge() if -mat_composite_merge 1 */
100: MatMult(B, v, y);
101: MatDestroy(&B);
103: /* Diff y and z */
104: VecAXPY(y, -1.0, z);
105: VecNorm(y, NORM_2, &rnorm);
106: if (rnorm > 10000.0 * PETSC_MACHINE_EPSILON) PetscPrintf(PETSC_COMM_WORLD, "Error with composite multiplicative %g\n", (double)rnorm);
108: /*
109: Test n/2 x n multiplicative composite B made up of A[2], A[3], A[4] without separate scalings
110: */
111: VecSet(x, 1.0);
112: MatMult(A[2], x, z);
113: for (i = 3; i < nmat + 1; i++) {
114: MatMult(A[i], z, y);
115: VecCopy(y, z);
116: }
117: MatMult(A[nmat + 1], z, v);
119: MatCreateComposite(PETSC_COMM_WORLD, nmat, A + 2, &B);
120: MatCompositeSetType(B, MAT_COMPOSITE_MULTIPLICATIVE);
121: MatSetFromOptions(B);
122: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
123: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY); /* do MatCompositeMerge() if -mat_composite_merge 1 */
124: MatMult(B, x, v2);
125: MatDestroy(&B);
127: VecAXPY(v2, -1.0, v);
128: VecNorm(v2, NORM_2, &rnorm);
129: if (rnorm > 10000.0 * PETSC_MACHINE_EPSILON) PetscPrintf(PETSC_COMM_WORLD, "Error with composite multiplicative %g\n", (double)rnorm);
131: /*
132: Test get functions
133: */
134: MatCreateComposite(PETSC_COMM_WORLD, nmat, A, &B);
135: MatCompositeGetNumberMat(B, &n);
136: if (nmat != n) PetscPrintf(PETSC_COMM_WORLD, "Error with GetNumberMat %" PetscInt_FMT " != %" PetscInt_FMT "\n", nmat, n);
137: MatCompositeGetMat(B, 0, &A[nmat + 2]);
138: if (A[0] != A[nmat + 2]) PetscPrintf(PETSC_COMM_WORLD, "Error with GetMat\n");
139: MatCompositeGetType(B, &type);
140: if (type != MAT_COMPOSITE_ADDITIVE) PetscPrintf(PETSC_COMM_WORLD, "Error with GetType\n");
141: MatDestroy(&B);
143: /*
144: Free work space. All PETSc objects should be destroyed when they
145: are no longer needed.
146: */
147: VecDestroy(&x);
148: VecDestroy(&y);
149: VecDestroy(&v);
150: VecDestroy(&v2);
151: VecDestroy(&z);
152: VecDestroy(&z2);
153: PetscRandomDestroy(&rctx);
154: for (i = 0; i < nmat + 2; i++) MatDestroy(&A[i]);
155: PetscFree(A);
157: PetscFinalize();
158: return 0;
159: }
161: /*TEST
163: test:
164: nsize: 2
165: requires: double
166: args: -mat_composite_merge {{0 1}shared output} -mat_composite_merge_mvctx {{0 1}shared output}
168: TEST*/