Actual source code: ex88.c
1: static char help[] = "Tests MatShift(), MatScale(), and MatDiagonalScale() for SHELL and NEST matrices\n\n";
3: #include <petscmat.h>
5: typedef struct _n_User *User;
6: struct _n_User {
7: Mat B;
8: };
10: static PetscErrorCode MatView_User(Mat A, PetscViewer viewer)
11: {
12: User user;
14: PetscFunctionBegin;
15: PetscCall(MatShellGetContext(A, &user));
16: PetscCall(MatView(user->B, viewer));
17: PetscFunctionReturn(PETSC_SUCCESS);
18: }
20: static PetscErrorCode MatMult_User(Mat A, Vec X, Vec Y)
21: {
22: User user;
24: PetscFunctionBegin;
25: PetscCall(MatShellGetContext(A, &user));
26: PetscCall(MatMult(user->B, X, Y));
27: PetscFunctionReturn(PETSC_SUCCESS);
28: }
30: static PetscErrorCode MatMultTranspose_User(Mat A, Vec X, Vec Y)
31: {
32: User user;
34: PetscFunctionBegin;
35: PetscCall(MatShellGetContext(A, &user));
36: PetscCall(MatMultTranspose(user->B, X, Y));
37: PetscFunctionReturn(PETSC_SUCCESS);
38: }
40: static PetscErrorCode MatGetDiagonal_User(Mat A, Vec X)
41: {
42: User user;
44: PetscFunctionBegin;
45: PetscCall(MatShellGetContext(A, &user));
46: PetscCall(MatGetDiagonal(user->B, X));
47: PetscFunctionReturn(PETSC_SUCCESS);
48: }
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
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));
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));
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
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));
112: PetscCall(PetscViewerASCIIPrintf(viewer, "Testing MatMultTransposeAdd\n"));
113: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, 2, multtadd, &diff));
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
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));
134: PetscCall(PetscViewerASCIIPrintf(viewer, "Testing MatMultHermitianTransposeAdd\n"));
135: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, 2, multtadd, &diff));
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
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));
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));
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: }
178: PetscCall(MatDestroy(&E));
179: PetscCall(VecDestroy(&W1));
180: PetscCall(VecDestroy(&W2));
181: PetscCall(VecDestroy(&W3));
182: PetscFunctionReturn(PETSC_SUCCESS);
183: }
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;
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));
220: PetscCall(PetscNew(&user));
221: user->B = A;
223: PetscCall(MatCreateShell(PETSC_COMM_WORLD, 2, 2, 2, 2, user, &S));
224: PetscCall(MatSetUp(S));
225: PetscCall(MatShellSetOperation(S, MATOP_VIEW, (void (*)(void))MatView_User));
226: PetscCall(MatShellSetOperation(S, MATOP_MULT, (void (*)(void))MatMult_User));
227: PetscCall(MatShellSetOperation(S, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultTranspose_User));
228: PetscCall(MatShellSetOperation(S, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_User));
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));
234: PetscCall(TestMatrix(S, X, Y, Z));
235: PetscCall(TestMatrix(A, X, Y, Z));
236: PetscCall(TestMatrix(N, X, Y, Z));
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: }
250: /*TEST
252: testset:
253: test:
254: suffix: 1
255: requires:!complex
256: test:
257: suffix: 2
258: requires: complex
260: TEST*/