Actual source code: ex78.c


  2: static char help[] = "Reads in a matrix in ASCII MATLAB format (I,J,A), read in vectors rhs and exact_solu in ASCII format.\n\
  3: Writes them using the PETSc sparse format.\n\
  4: Note: I and J start at 1, not 0, use -noshift if indices in file start with zero!\n\
  5: Input parameters are:\n\
  6:   -Ain  <filename> : input matrix in ascii format\n\
  7:   -rhs  <filename> : input rhs in ascii format\n\
  8:   -solu  <filename> : input true solution in ascii format\n\\n";

 10: /*
 11: Example: ./ex78 -Ain Ain -rhs rhs -solu solu -noshift -mat_view
 12:  with the datafiles in the following format:
 13: Ain (I and J start at 0):
 14: ------------------------
 15: 3 3 6
 16: 0 0 1.0
 17: 0 1 2.0
 18: 1 0 3.0
 19: 1 1 4.0
 20: 1 2 5.0
 21: 2 2 6.0

 23: rhs
 24: ---
 25: 0 3.0
 26: 1 12.0
 27: 2 6.0

 29: solu
 30: ----
 31: 0 1.0
 32: 0 1.0
 33: 0 1.0
 34: */

 36: #include <petscmat.h>

 38: int main(int argc,char **args)
 39: {
 40:   Mat            A = NULL;
 41:   Vec            b,u = NULL,u_tmp;
 42:   char           Ain[PETSC_MAX_PATH_LEN],rhs[PETSC_MAX_PATH_LEN],solu[PETSC_MAX_PATH_LEN];
 43:   int            m,n = 0,nz,dummy; /* these are fscaned so kept as int */
 44:   PetscInt       i,col,row,shift = 1,sizes[3],nsizes;
 45:   PetscScalar    val;
 46:   PetscReal      res_norm;
 47:   FILE           *Afile,*bfile,*ufile;
 48:   PetscViewer    view;
 49:   PetscBool      flg_A,flg_b,flg_u,flg;
 50:   PetscMPIInt    size;

 52:   PetscInitialize(&argc,&args,(char*)0,help);
 53:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

 56:   /* Read in matrix, rhs and exact solution from ascii files */
 57:   PetscOptionsGetString(NULL,NULL,"-Ain",Ain,sizeof(Ain),&flg_A);
 58:   PetscOptionsHasName(NULL,NULL,"-noshift",&flg);
 59:   if (flg) shift = 0;
 60:   if (flg_A) {
 61:     PetscPrintf(PETSC_COMM_SELF,"\n Read matrix in ascii format ...\n");
 62:     PetscFOpen(PETSC_COMM_SELF,Ain,"r",&Afile);
 63:     nsizes = 3;
 64:     PetscOptionsGetIntArray(NULL,NULL,"-nosizesinfile",sizes,&nsizes,&flg);
 65:     if (flg) {
 67:       m  = sizes[0];
 68:       n  = sizes[1];
 69:       nz = sizes[2];
 70:     } else {
 72:     }
 73:     PetscPrintf(PETSC_COMM_SELF,"m: %d, n: %d, nz: %d \n", m,n,nz);
 75:     MatCreate(PETSC_COMM_SELF,&A);
 76:     MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,n);
 77:     MatSetFromOptions(A);
 78:     MatSeqAIJSetPreallocation(A,nz/m,NULL);
 79:     MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);

 81:     for (i=0; i<nz; i++) {
 83:       row -= shift; col -= shift;  /* set index set starts at 0 */
 84:       MatSetValues(A,1,&row,1,&col,&val,INSERT_VALUES);
 85:     }
 86:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 87:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 88:     fclose(Afile);
 89:   }

 91:   PetscOptionsGetString(NULL,NULL,"-rhs",rhs,sizeof(rhs),&flg_b);
 92:   if (flg_b) {
 93:     VecCreate(PETSC_COMM_SELF,&b);
 94:     VecSetSizes(b,PETSC_DECIDE,n);
 95:     VecSetFromOptions(b);
 96:     PetscPrintf(PETSC_COMM_SELF,"\n Read rhs in ascii format ...\n");
 97:     PetscFOpen(PETSC_COMM_SELF,rhs,"r",&bfile);
 98:     for (i=0; i<n; i++) {
100:       VecSetValues(b,1,&i,&val,INSERT_VALUES);
101:     }
102:     VecAssemblyBegin(b);
103:     VecAssemblyEnd(b);
104:     fclose(bfile);
105:   }

107:   PetscOptionsGetString(NULL,NULL,"-solu",solu,sizeof(solu),&flg_u);
108:   if (flg_u) {
109:     VecCreate(PETSC_COMM_SELF,&u);
110:     VecSetSizes(u,PETSC_DECIDE,n);
111:     VecSetFromOptions(u);
112:     PetscPrintf(PETSC_COMM_SELF,"\n Read exact solution in ascii format ...\n");
113:     PetscFOpen(PETSC_COMM_SELF,solu,"r",&ufile);
114:     for (i=0; i<n; i++) {
116:       VecSetValues(u,1,&i,&val,INSERT_VALUES);
117:     }
118:     VecAssemblyBegin(u);
119:     VecAssemblyEnd(u);
120:     fclose(ufile);
121:   }

123:   /* Write matrix, rhs and exact solution in Petsc binary file */
124:   PetscPrintf(PETSC_COMM_SELF,"\n Write matrix in binary to 'matrix.dat' ...\n");
125:   PetscViewerBinaryOpen(PETSC_COMM_SELF,"matrix.dat",FILE_MODE_WRITE,&view);
126:   MatView(A,view);

128:   if (flg_b) { /* Write rhs in Petsc binary file */
129:     PetscPrintf(PETSC_COMM_SELF,"\n Write rhs in binary to 'matrix.dat' ...\n");
130:     VecView(b,view);
131:   }
132:   if (flg_u) {
133:     PetscPrintf(PETSC_COMM_SELF,"\n Write exact solution in binary to 'matrix.dat' ...\n");
134:     VecView(u,view);
135:   }

137:   /* Check accuracy of the data */
138:   if (flg_A & flg_b & flg_u) {
139:     VecDuplicate(u,&u_tmp);
140:     MatMult(A,u,u_tmp);
141:     VecAXPY(u_tmp,-1.0,b);
142:     VecNorm(u_tmp,NORM_2,&res_norm);
143:     PetscPrintf(PETSC_COMM_SELF,"\n Accuracy of the reading data: | b - A*u |_2 : %g \n",res_norm);
144:     VecDestroy(&u_tmp);
145:   }

147:   MatDestroy(&A);
148:   if (flg_b) VecDestroy(&b);
149:   if (flg_u) VecDestroy(&u);
150:   PetscViewerDestroy(&view);
151:   PetscFinalize();
152:   return 0;
153: }

155: /*TEST

157:    build:
158:       requires:  !defined(PETSC_USE_64BIT_INDICES) double !complex

160:    test:
161:       requires: datafilespath
162:       args: -Ain ${DATAFILESPATH}/matrices/indefinite/afiro_A.dat

164: TEST*/