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