Actual source code: ex62.c
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: MatGetFactor(C,MAT_SOLVER_PETSC,MAT_FACTOR_LU,&A);
57: MatLUFactorSymbolic(A,C,row,col,PETSC_NULL);
58: MatLUFactorNumeric(A,C,PETSC_NULL);
59: MatSolveTranspose(A,b,x);
61: VecAXPY(x,-1.0,u);
62: VecNorm(x,NORM_2,&norm);
63: PetscPrintf(PETSC_COMM_SELF,"Norm of error %G\n",norm);
65: ISDestroy(row);
66: ISDestroy(col);
67: VecDestroy(u);
68: VecDestroy(x);
69: VecDestroy(b);
70: MatDestroy(C);
71: MatDestroy(A);
72: PetscFinalize();
73: return 0;
74: }