Actual source code: ex47.c
2: static char help[] = "Tests the various routines in MatBAIJ format.\n\
3: Input arguments are:\n\
4: -f <input_file> : file to load. For a 5X5 example of the 5-pt. stencil,\n\
5: use the file petsc/src/mat/examples/matbinary.ex\n\n";
7: #include <petscmat.h>
11: int main(int argc,char **args)
12: {
13: Mat A,B,C;
14: PetscViewer va,vb,vc;
15: Vec x,y;
16: PetscErrorCode ierr;
17: PetscInt i,j,row,m,n,ncols1,ncols2,ct,m2,n2;
18: const PetscInt *cols1,*cols2;
19: char file[PETSC_MAX_PATH_LEN];
20: PetscBool tflg;
21: PetscScalar rval;
22: const PetscScalar *vals1,*vals2;
23: PetscReal norm1,norm2,rnorm;
24: PetscRandom r;
27: PetscInitialize(&argc,&args,(char *)0,help);
28: #if defined(PETSC_USE_COMPLEX)
29: SETERRQ(PETSC_COMM_WORLD,1,"This example does not work with complex numbers");
30: #else
31:
32: PetscOptionsGetString(PETSC_NULL,"-f",file,PETSC_MAX_PATH_LEN,PETSC_NULL);
34: /* Load the matrix as AIJ format */
35: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&va);
36: MatCreate(PETSC_COMM_WORLD,&A);
37: MatSetType(A,MATSEQAIJ);
38: MatLoad(A,va);
39: PetscViewerDestroy(&va);
41: /* Load the matrix as BAIJ format */
42: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&vb);
43: MatCreate(PETSC_COMM_WORLD,&B);
44: MatSetType(B,MATSEQBAIJ);
45: MatLoad(B,vb);
46: PetscViewerDestroy(&vb);
48: /* Load the matrix as BAIJ format */
49: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&vc);
50: MatCreate(PETSC_COMM_WORLD,&C);
51: MatSetType(C,MATSEQBAIJ);
52: MatLoad(C,vc);
53: PetscViewerDestroy(&vc);
55: MatGetSize(A,&m,&n);
56: MatGetSize(B,&m2,&n2);
57: if (m!=m2) SETERRQ(PETSC_COMM_SELF,1,"Matrices are of different size. Cannot run this example");
58:
59: /* Test MatEqual() */
60: MatEqual(B,C,&tflg);
61: if (!tflg) SETERRQ(PETSC_COMM_SELF,1,"MatEqual() failed");
63: /* Test MatGetDiagonal() */
64: VecCreateSeq(PETSC_COMM_SELF,m,&x);
65: VecCreateSeq(PETSC_COMM_SELF,m,&y);
67: MatGetDiagonal(A,x);
68: MatGetDiagonal(B,y);
69:
70: VecEqual(x,y,&tflg);
71: if (!tflg) SETERRQ(PETSC_COMM_SELF,1,"MatGetDiagonal() failed");
73: /* Test MatDiagonalScale() */
74: PetscRandomCreate(PETSC_COMM_SELF,&r);
75: PetscRandomSetFromOptions(r);
76: VecSetRandom(x,r);
77: VecSetRandom(y,r);
79: MatDiagonalScale(A,x,y);
80: MatDiagonalScale(B,x,y);
81: MatMult(A,x,y);
82: VecNorm(y,NORM_2,&norm1);
83: MatMult(B,x,y);
84: VecNorm(y,NORM_2,&norm2);
85: rnorm = ((norm1-norm2)*100)/norm1;
86: if (rnorm<-0.1 || rnorm>0.01) {
87: PetscPrintf(PETSC_COMM_SELF,"Norm1=%e Norm2=%e\n",norm1,norm2);
88: SETERRQ(PETSC_COMM_SELF,1,"MatDiagonalScale() failed");
89: }
91: /* Test MatGetRow()/ MatRestoreRow() */
92: for (ct=0; ct<100; ct++) {
93: PetscRandomGetValue(r,&rval);
94: row = (int)(rval*m);
95: MatGetRow(A,row,&ncols1,&cols1,&vals1);
96: MatGetRow(B,row,&ncols2,&cols2,&vals2);
97:
98: for (i=0,j=0; i<ncols1 && j<ncols2; i++) {
99: while (cols2[j] != cols1[i]) j++;
100: if (vals1[i] != vals2[j]) SETERRQ(PETSC_COMM_SELF,1,"MatGetRow() failed - vals incorrect.");
101: }
102: if (i<ncols1) SETERRQ(PETSC_COMM_SELF,1,"MatGetRow() failed - cols incorrect");
103:
104: MatRestoreRow(A,row,&ncols1,&cols1,&vals1);
105: MatRestoreRow(B,row,&ncols2,&cols2,&vals2);
106: }
107:
108: MatDestroy(&A);
109: MatDestroy(&B);
110: MatDestroy(&C);
111: VecDestroy(&x);
112: VecDestroy(&y);
113: PetscRandomDestroy(&r);
114: PetscFinalize();
115: #endif
116: return 0;
117: }