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