Actual source code: ex99.c
1: static const char help[] = "Tests MatCreateSubMatrix with MatSubMatrix versus MatAIJ, square, shifted (copied from ex97)\n";
3: #include <petscmat.h>
5: static PetscErrorCode AssembleMatrix(MPI_Comm comm, Mat *A)
6: {
7: Mat B;
8: PetscInt i, ms, me;
10: MatCreate(comm, &B);
11: MatSetSizes(B, 6, 6, PETSC_DETERMINE, PETSC_DETERMINE);
12: MatSetFromOptions(B);
13: MatSetUp(B);
14: MatGetOwnershipRange(B, &ms, &me);
15: for (i = ms; i < me; i++) MatSetValue(B, i, i, 1.0 * i, INSERT_VALUES);
16: MatSetValue(B, me - 1, me - 1, me * me, INSERT_VALUES);
17: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
18: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
19: *A = B;
20: return 0;
21: }
23: static PetscErrorCode Compare2(Vec *X, const char *test)
24: {
25: PetscReal norm;
26: Vec Y;
27: PetscInt verbose = 0;
29: VecDuplicate(X[0], &Y);
30: VecCopy(X[0], Y);
31: VecAYPX(Y, -1.0, X[1]);
32: VecNorm(Y, NORM_INFINITY, &norm);
34: PetscOptionsGetInt(NULL, NULL, "-verbose", &verbose, NULL);
35: if (norm < PETSC_SQRT_MACHINE_EPSILON && verbose < 1) {
36: PetscPrintf(PETSC_COMM_WORLD, "%30s: norm difference < sqrt(eps_machine)\n", test);
37: } else {
38: PetscPrintf(PETSC_COMM_WORLD, "%30s: norm difference %g\n", test, (double)norm);
39: }
40: if (verbose > 1) {
41: VecView(X[0], PETSC_VIEWER_STDOUT_WORLD);
42: VecView(X[1], PETSC_VIEWER_STDOUT_WORLD);
43: VecView(Y, PETSC_VIEWER_STDOUT_WORLD);
44: }
45: VecDestroy(&Y);
46: return 0;
47: }
49: static PetscErrorCode CheckMatrices(Mat A, Mat B, Vec left, Vec right, Vec X, Vec Y, Vec X1, Vec Y1)
50: {
51: Vec *ltmp, *rtmp;
53: VecDuplicateVecs(right, 2, &rtmp);
54: VecDuplicateVecs(left, 2, <mp);
55: MatScale(A, PETSC_PI);
56: MatScale(B, PETSC_PI);
57: MatDiagonalScale(A, left, right);
58: MatDiagonalScale(B, left, right);
59: MatShift(A, PETSC_PI);
60: MatShift(B, PETSC_PI);
62: MatMult(A, X, ltmp[0]);
63: MatMult(B, X, ltmp[1]);
64: Compare2(ltmp, "MatMult");
66: MatMultTranspose(A, Y, rtmp[0]);
67: MatMultTranspose(B, Y, rtmp[1]);
68: Compare2(rtmp, "MatMultTranspose");
70: VecCopy(Y1, ltmp[0]);
71: VecCopy(Y1, ltmp[1]);
72: MatMultAdd(A, X, ltmp[0], ltmp[0]);
73: MatMultAdd(B, X, ltmp[1], ltmp[1]);
74: Compare2(ltmp, "MatMultAdd v2==v3");
76: MatMultAdd(A, X, Y1, ltmp[0]);
77: MatMultAdd(B, X, Y1, ltmp[1]);
78: Compare2(ltmp, "MatMultAdd v2!=v3");
80: VecCopy(X1, rtmp[0]);
81: VecCopy(X1, rtmp[1]);
82: MatMultTransposeAdd(A, Y, rtmp[0], rtmp[0]);
83: MatMultTransposeAdd(B, Y, rtmp[1], rtmp[1]);
84: Compare2(rtmp, "MatMultTransposeAdd v2==v3");
86: MatMultTransposeAdd(A, Y, X1, rtmp[0]);
87: MatMultTransposeAdd(B, Y, X1, rtmp[1]);
88: Compare2(rtmp, "MatMultTransposeAdd v2!=v3");
90: VecDestroyVecs(2, <mp);
91: VecDestroyVecs(2, &rtmp);
92: return 0;
93: }
95: int main(int argc, char *argv[])
96: {
97: Mat A, B, Asub, Bsub;
98: PetscInt ms, idxrow[3], idxcol[3];
99: Vec left, right, X, Y, X1, Y1;
100: IS isrow, iscol;
101: PetscBool random = PETSC_TRUE;
104: PetscInitialize(&argc, &argv, NULL, help);
105: AssembleMatrix(PETSC_COMM_WORLD, &A);
106: AssembleMatrix(PETSC_COMM_WORLD, &B);
107: MatSetOperation(B, MATOP_CREATE_SUBMATRIX, NULL);
108: MatSetOperation(B, MATOP_CREATE_SUBMATRICES, NULL);
109: MatGetOwnershipRange(A, &ms, NULL);
111: idxrow[0] = ms + 1;
112: idxrow[1] = ms + 2;
113: idxrow[2] = ms + 4;
114: ISCreateGeneral(PETSC_COMM_WORLD, 3, idxrow, PETSC_USE_POINTER, &isrow);
116: idxcol[0] = ms + 1;
117: idxcol[1] = ms + 2;
118: idxcol[2] = ms + 4;
119: ISCreateGeneral(PETSC_COMM_WORLD, 3, idxcol, PETSC_USE_POINTER, &iscol);
121: MatCreateSubMatrix(A, isrow, iscol, MAT_INITIAL_MATRIX, &Asub);
122: MatCreateSubMatrix(B, isrow, iscol, MAT_INITIAL_MATRIX, &Bsub);
124: MatCreateVecs(Asub, &right, &left);
125: VecDuplicate(right, &X);
126: VecDuplicate(right, &X1);
127: VecDuplicate(left, &Y);
128: VecDuplicate(left, &Y1);
130: PetscOptionsGetBool(NULL, NULL, "-random", &random, NULL);
131: if (random) {
132: VecSetRandom(right, NULL);
133: VecSetRandom(left, NULL);
134: VecSetRandom(X, NULL);
135: VecSetRandom(Y, NULL);
136: VecSetRandom(X1, NULL);
137: VecSetRandom(Y1, NULL);
138: } else {
139: VecSet(right, 1.0);
140: VecSet(left, 2.0);
141: VecSet(X, 3.0);
142: VecSet(Y, 4.0);
143: VecSet(X1, 3.0);
144: VecSet(Y1, 4.0);
145: }
146: CheckMatrices(Asub, Bsub, left, right, X, Y, X1, Y1);
147: ISDestroy(&isrow);
148: ISDestroy(&iscol);
149: MatDestroy(&A);
150: MatDestroy(&B);
151: MatDestroy(&Asub);
152: MatDestroy(&Bsub);
153: VecDestroy(&left);
154: VecDestroy(&right);
155: VecDestroy(&X);
156: VecDestroy(&Y);
157: VecDestroy(&X1);
158: VecDestroy(&Y1);
159: PetscFinalize();
160: return 0;
161: }
163: /*TEST
165: test:
166: nsize: 3
168: TEST*/