Actual source code: ex53.c
petsc-3.6.4 2016-04-12
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: PetscBool 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(PETSC_COMM_WORLD,1,"This example does not work with complex numbers");
31: #else
33: /* Check out if MatLoad() works */
34: PetscOptionsGetString(NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg);
35: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Input file not specified");
36: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
37: MatCreate(PETSC_COMM_WORLD,&A);
38: MatSetType(A,MATBAIJ);
39: MatLoad(A,fd);
41: MatConvert(A,MATAIJ,MAT_INITIAL_MATRIX,&B);
42: PetscViewerDestroy(&fd);
44: PetscRandomCreate(PETSC_COMM_WORLD,&rand);
45: PetscRandomSetFromOptions(rand);
46: MatGetLocalSize(A,&m,&n);
47: VecCreate(PETSC_COMM_WORLD,&xx);
48: VecSetSizes(xx,m,PETSC_DECIDE);
49: VecSetFromOptions(xx);
50: VecDuplicate(xx,&s1);
51: VecDuplicate(xx,&s2);
52: VecDuplicate(xx,&yy);
53: MatGetBlockSize(A,&bs);
55: /* Test MatNorm() */
56: MatNorm(A,NORM_FROBENIUS,&s1norm);
57: MatNorm(B,NORM_FROBENIUS,&s2norm);
58: rnorm = PetscAbsScalar(s2norm-s1norm)/s2norm;
59: if (rnorm>tol) {
60: PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatNorm_FROBENIUS()- NormA=%16.14e NormB=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
61: }
62: MatNorm(A,NORM_INFINITY,&s1norm);
63: MatNorm(B,NORM_INFINITY,&s2norm);
64: rnorm = PetscAbsScalar(s2norm-s1norm)/s2norm;
65: if (rnorm>tol) {
66: PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatNorm_INFINITY()- NormA=%16.14e NormB=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
67: }
68: MatNorm(A,NORM_1,&s1norm);
69: MatNorm(B,NORM_1,&s2norm);
70: rnorm = PetscAbsScalar(s2norm-s1norm)/s2norm;
71: if (rnorm>tol) {
72: PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatNorm_NORM_1()- NormA=%16.14e NormB=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
73: }
75: /* Test MatMult() */
76: for (i=0; i<IMAX; i++) {
77: VecSetRandom(xx,rand);
78: MatMult(A,xx,s1);
79: MatMult(B,xx,s2);
80: VecAXPY(s2,-1.0,s1);
81: VecNorm(s2,NORM_2,&rnorm);
82: if (rnorm>tol) {
83: PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatMult - Norm2=%16.14e bs = %D\n",rank,rnorm,bs);
84: }
85: }
87: /* test MatMultAdd() */
88: for (i=0; i<IMAX; i++) {
89: VecSetRandom(xx,rand);
90: VecSetRandom(yy,rand);
91: MatMultAdd(A,xx,yy,s1);
92: MatMultAdd(B,xx,yy,s2);
93: VecAXPY(s2,-1.0,s1);
94: VecNorm(s2,NORM_2,&rnorm);
95: if (rnorm>tol) {
96: PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatMultAdd - Norm2=%16.14e bs = %D\n",rank,rnorm,bs);
97: }
98: }
100: /* Test MatMultTranspose() */
101: for (i=0; i<IMAX; i++) {
102: VecSetRandom(xx,rand);
103: MatMultTranspose(A,xx,s1);
104: MatMultTranspose(B,xx,s2);
105: VecNorm(s1,NORM_2,&s1norm);
106: VecNorm(s2,NORM_2,&s2norm);
107: rnorm = s2norm-s1norm;
108: if (rnorm<-tol || rnorm>tol) {
109: PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatMultTranspose - Norm1=%16.14e Norm2=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
110: }
111: }
112: /* Test MatMultTransposeAdd() */
113: for (i=0; i<IMAX; i++) {
114: VecSetRandom(xx,rand);
115: VecSetRandom(yy,rand);
116: MatMultTransposeAdd(A,xx,yy,s1);
117: MatMultTransposeAdd(B,xx,yy,s2);
118: VecNorm(s1,NORM_2,&s1norm);
119: VecNorm(s2,NORM_2,&s2norm);
120: rnorm = s2norm-s1norm;
121: if (rnorm<-tol || rnorm>tol) {
122: PetscPrintf(PETSC_COMM_SELF,"[%d] Error: MatMultTransposeAdd - Norm1=%16.14e Norm2=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
123: }
124: }
126: /* Check MatGetValues() */
127: MatGetOwnershipRange(A,&rstart,&rend);
128: MatGetSize(A,&M,&N);
131: for (i=0; i<IMAX; i++) {
132: /* Create random row numbers ad col numbers */
133: PetscRandomGetValue(rand,&v);
134: cols[0] = (int)(PetscRealPart(v)*N);
135: PetscRandomGetValue(rand,&v);
136: cols[1] = (int)(PetscRealPart(v)*N);
137: PetscRandomGetValue(rand,&v);
138: rows[0] = rstart + (int)(PetscRealPart(v)*m);
139: PetscRandomGetValue(rand,&v);
140: rows[1] = rstart + (int)(PetscRealPart(v)*m);
142: MatGetValues(A,2,rows,2,cols,vals1);
143: MatGetValues(B,2,rows,2,cols,vals2);
146: for (j=0; j<4; j++) {
147: if (vals1[j] != vals2[j]) {
148: 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);
149: }
150: }
151: }
153: /* Test MatGetRow()/ MatRestoreRow() */
154: for (ct=0; ct<100; ct++) {
155: PetscRandomGetValue(rand,&v);
156: row = rstart + (int)(PetscRealPart(v)*m);
157: MatGetRow(A,row,&ncols1,&cols1,&v1);
158: MatGetRow(B,row,&ncols2,&cols2,&v2);
160: for (i=0,j=0; i<ncols1 && j<ncols2; j++) {
161: while (cols2[j] != cols1[i]) i++;
162: if (v1[i] != v2[j]) SETERRQ(PETSC_COMM_SELF,1,"MatGetRow() failed - vals incorrect.");
163: }
164: if (j<ncols2) SETERRQ(PETSC_COMM_SELF,1,"MatGetRow() failed - cols incorrect");
166: MatRestoreRow(A,row,&ncols1,&cols1,&v1);
167: MatRestoreRow(B,row,&ncols2,&cols2,&v2);
168: }
170: /* Test MatConvert() */
171: MatConvert(A,MATSAME,MAT_INITIAL_MATRIX,&C);
173: /* See if MatMult Says both are same */
174: for (i=0; i<IMAX; i++) {
175: VecSetRandom(xx,rand);
176: MatMult(A,xx,s1);
177: MatMult(C,xx,s2);
178: VecNorm(s1,NORM_2,&s1norm);
179: VecNorm(s2,NORM_2,&s2norm);
180: rnorm = s2norm-s1norm;
181: if (rnorm<-tol || rnorm>tol) {
182: PetscPrintf(PETSC_COMM_SELF,"[%d] Error in MatConvert: MatMult - Norm1=%16.14e Norm2=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
183: }
184: }
185: MatDestroy(&C);
187: /* Test MatTranspose() */
188: MatTranspose(A,MAT_INITIAL_MATRIX,&At);
189: MatTranspose(B,MAT_INITIAL_MATRIX,&Bt);
190: for (i=0; i<IMAX; i++) {
191: VecSetRandom(xx,rand);
192: MatMult(At,xx,s1);
193: MatMult(Bt,xx,s2);
194: VecNorm(s1,NORM_2,&s1norm);
195: VecNorm(s2,NORM_2,&s2norm);
196: rnorm = s2norm-s1norm;
197: if (rnorm<-tol || rnorm>tol) {
198: PetscPrintf(PETSC_COMM_SELF,"[%d] Error in MatConvert:MatMult - Norm1=%16.14e Norm2=%16.14e bs = %D\n",rank,s1norm,s2norm,bs);
199: }
200: }
201: MatDestroy(&At);
202: MatDestroy(&Bt);
204: MatDestroy(&A);
205: MatDestroy(&B);
206: VecDestroy(&xx);
207: VecDestroy(&yy);
208: VecDestroy(&s1);
209: VecDestroy(&s2);
210: PetscRandomDestroy(&rand);
211: PetscFinalize();
212: #endif
213: return 0;
214: }