Actual source code: ex55.c

  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:       MatLoad(C,fd);
 57:     } else {
 58:       MatSetType(C,MATMPIBAIJ);
 59:       MatLoad(C,fd);
 60:     }
 61:     PetscViewerDestroy(&fd);
 62:   } else { /* Create a baij mat with bs>1  */
 63:     bs = 2; mbs=8;
 64:     PetscOptionsGetInt(PETSC_NULL,"-mbs",&mbs,PETSC_NULL);
 65:     PetscOptionsGetInt(PETSC_NULL,"-bs",&bs,PETSC_NULL);
 66:     if (bs <= 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG," bs must be >1 in this case");
 67:     m = mbs*bs;
 68:     if (size == 1){
 69:       MatCreateSeqBAIJ(PETSC_COMM_SELF,bs,m,m,d_nz,PETSC_NULL,&C);
 70:     } else {
 71:       MatCreateMPIBAIJ(PETSC_COMM_WORLD,bs,PETSC_DECIDE,PETSC_DECIDE,m,m,d_nz,PETSC_NULL,o_nz,PETSC_NULL,&C);
 72:     }
 73:     for (block=0; block<mbs; block++){
 74:       /* diagonal blocks */
 75:       value[0] = -1.0; value[1] = 4.0; value[2] = -1.0;
 76:       for (i=1+block*bs; i<bs-1+block*bs; i++) {
 77:         col[0] = i-1; col[1] = i; col[2] = i+1;
 78:         MatSetValues(C,1,&i,3,col,value,INSERT_VALUES);
 79:       }
 80:       i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs;
 81:       value[0]=-1.0; value[1]=4.0;
 82:       MatSetValues(C,1,&i,2,col,value,INSERT_VALUES);

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

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

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

143:         MatDestroy(&D);
144:       }
145:       MatDestroy(&B);
146:       MatDestroy(&D);
147:     }

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

156:     MatDestroy(&A);
157:   }
158:   MatDestroy(&C);

160:   PetscFinalize();
161:   return 0;
162: }