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

  4: #include <petscmat.h>

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

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

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

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

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

 50:   VecDuplicate(u,&x);
 51:   VecDuplicate(u,&b);

 53:   MatMultTranspose(C,u,b);

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

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

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

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