Actual source code: ex17.c
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;
9: PetscInt i,j,m = 5,n = 5,Ii,J;
10: PetscScalar v,five = 5.0,one = 1.0;
11: IS isrow,row,col;
12: Vec x,u,b;
13: PetscReal norm;
14: MatFactorInfo info;
16: PetscInitialize(&argc,&args,(char*)0,help);
17: PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
18: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
20: MatCreateSeqAIJ(PETSC_COMM_SELF,m*n,m*n,5,NULL,&C);
21: MatSetUp(C);
23: /* create the matrix for the five point stencil, YET AGAIN*/
24: for (i=0; i<m; i++) {
25: for (j=0; j<n; j++) {
26: v = -1.0; Ii = j + n*i;
27: if (i>0) {J = Ii - n; MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
28: if (i<m-1) {J = Ii + n; MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
29: if (j>0) {J = Ii - 1; MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
30: if (j<n-1) {J = Ii + 1; MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);}
31: v = 4.0; MatSetValues(C,1,&Ii,1,&Ii,&v,INSERT_VALUES);
32: }
33: }
34: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
35: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
37: ISCreateStride(PETSC_COMM_SELF,(m*n)/2,0,2,&isrow);
38: MatZeroRowsIS(C,isrow,five,0,0);
40: VecCreateSeq(PETSC_COMM_SELF,m*n,&u);
41: VecDuplicate(u,&x);
42: VecDuplicate(u,&b);
43: VecSet(u,one);
45: MatMultTranspose(C,u,b);
47: /* Set default ordering to be Quotient Minimum Degree; also read
48: orderings from the options database */
49: MatGetOrdering(C,MATORDERINGQMD,&row,&col);
51: MatFactorInfoInitialize(&info);
52: MatGetFactor(C,MATSOLVERPETSC,MAT_FACTOR_LU,&A);
53: MatLUFactorSymbolic(A,C,row,col,&info);
54: MatLUFactorNumeric(A,C,&info);
55: MatSolveTranspose(A,b,x);
57: ISView(row,PETSC_VIEWER_STDOUT_SELF);
58: VecAXPY(x,-1.0,u);
59: VecNorm(x,NORM_2,&norm);
60: PetscPrintf(PETSC_COMM_SELF,"Norm of error %g\n",(double)norm);
62: ISDestroy(&row);
63: ISDestroy(&col);
64: ISDestroy(&isrow);
65: VecDestroy(&u);
66: VecDestroy(&x);
67: VecDestroy(&b);
68: MatDestroy(&C);
69: MatDestroy(&A);
70: PetscFinalize();
71: return 0;
72: }
74: /*TEST
76: test:
78: TEST*/