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: PetscErrorCode 0;
17: const PetscInt *cols1,*cols2;
18: PetscScalar vals1[4],vals2[4],v;
19: const PetscScalar *v1,*v2;
20: PetscBool flg;
22: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
23: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
24: MPI_Comm_size(PETSC_COMM_WORLD,&size);
26: /* Check out if MatLoad() works */
27: PetscOptionsGetString(NULL,NULL,"-f",file,sizeof(file),&flg);
28: if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Input file not specified");
29: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
30: MatCreate(PETSC_COMM_WORLD,&A);
31: MatSetType(A,MATBAIJ);
32: MatLoad(A,fd);
34: MatConvert(A,MATAIJ,MAT_INITIAL_MATRIX,&B);
35: PetscViewerDestroy(&fd);
37: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
38: PetscRandomSetFromOptions(rand);
39: MatGetLocalSize(A,&m,&n);
40: VecCreate(PETSC_COMM_WORLD,&xx);
41: VecSetSizes(xx,m,PETSC_DECIDE);
42: VecSetFromOptions(xx);
43: VecDuplicate(xx,&s1);
44: VecDuplicate(xx,&s2);
45: VecDuplicate(xx,&yy);
46: MatGetBlockSize(A,&bs);
48: /* Test MatNorm() */
49: MatNorm(A,NORM_FROBENIUS,&s1norm);
50: MatNorm(B,NORM_FROBENIUS,&s2norm);
51: rnorm = PetscAbsScalar(s2norm-s1norm)/s2norm;
52: if (rnorm>tol) {
53: PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatNorm_FROBENIUS()- NormA=%16.14e NormB=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
54: }
55: MatNorm(A,NORM_INFINITY,&s1norm);
56: MatNorm(B,NORM_INFINITY,&s2norm);
57: rnorm = PetscAbsScalar(s2norm-s1norm)/s2norm;
58: if (rnorm>tol) {
59: PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatNorm_INFINITY()- NormA=%16.14e NormB=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
60: }
61: MatNorm(A,NORM_1,&s1norm);
62: MatNorm(B,NORM_1,&s2norm);
63: rnorm = PetscAbsScalar(s2norm-s1norm)/s2norm;
64: if (rnorm>tol) {
65: PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatNorm_NORM_1()- NormA=%16.14e NormB=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
66: }
68: /* Test MatMult() */
69: for (i=0; i<IMAX; i++) {
70: VecSetRandom(xx,rand);
71: MatMult(A,xx,s1);
72: MatMult(B,xx,s2);
73: VecAXPY(s2,-1.0,s1);
74: VecNorm(s2,NORM_2,&rnorm);
75: if (rnorm>tol) {
76: PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatMult - Norm2=%16.14e bs = %D\n",rank,rnorm,bs);
77: }
78: }
80: /* test MatMultAdd() */
81: for (i=0; i<IMAX; i++) {
82: VecSetRandom(xx,rand);
83: VecSetRandom(yy,rand);
84: MatMultAdd(A,xx,yy,s1);
85: MatMultAdd(B,xx,yy,s2);
86: VecAXPY(s2,-1.0,s1);
87: VecNorm(s2,NORM_2,&rnorm);
88: if (rnorm>tol) {
89: PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatMultAdd - Norm2=%16.14e bs = %D\n",rank,rnorm,bs);
90: }
91: }
93: /* Test MatMultTranspose() */
94: for (i=0; i<IMAX; i++) {
95: VecSetRandom(xx,rand);
96: MatMultTranspose(A,xx,s1);
97: MatMultTranspose(B,xx,s2);
98: VecNorm(s1,NORM_2,&s1norm);
99: VecNorm(s2,NORM_2,&s2norm);
100: rnorm = s2norm-s1norm;
101: if (rnorm<-tol || rnorm>tol) {
102: PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatMultTranspose - Norm1=%16.14e Norm2=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
103: }
104: }
105: /* Test MatMultTransposeAdd() */
106: for (i=0; i<IMAX; i++) {
107: VecSetRandom(xx,rand);
108: VecSetRandom(yy,rand);
109: MatMultTransposeAdd(A,xx,yy,s1);
110: MatMultTransposeAdd(B,xx,yy,s2);
111: VecNorm(s1,NORM_2,&s1norm);
112: VecNorm(s2,NORM_2,&s2norm);
113: rnorm = s2norm-s1norm;
114: if (rnorm<-tol || rnorm>tol) {
115: PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatMultTransposeAdd - Norm1=%16.14e Norm2=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
116: }
117: }
119: /* Check MatGetValues() */
120: MatGetOwnershipRange(A,&rstart,&rend);
121: MatGetSize(A,&M,&N);
123: for (i=0; i<IMAX; i++) {
124: /* Create random row numbers ad col numbers */
125: PetscRandomGetValue(rand,&v);
126: cols[0] = (int)(PetscRealPart(v)*N);
127: PetscRandomGetValue(rand,&v);
128: cols[1] = (int)(PetscRealPart(v)*N);
129: PetscRandomGetValue(rand,&v);
130: rows[0] = rstart + (int)(PetscRealPart(v)*m);
131: PetscRandomGetValue(rand,&v);
132: rows[1] = rstart + (int)(PetscRealPart(v)*m);
134: MatGetValues(A,2,rows,2,cols,vals1);
135: MatGetValues(B,2,rows,2,cols,vals2);
137: for (j=0; j<4; j++) {
138: if (vals1[j] != vals2[j]) {
139: 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);
140: }
141: }
142: }
144: /* Test MatGetRow()/ MatRestoreRow() */
145: for (ct=0; ct<100; ct++) {
146: PetscRandomGetValue(rand,&v);
147: row = rstart + (int)(PetscRealPart(v)*m);
148: MatGetRow(A,row,&ncols1,&cols1,&v1);
149: MatGetRow(B,row,&ncols2,&cols2,&v2);
151: for (i=0,j=0; i<ncols1 && j<ncols2; j++) {
152: while (cols2[j] != cols1[i]) i++;
153: if (v1[i] != v2[j]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"MatGetRow() failed - vals incorrect.");
154: }
155: if (j<ncols2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"MatGetRow() failed - cols incorrect");
157: MatRestoreRow(A,row,&ncols1,&cols1,&v1);
158: MatRestoreRow(B,row,&ncols2,&cols2,&v2);
159: }
161: /* Test MatConvert() */
162: MatConvert(A,MATSAME,MAT_INITIAL_MATRIX,&C);
164: /* See if MatMult Says both are same */
165: for (i=0; i<IMAX; i++) {
166: VecSetRandom(xx,rand);
167: MatMult(A,xx,s1);
168: MatMult(C,xx,s2);
169: VecNorm(s1,NORM_2,&s1norm);
170: VecNorm(s2,NORM_2,&s2norm);
171: rnorm = s2norm-s1norm;
172: if (rnorm<-tol || rnorm>tol) {
173: PetscPrintf(PETSC_COMM_SELF,"[%d] Error in MatConvert: MatMult - Norm1=%16.14e Norm2=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
174: }
175: }
176: MatDestroy(&C);
178: /* Test MatTranspose() */
179: MatTranspose(A,MAT_INITIAL_MATRIX,&At);
180: MatTranspose(B,MAT_INITIAL_MATRIX,&Bt);
181: for (i=0; i<IMAX; i++) {
182: VecSetRandom(xx,rand);
183: MatMult(At,xx,s1);
184: MatMult(Bt,xx,s2);
185: VecNorm(s1,NORM_2,&s1norm);
186: VecNorm(s2,NORM_2,&s2norm);
187: rnorm = s2norm-s1norm;
188: if (rnorm<-tol || rnorm>tol) {
189: PetscPrintf(PETSC_COMM_SELF,"[%d] Error in MatConvert:MatMult - Norm1=%16.14e Norm2=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
190: }
191: }
192: MatDestroy(&At);
193: MatDestroy(&Bt);
195: MatDestroy(&A);
196: MatDestroy(&B);
197: VecDestroy(&xx);
198: VecDestroy(&yy);
199: VecDestroy(&s1);
200: VecDestroy(&s2);
201: PetscRandomDestroy(&rand);
202: PetscFinalize();
203: return ierr;
204: }
206: /*TEST
208: build:
209: requires: !complex
211: test:
212: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
213: nsize: 3
214: args: -matload_block_size 1 -f ${DATAFILESPATH}/matrices/small
216: test:
217: suffix: 2
218: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
219: nsize: 3
220: args: -matload_block_size 2 -f ${DATAFILESPATH}/matrices/small
221: output_file: output/ex53_1.out
223: test:
224: suffix: 3
225: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
226: nsize: 3
227: args: -matload_block_size 4 -f ${DATAFILESPATH}/matrices/small
228: output_file: output/ex53_1.out
230: test:
231: TODO: Matrix row/column sizes are not compatible with block size
232: suffix: 4
233: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
234: nsize: 3
235: args: -matload_block_size 5 -f ${DATAFILESPATH}/matrices/small
236: output_file: output/ex53_1.out
238: test:
239: suffix: 5
240: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
241: nsize: 3
242: args: -matload_block_size 6 -f ${DATAFILESPATH}/matrices/small
243: output_file: output/ex53_1.out
245: test:
246: TODO: Matrix row/column sizes are not compatible with block size
247: suffix: 6
248: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
249: nsize: 3
250: args: -matload_block_size 7 -f ${DATAFILESPATH}/matrices/small
251: output_file: output/ex53_1.out
253: test:
254: TODO: Matrix row/column sizes are not compatible with block size
255: suffix: 7
256: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
257: nsize: 3
258: args: -matload_block_size 8 -f ${DATAFILESPATH}/matrices/small
259: output_file: output/ex53_1.out
261: test:
262: suffix: 8
263: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
264: nsize: 4
265: args: -matload_block_size 3 -f ${DATAFILESPATH}/matrices/small
266: output_file: output/ex53_1.out
268: TEST*/