Actual source code: ex53.c
petsc-3.14.6 2021-03-30
2: static char help[] = "Tests various routines in MatMPIBAIJ format.\n";
5: #include <petscmat.h>
6: #define IMAX 15
7: int main(int argc,char **args)
8: {
9: Mat A,B,C,At,Bt;
10: PetscViewer fd;
11: char file[PETSC_MAX_PATH_LEN];
12: PetscRandom rand;
13: Vec xx,yy,s1,s2;
14: PetscReal s1norm,s2norm,rnorm,tol=1.e-10;
15: PetscInt rstart,rend,rows[2],cols[2],m,n,i,j,M,N,ct,row,ncols1,ncols2,bs;
16: PetscMPIInt rank,size;
17: PetscErrorCode 0;
18: const PetscInt *cols1,*cols2;
19: PetscScalar vals1[4],vals2[4],v;
20: const PetscScalar *v1,*v2;
21: PetscBool flg;
23: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
24: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
25: MPI_Comm_size(PETSC_COMM_WORLD,&size);
27: /* Check out if MatLoad() works */
28: PetscOptionsGetString(NULL,NULL,"-f",file,sizeof(file),&flg);
29: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Input file not specified");
30: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
31: MatCreate(PETSC_COMM_WORLD,&A);
32: MatSetType(A,MATBAIJ);
33: MatLoad(A,fd);
35: MatConvert(A,MATAIJ,MAT_INITIAL_MATRIX,&B);
36: PetscViewerDestroy(&fd);
38: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
39: PetscRandomSetFromOptions(rand);
40: MatGetLocalSize(A,&m,&n);
41: VecCreate(PETSC_COMM_WORLD,&xx);
42: VecSetSizes(xx,m,PETSC_DECIDE);
43: VecSetFromOptions(xx);
44: VecDuplicate(xx,&s1);
45: VecDuplicate(xx,&s2);
46: VecDuplicate(xx,&yy);
47: MatGetBlockSize(A,&bs);
49: /* Test MatNorm() */
50: MatNorm(A,NORM_FROBENIUS,&s1norm);
51: MatNorm(B,NORM_FROBENIUS,&s2norm);
52: rnorm = PetscAbsScalar(s2norm-s1norm)/s2norm;
53: if (rnorm>tol) {
54: PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatNorm_FROBENIUS()- NormA=%16.14e NormB=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
55: }
56: MatNorm(A,NORM_INFINITY,&s1norm);
57: MatNorm(B,NORM_INFINITY,&s2norm);
58: rnorm = PetscAbsScalar(s2norm-s1norm)/s2norm;
59: if (rnorm>tol) {
60: PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatNorm_INFINITY()- NormA=%16.14e NormB=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
61: }
62: MatNorm(A,NORM_1,&s1norm);
63: MatNorm(B,NORM_1,&s2norm);
64: rnorm = PetscAbsScalar(s2norm-s1norm)/s2norm;
65: if (rnorm>tol) {
66: PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatNorm_NORM_1()- NormA=%16.14e NormB=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
67: }
69: /* Test MatMult() */
70: for (i=0; i<IMAX; i++) {
71: VecSetRandom(xx,rand);
72: MatMult(A,xx,s1);
73: MatMult(B,xx,s2);
74: VecAXPY(s2,-1.0,s1);
75: VecNorm(s2,NORM_2,&rnorm);
76: if (rnorm>tol) {
77: PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatMult - Norm2=%16.14e bs = %D\n",rank,rnorm,bs);
78: }
79: }
81: /* test MatMultAdd() */
82: for (i=0; i<IMAX; i++) {
83: VecSetRandom(xx,rand);
84: VecSetRandom(yy,rand);
85: MatMultAdd(A,xx,yy,s1);
86: MatMultAdd(B,xx,yy,s2);
87: VecAXPY(s2,-1.0,s1);
88: VecNorm(s2,NORM_2,&rnorm);
89: if (rnorm>tol) {
90: PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatMultAdd - Norm2=%16.14e bs = %D\n",rank,rnorm,bs);
91: }
92: }
94: /* Test MatMultTranspose() */
95: for (i=0; i<IMAX; i++) {
96: VecSetRandom(xx,rand);
97: MatMultTranspose(A,xx,s1);
98: MatMultTranspose(B,xx,s2);
99: VecNorm(s1,NORM_2,&s1norm);
100: VecNorm(s2,NORM_2,&s2norm);
101: rnorm = s2norm-s1norm;
102: if (rnorm<-tol || rnorm>tol) {
103: PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatMultTranspose - Norm1=%16.14e Norm2=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
104: }
105: }
106: /* Test MatMultTransposeAdd() */
107: for (i=0; i<IMAX; i++) {
108: VecSetRandom(xx,rand);
109: VecSetRandom(yy,rand);
110: MatMultTransposeAdd(A,xx,yy,s1);
111: MatMultTransposeAdd(B,xx,yy,s2);
112: VecNorm(s1,NORM_2,&s1norm);
113: VecNorm(s2,NORM_2,&s2norm);
114: rnorm = s2norm-s1norm;
115: if (rnorm<-tol || rnorm>tol) {
116: PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatMultTransposeAdd - Norm1=%16.14e Norm2=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
117: }
118: }
120: /* Check MatGetValues() */
121: MatGetOwnershipRange(A,&rstart,&rend);
122: MatGetSize(A,&M,&N);
125: for (i=0; i<IMAX; i++) {
126: /* Create random row numbers ad col numbers */
127: PetscRandomGetValue(rand,&v);
128: cols[0] = (int)(PetscRealPart(v)*N);
129: PetscRandomGetValue(rand,&v);
130: cols[1] = (int)(PetscRealPart(v)*N);
131: PetscRandomGetValue(rand,&v);
132: rows[0] = rstart + (int)(PetscRealPart(v)*m);
133: PetscRandomGetValue(rand,&v);
134: rows[1] = rstart + (int)(PetscRealPart(v)*m);
136: MatGetValues(A,2,rows,2,cols,vals1);
137: MatGetValues(B,2,rows,2,cols,vals2);
140: for (j=0; j<4; j++) {
141: if (vals1[j] != vals2[j]) {
142: PetscPrintf(PETSC_COMM_SELF,"[%d]: Error: MatGetValues rstart = %2d row = %2d col = %2d val1 = %e val2 = %e bs = %D\n",rank,rstart,rows[j/2],cols[j%2],PetscRealPart(vals1[j]),PetscRealPart(vals2[j]),bs);
143: }
144: }
145: }
147: /* Test MatGetRow()/ MatRestoreRow() */
148: for (ct=0; ct<100; ct++) {
149: PetscRandomGetValue(rand,&v);
150: row = rstart + (int)(PetscRealPart(v)*m);
151: MatGetRow(A,row,&ncols1,&cols1,&v1);
152: MatGetRow(B,row,&ncols2,&cols2,&v2);
154: for (i=0,j=0; i<ncols1 && j<ncols2; j++) {
155: while (cols2[j] != cols1[i]) i++;
156: if (v1[i] != v2[j]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"MatGetRow() failed - vals incorrect.");
157: }
158: if (j<ncols2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"MatGetRow() failed - cols incorrect");
160: MatRestoreRow(A,row,&ncols1,&cols1,&v1);
161: MatRestoreRow(B,row,&ncols2,&cols2,&v2);
162: }
164: /* Test MatConvert() */
165: MatConvert(A,MATSAME,MAT_INITIAL_MATRIX,&C);
167: /* See if MatMult Says both are same */
168: for (i=0; i<IMAX; i++) {
169: VecSetRandom(xx,rand);
170: MatMult(A,xx,s1);
171: MatMult(C,xx,s2);
172: VecNorm(s1,NORM_2,&s1norm);
173: VecNorm(s2,NORM_2,&s2norm);
174: rnorm = s2norm-s1norm;
175: if (rnorm<-tol || rnorm>tol) {
176: PetscPrintf(PETSC_COMM_SELF,"[%d] Error in MatConvert: MatMult - Norm1=%16.14e Norm2=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
177: }
178: }
179: MatDestroy(&C);
181: /* Test MatTranspose() */
182: MatTranspose(A,MAT_INITIAL_MATRIX,&At);
183: MatTranspose(B,MAT_INITIAL_MATRIX,&Bt);
184: for (i=0; i<IMAX; i++) {
185: VecSetRandom(xx,rand);
186: MatMult(At,xx,s1);
187: MatMult(Bt,xx,s2);
188: VecNorm(s1,NORM_2,&s1norm);
189: VecNorm(s2,NORM_2,&s2norm);
190: rnorm = s2norm-s1norm;
191: if (rnorm<-tol || rnorm>tol) {
192: PetscPrintf(PETSC_COMM_SELF,"[%d] Error in MatConvert:MatMult - Norm1=%16.14e Norm2=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
193: }
194: }
195: MatDestroy(&At);
196: MatDestroy(&Bt);
198: MatDestroy(&A);
199: MatDestroy(&B);
200: VecDestroy(&xx);
201: VecDestroy(&yy);
202: VecDestroy(&s1);
203: VecDestroy(&s2);
204: PetscRandomDestroy(&rand);
205: PetscFinalize();
206: return ierr;
207: }
210: /*TEST
212: build:
213: requires: !complex
215: test:
216: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
217: nsize: 3
218: args: -matload_block_size 1 -f ${DATAFILESPATH}/matrices/small
220: test:
221: suffix: 2
222: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
223: nsize: 3
224: args: -matload_block_size 2 -f ${DATAFILESPATH}/matrices/small
225: output_file: output/ex53_1.out
227: test:
228: suffix: 3
229: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
230: nsize: 3
231: args: -matload_block_size 4 -f ${DATAFILESPATH}/matrices/small
232: output_file: output/ex53_1.out
234: test:
235: TODO: Matrix row/column sizes are not compatible with block size
236: suffix: 4
237: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
238: nsize: 3
239: args: -matload_block_size 5 -f ${DATAFILESPATH}/matrices/small
240: output_file: output/ex53_1.out
242: test:
243: suffix: 5
244: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
245: nsize: 3
246: args: -matload_block_size 6 -f ${DATAFILESPATH}/matrices/small
247: output_file: output/ex53_1.out
249: test:
250: TODO: Matrix row/column sizes are not compatible with block size
251: suffix: 6
252: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
253: nsize: 3
254: args: -matload_block_size 7 -f ${DATAFILESPATH}/matrices/small
255: output_file: output/ex53_1.out
257: test:
258: TODO: Matrix row/column sizes are not compatible with block size
259: suffix: 7
260: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
261: nsize: 3
262: args: -matload_block_size 8 -f ${DATAFILESPATH}/matrices/small
263: output_file: output/ex53_1.out
265: test:
266: suffix: 8
267: requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES)
268: nsize: 4
269: args: -matload_block_size 3 -f ${DATAFILESPATH}/matrices/small
270: output_file: output/ex53_1.out
272: TEST*/