Actual source code: ex50.c
2: static char help[] = "Tests SeqBAIJ point block Jacobi for different block sizes\n\n";
4: #include <petscksp.h>
6: int main(int argc,char **args)
7: {
8: Vec x,b,u;
9: Mat A; /* linear system matrix */
10: KSP ksp; /* linear solver context */
11: PetscRandom rctx; /* random number generator context */
12: PetscReal norm; /* norm of solution error */
13: PetscInt i,j,k,l,n = 27,its,bs = 2,Ii,J;
14: PetscScalar v;
16: PetscInitialize(&argc,&args,(char*)0,help);
17: PetscOptionsGetInt(NULL,NULL,"-bs",&bs,NULL);
18: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
20: MatCreate(PETSC_COMM_WORLD,&A);
21: MatSetSizes(A,n*bs,n*bs,PETSC_DETERMINE,PETSC_DETERMINE);
22: MatSetBlockSize(A,bs);
23: MatSetType(A,MATSEQBAIJ);
24: MatSetFromOptions(A);
25: MatSetUp(A);
27: /*
28: Don't bother to preallocate matrix
29: */
30: PetscRandomCreate(PETSC_COMM_SELF,&rctx);
31: for (i=0; i<n; i++) {
32: for (j=0; j<n; j++) {
33: PetscRandomGetValue(rctx,&v);
34: if (PetscRealPart(v) < .25 || i == j) {
35: for (k=0; k<bs; k++) {
36: for (l=0; l<bs; l++) {
37: PetscRandomGetValue(rctx,&v);
38: Ii = i*bs + k;
39: J = j*bs + l;
40: if (Ii == J) v += 10.;
41: MatSetValue(A,Ii,J,v,INSERT_VALUES);
42: }
43: }
44: }
45: }
46: }
48: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
49: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
51: VecCreate(PETSC_COMM_WORLD,&u);
52: VecSetSizes(u,PETSC_DECIDE,n*bs);
53: VecSetFromOptions(u);
54: VecDuplicate(u,&b);
55: VecDuplicate(b,&x);
57: VecSet(u,1.0);
58: MatMult(A,u,b);
60: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
61: Create the linear solver and set various options
62: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
64: /*
65: Create linear solver context
66: */
67: KSPCreate(PETSC_COMM_WORLD,&ksp);
69: /*
70: Set operators. Here the matrix that defines the linear system
71: also serves as the preconditioning matrix.
72: */
73: KSPSetOperators(ksp,A,A);
75: KSPSetFromOptions(ksp);
77: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
78: Solve the linear system
79: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
81: KSPSolve(ksp,b,x);
83: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
84: Check solution and clean up
85: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
87: /*
88: Check the error
89: */
90: VecAXPY(x,-1.0,u);
91: VecNorm(x,NORM_2,&norm);
92: KSPGetIterationNumber(ksp,&its);
94: /*
95: Print convergence information. PetscPrintf() produces a single
96: print statement from all processes that share a communicator.
97: An alternative is PetscFPrintf(), which prints to a file.
98: */
99: if (norm > .1) {
100: PetscPrintf(PETSC_COMM_WORLD,"Norm of residual %g iterations %D bs %D\n",(double)norm,its,bs);
101: }
103: /*
104: Free work space. All PETSc objects should be destroyed when they
105: are no longer needed.
106: */
107: KSPDestroy(&ksp);
108: VecDestroy(&u);
109: VecDestroy(&x);
110: VecDestroy(&b);
111: MatDestroy(&A);
112: PetscRandomDestroy(&rctx);
114: /*
115: Always call PetscFinalize() before exiting a program. This routine
116: - finalizes the PETSc libraries as well as MPI
117: - provides summary and diagnostic information if certain runtime
118: options are chosen (e.g., -log_view).
119: */
120: PetscFinalize();
121: return 0;
122: }
124: /*TEST
126: test:
127: args: -bs {{1 2 3 4 5 6 7}} -pc_type pbjacobi
129: test:
130: suffix: 2
131: args: -bs {{8 9 10 11 12 13 14 15}} -pc_type ilu
133: TEST*/