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> */

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

 19:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
 20:   PetscOptionsGetInt(NULL,NULL,"-verbose",&verbose,NULL);
 21:   PetscOptionsGetString(NULL,NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg_loadmat);
 22:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 23:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

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

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

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

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

 93:   {
 94:     /* Check the symmetry of C because it will be converted to a sbaij matrix */
 95:     Mat Ctrans;
 96:     MatTranspose(C, MAT_INITIAL_MATRIX,&Ctrans);
 97:     MatEqual(C, Ctrans, &flg);
 98:     if (flg) {
 99:       MatSetOption(C,MAT_SYMMETRIC,PETSC_TRUE);
100:     } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"C must be symmetric for this example");
101:     MatDestroy(&Ctrans);
102:   }
103:   /*MatView(C,PETSC_VIEWER_STDOUT_WORLD);*/

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

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

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

134:         MatDestroy(&D);
135:       }
136:       MatDestroy(&B);
137:       MatDestroy(&D);
138:     }

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

147:     MatDestroy(&A);
148:   }
149:   MatDestroy(&C);

151:   PetscFinalize();
152:   return ierr;
153: }