Actual source code: ex55.c

petsc-3.3-p6 2013-02-11
  2: static char help[] = "Tests converting a matrix to another format with MatConvert().\n\n";

  4: #include <petscmat.h>
  5: /* Usage: mpiexec -n <np> ex55 -verbose <0 or 1> */

  9: int main(int argc,char **args)
 10: {
 11:   Mat            C,A,B,D;
 13:   PetscInt       i,j,ntypes,bs,mbs,m,block,d_nz=6, o_nz=3,col[3],row,verbose=0;
 14:   PetscMPIInt    size,rank;
 15:   const MatType  type[9];
 16:   char           file[PETSC_MAX_PATH_LEN];
 17:   PetscViewer    fd;
 18:   PetscBool      equal,flg_loadmat,flg;
 19:   PetscScalar    value[3];

 21:   PetscInitialize(&argc,&args,(char *)0,help);
 22:   PetscOptionsGetInt(PETSC_NULL,"-verbose",&verbose,PETSC_NULL);
 23:   PetscOptionsGetString(PETSC_NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg_loadmat);
 24:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 25:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

 27:   PetscOptionsHasName(PETSC_NULL,"-testseqaij",&flg);
 28:   if (flg){
 29:     if (size == 1){
 30:       type[0] = MATSEQAIJ;
 31:     } else {
 32:       type[0] = MATMPIAIJ;
 33:     }
 34:   } else {
 35:     type[0] = MATAIJ;
 36:   }
 37:   if (size == 1){
 38:     ntypes = 3;
 39:     type[1] = MATSEQBAIJ;
 40:     type[2] = MATSEQSBAIJ;
 41:   } else {
 42:     ntypes = 3;
 43:     type[1] = MATMPIBAIJ;
 44:     type[2] = MATMPISBAIJ;
 45:   }

 47:   /* input matrix C */
 48:   if (flg_loadmat){
 49:     /* Open binary file. */
 50:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);

 52:     /* Load a baij matrix, then destroy the viewer. */
 53:     MatCreate(PETSC_COMM_WORLD,&C);
 54:     if (size == 1){
 55:       MatSetType(C,MATSEQBAIJ);
 56:     } else {
 57:       MatSetType(C,MATMPIBAIJ);
 58:     }
 59:     MatLoad(C,fd);
 60:     PetscViewerDestroy(&fd);
 61:   } else { /* Create a baij mat with bs>1  */
 62:     bs = 2; mbs=8;
 63:     PetscOptionsGetInt(PETSC_NULL,"-mbs",&mbs,PETSC_NULL);
 64:     PetscOptionsGetInt(PETSC_NULL,"-bs",&bs,PETSC_NULL);
 65:     if (bs <= 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG," bs must be >1 in this case");
 66:     m = mbs*bs;
 67:     MatCreateBAIJ(PETSC_COMM_WORLD,bs,PETSC_DECIDE,PETSC_DECIDE,m,m,d_nz,PETSC_NULL,o_nz,PETSC_NULL,&C);
 68:     for (block=0; block<mbs; block++){
 69:       /* diagonal blocks */
 70:       value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
 71:       for (i=1+block*bs; i<bs-1+block*bs; i++) {
 72:         col[0] = i-1; col[1] = i; col[2] = i+1;
 73:         MatSetValues(C,1,&i,3,col,value,INSERT_VALUES);
 74:       }
 75:       i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
 76:       value[0]=-1.0; value[1]=4.0;
 77:       MatSetValues(C,1,&i,2,col,value,INSERT_VALUES);

 79:       i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs;
 80:       value[0]=4.0; value[1] = -1.0;
 81:       MatSetValues(C,1,&i,2,col,value,INSERT_VALUES);
 82:     }
 83:     /* off-diagonal blocks */
 84:     value[0]=-1.0;
 85:     for (i=0; i<(mbs-1)*bs; i++){
 86:       col[0]=i+bs;
 87:       MatSetValues(C,1,&i,1,col,value,INSERT_VALUES);
 88:       col[0]=i; row=i+bs;
 89:       MatSetValues(C,1,&row,1,col,value,INSERT_VALUES);
 90:     }
 91:     MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
 92:     MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
 93:   }

 95:   {
 96:     /* Check the symmetry of C because it will be converted to a sbaij matrix */
 97:     Mat Ctrans;
 98:     MatTranspose(C, MAT_INITIAL_MATRIX,&Ctrans);
 99:     MatEqual(C, Ctrans, &flg);
100:     if (flg) {
101:       MatSetOption(C,MAT_SYMMETRIC,PETSC_TRUE);
102:     } else {
103:       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"C must be symmetric for this example");
104:     }
105:     MatDestroy(&Ctrans);
106:   }
107:   //MatView(C,PETSC_VIEWER_STDOUT_WORLD);
108: 
109:   /* convert C to other formats */
110:   for (i=0; i<ntypes; i++) {
111:     MatConvert(C,type[i],MAT_INITIAL_MATRIX,&A);
112:     MatMultEqual(A,C,10,&equal);
113:     if (!equal) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"Error in conversion from BAIJ to %s",type[i]);
114:     for (j=i+1; j<ntypes; j++) {
115:       if (verbose>0) {
116:         PetscPrintf(PETSC_COMM_WORLD," \n[%d] test conversion between %s and %s\n",rank,type[i],type[j]);
117:       }
118: 
119:       if (!rank && verbose) printf("Convert %s A to %s B\n",type[i],type[j]);
120:       MatConvert(A,type[j],MAT_INITIAL_MATRIX,&B);
121:       /*
122:       if (j == 2){
123:         PetscPrintf(PETSC_COMM_SELF," A: %s\n",type[i]);
124:         MatView(A,PETSC_VIEWER_STDOUT_WORLD);
125:         PetscPrintf(PETSC_COMM_SELF," B: %s\n",type[j]);
126:         MatView(B,PETSC_VIEWER_STDOUT_WORLD);
127:       }
128:        */
129:       MatMultEqual(A,B,10,&equal);
130:       if (!equal) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"Error in conversion from %s to %s",type[i],type[j]);

132:       if (size == 1 || j != 2){ /* Matconvert from mpisbaij mat to other formats are not supported */
133:         if (!rank && verbose) printf("Convert %s B to %s D\n",type[j],type[i]);
134:         MatConvert(B,type[i],MAT_INITIAL_MATRIX,&D);
135:         MatMultEqual(B,D,10,&equal);
136:         if (!equal) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"Error in conversion from %s to %s",type[j],type[i]);

138:         MatDestroy(&D);
139:       }
140:       MatDestroy(&B);
141:       MatDestroy(&D);
142:     }

144:     /* Test in-place convert */
145:     if (size == 1){ /* size > 1 is not working yet! */
146:       j = (i+1)%ntypes;
147:       /* printf("[%d] i: %d, j: %d\n",rank,i,j); */
148:       MatConvert(A,type[j],MAT_REUSE_MATRIX,&A);
149:     }

151:     MatDestroy(&A);
152:   }
153:   MatDestroy(&C);

155:   PetscFinalize();
156:   return 0;
157: }