Actual source code: ex55.c
1: /*$Id: ex55.c,v 1.18 2001/04/10 19:35:44 bsmith Exp $*/
3: static char help[] = "Tests converting a matrix to another format with MatConvert().\n\n";
5: #include petscmat.h
9: int main(int argc,char **args)
10: {
11: Mat C,A,B,D;
12: int ierr,i,j,k,ntypes = 3,size;
13: MatType type[9] = {MATSEQAIJ,MATSEQBAIJ,MATSEQSBAIJ,MATMPIROWBS};
14: char file[128];
15: PetscViewer fd;
16: PetscTruth equal,flg_loadmat;
17: int bs,mbs,m,n,block,d_nz=3, o_nz=3,col[3],row,msglvl=0;
18: PetscScalar one=1.0, neg_one=-1.0, value[3], four=4.0;
20: PetscInitialize(&argc,&args,(char *)0,help);
21: PetscOptionsGetString(PETSC_NULL,"-f",file,127,&flg_loadmat);
22: MPI_Comm_size(PETSC_COMM_WORLD,&size);
23: if (size > 1) ntypes = 5;
25: /* input matrix C */
26: if (flg_loadmat){
27: /* Open binary file. Note that we use PETSC_FILE_RDONLY to indicate
28: reading from this file. */
29: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,PETSC_FILE_RDONLY,&fd);
31: /* Load the matrix, then destroy the viewer. */
32: MatLoad(fd,MATMPIAIJ,&C);
33: PetscViewerDestroy(fd);
34: bs = 1; /* assume the loaded matrix has block size 1 */
35: } else { /* Convert a baij mat with bs>1 to other formats */
36: bs = 2; mbs=8;
37: PetscOptionsGetInt(PETSC_NULL,"-mbs",&mbs,PETSC_NULL);
38: PetscOptionsGetInt(PETSC_NULL,"-bs",&bs,PETSC_NULL);
39: if (bs <= 1) SETERRQ(PETSC_ERR_ARG_WRONG," bs must be >1 in this case");
40: m = mbs*bs;
41: MatCreateMPIBAIJ(PETSC_COMM_WORLD,bs,PETSC_DECIDE,PETSC_DECIDE,m,m,d_nz,PETSC_NULL,o_nz,PETSC_NULL,&C);
42: for (block=0; block<mbs; block++){
43: /* diagonal blocks */
44: value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
45: for (i=1+block*bs; i<bs-1+block*bs; i++) {
46: col[0] = i-1; col[1] = i; col[2] = i+1;
47: MatSetValues(C,1,&i,3,col,value,INSERT_VALUES);
48: }
49: i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
50: value[0]=-1.0; value[1]=4.0;
51: MatSetValues(C,1,&i,2,col,value,INSERT_VALUES);
53: i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
54: value[0]=4.0; value[1] = -1.0;
55: MatSetValues(C,1,&i,2,col,value,INSERT_VALUES);
56: }
57: /* off-diagonal blocks */
58: value[0]=-1.0;
59: for (i=0; i<(mbs-1)*bs; i++){
60: col[0]=i+bs;
61: MatSetValues(C,1,&i,1,col,value,INSERT_VALUES);
62: col[0]=i; row=i+bs;
63: MatSetValues(C,1,&row,1,col,value,INSERT_VALUES);
64: }
65:
66: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
67: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
68: }
69:
70: /* convert C to other formats */
71: for (i=0; i<ntypes; i++) {
72: MatConvert(C,type[i],&A);
73: for (j=0; j<ntypes; j++) {
74: if (j==i) continue;
75: if (msglvl>0)
76: PetscPrintf(PETSC_COMM_SELF," test conversion between %s and %s\n",type[i],type[j]);
78: MatConvert(A,type[j],&B);
79: MatConvert(B,type[i],&D);
81: if (bs == 1){
82: MatEqual(A,D,&equal);
83: if (!equal){
84: PetscPrintf(PETSC_COMM_SELF," A: %s\n",type[i]);
85: MatView(A,PETSC_VIEWER_STDOUT_WORLD);
86: PetscPrintf(PETSC_COMM_SELF," B: %s\n",type[j]);
87: MatView(B,PETSC_VIEWER_STDOUT_WORLD);
88: PetscPrintf(PETSC_COMM_SELF," D: %s\n",type[i]);
89: MatView(D,PETSC_VIEWER_STDOUT_WORLD);
90: SETERRQ2(1,"Error in conversion from %s to %s",type[i],type[j]);
91: }
92: } else { /* bs > 1 */
93: Vec x,s1,s2;
94: PetscRandom rctx;
95: PetscReal r1,r2,tol=1.e-10;
97: MatGetSize(A,&m,&n);
98: PetscRandomCreate(PETSC_COMM_WORLD,RANDOM_DEFAULT,&rctx);
99: VecCreate(PETSC_COMM_WORLD,&x);
100: VecSetSizes(x,m,PETSC_DECIDE);
101: VecSetFromOptions(x);
102: VecDuplicate(x,&s1);
103: VecDuplicate(x,&s2);
105: for (k=0; k<10; k++) {
106: VecSetRandom(rctx,x);
107: MatMult(A,x,s1);
108: MatMult(B,x,s2);
109: /* MatMult(D,x,s2); */
110: VecNorm(s1,NORM_1,&r1);
111: VecNorm(s2,NORM_1,&r2);
112: r1 -= r2;
113: if (r1<-tol || r1>tol) {
114: SETERRQ2(1,"Error in conversion from %s to %s",type[i],type[j]);
115: }
116: }
117: PetscRandomDestroy(rctx);
118: VecDestroy(x);
119: VecDestroy(s1);
120: VecDestroy(s2);
121: }
122: MatDestroy(B);
123: MatDestroy(D);
124: }
125: MatDestroy(A);
126: }
127: MatDestroy(C);
129: PetscFinalize();
130: return 0;
131: }