Actual source code: ex72.c

petsc-3.8.4 2018-03-24
Report Typos and Errors

  2:  #include <petscmat.h>

  4: static char help[] = "Read in a Symmetric matrix in MatrixMarket format (only the lower triangle). \n\
  5:   Assemble it to a PETSc sparse SBAIJ (upper triangle) matrix. \n\
  6:   Write it in a AIJ matrix (entire matrix) to a file. \n\
  7:   Input parameters are:            \n\
  8:     -fin <filename> : input file   \n\
  9:     -fout <filename> : output file \n\n";

 11: int main(int argc,char **args)
 12: {
 13:   Mat            A;
 14:   char           filein[PETSC_MAX_PATH_LEN],fileout[PETSC_MAX_PATH_LEN],buf[PETSC_MAX_PATH_LEN];
 15:   PetscInt       i,m,n,nnz;
 17:   PetscMPIInt    size;
 18:   PetscScalar    *val,zero=0.0;
 19:   FILE           *file;
 20:   PetscViewer    view;
 21:   int            *row,*col,*rownz;
 22:   PetscBool      flg;

 24:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
 25: #if defined(PETSC_USE_COMPLEX)
 26:   SETERRQ(PETSC_COMM_WORLD,1,"This example does not work with complex numbers");
 27:   PetscFinalize();
 28:   return 0;
 29: #endif

 31:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 32:   if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Uniprocessor Example only\n");

 34:   /* Read in matrix and RHS */
 35:   PetscOptionsGetString(NULL,NULL,"-fin",filein,PETSC_MAX_PATH_LEN,&flg);
 36:   if (!flg) SETERRQ(PETSC_COMM_SELF,1,"Must indicate input file with -fin option");
 37:   PetscFOpen(PETSC_COMM_SELF,filein,"r",&file);

 39:   /* process header with comments */
 40:   do {
 41:     char *str = fgets(buf,PETSC_MAX_PATH_LEN-1,file);
 42:     if (!str) SETERRQ(PETSC_COMM_SELF,1,"Incorrect format in file");
 43:   }while (buf[0] == '%');

 45:   /* The first non-comment line has the matrix dimensions */
 46:   sscanf(buf,"%d %d %d\n",&m,&n,&nnz);
 47:   PetscPrintf (PETSC_COMM_SELF,"m = %d, n = %d, nnz = %d\n",m,n,nnz);

 49:   /* reseve memory for matrices */
 50:   PetscMalloc4(nnz,&row,nnz,&col,nnz,&val,m,&rownz);
 51:   for (i=0; i<m; i++) rownz[i] = 1; /* add 0.0 to diagonal entries */

 53:   for (i=0; i<nnz; i++) {
 54:     fscanf(file,"%d %d %le\n",&row[i],&col[i],(double*)&val[i]);
 55:     if (ierr == EOF) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"i=%d, reach EOF\n",i);
 56:     row[i]--; col[i]--;    /* adjust from 1-based to 0-based */
 57:     rownz[col[i]]++;
 58:   }
 59:   fclose(file);
 60:   PetscPrintf(PETSC_COMM_SELF,"Read file completes.\n");

 62:   /* Creat and asseble SBAIJ matrix */
 63:   MatCreate(PETSC_COMM_SELF,&A);
 64:   MatSetType(A,MATSBAIJ);
 65:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,n);
 66:   MatSetFromOptions(A);
 67:   MatSeqSBAIJSetPreallocation(A,1,0,rownz);

 69:   /* Add zero to diagonals, in case the matrix missing diagonals */
 70:   for (i=0; i<m; i++){
 71:     MatSetValues(A,1,&i,1,&i,&zero,INSERT_VALUES);
 72:   }
 73:   for (i=0; i<nnz; i++) {
 74:     MatSetValues(A,1,&col[i],1,&row[i],&val[i],INSERT_VALUES);
 75:   }
 76:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 77:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 78:   PetscPrintf(PETSC_COMM_SELF,"Assemble SBAIJ matrix completes.\n");

 80:   /* Write the entire matrix in AIJ format to a file */
 81:   PetscOptionsGetString(NULL,NULL,"-fout",fileout,PETSC_MAX_PATH_LEN,&flg);
 82:   if (flg) {
 83:     PetscPrintf(PETSC_COMM_SELF,"Write the entire matrix in AIJ format to file %s\n",fileout);
 84:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,fileout,FILE_MODE_WRITE,&view);
 85:     MatView(A,view);
 86:     PetscViewerDestroy(&view);
 87:   }

 89:   PetscFree4(row,col,val,rownz);
 90:   MatDestroy(&A);
 91:   PetscFinalize();
 92:   return ierr;
 93: }