Actual source code: ex81.c

  2: static char help[] = "Reads in a PETSc binary matrix and saves in Harwell-Boeing format.\n\
  3:   -fout <output_file> : file to load.\n\
  4:   -fin <input_file> : For a 5X5 example of the 5-pt. stencil,\n\
  5:                        use the file petsc/src/mat/examples/matbinary.ex\n\n";

  7: /*
  8:   Include the private file (not included by most applications) so we have direct
  9:   access to the matrix data structure.

 11:   This code is buggy! What is it doing here?
 12: */
 13: #include <../src/mat/impls/aij/seq/aij.h>

 17: int main(int argc,char **args)
 18: {
 20:   PetscInt       n,m,i,*ai,*aj,nz;
 21:   PetscMPIInt    size;
 22:   Mat            A;
 23:   Vec            x;
 24:   char           bfile[PETSC_MAX_PATH_LEN],hbfile[PETSC_MAX_PATH_LEN];
 25:   PetscViewer    fd;
 26:   Mat_SeqAIJ     *a;
 27:   PetscScalar    *aa,*xx;
 28:   FILE           *file;
 29:   char           head[81];

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

 33: #if defined(PETSC_USE_COMPLEX)
 34:   SETERRQ(PETSC_COMM_WORLD,1,"This example does not work with complex numbers");
 35: #endif
 36:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 37:   if (size > 1) SETERRQ(PETSC_COMM_WORLD,1,"Only runs on one processor");

 39:   PetscOptionsGetString(PETSC_NULL,"-fin",bfile,PETSC_MAX_PATH_LEN,PETSC_NULL);
 40:   PetscOptionsGetString(PETSC_NULL,"-fout",hbfile,PETSC_MAX_PATH_LEN,PETSC_NULL);

 42:   /* Read matrix and RHS */
 43:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,bfile,FILE_MODE_READ,&fd);
 44:   MatCreate(PETSC_COMM_WORLD,&A);
 45:   MatSetType(A,MATSEQAIJ);
 46:   MatLoad(A,fd);
 47:   VecCreate(PETSC_COMM_WORLD,&x);
 48:   VecLoad(x,fd);
 49:   PetscViewerDestroy(&fd);

 51:   /* Format is in column storage so we print transpose matrix */
 52:   MatTranspose(A,MAT_REUSE_MATRIX,&A);

 54:   m = A->rmap->n;
 55:   n = A->cmap->n;
 56:   if (n != m) SETERRQ(PETSC_COMM_SELF,1,"Only for square matrices");

 58:   /* charrage returns \n may not belong below
 59:     depends on what 80 character fixed format means to Fortran */

 61:   file = fopen(hbfile,"w"); if (!file) SETERRQ(PETSC_COMM_SELF,1,"Cannot open HB file");
 62:   sprintf(head,"%-72s%-8s\n","Title","Key");
 63:   fprintf(file,head);
 64:   a  = (Mat_SeqAIJ*)A->data;
 65:   aa = a->a;
 66:   ai = a->i;
 67:   aj = a->j;
 68:   nz = a->nz;


 71:   sprintf(head,"%14d%14d%14d%14d%14d%10s\n",3*m+1,m+1,nz,nz," ");
 72:   fprintf(file,head);
 73:   sprintf(head,"RUA%14d%14d%14d%14d%13s\n",m,m,nz," ");
 74:   fprintf(file,head);

 76:   fprintf(file,"Formats I don't know\n");

 78:   for (i=0; i<m+1; i++) {
 79:     fprintf(file,"%10d%70s\n",ai[i]," ");
 80:   }
 81:   for (i=0; i<nz; i++) {
 82:     fprintf(file,"%10d%70s\n",aj[i]," ");
 83:   }

 85:   for (i=0; i<nz; i++) {
 86:     fprintf(file,"%16.14e,%64s\n",aa[i]," ");
 87:   }

 89:   /* print the vector to the file */
 90:   VecGetArray(x,&xx);
 91:   for (i=0; i<m; i++) {
 92:     fprintf(file,"%16.14e%64s\n",xx[i]," ");
 93:   }
 94:   VecRestoreArray(x,&xx);

 96:   fclose(file);
 97:   MatDestroy(&A);
 98:   VecDestroy(&x);

100:   PetscFinalize();
101:   return 0;
102: }