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 followig 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

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

 55:   PetscInitialize(&argc,&args,(char *)0,help);
 56:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 57:   if (size != 1) SETERRQ(PETSC_ERR_SUP,"This is a uniprocessor example only!");

 59:   /* Read in matrix, rhs and exact solution from ascii files */
 60:   PetscOptionsGetString(PETSC_NULL,"-Ain",Ain,PETSC_MAX_PATH_LEN-1,&flg_A);
 61:   PetscOptionsHasName(PETSC_NULL,"-noshift",&flg);
 62:   if (flg) shift = 0;
 63:   if (flg_A){
 64:     PetscPrintf(PETSC_COMM_SELF,"\n Read matrix in ascii format ...\n");
 65:     PetscFOpen(PETSC_COMM_SELF,Ain,"r",&Afile);
 66:     nsizes = 3;
 67:     PetscOptionsGetIntArray(PETSC_NULL,"-nosizesinfile",sizes,&nsizes,&flg);
 68:     if (flg) {
 69:       if (nsizes != 3) SETERRQ(1,"Must pass in three m,n,nz as arguments for -nosizesinfile");
 70:       m = sizes[0];
 71:       n = sizes[1];
 72:       nz = sizes[2];
 73:     } else {
 74:       fscanf(Afile,"%d %d %d\n",&m,&n,&nz);
 75:     }
 76:     printf("m: %d, n: %d, nz: %d \n", m,n,nz);
 77:     if (m != n) SETERRQ(PETSC_ERR_ARG_SIZ, "Number of rows, cols must be same for this example\n");
 78:     MatCreate(PETSC_COMM_SELF,&A);
 79:     MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,n);
 80:     MatSetFromOptions(A);
 81:     MatSeqAIJSetPreallocation(A,nz/m,PETSC_NULL);

 83:     for (i=0; i<nz; i++) {
 84:       fscanf(Afile,"%d %d %le\n",&row,&col,(double*)&val);
 85:       row -= shift; col -= shift;  /* set index set starts at 0 */
 86:       MatSetValues(A,1,&row,1,&col,&val,INSERT_VALUES);
 87:     }
 88:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 89:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 90:     fflush(stdout);
 91:     fclose(Afile);
 92:   }

 94:   PetscOptionsGetString(PETSC_NULL,"-rhs",rhs,PETSC_MAX_PATH_LEN-1,&flg_b);
 95:   if (flg_b){
 96:     VecCreate(PETSC_COMM_SELF,&b);
 97:     VecSetSizes(b,PETSC_DECIDE,n);
 98:     VecSetFromOptions(b);
 99:     PetscPrintf(PETSC_COMM_SELF,"\n Read rhs in ascii format ...\n");
100:     PetscFOpen(PETSC_COMM_SELF,rhs,"r",&bfile);
101:     for (i=0; i<n; i++) {
102:       fscanf(bfile,"%d %le\n",&dummy,(double*)&val);
103:       VecSetValues(b,1,&i,&val,INSERT_VALUES);
104:     }
105:     VecAssemblyBegin(b);
106:     VecAssemblyEnd(b);
107:     fflush(stdout);
108:     fclose(bfile);
109:   }

111:   PetscOptionsGetString(PETSC_NULL,"-solu",solu,PETSC_MAX_PATH_LEN-1,&flg_u);
112:   if (flg_u){
113:     VecCreate(PETSC_COMM_SELF,&u);
114:     VecSetSizes(u,PETSC_DECIDE,n);
115:     VecSetFromOptions(u);
116:     PetscPrintf(PETSC_COMM_SELF,"\n Read exact solution in ascii format ...\n");
117:     PetscFOpen(PETSC_COMM_SELF,solu,"r",&ufile);
118:     for (i=0; i<n; i++) {
119:       fscanf(ufile,"%d  %le\n",&dummy,(double*)&val);
120:       VecSetValues(u,1,&i,&val,INSERT_VALUES);
121:     }
122:     VecAssemblyBegin(u);
123:     VecAssemblyEnd(u);
124:     fflush(stdout);
125:     fclose(ufile);
126:   }
127: 
128:   /* Check accuracy of the data */
129:   if (flg_A & flg_b & flg_u){
130:     VecDuplicate(u,&u_tmp);
131:     MatMult(A,u,u_tmp);
132:     VecAXPY(u_tmp,-1.0,b);
133:     VecNorm(u_tmp,NORM_2,&res_norm);
134:     printf("\n Accuracy of the reading data: | b - A*u |_2 : %g \n",res_norm);

136:     /* Write the matrix, rhs and exact solution in Petsc binary file */
137:     PetscPrintf(PETSC_COMM_SELF,"\n Write matrix and rhs in binary to 'matrix.dat' ...\n");
138:     PetscViewerBinaryOpen(PETSC_COMM_SELF,"matrix.dat",FILE_MODE_WRITE,&view);
139:     MatView(A,view);
140:     VecView(b,view);
141:     VecView(u,view);
142:     PetscViewerDestroy(view);

144:     VecDestroy(b);
145:     VecDestroy(u);
146:     VecDestroy(u_tmp);
147:     MatDestroy(A);
148:   }
149:   PetscFinalize();
150:   return 0;
151: }