Actual source code: ex97.c
1: static const char help[] = "Tests MatCreateSubMatrix with MatSubMatrix versus MatAIJ, non-square\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, 5, 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, 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);
60: MatMult(A, X, ltmp[0]);
61: MatMult(B, X, ltmp[1]);
62: Compare2(ltmp, "MatMult");
64: MatMultTranspose(A, Y, rtmp[0]);
65: MatMultTranspose(B, Y, rtmp[1]);
66: Compare2(rtmp, "MatMultTranspose");
68: VecCopy(Y1, ltmp[0]);
69: VecCopy(Y1, ltmp[1]);
70: MatMultAdd(A, X, ltmp[0], ltmp[0]);
71: MatMultAdd(B, X, ltmp[1], ltmp[1]);
72: Compare2(ltmp, "MatMultAdd v2==v3");
74: MatMultAdd(A, X, Y1, ltmp[0]);
75: MatMultAdd(B, X, Y1, ltmp[1]);
76: Compare2(ltmp, "MatMultAdd v2!=v3");
78: VecCopy(X1, rtmp[0]);
79: VecCopy(X1, rtmp[1]);
80: MatMultTransposeAdd(A, Y, rtmp[0], rtmp[0]);
81: MatMultTransposeAdd(B, Y, rtmp[1], rtmp[1]);
82: Compare2(rtmp, "MatMultTransposeAdd v2==v3");
84: MatMultTransposeAdd(A, Y, X1, rtmp[0]);
85: MatMultTransposeAdd(B, Y, X1, rtmp[1]);
86: Compare2(rtmp, "MatMultTransposeAdd v2!=v3");
88: VecDestroyVecs(2, <mp);
89: VecDestroyVecs(2, &rtmp);
90: return 0;
91: }
93: int main(int argc, char *argv[])
94: {
95: Mat A, B, Asub, Bsub;
96: PetscInt ms, idxrow[3], idxcol[4];
97: Vec left, right, X, Y, X1, Y1;
98: IS isrow, iscol;
99: PetscBool random = PETSC_TRUE;
102: PetscInitialize(&argc, &argv, NULL, help);
103: AssembleMatrix(PETSC_COMM_WORLD, &A);
104: AssembleMatrix(PETSC_COMM_WORLD, &B);
105: MatSetOperation(B, MATOP_CREATE_SUBMATRIX, NULL);
106: MatSetOperation(B, MATOP_CREATE_SUBMATRICES, NULL);
107: MatGetOwnershipRange(A, &ms, NULL);
109: idxrow[0] = ms + 1;
110: idxrow[1] = ms + 2;
111: idxrow[2] = ms + 4;
112: ISCreateGeneral(PETSC_COMM_WORLD, 3, idxrow, PETSC_USE_POINTER, &isrow);
114: idxcol[0] = ms + 1;
115: idxcol[1] = ms + 2;
116: idxcol[2] = ms + 4;
117: idxcol[3] = ms + 5;
118: ISCreateGeneral(PETSC_COMM_WORLD, 4, idxcol, PETSC_USE_POINTER, &iscol);
120: MatCreateSubMatrix(A, isrow, iscol, MAT_INITIAL_MATRIX, &Asub);
121: MatCreateSubMatrix(B, isrow, iscol, MAT_INITIAL_MATRIX, &Bsub);
123: MatCreateVecs(Asub, &right, &left);
124: VecDuplicate(right, &X);
125: VecDuplicate(right, &X1);
126: VecDuplicate(left, &Y);
127: VecDuplicate(left, &Y1);
129: PetscOptionsGetBool(NULL, NULL, "-random", &random, NULL);
130: if (random) {
131: VecSetRandom(right, NULL);
132: VecSetRandom(left, NULL);
133: VecSetRandom(X, NULL);
134: VecSetRandom(Y, NULL);
135: VecSetRandom(X1, NULL);
136: VecSetRandom(Y1, NULL);
137: } else {
138: VecSet(right, 1.0);
139: VecSet(left, 2.0);
140: VecSet(X, 3.0);
141: VecSet(Y, 4.0);
142: VecSet(X1, 3.0);
143: VecSet(Y1, 4.0);
144: }
145: CheckMatrices(Asub, Bsub, left, right, X, Y, X1, Y1);
146: ISDestroy(&isrow);
147: ISDestroy(&iscol);
148: MatDestroy(&A);
149: MatDestroy(&B);
150: MatDestroy(&Asub);
151: MatDestroy(&Bsub);
152: VecDestroy(&left);
153: VecDestroy(&right);
154: VecDestroy(&X);
155: VecDestroy(&Y);
156: VecDestroy(&X1);
157: VecDestroy(&Y1);
158: PetscFinalize();
159: return 0;
160: }
162: /*TEST
164: test:
165: nsize: 3
167: TEST*/