Actual source code: ex88.c
2: static char help[] = "Tests MatShift(), MatScale(), and MatDiagonalScale() for SHELL and NEST matrices\n\n";
4: #include <petscmat.h>
6: typedef struct _n_User *User;
7: struct _n_User {
8: Mat B;
9: };
11: static PetscErrorCode MatView_User(Mat A, PetscViewer viewer)
12: {
13: User user;
15: PetscFunctionBegin;
16: PetscCall(MatShellGetContext(A, &user));
17: PetscCall(MatView(user->B, viewer));
18: PetscFunctionReturn(PETSC_SUCCESS);
19: }
21: static PetscErrorCode MatMult_User(Mat A, Vec X, Vec Y)
22: {
23: User user;
25: PetscFunctionBegin;
26: PetscCall(MatShellGetContext(A, &user));
27: PetscCall(MatMult(user->B, X, Y));
28: PetscFunctionReturn(PETSC_SUCCESS);
29: }
31: static PetscErrorCode MatMultTranspose_User(Mat A, Vec X, Vec Y)
32: {
33: User user;
35: PetscFunctionBegin;
36: PetscCall(MatShellGetContext(A, &user));
37: PetscCall(MatMultTranspose(user->B, X, Y));
38: PetscFunctionReturn(PETSC_SUCCESS);
39: }
41: static PetscErrorCode MatGetDiagonal_User(Mat A, Vec X)
42: {
43: User user;
45: PetscFunctionBegin;
46: PetscCall(MatShellGetContext(A, &user));
47: PetscCall(MatGetDiagonal(user->B, X));
48: PetscFunctionReturn(PETSC_SUCCESS);
49: }
51: static PetscErrorCode TestMatrix(Mat A, Vec X, Vec Y, Vec Z)
52: {
53: Vec W1, W2, diff;
54: Mat E;
55: const char *mattypename;
56: PetscViewer viewer = PETSC_VIEWER_STDOUT_WORLD;
57: PetscScalar diag[2] = {2.9678190300000000e+08, 1.4173141580000000e+09};
58: PetscScalar multadd[2] = {-6.8966198500000000e+08, -2.0310609940000000e+09};
59: PetscScalar multtadd[2] = {-9.1052873900000000e+08, -1.8101942400000000e+09};
60: PetscReal nrm;
62: PetscFunctionBegin;
63: PetscCall(PetscObjectGetType((PetscObject)A, &mattypename));
64: PetscCall(PetscViewerASCIIPrintf(viewer, "\nMatrix of type: %s\n", mattypename));
65: PetscCall(VecDuplicate(X, &W1));
66: PetscCall(VecDuplicate(X, &W2));
67: PetscCall(MatScale(A, 31));
68: PetscCall(MatShift(A, 37));
69: PetscCall(MatDiagonalScale(A, X, Y));
70: PetscCall(MatScale(A, 41));
71: PetscCall(MatDiagonalScale(A, Y, Z));
72: PetscCall(MatComputeOperator(A, MATDENSE, &E));
74: PetscCall(MatView(E, viewer));
75: PetscCall(PetscViewerASCIIPrintf(viewer, "Testing MatMult + MatMultTranspose\n"));
76: PetscCall(MatMult(A, Z, W1));
77: PetscCall(MatMultTranspose(A, W1, W2));
78: PetscCall(VecView(W2, viewer));
80: PetscCall(PetscViewerASCIIPrintf(viewer, "Testing MatMultAdd\n"));
81: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, 2, multadd, &diff));
82: PetscCall(VecSet(W1, -1.0));
83: PetscCall(MatMultAdd(A, W1, W1, W2));
84: PetscCall(VecView(W2, viewer));
85: PetscCall(VecAXPY(W2, -1.0, diff));
86: PetscCall(VecNorm(W2, NORM_2, &nrm));
87: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
88: PetscCheck(nrm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatMultAdd(A,x,x,y) produces incorrect result");
89: #endif
91: PetscCall(VecSet(W2, -1.0));
92: PetscCall(MatMultAdd(A, W1, W2, W2));
93: PetscCall(VecView(W2, viewer));
94: PetscCall(VecAXPY(W2, -1.0, diff));
95: PetscCall(VecNorm(W2, NORM_2, &nrm));
96: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
97: PetscCheck(nrm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatMultAdd(A,x,y,y) produces incorrect result");
98: #endif
99: PetscCall(VecDestroy(&diff));
101: PetscCall(PetscViewerASCIIPrintf(viewer, "Testing MatMultTransposeAdd\n"));
102: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, 2, multtadd, &diff));
104: PetscCall(VecSet(W1, -1.0));
105: PetscCall(MatMultTransposeAdd(A, W1, W1, W2));
106: PetscCall(VecView(W2, viewer));
107: PetscCall(VecAXPY(W2, -1.0, diff));
108: PetscCall(VecNorm(W2, NORM_2, &nrm));
109: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
110: PetscCheck(nrm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatMultTransposeAdd(A,x,x,y) produces incorrect result");
111: #endif
113: PetscCall(VecSet(W2, -1.0));
114: PetscCall(MatMultTransposeAdd(A, W1, W2, W2));
115: PetscCall(VecView(W2, viewer));
116: PetscCall(VecAXPY(W2, -1.0, diff));
117: PetscCall(VecNorm(W2, NORM_2, &nrm));
118: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
119: PetscCheck(nrm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatMultTransposeAdd(A,x,y,y) produces incorrect result");
120: #endif
121: PetscCall(VecDestroy(&diff));
123: PetscCall(PetscViewerASCIIPrintf(viewer, "Testing MatGetDiagonal\n"));
124: PetscCall(MatGetDiagonal(A, W2));
125: PetscCall(VecView(W2, viewer));
126: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, 2, diag, &diff));
127: PetscCall(VecAXPY(diff, -1.0, W2));
128: PetscCall(VecNorm(diff, NORM_2, &nrm));
129: PetscCheck(nrm <= PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_PLIB, "MatGetDiagonal() produces incorrect result");
130: PetscCall(VecDestroy(&diff));
132: /* MATSHELL does not support MatDiagonalSet after MatScale */
133: if (strncmp(mattypename, "shell", 5)) {
134: PetscCall(MatDiagonalSet(A, X, INSERT_VALUES));
135: PetscCall(MatGetDiagonal(A, W1));
136: PetscCall(VecView(W1, viewer));
137: } else {
138: PetscCall(PetscViewerASCIIPrintf(viewer, "MatDiagonalSet not tested on MATSHELL\n"));
139: }
141: PetscCall(MatDestroy(&E));
142: PetscCall(VecDestroy(&W1));
143: PetscCall(VecDestroy(&W2));
144: PetscFunctionReturn(PETSC_SUCCESS);
145: }
147: int main(int argc, char **args)
148: {
149: const PetscScalar xvals[] = {11, 13}, yvals[] = {17, 19}, zvals[] = {23, 29};
150: const PetscInt inds[] = {0, 1};
151: PetscScalar avals[] = {2, 3, 5, 7};
152: Mat A, S, D[4], N;
153: Vec X, Y, Z;
154: User user;
155: PetscInt i;
157: PetscFunctionBeginUser;
158: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
159: PetscCall(MatCreateSeqAIJ(PETSC_COMM_WORLD, 2, 2, 2, NULL, &A));
160: PetscCall(MatSetUp(A));
161: PetscCall(MatSetValues(A, 2, inds, 2, inds, avals, INSERT_VALUES));
162: PetscCall(VecCreateSeq(PETSC_COMM_WORLD, 2, &X));
163: PetscCall(VecDuplicate(X, &Y));
164: PetscCall(VecDuplicate(X, &Z));
165: PetscCall(VecSetValues(X, 2, inds, xvals, INSERT_VALUES));
166: PetscCall(VecSetValues(Y, 2, inds, yvals, INSERT_VALUES));
167: PetscCall(VecSetValues(Z, 2, inds, zvals, INSERT_VALUES));
168: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
169: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
170: PetscCall(VecAssemblyBegin(X));
171: PetscCall(VecAssemblyBegin(Y));
172: PetscCall(VecAssemblyBegin(Z));
173: PetscCall(VecAssemblyEnd(X));
174: PetscCall(VecAssemblyEnd(Y));
175: PetscCall(VecAssemblyEnd(Z));
177: PetscCall(PetscNew(&user));
178: user->B = A;
180: PetscCall(MatCreateShell(PETSC_COMM_WORLD, 2, 2, 2, 2, user, &S));
181: PetscCall(MatSetUp(S));
182: PetscCall(MatShellSetOperation(S, MATOP_VIEW, (void (*)(void))MatView_User));
183: PetscCall(MatShellSetOperation(S, MATOP_MULT, (void (*)(void))MatMult_User));
184: PetscCall(MatShellSetOperation(S, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultTranspose_User));
185: PetscCall(MatShellSetOperation(S, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_User));
187: for (i = 0; i < 4; i++) PetscCall(MatCreateSeqDense(PETSC_COMM_WORLD, 1, 1, &avals[i], &D[i]));
188: PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, D, &N));
189: PetscCall(MatSetUp(N));
191: PetscCall(TestMatrix(S, X, Y, Z));
192: PetscCall(TestMatrix(A, X, Y, Z));
193: PetscCall(TestMatrix(N, X, Y, Z));
195: for (i = 0; i < 4; i++) PetscCall(MatDestroy(&D[i]));
196: PetscCall(MatDestroy(&A));
197: PetscCall(MatDestroy(&S));
198: PetscCall(MatDestroy(&N));
199: PetscCall(VecDestroy(&X));
200: PetscCall(VecDestroy(&Y));
201: PetscCall(VecDestroy(&Z));
202: PetscCall(PetscFree(user));
203: PetscCall(PetscFinalize());
204: return 0;
205: }
207: /*TEST
209: test:
211: TEST*/