2: static char help[] = "Tests the use of MatSolveTranspose().\n\n";

  4:  #include <petscmat.h>

  6: int main(int argc,char **args)
  7: {
  8:   Mat            C,A;
 10:   PetscMPIInt    size;
 11:   IS             row,col;
 12:   Vec            x,u,b;
 13:   PetscReal      norm;
 14:   PetscViewer    fd;
 15:   char           type[256],file[PETSC_MAX_PATH_LEN];
 16:   PetscScalar    mone = -1.0;
 17:   PetscBool      flg;

 19:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
 20:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 21:   if (size > 1) SETERRQ(PETSC_COMM_WORLD,1,"Can only run on one processor");

 23:   PetscOptionsGetString(NULL,NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg);
 24:   if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f option");
 25:   /*
 26:      Open binary file.  Note that we use FILE_MODE_READ to indicate
 27:      reading from this file.
 28:   */
 29:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);

 31:   /*
 32:      Determine matrix format to be used (specified at runtime).
 33:      See the manpage for MatLoad() for available formats.
 34:   */
 35:   PetscStrcpy(type,MATSEQAIJ);
 36:   PetscOptionsGetString(NULL,NULL,"-mat_type",type,256,NULL);

 38:   /*
 39:      Load the matrix and vector; then destroy the viewer.
 40:   */
 41:   MatCreate(PETSC_COMM_WORLD,&C);
 42:   MatSetType(C,type);
 43:   MatLoad(C,fd);
 44:   VecCreate(PETSC_COMM_WORLD,&u);
 45:   VecLoad(u,fd);
 46:   PetscViewerDestroy(&fd);

 48:   VecDuplicate(u,&x);
 49:   VecDuplicate(u,&b);

 51:   MatMultTranspose(C,u,b);

 53:   /* Set default ordering to be Quotient Minimum Degree; also read
 54:      orderings from the options database */
 55:   MatGetOrdering(C,MATORDERINGQMD,&row,&col);

 57:   MatGetFactor(C,MATSOLVERPETSC,MAT_FACTOR_LU,&A);
 58:   MatLUFactorSymbolic(A,C,row,col,NULL);
 59:   MatLUFactorNumeric(A,C,NULL);
 60:   MatSolveTranspose(A,b,x);

 62:   VecAXPY(x,-1.0,u);
 63:   VecNorm(x,NORM_2,&norm);
 64:   PetscPrintf(PETSC_COMM_SELF,"Norm of error %g\n",(double)norm);

 66:   ISDestroy(&row);
 67:   ISDestroy(&col);
 68:   VecDestroy(&u);
 69:   VecDestroy(&x);
 70:   VecDestroy(&b);
 71:   MatDestroy(&C);
 72:   MatDestroy(&A);
 73:   PetscFinalize();
 74:   return ierr;
 75: }