2: static char help[] = "Test MatFwkAIJ: a block matrix with an AIJ-like datastructure keeping track of nonzero blocks.\n\
  3: Each block is a matrix of (generally) any type.\n\n";

  5: /* 
  6:   Include "petscmat.h" so that we can use matrices.
  7:   automatically includes:
  8:      petscsys.h       - base PETSc routines   petscvec.h    - vectors
  9:      petscmat.h    - matrices
 10:      petscis.h     - index sets            petscviewer.h - viewers               
 11: */
 12: #include <petscmat.h>
 13: #include <petsc-private/matimpl.h>
 14: extern PetscErrorCode MatSolveTranspose_SeqBAIJ_N(Mat,Vec,Vec);

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

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

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

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


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

 47:   /*
 48:     Load the matrix; then destroy the viewer.
 49:   */
 50:   MatCreate(PETSC_COMM_WORLD,&A);
 51:   MatSetType(A,MATSEQBAIJ);
 52:   MatLoad(A,fd);
 53:   VecCreate(PETSC_COMM_WORLD,&x);
 54:   VecLoad(x,fd);
 55:   PetscViewerDestroy(&fd);
 56:   VecDuplicate(x,&y);
 57:   VecDuplicate(x,&w);

 59:   MatGetFactor(A,"petsc",MAT_FACTOR_ILU,&F);
 60:   iluinfo.fill = 1.0;
 61:   MatGetSize(A,&m,0);
 62:   ISCreateStride(PETSC_COMM_WORLD,m,0,1,&perm);

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

 82:   PetscFinalize();
 83:   return 0;
 84: }