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:   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(NULL,NULL,"-verbose",&verbose,NULL);
 23:   PetscOptionsGetString(NULL,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(NULL,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(NULL,NULL,"-mbs",&mbs,NULL);
 64:     PetscOptionsGetInt(NULL,NULL,"-bs",&bs,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,NULL,o_nz,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 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"C must be symmetric for this example");
103:     MatDestroy(&Ctrans);
104:   }
105:   /*MatView(C,PETSC_VIEWER_STDOUT_WORLD);*/

107:   /* convert C to other formats */
108:   for (i=0; i<ntypes; i++) {
109:     MatConvert(C,type[i],MAT_INITIAL_MATRIX,&A);
110:     MatMultEqual(A,C,10,&equal);
111:     if (!equal) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"Error in conversion from BAIJ to %s",type[i]);
112:     for (j=i+1; j<ntypes; j++) {
113:       if (verbose>0) {
114:         PetscPrintf(PETSC_COMM_WORLD," \n[%d] test conversion between %s and %s\n",rank,type[i],type[j]);
115:       }

117:       if (!rank && verbose) printf("Convert %s A to %s B\n",type[i],type[j]);
118:       MatConvert(A,type[j],MAT_INITIAL_MATRIX,&B);
119:       /*
120:       if (j == 2) {
121:         PetscPrintf(PETSC_COMM_SELF," A: %s\n",type[i]);
122:         MatView(A,PETSC_VIEWER_STDOUT_WORLD);
123:         PetscPrintf(PETSC_COMM_SELF," B: %s\n",type[j]);
124:         MatView(B,PETSC_VIEWER_STDOUT_WORLD);
125:       }
126:        */
127:       MatMultEqual(A,B,10,&equal);
128:       if (!equal) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"Error in conversion from %s to %s",type[i],type[j]);

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

136:         MatDestroy(&D);
137:       }
138:       MatDestroy(&B);
139:       MatDestroy(&D);
140:     }

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

149:     MatDestroy(&A);
150:   }
151:   MatDestroy(&C);

153:   PetscFinalize();
154:   return 0;
155: }