1 static char help[] = "Tests MatCreateComposite()\n\n";
2
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>
11
main(int argc,char ** args)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 PetscInt i;
20 PetscRandom rctx;
21 MatCompositeType type;
22 PetscScalar scalings[5] = {2, 3, 4, 5, 6};
23
24 PetscFunctionBeginUser;
25 PetscCall(PetscInitialize(&argc, &args, NULL, help));
26 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
27 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nmat", &nmat, NULL));
28
29 /*
30 Create random matrices
31 */
32 PetscCall(PetscMalloc1(nmat + 3, &A));
33 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx));
34 PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, n, n / 2, 3, NULL, 3, NULL, &A[0]));
35 for (i = 1; i < nmat + 1; i++) PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, n, n, 3, NULL, 3, NULL, &A[i]));
36 PetscCall(MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, n / 2, n, 3, NULL, 3, NULL, &A[nmat + 1]));
37 for (i = 0; i < nmat + 2; i++) PetscCall(MatSetRandom(A[i], rctx));
38
39 PetscCall(MatCreateVecs(A[1], &x, &y));
40 PetscCall(VecDuplicate(y, &z));
41 PetscCall(VecDuplicate(z, &z2));
42 PetscCall(MatCreateVecs(A[0], &v, NULL));
43 PetscCall(VecDuplicate(v, &v2));
44
45 /* Test MatMult of an ADDITIVE MatComposite B made up of A[1],A[2],A[3] with separate scalings */
46
47 /* Do MatMult with A[1],A[2],A[3] by hand and store the result in z */
48 PetscCall(VecSet(x, 1.0));
49 PetscCall(MatMult(A[1], x, z));
50 PetscCall(VecScale(z, scalings[1]));
51 for (i = 2; i < nmat + 1; i++) {
52 PetscCall(MatMult(A[i], x, z2));
53 PetscCall(VecAXPY(z, scalings[i], z2));
54 }
55
56 /* Do MatMult using MatComposite and store the result in y */
57 PetscCall(VecSet(y, 0.0));
58 PetscCall(MatCreateComposite(PETSC_COMM_WORLD, nmat, A + 1, &B));
59 PetscCall(MatSetFromOptions(B));
60 PetscCall(MatCompositeSetScalings(B, &scalings[1]));
61 PetscCall(MatMultAdd(B, x, y, y));
62
63 /* Diff y and z */
64 PetscCall(VecAXPY(y, -1.0, z));
65 PetscCall(VecNorm(y, NORM_2, &rnorm));
66 if (rnorm > 10000.0 * PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with composite add %g\n", (double)rnorm));
67
68 /* Test MatCompositeMerge on ADDITIVE MatComposite */
69 PetscCall(MatCompositeSetMatStructure(B, DIFFERENT_NONZERO_PATTERN)); /* default */
70 PetscCall(MatCompositeMerge(B));
71 PetscCall(MatMult(B, x, y));
72 PetscCall(MatDestroy(&B));
73 PetscCall(VecAXPY(y, -1.0, z));
74 PetscCall(VecNorm(y, NORM_2, &rnorm));
75 if (rnorm > 10000.0 * PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with composite add after merge %g\n", (double)rnorm));
76
77 /*
78 Test n x n/2 multiplicative composite B made up of A[0],A[1],A[2] with separate scalings
79 */
80
81 /* Do MatMult with A[0],A[1],A[2] by hand and store the result in z */
82 PetscCall(VecSet(v, 1.0));
83 PetscCall(MatMult(A[0], v, z));
84 PetscCall(VecScale(z, scalings[0]));
85 for (i = 1; i < nmat; i++) {
86 PetscCall(MatMult(A[i], z, y));
87 PetscCall(VecScale(y, scalings[i]));
88 PetscCall(VecCopy(y, z));
89 }
90
91 /* Do MatMult using MatComposite and store the result in y */
92 PetscCall(MatCreateComposite(PETSC_COMM_WORLD, nmat, A, &B));
93 PetscCall(MatCompositeSetType(B, MAT_COMPOSITE_MULTIPLICATIVE));
94 PetscCall(MatCompositeSetMergeType(B, MAT_COMPOSITE_MERGE_LEFT));
95 PetscCall(MatSetFromOptions(B));
96 PetscCall(MatCompositeSetScalings(B, &scalings[0]));
97 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
98 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); /* do MatCompositeMerge() if -mat_composite_merge 1 */
99 PetscCall(MatMult(B, v, y));
100 PetscCall(MatDestroy(&B));
101
102 /* Diff y and z */
103 PetscCall(VecAXPY(y, -1.0, z));
104 PetscCall(VecNorm(y, NORM_2, &rnorm));
105 if (rnorm > 10000.0 * PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with composite multiplicative %g\n", (double)rnorm));
106
107 /*
108 Test n/2 x n multiplicative composite B made up of A[2], A[3], A[4] without separate scalings
109 */
110 PetscCall(VecSet(x, 1.0));
111 PetscCall(MatMult(A[2], x, z));
112 for (i = 3; i < nmat + 1; i++) {
113 PetscCall(MatMult(A[i], z, y));
114 PetscCall(VecCopy(y, z));
115 }
116 PetscCall(MatMult(A[nmat + 1], z, v));
117
118 PetscCall(MatCreateComposite(PETSC_COMM_WORLD, nmat, A + 2, &B));
119 PetscCall(MatCompositeSetType(B, MAT_COMPOSITE_MULTIPLICATIVE));
120 PetscCall(MatSetFromOptions(B));
121 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
122 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); /* do MatCompositeMerge() if -mat_composite_merge 1 */
123 PetscCall(MatMult(B, x, v2));
124 PetscCall(MatDestroy(&B));
125
126 PetscCall(VecAXPY(v2, -1.0, v));
127 PetscCall(VecNorm(v2, NORM_2, &rnorm));
128 if (rnorm > 10000.0 * PETSC_MACHINE_EPSILON) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with composite multiplicative %g\n", (double)rnorm));
129
130 /*
131 Test get functions
132 */
133 PetscCall(MatCreateComposite(PETSC_COMM_WORLD, nmat, A, &B));
134 PetscCall(MatCompositeGetNumberMat(B, &n));
135 if (nmat != n) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with GetNumberMat %" PetscInt_FMT " != %" PetscInt_FMT "\n", nmat, n));
136 PetscCall(MatCompositeGetMat(B, 0, &A[nmat + 2]));
137 if (A[0] != A[nmat + 2]) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with GetMat\n"));
138 PetscCall(MatCompositeGetType(B, &type));
139 if (type != MAT_COMPOSITE_ADDITIVE) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error with GetType\n"));
140 PetscCall(MatDestroy(&B));
141
142 /*
143 Free work space. All PETSc objects should be destroyed when they
144 are no longer needed.
145 */
146 PetscCall(VecDestroy(&x));
147 PetscCall(VecDestroy(&y));
148 PetscCall(VecDestroy(&v));
149 PetscCall(VecDestroy(&v2));
150 PetscCall(VecDestroy(&z));
151 PetscCall(VecDestroy(&z2));
152 PetscCall(PetscRandomDestroy(&rctx));
153 for (i = 0; i < nmat + 2; i++) PetscCall(MatDestroy(&A[i]));
154 PetscCall(PetscFree(A));
155
156 PetscCall(PetscFinalize());
157 return 0;
158 }
159
160 /*TEST
161
162 test:
163 nsize: 2
164 requires: double
165 args: -mat_composite_merge {{0 1}shared output} -mat_composite_merge_mvctx {{0 1}shared output}
166 output_file: output/empty.out
167
168 TEST*/
169