Actual source code: ex13.c

  2: static char help[] = "Test MatSolveTranspose() for BAIJ matrix  -f <input_file> : file to load \n\n";

  4: /* 
  5:   Include "petscmat.h" so that we can use matrices.
  6:   automatically includes:
  7:      petscsys.h       - base PETSc routines   petscvec.h    - vectors
  8:      petscmat.h    - matrices
  9:      petscis.h     - index sets            petscviewer.h - viewers               
 10: */
 11:  #include petscmat.h
 12:  #include private/matimpl.h

 17: int main(int argc,char **args)
 18: {
 19:   Mat                   A,F;
 20:   PetscViewer           fd;               /* viewer */
 21:   char                  file[PETSC_MAX_PATH_LEN];     /* input file name */
 22:   PetscErrorCode        ierr;
 23:   PetscTruth            flg;
 24:   Vec                   x,y,w;
 25:   MatFactorInfo         iluinfo;
 26:   IS                    perm;
 27:   PetscInt              m;
 28:   PetscReal             norm;

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

 32:   /* 
 33:      Determine file from which we read the matrix

 35:   */
 36:   PetscOptionsGetString(PETSC_NULL,"-f",file,PETSC_MAX_PATH_LEN-1,&flg);
 37:   if (!flg) SETERRQ(1,"Must indicate binary file with the -f option");


 40:   /* 
 41:      Open binary file.  Note that we use FILE_MODE_READ to indicate
 42:      reading from this file.
 43:   */
 44:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);

 46:   /*
 47:     Load the matrix; then destroy the viewer.
 48:   */
 49:   MatLoad(fd,MATSEQBAIJ,&A);
 50:   VecLoad(fd,VECSEQ,&x);
 51:   PetscViewerDestroy(fd);
 52:   VecDuplicate(x,&y);
 53:   VecDuplicate(x,&w);

 55:   MatGetFactor(A,"petsc",MAT_FACTOR_ILU,&F);
 56:   iluinfo.fill = 1.0;
 57:   MatGetSize(A,&m,0);
 58:   ISCreateStride(PETSC_COMM_WORLD,m,0,1,&perm);

 60:   MatLUFactorSymbolic(F,A,perm,perm,&iluinfo);
 61:   MatLUFactorNumeric(F,A,&iluinfo);
 62:   MatSolveTranspose(F,x,y);
 63:   F->ops->solvetranspose = MatSolveTranspose_SeqBAIJ_N;
 64:   MatSolveTranspose(F,x,w);
 65:   //  VecView(w,0);VecView(y,0);
 66:   VecAXPY(w,-1.0,y);
 67:   VecNorm(w,NORM_2,&norm);
 68:   if (norm) {
 69:     PetscPrintf(PETSC_COMM_SELF,"Norm of difference is nonzero %g\n",norm);
 70:   }
 71:   ISDestroy(perm);
 72:   MatDestroy(A);
 73:   MatDestroy(F);
 74:   VecDestroy(x);
 75:   VecDestroy(y);
 76:   VecDestroy(w);

 78:   PetscFinalize();
 79:   return 0;
 80: }