1 static char help[] = "Tests MatShift(), MatScale(), and MatDiagonalScale() for SHELL and NEST matrices\n\n";
2
3 #include <petscmat.h>
4
5 typedef struct _n_User *User;
6 struct _n_User {
7 Mat B;
8 };
9
MatView_User(Mat A,PetscViewer viewer)10 static PetscErrorCode MatView_User(Mat A, PetscViewer viewer)
11 {
12 User user;
13
14 PetscFunctionBegin;
15 PetscCall(MatShellGetContext(A, &user));
16 PetscCall(MatView(user->B, viewer));
17 PetscFunctionReturn(PETSC_SUCCESS);
18 }
19
MatMult_User(Mat A,Vec X,Vec Y)20 static PetscErrorCode MatMult_User(Mat A, Vec X, Vec Y)
21 {
22 User user;
23
24 PetscFunctionBegin;
25 PetscCall(MatShellGetContext(A, &user));
26 PetscCall(MatMult(user->B, X, Y));
27 PetscFunctionReturn(PETSC_SUCCESS);
28 }
29
MatMultTranspose_User(Mat A,Vec X,Vec Y)30 static PetscErrorCode MatMultTranspose_User(Mat A, Vec X, Vec Y)
31 {
32 User user;
33
34 PetscFunctionBegin;
35 PetscCall(MatShellGetContext(A, &user));
36 PetscCall(MatMultTranspose(user->B, X, Y));
37 PetscFunctionReturn(PETSC_SUCCESS);
38 }
39
MatGetDiagonal_User(Mat A,Vec X)40 static PetscErrorCode MatGetDiagonal_User(Mat A, Vec X)
41 {
42 User user;
43
44 PetscFunctionBegin;
45 PetscCall(MatShellGetContext(A, &user));
46 PetscCall(MatGetDiagonal(user->B, X));
47 PetscFunctionReturn(PETSC_SUCCESS);
48 }
49
TestMatrix(Mat A,Vec X,Vec Y,Vec Z)50 static PetscErrorCode TestMatrix(Mat A, Vec X, Vec Y, Vec Z)
51 {
52 Vec W1, W2, W3, diff;
53 Mat E;
54 const char *mattypename;
55 PetscViewer viewer = PETSC_VIEWER_STDOUT_WORLD;
56 PetscReal nrm;
57 #if defined(PETSC_USE_COMPLEX)
58 const PetscScalar diag[2] = {PetscCMPLX(-6.2902938000000000e+07, 4.5741953400000000e+08), PetscCMPLX(1.0828994620000000e+09, 1.2955916360000000e+09)};
59 const PetscScalar multadd[2] = {PetscCMPLX(1.4926230300000000e+08, -1.2811063360000000e+09), PetscCMPLX(-1.2985220710000000e+09, -2.1029893020000000e+09)};
60 const PetscScalar multtadd[2] = {PetscCMPLX(-1.5271967100000000e+08, -1.2648172000000000e+09), PetscCMPLX(-9.9654009700000000e+08, -2.1192784380000000e+09)};
61 #else
62 const PetscScalar diag[2] = {2.9678190300000000e+08, 1.4173141580000000e+09};
63 const PetscScalar multadd[2] = {-6.8966198500000000e+08, -2.0310609940000000e+09};
64 const PetscScalar multtadd[2] = {-9.1052873900000000e+08, -1.8101942400000000e+09};
65 #endif
66
67 PetscFunctionBegin;
68 PetscCall(PetscObjectGetType((PetscObject)A, &mattypename));
69 PetscCall(PetscViewerASCIIPrintf(viewer, "\nMatrix of type: %s\n", mattypename));
70 PetscCall(VecDuplicate(X, &W1));
71 PetscCall(VecDuplicate(X, &W2));
72 PetscCall(VecDuplicate(X, &W3));
73 PetscCall(MatScale(A, 31));
74 PetscCall(MatShift(A, 37));
75 PetscCall(MatDiagonalScale(A, X, Y));
76 PetscCall(MatScale(A, 41));
77 PetscCall(MatDiagonalScale(A, Y, Z));
78 PetscCall(MatComputeOperator(A, MATDENSE, &E));
79
80 PetscCall(MatView(E, viewer));
81 PetscCall(PetscViewerASCIIPrintf(viewer, "Testing MatMult + MatMultTranspose\n"));
82 PetscCall(MatMult(A, Z, W1));
83 PetscCall(MatMultTranspose(A, W1, W2));
84 PetscCall(VecView(W2, viewer));
85 PetscCall(PetscViewerASCIIPrintf(viewer, "Testing MatMultHermitianTranspose\n"));
86 PetscCall(VecConjugate(W1));
87 PetscCall(MatMultHermitianTranspose(A, W1, W2));
88 PetscCall(VecConjugate(W2));
89 PetscCall(VecView(W2, viewer));
90
91 PetscCall(PetscViewerASCIIPrintf(viewer, "Testing MatMultAdd\n"));
92 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, 2, multadd, &diff));
93 PetscCall(VecSet(W1, -1.0));
94 PetscCall(MatMultAdd(A, W1, W1, W2));
95 PetscCall(VecView(W2, viewer));
96 PetscCall(VecAXPY(W2, -1.0, diff));
97 PetscCall(VecNorm(W2, NORM_2, &nrm));
98 #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
99 PetscCheck(nrm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatMultAdd(A,x,x,y) produces incorrect result");
100 #endif
101
102 PetscCall(VecSet(W2, -1.0));
103 PetscCall(MatMultAdd(A, W1, W2, W2));
104 PetscCall(VecView(W2, viewer));
105 PetscCall(VecAXPY(W2, -1.0, diff));
106 PetscCall(VecNorm(W2, NORM_2, &nrm));
107 #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
108 PetscCheck(nrm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatMultAdd(A,x,y,y) produces incorrect result");
109 #endif
110 PetscCall(VecDestroy(&diff));
111
112 PetscCall(PetscViewerASCIIPrintf(viewer, "Testing MatMultTransposeAdd\n"));
113 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, 2, multtadd, &diff));
114
115 PetscCall(VecSet(W1, -1.0));
116 PetscCall(MatMultTransposeAdd(A, W1, W1, W2));
117 PetscCall(VecView(W2, viewer));
118 PetscCall(VecAXPY(W2, -1.0, diff));
119 PetscCall(VecNorm(W2, NORM_2, &nrm));
120 #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
121 PetscCheck(nrm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatMultTransposeAdd(A,x,x,y) produces incorrect result");
122 #endif
123
124 PetscCall(VecSet(W2, -1.0));
125 PetscCall(MatMultTransposeAdd(A, W1, W2, W2));
126 PetscCall(VecView(W2, viewer));
127 PetscCall(VecAXPY(W2, -1.0, diff));
128 PetscCall(VecNorm(W2, NORM_2, &nrm));
129 #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
130 PetscCheck(nrm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatMultTransposeAdd(A,x,y,y) produces incorrect result");
131 #endif
132 PetscCall(VecDestroy(&diff));
133
134 PetscCall(PetscViewerASCIIPrintf(viewer, "Testing MatMultHermitianTransposeAdd\n"));
135 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, 2, multtadd, &diff));
136
137 PetscCall(VecSet(W1, -1.0));
138 PetscCall(MatMultHermitianTransposeAdd(A, W1, W1, W3));
139 PetscCall(VecConjugate(W3));
140 PetscCall(VecView(W3, viewer));
141 PetscCall(VecAXPY(W3, -1.0, diff));
142 PetscCall(VecNorm(W3, NORM_2, &nrm));
143 #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
144 PetscCheck(nrm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatMultHermitianTransposeAdd(A,x,x,y) produces incorrect result");
145 #endif
146
147 PetscCall(VecSet(W3, -1.0));
148 PetscCall(MatMultHermitianTransposeAdd(A, W1, W3, W3));
149 PetscCall(VecConjugate(W3));
150 PetscCall(VecView(W3, viewer));
151 PetscCall(VecAXPY(W3, -1.0, diff));
152 PetscCall(VecNorm(W3, NORM_2, &nrm));
153 #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
154 PetscCheck(nrm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatMultHermitianTransposeAdd(A,x,y,y) produces incorrect result");
155 #endif
156 PetscCall(VecDestroy(&diff));
157
158 PetscCall(PetscViewerASCIIPrintf(viewer, "Testing MatGetDiagonal\n"));
159 PetscCall(MatGetDiagonal(A, W2));
160 PetscCall(VecView(W2, viewer));
161 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, 2, diag, &diff));
162 PetscCall(VecAXPY(diff, -1.0, W2));
163 PetscCall(VecNorm(diff, NORM_2, &nrm));
164 #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
165 PetscCheck(nrm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatGetDiagonal() produces incorrect result");
166 #endif
167 PetscCall(VecDestroy(&diff));
168
169 /* MATSHELL does not support MatDiagonalSet after MatScale */
170 if (strncmp(mattypename, "shell", 5)) {
171 PetscCall(MatDiagonalSet(A, X, INSERT_VALUES));
172 PetscCall(MatGetDiagonal(A, W1));
173 PetscCall(VecView(W1, viewer));
174 } else {
175 PetscCall(PetscViewerASCIIPrintf(viewer, "MatDiagonalSet not tested on MATSHELL\n"));
176 }
177
178 PetscCall(MatDestroy(&E));
179 PetscCall(VecDestroy(&W1));
180 PetscCall(VecDestroy(&W2));
181 PetscCall(VecDestroy(&W3));
182 PetscFunctionReturn(PETSC_SUCCESS);
183 }
184
main(int argc,char ** args)185 int main(int argc, char **args)
186 {
187 const PetscInt inds[] = {0, 1};
188 #if defined(PETSC_USE_COMPLEX)
189 const PetscScalar xvals[] = {PetscCMPLX(11, 4), PetscCMPLX(13, 2)}, yvals[] = {PetscCMPLX(17, 3), PetscCMPLX(19, 1)}, zvals[] = {PetscCMPLX(23, 6), PetscCMPLX(29, 2)};
190 PetscScalar avals[] = {PetscCMPLX(2, 3), PetscCMPLX(3, 5), PetscCMPLX(5, 4), PetscCMPLX(7, 5)};
191 #else
192 const PetscScalar xvals[] = {11, 13}, yvals[] = {17, 19}, zvals[] = {23, 29};
193 PetscScalar avals[] = {2, 3, 5, 7};
194 #endif
195 Mat A, S, D[4], N;
196 Vec X, Y, Z;
197 User user;
198 PetscInt i;
199
200 PetscFunctionBeginUser;
201 PetscCall(PetscInitialize(&argc, &args, NULL, help));
202 PetscCall(MatCreateSeqAIJ(PETSC_COMM_WORLD, 2, 2, 2, NULL, &A));
203 PetscCall(MatSetUp(A));
204 PetscCall(VecCreateSeq(PETSC_COMM_WORLD, 2, &X));
205 PetscCall(VecDuplicate(X, &Y));
206 PetscCall(VecDuplicate(X, &Z));
207 PetscCall(MatSetValues(A, 2, inds, 2, inds, avals, INSERT_VALUES));
208 PetscCall(VecSetValues(X, 2, inds, xvals, INSERT_VALUES));
209 PetscCall(VecSetValues(Y, 2, inds, yvals, INSERT_VALUES));
210 PetscCall(VecSetValues(Z, 2, inds, zvals, INSERT_VALUES));
211 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
212 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
213 PetscCall(VecAssemblyBegin(X));
214 PetscCall(VecAssemblyBegin(Y));
215 PetscCall(VecAssemblyBegin(Z));
216 PetscCall(VecAssemblyEnd(X));
217 PetscCall(VecAssemblyEnd(Y));
218 PetscCall(VecAssemblyEnd(Z));
219
220 PetscCall(PetscNew(&user));
221 user->B = A;
222
223 PetscCall(MatCreateShell(PETSC_COMM_WORLD, 2, 2, 2, 2, user, &S));
224 PetscCall(MatSetUp(S));
225 PetscCall(MatShellSetOperation(S, MATOP_VIEW, (PetscErrorCodeFn *)MatView_User));
226 PetscCall(MatShellSetOperation(S, MATOP_MULT, (PetscErrorCodeFn *)MatMult_User));
227 PetscCall(MatShellSetOperation(S, MATOP_MULT_TRANSPOSE, (PetscErrorCodeFn *)MatMultTranspose_User));
228 PetscCall(MatShellSetOperation(S, MATOP_GET_DIAGONAL, (PetscErrorCodeFn *)MatGetDiagonal_User));
229
230 for (i = 0; i < 4; i++) PetscCall(MatCreateSeqDense(PETSC_COMM_WORLD, 1, 1, &avals[i], &D[i]));
231 PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, D, &N));
232 PetscCall(MatSetUp(N));
233
234 PetscCall(TestMatrix(S, X, Y, Z));
235 PetscCall(TestMatrix(A, X, Y, Z));
236 PetscCall(TestMatrix(N, X, Y, Z));
237
238 for (i = 0; i < 4; i++) PetscCall(MatDestroy(&D[i]));
239 PetscCall(MatDestroy(&A));
240 PetscCall(MatDestroy(&S));
241 PetscCall(MatDestroy(&N));
242 PetscCall(VecDestroy(&X));
243 PetscCall(VecDestroy(&Y));
244 PetscCall(VecDestroy(&Z));
245 PetscCall(PetscFree(user));
246 PetscCall(PetscFinalize());
247 return 0;
248 }
249
250 /*TEST
251
252 testset:
253 test:
254 suffix: 1
255 requires:!complex
256 test:
257 suffix: 2
258 requires: complex
259
260 TEST*/
261