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: }