1 static const char help[] = "Tests MatCreateSubMatrix with MatSubMatrix versus MatAIJ, non-square\n";
2
3 #include <petscmat.h>
4
AssembleMatrix(MPI_Comm comm,Mat * A)5 static PetscErrorCode AssembleMatrix(MPI_Comm comm, Mat *A)
6 {
7 Mat B;
8 PetscInt i, ms, me;
9
10 PetscFunctionBegin;
11 PetscCall(MatCreate(comm, &B));
12 PetscCall(MatSetSizes(B, 5, 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, 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 }
23
Compare2(Vec * X,const char * test)24 static PetscErrorCode Compare2(Vec *X, const char *test)
25 {
26 PetscReal norm;
27 Vec Y;
28 PetscInt verbose = 0;
29
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));
35
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 }
50
CheckMatrices(Mat A,Mat B,Vec left,Vec right,Vec X,Vec Y,Vec X1,Vec Y1)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;
54
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
63 PetscCall(MatMult(A, X, ltmp[0]));
64 PetscCall(MatMult(B, X, ltmp[1]));
65 PetscCall(Compare2(ltmp, "MatMult"));
66
67 PetscCall(MatMultTranspose(A, Y, rtmp[0]));
68 PetscCall(MatMultTranspose(B, Y, rtmp[1]));
69 PetscCall(Compare2(rtmp, "MatMultTranspose"));
70
71 PetscCall(VecCopy(Y1, ltmp[0]));
72 PetscCall(VecCopy(Y1, ltmp[1]));
73 PetscCall(MatMultAdd(A, X, ltmp[0], ltmp[0]));
74 PetscCall(MatMultAdd(B, X, ltmp[1], ltmp[1]));
75 PetscCall(Compare2(ltmp, "MatMultAdd v2==v3"));
76
77 PetscCall(MatMultAdd(A, X, Y1, ltmp[0]));
78 PetscCall(MatMultAdd(B, X, Y1, ltmp[1]));
79 PetscCall(Compare2(ltmp, "MatMultAdd v2!=v3"));
80
81 PetscCall(VecCopy(X1, rtmp[0]));
82 PetscCall(VecCopy(X1, rtmp[1]));
83 PetscCall(MatMultTransposeAdd(A, Y, rtmp[0], rtmp[0]));
84 PetscCall(MatMultTransposeAdd(B, Y, rtmp[1], rtmp[1]));
85 PetscCall(Compare2(rtmp, "MatMultTransposeAdd v2==v3"));
86
87 PetscCall(MatMultTransposeAdd(A, Y, X1, rtmp[0]));
88 PetscCall(MatMultTransposeAdd(B, Y, X1, rtmp[1]));
89 PetscCall(Compare2(rtmp, "MatMultTransposeAdd v2!=v3"));
90
91 PetscCall(VecDestroyVecs(2, <mp));
92 PetscCall(VecDestroyVecs(2, &rtmp));
93 PetscFunctionReturn(PETSC_SUCCESS);
94 }
95
main(int argc,char * argv[])96 int main(int argc, char *argv[])
97 {
98 Mat A, B, Asub, Bsub;
99 PetscInt ms, idxrow[3], idxcol[4];
100 Vec left, right, X, Y, X1, Y1;
101 IS isrow, iscol;
102 PetscBool random = PETSC_TRUE;
103
104 PetscFunctionBeginUser;
105 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
106 PetscCall(AssembleMatrix(PETSC_COMM_WORLD, &A));
107 PetscCall(AssembleMatrix(PETSC_COMM_WORLD, &B));
108 PetscCall(MatSetOperation(B, MATOP_CREATE_SUBMATRIX, NULL));
109 PetscCall(MatSetOperation(B, MATOP_CREATE_SUBMATRICES, NULL));
110 PetscCall(MatGetOwnershipRange(A, &ms, NULL));
111
112 idxrow[0] = ms + 1;
113 idxrow[1] = ms + 2;
114 idxrow[2] = ms + 4;
115 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, 3, idxrow, PETSC_USE_POINTER, &isrow));
116
117 idxcol[0] = ms + 1;
118 idxcol[1] = ms + 2;
119 idxcol[2] = ms + 4;
120 idxcol[3] = ms + 5;
121 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, 4, idxcol, PETSC_USE_POINTER, &iscol));
122
123 PetscCall(MatCreateSubMatrix(A, isrow, iscol, MAT_INITIAL_MATRIX, &Asub));
124 PetscCall(MatCreateSubMatrix(B, isrow, iscol, MAT_INITIAL_MATRIX, &Bsub));
125
126 PetscCall(MatCreateVecs(Asub, &right, &left));
127 PetscCall(VecDuplicate(right, &X));
128 PetscCall(VecDuplicate(right, &X1));
129 PetscCall(VecDuplicate(left, &Y));
130 PetscCall(VecDuplicate(left, &Y1));
131
132 PetscCall(PetscOptionsGetBool(NULL, NULL, "-random", &random, NULL));
133 if (random) {
134 PetscCall(VecSetRandom(right, NULL));
135 PetscCall(VecSetRandom(left, NULL));
136 PetscCall(VecSetRandom(X, NULL));
137 PetscCall(VecSetRandom(Y, NULL));
138 PetscCall(VecSetRandom(X1, NULL));
139 PetscCall(VecSetRandom(Y1, NULL));
140 } else {
141 PetscCall(VecSet(right, 1.0));
142 PetscCall(VecSet(left, 2.0));
143 PetscCall(VecSet(X, 3.0));
144 PetscCall(VecSet(Y, 4.0));
145 PetscCall(VecSet(X1, 3.0));
146 PetscCall(VecSet(Y1, 4.0));
147 }
148 PetscCall(CheckMatrices(Asub, Bsub, left, right, X, Y, X1, Y1));
149 PetscCall(ISDestroy(&isrow));
150 PetscCall(ISDestroy(&iscol));
151 PetscCall(MatDestroy(&A));
152 PetscCall(MatDestroy(&B));
153 PetscCall(MatDestroy(&Asub));
154 PetscCall(MatDestroy(&Bsub));
155 PetscCall(VecDestroy(&left));
156 PetscCall(VecDestroy(&right));
157 PetscCall(VecDestroy(&X));
158 PetscCall(VecDestroy(&Y));
159 PetscCall(VecDestroy(&X1));
160 PetscCall(VecDestroy(&Y1));
161 PetscCall(PetscFinalize());
162 return 0;
163 }
164
165 /*TEST
166
167 test:
168 nsize: 3
169
170 TEST*/
171