Actual source code: ex78.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  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:   PetscBool      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_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");

 59:   /* Read in matrix, rhs and exact solution from ascii files */
 60:   PetscOptionsGetString(NULL,NULL,"-Ain",Ain,PETSC_MAX_PATH_LEN,&flg_A);
 61:   PetscOptionsHasName(NULL,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(NULL,NULL,"-nosizesinfile",sizes,&nsizes,&flg);
 68:     if (flg) {
 69:       if (nsizes != 3) SETERRQ(PETSC_COMM_WORLD,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:     PetscPrintf(PETSC_COMM_SELF,"m: %d, n: %d, nz: %d \n", m,n,nz);
 77:     if (m != n) SETERRQ(PETSC_COMM_SELF,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,NULL);
 82:     MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);

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

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

112:   PetscOptionsGetString(NULL,NULL,"-solu",solu,PETSC_MAX_PATH_LEN,&flg_u);
113:   if (flg_u) {
114:     VecCreate(PETSC_COMM_SELF,&u);
115:     VecSetSizes(u,PETSC_DECIDE,n);
116:     VecSetFromOptions(u);
117:     PetscPrintf(PETSC_COMM_SELF,"\n Read exact solution in ascii format ...\n");
118:     PetscFOpen(PETSC_COMM_SELF,solu,"r",&ufile);
119:     for (i=0; i<n; i++) {
120:       fscanf(ufile,"%d  %le\n",&dummy,(double*)&val);
121:       VecSetValues(u,1,&i,&val,INSERT_VALUES);
122:     }
123:     VecAssemblyBegin(u);
124:     VecAssemblyEnd(u);
125:     fflush(stdout);
126:     fclose(ufile);
127:   }

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

134:   if (flg_b) { /* Write rhs in Petsc binary file */
135:     PetscPrintf(PETSC_COMM_SELF,"\n Write rhs in binary to 'matrix.dat' ...\n");
136:     VecView(b,view);
137:   }
138:   if (flg_u) {
139:     PetscPrintf(PETSC_COMM_SELF,"\n Write exact solution in binary to 'matrix.dat' ...\n");
140:     VecView(u,view);
141:   }

143:   /* Check accuracy of the data */
144:   if (flg_A & flg_b & flg_u) {
145:     VecDuplicate(u,&u_tmp);
146:     MatMult(A,u,u_tmp);
147:     VecAXPY(u_tmp,-1.0,b);
148:     VecNorm(u_tmp,NORM_2,&res_norm);
149:     PetscPrintf(PETSC_COMM_SELF,"\n Accuracy of the reading data: | b - A*u |_2 : %g \n",res_norm);
150:     VecDestroy(&u_tmp);
151:   }

153:   MatDestroy(&A);
154:   if (flg_b) {VecDestroy(&b);}
155:   if (flg_u) {VecDestroy(&u);}
156:   PetscViewerDestroy(&view);
157:   PetscFinalize();
158:   return 0;
159: }