Actual source code: ex49.c
2: static char help[] = "Tests SeqSBAIJ factorizations 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,A2;
10: KSP ksp;
11: PetscRandom rctx;
12: PetscReal norm;
13: PetscInt i,j,k,l,n = 27,its,bs = 2,Ii,J;
14: PetscBool test_hermitian = PETSC_FALSE, convert = PETSC_FALSE;
15: PetscScalar v;
17: PetscInitialize(&argc,&args,(char*)0,help);
18: PetscOptionsGetInt(NULL,NULL,"-bs",&bs,NULL);
19: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
20: PetscOptionsGetBool(NULL,NULL,"-herm",&test_hermitian,NULL);
21: PetscOptionsGetBool(NULL,NULL,"-conv",&convert,NULL);
23: MatCreate(PETSC_COMM_SELF,&A);
24: MatSetSizes(A,n*bs,n*bs,PETSC_DETERMINE,PETSC_DETERMINE);
25: MatSetBlockSize(A,bs);
26: MatSetType(A,MATSEQSBAIJ);
27: MatSetFromOptions(A);
28: MatSeqSBAIJSetPreallocation(A,bs,n,NULL);
29: MatSeqBAIJSetPreallocation(A,bs,n,NULL);
30: MatSeqAIJSetPreallocation(A,n*bs,NULL);
31: MatMPIAIJSetPreallocation(A,n*bs,NULL,n*bs,NULL);
33: PetscRandomCreate(PETSC_COMM_SELF,&rctx);
34: for (i=0; i<n; i++) {
35: for (j=i; j<n; j++) {
36: PetscRandomGetValue(rctx,&v);
37: if (PetscRealPart(v) < .1 || i == j) {
38: for (k=0; k<bs; k++) {
39: for (l=0; l<bs; l++) {
40: Ii = i*bs + k;
41: J = j*bs + l;
42: PetscRandomGetValue(rctx,&v);
43: if (Ii == J) v = PetscRealPart(v+3*n*bs);
44: MatSetValue(A,Ii,J,v,INSERT_VALUES);
45: if (test_hermitian) {
46: MatSetValue(A,J,Ii,PetscConj(v),INSERT_VALUES);
47: } else {
48: MatSetValue(A,J,Ii,v,INSERT_VALUES);
49: }
50: }
51: }
52: }
53: }
54: }
55: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
56: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
58: /* With complex numbers:
59: - PETSc cholesky does not support hermitian matrices
60: - CHOLMOD only supports hermitian matrices
61: - SUPERLU_DIST seems supporting both
62: */
63: if (test_hermitian) {
64: MatSetOption(A,MAT_HERMITIAN,PETSC_TRUE);
65: }
67: {
68: Mat M;
69: MatComputeOperator(A,MATAIJ,&M);
70: MatViewFromOptions(M,NULL,"-expl_view");
71: MatDestroy(&M);
72: }
74: A2 = NULL;
75: if (convert) {
76: MatConvert(A,MATAIJ,MAT_INITIAL_MATRIX,&A2);
77: }
79: VecCreate(PETSC_COMM_SELF,&u);
80: VecSetSizes(u,PETSC_DECIDE,n*bs);
81: VecSetFromOptions(u);
82: VecDuplicate(u,&b);
83: VecDuplicate(b,&x);
85: VecSet(u,1.0);
86: MatMult(A,u,b);
88: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
89: Create the linear solver and set various options
90: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
92: /*
93: Create linear solver context
94: */
95: KSPCreate(PETSC_COMM_SELF,&ksp);
97: /*
98: Set operators.
99: */
100: KSPSetOperators(ksp,A2 ? A2 : A,A);
102: KSPSetFromOptions(ksp);
104: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
105: Solve the linear system
106: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
108: KSPSolve(ksp,b,x);
110: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111: Check solution and clean up
112: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
114: /*
115: Check the error
116: */
117: VecAXPY(x,-1.0,u);
118: VecNorm(x,NORM_2,&norm);
119: KSPGetIterationNumber(ksp,&its);
121: /*
122: Print convergence information. PetscPrintf() produces a single
123: print statement from all processes that share a communicator.
124: An alternative is PetscFPrintf(), which prints to a file.
125: */
126: if (norm > 100*PETSC_SMALL) {
127: PetscPrintf(PETSC_COMM_SELF,"Norm of residual %g iterations %D bs %D\n",(double)norm,its,bs);
128: }
130: /*
131: Free work space. All PETSc objects should be destroyed when they
132: are no longer needed.
133: */
134: KSPDestroy(&ksp);
135: VecDestroy(&u);
136: VecDestroy(&x);
137: VecDestroy(&b);
138: MatDestroy(&A);
139: MatDestroy(&A2);
140: PetscRandomDestroy(&rctx);
142: /*
143: Always call PetscFinalize() before exiting a program. This routine
144: - finalizes the PETSc libraries as well as MPI
145: - provides summary and diagnostic information if certain runtime
146: options are chosen (e.g., -log_view).
147: */
148: PetscFinalize();
149: return 0;
150: }
152: /*TEST
154: test:
155: args: -mat_type {{aij baij sbaij}} -bs {{1 2 3 4 5 6 7 8 9 10 11 12}} -pc_type cholesky -herm 0 -conv {{0 1}}
157: test:
158: nsize: {{1 4}}
159: suffix: cholmod
160: requires: suitesparse
161: args: -mat_type {{aij sbaij}} -bs 1 -pc_type cholesky -pc_factor_mat_solver_type cholmod -herm -conv {{0 1}}
163: test:
164: nsize: {{1 4}}
165: suffix: superlu_dist
166: requires: superlu_dist
167: output_file: output/ex49_cholmod.out
168: args: -mat_type mpiaij -bs 3 -pc_type cholesky -pc_factor_mat_solver_type superlu_dist -herm {{0 1}} -conv {{0 1}}
170: test:
171: suffix: mkl_pardiso
172: requires: mkl_pardiso
173: output_file: output/ex49_1.out
174: args: -bs {{1 3}} -pc_type cholesky -pc_factor_mat_solver_type mkl_pardiso
176: test:
177: suffix: cg
178: requires: complex
179: output_file: output/ex49_cg.out
180: args: -herm -ksp_cg_type hermitian -mat_type aij -ksp_type cg -pc_type jacobi -ksp_rtol 4e-07
182: test:
183: suffix: pipecg2
184: requires: complex
185: output_file: output/ex49_pipecg2.out
186: args: -herm -mat_type aij -ksp_type pipecg2 -pc_type jacobi -ksp_rtol 4e-07 -ksp_norm_type {{preconditioned unpreconditioned natural}}
188: TEST*/