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: PetscFunctionBegin;
11: PetscCall(MatCreate(comm, &B));
12: PetscCall(MatSetSizes(B, 6, 6, PETSC_DETERMINE, PETSC_DETERMINE));
13: PetscCall(MatSetFromOptions(B));
14: PetscCall(MatSetUp(B));
15: PetscCall(MatGetOwnershipRange(B, &ms, &me));
16: for (i = ms; i < me; i++) PetscCall(MatSetValue(B, i, i, 1.0 * i, INSERT_VALUES));
17: PetscCall(MatSetValue(B, me - 1, me - 1, me * me, INSERT_VALUES));
18: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
19: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
20: *A = B;
21: PetscFunctionReturn(PETSC_SUCCESS);
22: }
24: static PetscErrorCode Compare2(Vec *X, const char *test)
25: {
26: PetscReal norm;
27: Vec Y;
28: PetscInt verbose = 0;
30: PetscFunctionBegin;
31: PetscCall(VecDuplicate(X[0], &Y));
32: PetscCall(VecCopy(X[0], Y));
33: PetscCall(VecAYPX(Y, -1.0, X[1]));
34: PetscCall(VecNorm(Y, NORM_INFINITY, &norm));
36: PetscCall(PetscOptionsGetInt(NULL, NULL, "-verbose", &verbose, NULL));
37: if (norm < PETSC_SQRT_MACHINE_EPSILON && verbose < 1) {
38: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%30s: norm difference < sqrt(eps_machine)\n", test));
39: } else {
40: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%30s: norm difference %g\n", test, (double)norm));
41: }
42: if (verbose > 1) {
43: PetscCall(VecView(X[0], PETSC_VIEWER_STDOUT_WORLD));
44: PetscCall(VecView(X[1], PETSC_VIEWER_STDOUT_WORLD));
45: PetscCall(VecView(Y, PETSC_VIEWER_STDOUT_WORLD));
46: }
47: PetscCall(VecDestroy(&Y));
48: PetscFunctionReturn(PETSC_SUCCESS);
49: }
51: static PetscErrorCode CheckMatrices(Mat A, Mat B, Vec left, Vec right, Vec X, Vec Y, Vec X1, Vec Y1)
52: {
53: Vec *ltmp, *rtmp;
55: PetscFunctionBegin;
56: PetscCall(VecDuplicateVecs(right, 2, &rtmp));
57: PetscCall(VecDuplicateVecs(left, 2, <mp));
58: PetscCall(MatScale(A, PETSC_PI));
59: PetscCall(MatScale(B, PETSC_PI));
60: PetscCall(MatDiagonalScale(A, left, right));
61: PetscCall(MatDiagonalScale(B, left, right));
62: PetscCall(MatShift(A, PETSC_PI));
63: PetscCall(MatShift(B, PETSC_PI));
65: PetscCall(MatMult(A, X, ltmp[0]));
66: PetscCall(MatMult(B, X, ltmp[1]));
67: PetscCall(Compare2(ltmp, "MatMult"));
69: PetscCall(MatMultTranspose(A, Y, rtmp[0]));
70: PetscCall(MatMultTranspose(B, Y, rtmp[1]));
71: PetscCall(Compare2(rtmp, "MatMultTranspose"));
73: PetscCall(VecCopy(Y1, ltmp[0]));
74: PetscCall(VecCopy(Y1, ltmp[1]));
75: PetscCall(MatMultAdd(A, X, ltmp[0], ltmp[0]));
76: PetscCall(MatMultAdd(B, X, ltmp[1], ltmp[1]));
77: PetscCall(Compare2(ltmp, "MatMultAdd v2==v3"));
79: PetscCall(MatMultAdd(A, X, Y1, ltmp[0]));
80: PetscCall(MatMultAdd(B, X, Y1, ltmp[1]));
81: PetscCall(Compare2(ltmp, "MatMultAdd v2!=v3"));
83: PetscCall(VecCopy(X1, rtmp[0]));
84: PetscCall(VecCopy(X1, rtmp[1]));
85: PetscCall(MatMultTransposeAdd(A, Y, rtmp[0], rtmp[0]));
86: PetscCall(MatMultTransposeAdd(B, Y, rtmp[1], rtmp[1]));
87: PetscCall(Compare2(rtmp, "MatMultTransposeAdd v2==v3"));
89: PetscCall(MatMultTransposeAdd(A, Y, X1, rtmp[0]));
90: PetscCall(MatMultTransposeAdd(B, Y, X1, rtmp[1]));
91: PetscCall(Compare2(rtmp, "MatMultTransposeAdd v2!=v3"));
93: PetscCall(VecDestroyVecs(2, <mp));
94: PetscCall(VecDestroyVecs(2, &rtmp));
95: PetscFunctionReturn(PETSC_SUCCESS);
96: }
98: int main(int argc, char *argv[])
99: {
100: Mat A, B, Asub, Bsub;
101: PetscInt ms, idxrow[3], idxcol[3];
102: Vec left, right, X, Y, X1, Y1;
103: IS isrow, iscol;
104: PetscBool random = PETSC_TRUE;
106: PetscFunctionBeginUser;
107: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
108: PetscCall(AssembleMatrix(PETSC_COMM_WORLD, &A));
109: PetscCall(AssembleMatrix(PETSC_COMM_WORLD, &B));
110: PetscCall(MatSetOperation(B, MATOP_CREATE_SUBMATRIX, NULL));
111: PetscCall(MatSetOperation(B, MATOP_CREATE_SUBMATRICES, NULL));
112: PetscCall(MatGetOwnershipRange(A, &ms, NULL));
114: idxrow[0] = ms + 1;
115: idxrow[1] = ms + 2;
116: idxrow[2] = ms + 4;
117: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, 3, idxrow, PETSC_USE_POINTER, &isrow));
119: idxcol[0] = ms + 1;
120: idxcol[1] = ms + 2;
121: idxcol[2] = ms + 4;
122: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, 3, idxcol, PETSC_USE_POINTER, &iscol));
124: PetscCall(MatCreateSubMatrix(A, isrow, iscol, MAT_INITIAL_MATRIX, &Asub));
125: PetscCall(MatCreateSubMatrix(B, isrow, iscol, MAT_INITIAL_MATRIX, &Bsub));
127: PetscCall(MatCreateVecs(Asub, &right, &left));
128: PetscCall(VecDuplicate(right, &X));
129: PetscCall(VecDuplicate(right, &X1));
130: PetscCall(VecDuplicate(left, &Y));
131: PetscCall(VecDuplicate(left, &Y1));
133: PetscCall(PetscOptionsGetBool(NULL, NULL, "-random", &random, NULL));
134: if (random) {
135: PetscCall(VecSetRandom(right, NULL));
136: PetscCall(VecSetRandom(left, NULL));
137: PetscCall(VecSetRandom(X, NULL));
138: PetscCall(VecSetRandom(Y, NULL));
139: PetscCall(VecSetRandom(X1, NULL));
140: PetscCall(VecSetRandom(Y1, NULL));
141: } else {
142: PetscCall(VecSet(right, 1.0));
143: PetscCall(VecSet(left, 2.0));
144: PetscCall(VecSet(X, 3.0));
145: PetscCall(VecSet(Y, 4.0));
146: PetscCall(VecSet(X1, 3.0));
147: PetscCall(VecSet(Y1, 4.0));
148: }
149: PetscCall(CheckMatrices(Asub, Bsub, left, right, X, Y, X1, Y1));
150: PetscCall(ISDestroy(&isrow));
151: PetscCall(ISDestroy(&iscol));
152: PetscCall(MatDestroy(&A));
153: PetscCall(MatDestroy(&B));
154: PetscCall(MatDestroy(&Asub));
155: PetscCall(MatDestroy(&Bsub));
156: PetscCall(VecDestroy(&left));
157: PetscCall(VecDestroy(&right));
158: PetscCall(VecDestroy(&X));
159: PetscCall(VecDestroy(&Y));
160: PetscCall(VecDestroy(&X1));
161: PetscCall(VecDestroy(&Y1));
162: PetscCall(PetscFinalize());
163: return 0;
164: }
166: /*TEST
168: test:
169: nsize: 3
171: TEST*/