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: MatShellGetContext(A, &user);
16: MatView(user->B, viewer);
17: return 0;
18: }
20: static PetscErrorCode MatMult_User(Mat A, Vec X, Vec Y)
21: {
22: User user;
24: MatShellGetContext(A, &user);
25: MatMult(user->B, X, Y);
26: return 0;
27: }
29: static PetscErrorCode MatMultTranspose_User(Mat A, Vec X, Vec Y)
30: {
31: User user;
33: MatShellGetContext(A, &user);
34: MatMultTranspose(user->B, X, Y);
35: return 0;
36: }
38: static PetscErrorCode MatGetDiagonal_User(Mat A, Vec X)
39: {
40: User user;
42: MatShellGetContext(A, &user);
43: MatGetDiagonal(user->B, X);
44: return 0;
45: }
47: static PetscErrorCode TestMatrix(Mat A, Vec X, Vec Y, Vec Z)
48: {
49: Vec W1, W2, diff;
50: Mat E;
51: const char *mattypename;
52: PetscViewer viewer = PETSC_VIEWER_STDOUT_WORLD;
53: PetscScalar diag[2] = {2.9678190300000000e+08, 1.4173141580000000e+09};
54: PetscScalar multadd[2] = {-6.8966198500000000e+08, -2.0310609940000000e+09};
55: PetscScalar multtadd[2] = {-9.1052873900000000e+08, -1.8101942400000000e+09};
56: PetscReal nrm;
58: PetscObjectGetType((PetscObject)A, &mattypename);
59: PetscViewerASCIIPrintf(viewer, "\nMatrix of type: %s\n", mattypename);
60: VecDuplicate(X, &W1);
61: VecDuplicate(X, &W2);
62: MatScale(A, 31);
63: MatShift(A, 37);
64: MatDiagonalScale(A, X, Y);
65: MatScale(A, 41);
66: MatDiagonalScale(A, Y, Z);
67: MatComputeOperator(A, MATDENSE, &E);
69: MatView(E, viewer);
70: PetscViewerASCIIPrintf(viewer, "Testing MatMult + MatMultTranspose\n");
71: MatMult(A, Z, W1);
72: MatMultTranspose(A, W1, W2);
73: VecView(W2, viewer);
75: PetscViewerASCIIPrintf(viewer, "Testing MatMultAdd\n");
76: VecCreateSeqWithArray(PETSC_COMM_SELF, 1, 2, multadd, &diff);
77: VecSet(W1, -1.0);
78: MatMultAdd(A, W1, W1, W2);
79: VecView(W2, viewer);
80: VecAXPY(W2, -1.0, diff);
81: VecNorm(W2, NORM_2, &nrm);
82: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
84: #endif
86: VecSet(W2, -1.0);
87: MatMultAdd(A, W1, W2, W2);
88: VecView(W2, viewer);
89: VecAXPY(W2, -1.0, diff);
90: VecNorm(W2, NORM_2, &nrm);
91: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
93: #endif
94: VecDestroy(&diff);
96: PetscViewerASCIIPrintf(viewer, "Testing MatMultTransposeAdd\n");
97: VecCreateSeqWithArray(PETSC_COMM_SELF, 1, 2, multtadd, &diff);
99: VecSet(W1, -1.0);
100: MatMultTransposeAdd(A, W1, W1, W2);
101: VecView(W2, viewer);
102: VecAXPY(W2, -1.0, diff);
103: VecNorm(W2, NORM_2, &nrm);
104: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
106: #endif
108: VecSet(W2, -1.0);
109: MatMultTransposeAdd(A, W1, W2, W2);
110: VecView(W2, viewer);
111: VecAXPY(W2, -1.0, diff);
112: VecNorm(W2, NORM_2, &nrm);
113: #if defined(PETSC_USE_REAL_DOUBLE) || defined(PETSC_USE_REAL___FLOAT128)
115: #endif
116: VecDestroy(&diff);
118: PetscViewerASCIIPrintf(viewer, "Testing MatGetDiagonal\n");
119: MatGetDiagonal(A, W2);
120: VecView(W2, viewer);
121: VecCreateSeqWithArray(PETSC_COMM_SELF, 1, 2, diag, &diff);
122: VecAXPY(diff, -1.0, W2);
123: VecNorm(diff, NORM_2, &nrm);
125: VecDestroy(&diff);
127: /* MATSHELL does not support MatDiagonalSet after MatScale */
128: if (strncmp(mattypename, "shell", 5)) {
129: MatDiagonalSet(A, X, INSERT_VALUES);
130: MatGetDiagonal(A, W1);
131: VecView(W1, viewer);
132: } else {
133: PetscViewerASCIIPrintf(viewer, "MatDiagonalSet not tested on MATSHELL\n");
134: }
136: MatDestroy(&E);
137: VecDestroy(&W1);
138: VecDestroy(&W2);
139: return 0;
140: }
142: int main(int argc, char **args)
143: {
144: const PetscScalar xvals[] = {11, 13}, yvals[] = {17, 19}, zvals[] = {23, 29};
145: const PetscInt inds[] = {0, 1};
146: PetscScalar avals[] = {2, 3, 5, 7};
147: Mat A, S, D[4], N;
148: Vec X, Y, Z;
149: User user;
150: PetscInt i;
153: PetscInitialize(&argc, &args, (char *)0, help);
154: MatCreateSeqAIJ(PETSC_COMM_WORLD, 2, 2, 2, NULL, &A);
155: MatSetUp(A);
156: MatSetValues(A, 2, inds, 2, inds, avals, INSERT_VALUES);
157: VecCreateSeq(PETSC_COMM_WORLD, 2, &X);
158: VecDuplicate(X, &Y);
159: VecDuplicate(X, &Z);
160: VecSetValues(X, 2, inds, xvals, INSERT_VALUES);
161: VecSetValues(Y, 2, inds, yvals, INSERT_VALUES);
162: VecSetValues(Z, 2, inds, zvals, INSERT_VALUES);
163: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
164: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
165: VecAssemblyBegin(X);
166: VecAssemblyBegin(Y);
167: VecAssemblyBegin(Z);
168: VecAssemblyEnd(X);
169: VecAssemblyEnd(Y);
170: VecAssemblyEnd(Z);
172: PetscNew(&user);
173: user->B = A;
175: MatCreateShell(PETSC_COMM_WORLD, 2, 2, 2, 2, user, &S);
176: MatSetUp(S);
177: MatShellSetOperation(S, MATOP_VIEW, (void (*)(void))MatView_User);
178: MatShellSetOperation(S, MATOP_MULT, (void (*)(void))MatMult_User);
179: MatShellSetOperation(S, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultTranspose_User);
180: MatShellSetOperation(S, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_User);
182: for (i = 0; i < 4; i++) MatCreateSeqDense(PETSC_COMM_WORLD, 1, 1, &avals[i], &D[i]);
183: MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, D, &N);
184: MatSetUp(N);
186: TestMatrix(S, X, Y, Z);
187: TestMatrix(A, X, Y, Z);
188: TestMatrix(N, X, Y, Z);
190: for (i = 0; i < 4; i++) MatDestroy(&D[i]);
191: MatDestroy(&A);
192: MatDestroy(&S);
193: MatDestroy(&N);
194: VecDestroy(&X);
195: VecDestroy(&Y);
196: VecDestroy(&Z);
197: PetscFree(user);
198: PetscFinalize();
199: return 0;
200: }
202: /*TEST
204: test:
206: TEST*/