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: }