Actual source code: ex78.c

  1: /*$Id: ex78.c,v 1.14 2001/08/07 21:30:08 bsmith Exp $*/

  3: 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\
  4: Writes them using the PETSc sparse format.\n\
  5: Note: I and J start at 1, not 0!\n\
  6: Input parameters are:\n\
  7:   -Ain  <filename> : input matrix in ascii format\n\
  8:   -bin  <filename> : input rhs in ascii format\n\
  9:   -uin  <filename> : input true solution in ascii format\n\
 10: Run this program: ex33h -Ain Ain -bin bin -uin uin\n\n";

 12:  #include petscmat.h

 16: int main(int argc,char **args)
 17: {
 18:   Mat         A;
 19:   Vec         b,u,u_tmp;
 20:   char        Ain[128],bin[128],uin[128];
 21:   int         i,m,n,nz,ierr,*ib=0,col_i,row_i;
 22:   PetscScalar val_i,*work=0,mone=-1.0;
 23:   PetscReal   *col=0,*row=0,res_norm,*val=0,*bval=0,*uval=0;
 24:   FILE        *Afile,*bfile,*ufile;
 25:   PetscViewer view;
 26:   PetscTruth  flg_A,flg_b,flg_u;

 28:   PetscInitialize(&argc,&args,(char *)0,help);

 30:   /* Read in matrix, rhs and exact solution from an ascii file */
 31:   PetscOptionsGetString(PETSC_NULL,"-Ain",Ain,127,&flg_A);
 32:   if (flg_A){
 33:     PetscPrintf(PETSC_COMM_SELF,"\n Read matrix in ascii format ...\n");
 34:     PetscFOpen(PETSC_COMM_SELF,Ain,"r",&Afile);
 35:     fscanf(Afile,"%d %d %d\n",&m,&n,&nz);
 36:     printf("m: %d, n: %d, nz: %d \n", m,n,nz);
 37:     if (m != n) SETERRQ(PETSC_ERR_ARG_SIZ, "Number of rows, cols must be same for SBAIJ format\n");

 39:     MatCreateSeqBAIJ(PETSC_COMM_SELF,1,n,n,20,0,&A);
 40:     VecCreateSeq(PETSC_COMM_SELF,n,&b);
 41:     VecCreateSeq(PETSC_COMM_SELF,n,&u);

 43:     PetscMalloc(nz*sizeof(PetscReal),&val);
 44:     PetscMalloc(nz*sizeof(PetscReal),&row);
 45:     PetscMalloc(nz*sizeof(PetscReal),&col);
 46:     for (i=0; i<nz; i++) {
 47:       fscanf(Afile,"%le %le %le\n",row+i,col+i,val+i); /* modify according to data file! */
 48:       row[i]--; col[i]--;  /* set index set starts at 0 */
 49:     }
 50:     printf("row[0]: %g, col[0]: %g, val: %g\n",row[0],col[0],val[0]);
 51:     printf("row[last]: %g, col: %g, val: %g\n",row[nz-1],col[nz-1],val[nz-1]);
 52:     fclose(Afile);
 53:   }

 55:   PetscOptionsGetString(PETSC_NULL,"-bin",bin,127,&flg_b);
 56:   if (flg_b){
 57:     PetscPrintf(PETSC_COMM_SELF,"\n Read rhs in ascii format ...\n");
 58:     PetscFOpen(PETSC_COMM_SELF,bin,"r",&bfile);
 59:     PetscMalloc(n*sizeof(PetscReal),&bval);
 60:     PetscMalloc(n*sizeof(PetscScalar),&work);
 61:     PetscMalloc(n*sizeof(int),&ib);
 62:     for (i=0; i<n; i++) {
 63:       /* fscanf(bfile,"%d %le\n",ib+i,bval+i); ib[i]--;  */  /* modify according to data file! */
 64:       fscanf(bfile,"%le\n",bval+i); ib[i] = i;         /* modify according to data file! */
 65:     }
 66:     printf("bval[0]: %g, bval[%d]: %g\n",bval[0],ib[n-1],bval[n-1]);
 67:     fclose(bfile);
 68:   }

 70:   PetscOptionsGetString(PETSC_NULL,"-uin",uin,127,&flg_u);
 71:   if (flg_u){
 72:     PetscPrintf(PETSC_COMM_SELF,"\n Read exact solution in ascii format ...\n");
 73:     PetscFOpen(PETSC_COMM_SELF,uin,"r",&ufile);
 74:     PetscMalloc(n*sizeof(PetscReal),&uval);
 75:     for (i=0; i<n; i++) {
 76:       fscanf(ufile,"  %le\n",uval+i);  /* modify according to data file! */
 77:     }
 78:     printf("uval[0]: %g, uval[%d]: %g\n",uval[0], n-1, uval[n-1]);
 79:     fclose(ufile);
 80:   }

 82:   if(flg_A){
 83:     /*
 84:     for (i=0; i<n; i++){ 
 85:       MatSetValues(A,1,&i,1,&i,&zero,INSERT_VALUES); 
 86:     }
 87:     */
 88:     for (i=0; i<nz; i++) {
 89:       row_i =(int)row[i]; col_i =(int)col[i]; val_i = (PetscScalar)val[i];
 90:       MatSetValues(A,1,&row_i,1,&col_i,&val_i,INSERT_VALUES);
 91:     }
 92:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 93:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 94:     PetscFree(col);
 95:     PetscFree(val);
 96:     PetscFree(row);
 97:   }
 98:   if(flg_b){
 99:     for (i=0; i<n; i++) work[i]=(PetscScalar)bval[i];
100:     VecSetValues(b,n,ib,work,INSERT_VALUES);
101:     VecAssemblyBegin(b);
102:     VecAssemblyEnd(b);
103:     /* printf("b: \n"); VecView(b,PETSC_VIEWER_STDOUT_SELF); */
104:     PetscFree(bval);
105:   }

107:   if(flg_u & flg_b){
108:     for (i=0; i<n; i++) work[i]=(PetscScalar)uval[i];
109:     VecSetValues(u,n,ib,work,INSERT_VALUES);
110:     VecAssemblyBegin(u);
111:     VecAssemblyEnd(u);
112:     /* printf("u: \n"); VecView(u,PETSC_VIEWER_STDOUT_SELF); */
113:     PetscFree(uval);
114:   }
115: 
116:   if(flg_b) {
117:     PetscFree(ib);
118:     PetscFree(work);
119:   }
120:   /* Check accuracy of the data */
121:   if (flg_A & flg_b & flg_u){
122:     VecDuplicate(u,&u_tmp);
123:     MatMult(A,u,u_tmp);
124:     VecAXPY(&mone,b,u_tmp);
125:     VecNorm(u_tmp,NORM_2,&res_norm);
126:     printf("\n Accuracy of the reading data: | b - A*u |_2 : %g \n",res_norm);

128:   /* Write the matrix, rhs and exact solution in Petsc binary file */
129:     PetscPrintf(PETSC_COMM_SELF,"\n Write matrix and rhs in binary to 'matrix.dat' ...\n");
130:     PetscViewerBinaryOpen(PETSC_COMM_SELF,"matrix.dat",PETSC_BINARY_CREATE,&view);
131:     MatView(A,view);
132:     VecView(b,view);
133:     VecView(u,view);
134:     PetscViewerDestroy(view);

136:     VecDestroy(b);
137:     VecDestroy(u);
138:     VecDestroy(u_tmp);
139: 
140:     MatDestroy(A);
141:   }

143:   PetscFinalize();
144:   return 0;
145: }