Actual source code: ex53.c
1: /*$Id: ex53.c,v 1.25 2001/09/11 16:32:50 bsmith Exp $*/
3: static char help[] = "Tests the vatious routines in MatMPIBAIJ format.\n";
6: #include petscmat.h
7: #define IMAX 15
10: int main(int argc,char **args)
11: {
12: Mat A,B,C,At,Bt;
13: PetscViewer fd;
14: char file[128];
15: PetscRandom rand;
16: Vec xx,yy,s1,s2;
17: PetscReal s1norm,s2norm,rnorm,tol = 1.e-10;
18: int rstart,rend,rows[2],cols[2],m,n,i,j,ierr,M,N,rank,ct,row,ncols1;
19: int *cols1,ncols2,*cols2,bs;
20: PetscScalar vals1[4],vals2[4],v,*v1,*v2;
21: PetscTruth flg;
23: PetscInitialize(&argc,&args,(char *)0,help);
24: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
26: #if defined(PETSC_USE_COMPLEX)
27: SETERRQ(1,"This example does not work with complex numbers");
28: #else
30: /* Check out if MatLoad() works */
31: PetscOptionsGetString(PETSC_NULL,"-f",file,127,&flg);
32: if (!flg) SETERRQ(1,"Input file not specified");
33: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,PETSC_BINARY_RDONLY,&fd);
34: MatLoad(fd,MATMPIBAIJ,&A);
35: PetscViewerDestroy(fd);
37: MatConvert(A,MATMPIAIJ,&B);
38:
39: PetscRandomCreate(PETSC_COMM_WORLD,RANDOM_DEFAULT,&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);
48: MatGetBlockSize(A,&bs);
49: /* Test MatMult() */
50: for (i=0; i<IMAX; i++) {
51: VecSetRandom(rand,xx);
52: MatMult(A,xx,s1);
53: MatMult(B,xx,s2);
54: VecNorm(s1,NORM_2,&s1norm);
55: VecNorm(s2,NORM_2,&s2norm);
56: rnorm = s2norm-s1norm;
57: if (rnorm<-tol || rnorm>tol) {
58: PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e bs = %d\n",
59: s1norm,s2norm,bs);
60: }
61: }
62: /* test MatMultAdd() */
63: for (i=0; i<IMAX; i++) {
64: VecSetRandom(rand,xx);
65: VecSetRandom(rand,yy);
66: MatMultAdd(A,xx,yy,s1);
67: MatMultAdd(B,xx,yy,s2);
68: VecNorm(s1,NORM_2,&s1norm);
69: VecNorm(s2,NORM_2,&s2norm);
70: rnorm = s2norm-s1norm;
71: if (rnorm<-tol || rnorm>tol) {
72: PetscPrintf(PETSC_COMM_SELF,"Error:MatMultAdd - Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
73: }
74: }
75: /* Test MatMultTranspose() */
76: for (i=0; i<IMAX; i++) {
77: VecSetRandom(rand,xx);
78: MatMultTranspose(A,xx,s1);
79: MatMultTranspose(B,xx,s2);
80: VecNorm(s1,NORM_2,&s1norm);
81: VecNorm(s2,NORM_2,&s2norm);
82: rnorm = s2norm-s1norm;
83: if (rnorm<-tol || rnorm>tol) {
84: PetscPrintf(PETSC_COMM_SELF,"Error:MatMultTranspose - Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
85: }
86: }
87: /* Test MatMultTransposeAdd() */
88: for (i=0; i<IMAX; i++) {
89: VecSetRandom(rand,xx);
90: VecSetRandom(rand,yy);
91: MatMultTransposeAdd(A,xx,yy,s1);
92: MatMultTransposeAdd(B,xx,yy,s2);
93: VecNorm(s1,NORM_2,&s1norm);
94: VecNorm(s2,NORM_2,&s2norm);
95: rnorm = s2norm-s1norm;
96: if (rnorm<-tol || rnorm>tol) {
97: PetscPrintf(PETSC_COMM_SELF,"Error:MatMultTransposeAdd - Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
98: }
99: }
101: /* Check MatGetValues() */
102: MatGetOwnershipRange(A,&rstart,&rend);
103: MatGetSize(A,&M,&N);
106: for (i=0; i<IMAX; i++) {
107: /* Create random row numbers ad col numbers */
108: PetscRandomGetValue(rand,&v);
109: cols[0] = (int)(PetscRealPart(v)*N);
110: PetscRandomGetValue(rand,&v);
111: cols[1] = (int)(PetscRealPart(v)*N);
112: PetscRandomGetValue(rand,&v);
113: rows[0] = rstart + (int)(PetscRealPart(v)*m);
114: PetscRandomGetValue(rand,&v);
115: rows[1] = rstart + (int)(PetscRealPart(v)*m);
116:
117: MatGetValues(A,2,rows,2,cols,vals1);
118: MatGetValues(B,2,rows,2,cols,vals2);
121: for (j=0; j<4; j++) {
122: if(vals1[j] != vals2[j])
123: 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);
124: }
125: }
127: /* Test MatGetRow()/ MatRestoreRow() */
128: for (ct=0; ct<100; ct++) {
129: PetscRandomGetValue(rand,&v);
130: row = rstart + (int)(PetscRealPart(v)*m);
131: MatGetRow(A,row,&ncols1,&cols1,&v1);
132: MatGetRow(B,row,&ncols2,&cols2,&v2);
133:
134: for (i=0,j=0; i<ncols1 && j<ncols2; j++) {
135: while (cols2[j] != cols1[i]) i++;
136: if (v1[i] != v2[j]) SETERRQ(1,"MatGetRow() failed - vals incorrect.");
137: }
138: if (j<ncols2) SETERRQ(1,"MatGetRow() failed - cols incorrect");
139:
140: MatRestoreRow(A,row,&ncols1,&cols1,&v1);
141: MatRestoreRow(B,row,&ncols2,&cols2,&v2);
142: }
143:
144: /* Test MatConvert() */
145: MatConvert(A,MATSAME,&C);
146:
147: /* See if MatMult Says both are same */
148: for (i=0; i<IMAX; i++) {
149: VecSetRandom(rand,xx);
150: MatMult(A,xx,s1);
151: MatMult(C,xx,s2);
152: VecNorm(s1,NORM_2,&s1norm);
153: VecNorm(s2,NORM_2,&s2norm);
154: rnorm = s2norm-s1norm;
155: if (rnorm<-tol || rnorm>tol) {
156: PetscPrintf(PETSC_COMM_SELF,"Error in MatConvert:MatMult - Norm1=%16.14e Norm2=%16.14e bs = %d\n",
157: s1norm,s2norm,bs);
158: }
159: }
160: MatDestroy(C);
162: /* Test MatTranspose() */
163: MatTranspose(A,&At);
164: MatTranspose(B,&Bt);
165: for (i=0; i<IMAX; i++) {
166: VecSetRandom(rand,xx);
167: MatMult(At,xx,s1);
168: MatMult(Bt,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,"Error in MatConvert:MatMult - Norm1=%16.14e Norm2=%16.14e bs = %d\n",
174: s1norm,s2norm,bs);
175: }
176: }
177: MatDestroy(At);
178: MatDestroy(Bt);
180: MatDestroy(A);
181: MatDestroy(B);
182: VecDestroy(xx);
183: VecDestroy(yy);
184: VecDestroy(s1);
185: VecDestroy(s2);
186: PetscRandomDestroy(rand);
187: PetscFinalize();
188: #endif
189: return 0;
190: }