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:   PetscTruth     flg;

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

 25:   PetscOptionsGetString(PETSC_NULL,"-f",file,PETSC_MAX_PATH_LEN-1,&flg);
 26:   if (!flg) SETERRQ(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(PETSC_NULL,"-mat_type",type,256,PETSC_NULL);

 40:   /*
 41:      Load the matrix and vector; then destroy the viewer.
 42:   */
 43:   MatLoad(fd,type,&C);
 44:   VecLoad(fd,PETSC_NULL,&u);
 45:   PetscViewerDestroy(fd);

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

 50:   MatMultTranspose(C,u,b);

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

 56:   MatLUFactorSymbolic(C,row,col,PETSC_NULL,&A);
 57:   MatLUFactorNumeric(C,PETSC_NULL,&A);
 58:   MatSolveTranspose(A,b,x);

 60:   VecAXPY(x,-1.0,u);
 61:   VecNorm(x,NORM_2,&norm);
 62:   PetscPrintf(PETSC_COMM_SELF,"Norm of error %G\n",norm);

 64:   ISDestroy(row);
 65:   ISDestroy(col);
 66:   VecDestroy(u);
 67:   VecDestroy(x);
 68:   VecDestroy(b);
 69:   MatDestroy(C);
 70:   MatDestroy(A);
 71:   PetscFinalize();
 72:   return 0;
 73: }