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