Actual source code: ex53.c
2: static char help[] = "Tests the vatious routines in MatMPIBAIJ format.\n";
5: #include petscmat.h
6: #define IMAX 15
9: int main(int argc,char **args)
10: {
11: Mat A,B,C,At,Bt;
12: PetscViewer fd;
13: char file[PETSC_MAX_PATH_LEN];
14: PetscRandom rand;
15: Vec xx,yy,s1,s2;
16: PetscReal s1norm,s2norm,rnorm,tol=1.e-10;
17: PetscInt rstart,rend,rows[2],cols[2],m,n,i,j,M,N,ct,row,ncols1,ncols2,bs;
18: PetscMPIInt rank,size;
19: PetscErrorCode ierr;
20: const PetscInt *cols1,*cols2;
21: PetscScalar vals1[4],vals2[4],v;
22: const PetscScalar *v1,*v2;
23: PetscTruth flg;
25: PetscInitialize(&argc,&args,(char *)0,help);
26: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
27: MPI_Comm_size(PETSC_COMM_WORLD,&size);
29: #if defined(PETSC_USE_COMPLEX)
30: SETERRQ(1,"This example does not work with complex numbers");
31: #else
33: /* Check out if MatLoad() works */
34: PetscOptionsGetString(PETSC_NULL,"-f",file,PETSC_MAX_PATH_LEN-1,&flg);
35: if (!flg) SETERRQ(1,"Input file not specified");
36: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
37: MatLoad(fd,MATBAIJ,&A);
38: /*
39: if (size == 1){
40: MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&B);
41: } else {
42: MatConvert(A,MATMPIAIJ,MAT_INITIAL_MATRIX,&B);
43: }
44: */
45: MatConvert(A,MATAIJ,MAT_INITIAL_MATRIX,&B);
46: PetscViewerDestroy(fd);
47:
48: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
49: PetscRandomSetFromOptions(rand);
50: MatGetLocalSize(A,&m,&n);
51: VecCreate(PETSC_COMM_WORLD,&xx);
52: VecSetSizes(xx,m,PETSC_DECIDE);
53: VecSetFromOptions(xx);
54: VecDuplicate(xx,&s1);
55: VecDuplicate(xx,&s2);
56: VecDuplicate(xx,&yy);
57: MatGetBlockSize(A,&bs);
59: /* Test MatNorm() */
60: MatNorm(A,NORM_FROBENIUS,&s1norm);
61: MatNorm(B,NORM_FROBENIUS,&s2norm);
62: rnorm = PetscAbsScalar(s2norm-s1norm)/s2norm;
63: if ( rnorm>tol ) {
64: PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatNorm_FROBENIUS()- NormA=%16.14e NormB=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
65: }
66: MatNorm(A,NORM_INFINITY,&s1norm);
67: MatNorm(B,NORM_INFINITY,&s2norm);
68: rnorm = PetscAbsScalar(s2norm-s1norm)/s2norm;
69: if ( rnorm>tol ) {
70: PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatNorm_INFINITY()- NormA=%16.14e NormB=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
71: }
72: MatNorm(A,NORM_1,&s1norm);
73: MatNorm(B,NORM_1,&s2norm);
74: rnorm = PetscAbsScalar(s2norm-s1norm)/s2norm;
75: if ( rnorm>tol ) {
76: PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatNorm_NORM_1()- NormA=%16.14e NormB=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
77: }
79: /* Test MatMult() */
80: for (i=0; i<IMAX; i++) {
81: VecSetRandom(xx,rand);
82: MatMult(A,xx,s1);
83: MatMult(B,xx,s2);
84: VecAXPY(s2,-1.0,s1);
85: VecNorm(s2,NORM_2,&rnorm);
86: if (rnorm>tol) {
87: PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatMult - Norm2=%16.14e bs = %D\n",rank,rnorm,bs);
88: }
89: }
91: /* test MatMultAdd() */
92: for (i=0; i<IMAX; i++) {
93: VecSetRandom(xx,rand);
94: VecSetRandom(yy,rand);
95: MatMultAdd(A,xx,yy,s1);
96: MatMultAdd(B,xx,yy,s2);
97: VecAXPY(s2,-1.0,s1);
98: VecNorm(s2,NORM_2,&rnorm);
99: if (rnorm>tol) {
100: PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatMultAdd - Norm2=%16.14e bs = %D\n",rank,rnorm,bs);
101: }
102: }
104: /* Test MatMultTranspose() */
105: for (i=0; i<IMAX; i++) {
106: VecSetRandom(xx,rand);
107: MatMultTranspose(A,xx,s1);
108: MatMultTranspose(B,xx,s2);
109: VecNorm(s1,NORM_2,&s1norm);
110: VecNorm(s2,NORM_2,&s2norm);
111: rnorm = s2norm-s1norm;
112: if (rnorm<-tol || rnorm>tol) {
113: PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatMultTranspose - Norm1=%16.14e Norm2=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
114: }
115: }
116: /* Test MatMultTransposeAdd() */
117: for (i=0; i<IMAX; i++) {
118: VecSetRandom(xx,rand);
119: VecSetRandom(yy,rand);
120: MatMultTransposeAdd(A,xx,yy,s1);
121: MatMultTransposeAdd(B,xx,yy,s2);
122: VecNorm(s1,NORM_2,&s1norm);
123: VecNorm(s2,NORM_2,&s2norm);
124: rnorm = s2norm-s1norm;
125: if (rnorm<-tol || rnorm>tol) {
126: PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatMultTransposeAdd - Norm1=%16.14e Norm2=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
127: }
128: }
130: /* Check MatGetValues() */
131: MatGetOwnershipRange(A,&rstart,&rend);
132: MatGetSize(A,&M,&N);
135: for (i=0; i<IMAX; i++) {
136: /* Create random row numbers ad col numbers */
137: PetscRandomGetValue(rand,&v);
138: cols[0] = (int)(PetscRealPart(v)*N);
139: PetscRandomGetValue(rand,&v);
140: cols[1] = (int)(PetscRealPart(v)*N);
141: PetscRandomGetValue(rand,&v);
142: rows[0] = rstart + (int)(PetscRealPart(v)*m);
143: PetscRandomGetValue(rand,&v);
144: rows[1] = rstart + (int)(PetscRealPart(v)*m);
145:
146: MatGetValues(A,2,rows,2,cols,vals1);
147: MatGetValues(B,2,rows,2,cols,vals2);
150: for (j=0; j<4; j++) {
151: if(vals1[j] != vals2[j])
152: 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);
153: }
154: }
156: /* Test MatGetRow()/ MatRestoreRow() */
157: for (ct=0; ct<100; ct++) {
158: PetscRandomGetValue(rand,&v);
159: row = rstart + (int)(PetscRealPart(v)*m);
160: MatGetRow(A,row,&ncols1,&cols1,&v1);
161: MatGetRow(B,row,&ncols2,&cols2,&v2);
162:
163: for (i=0,j=0; i<ncols1 && j<ncols2; j++) {
164: while (cols2[j] != cols1[i]) i++;
165: if (v1[i] != v2[j]) SETERRQ(1,"MatGetRow() failed - vals incorrect.");
166: }
167: if (j<ncols2) SETERRQ(1,"MatGetRow() failed - cols incorrect");
168:
169: MatRestoreRow(A,row,&ncols1,&cols1,&v1);
170: MatRestoreRow(B,row,&ncols2,&cols2,&v2);
171: }
172:
173: /* Test MatConvert() */
174: MatConvert(A,MATSAME,MAT_INITIAL_MATRIX,&C);
175:
176: /* See if MatMult Says both are same */
177: for (i=0; i<IMAX; i++) {
178: VecSetRandom(xx,rand);
179: MatMult(A,xx,s1);
180: MatMult(C,xx,s2);
181: VecNorm(s1,NORM_2,&s1norm);
182: VecNorm(s2,NORM_2,&s2norm);
183: rnorm = s2norm-s1norm;
184: if (rnorm<-tol || rnorm>tol) {
185: PetscPrintf(PETSC_COMM_SELF,"[%d] Error in MatConvert: MatMult - Norm1=%16.14e Norm2=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
186: }
187: }
188: MatDestroy(C);
190: /* Test MatTranspose() */
191: MatTranspose(A,&At);
192: MatTranspose(B,&Bt);
193: for (i=0; i<IMAX; i++) {
194: VecSetRandom(xx,rand);
195: MatMult(At,xx,s1);
196: MatMult(Bt,xx,s2);
197: VecNorm(s1,NORM_2,&s1norm);
198: VecNorm(s2,NORM_2,&s2norm);
199: rnorm = s2norm-s1norm;
200: if (rnorm<-tol || rnorm>tol) {
201: PetscPrintf(PETSC_COMM_SELF,"[%d] Error in MatConvert:MatMult - Norm1=%16.14e Norm2=%16.14e bs = %D\n",
202: rank,s1norm,s2norm,bs);
203: }
204: }
205: MatDestroy(At);
206: MatDestroy(Bt);
208: MatDestroy(A);
209: MatDestroy(B);
210: VecDestroy(xx);
211: VecDestroy(yy);
212: VecDestroy(s1);
213: VecDestroy(s2);
214: PetscRandomDestroy(rand);
215: PetscFinalize();
216: #endif
217: return 0;
218: }