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