Actual source code: ex49.c

  1: static char help[] = "Tests SeqSBAIJ factorizations for different block sizes\n\n";

  3: #include <petscksp.h>

  5: int main(int argc, char **args)
  6: {
  7:   Vec         x, b, u;
  8:   Mat         A, A2;
  9:   KSP         ksp;
 10:   PetscRandom rctx;
 11:   PetscReal   norm;
 12:   PetscInt    i, j, k, l, n = 27, its, bs = 2, Ii, J;
 13:   PetscBool   test_hermitian = PETSC_FALSE, convert = PETSC_FALSE;
 14:   PetscScalar v;

 16:   PetscFunctionBeginUser;
 17:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
 18:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-bs", &bs, NULL));
 19:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
 20:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-herm", &test_hermitian, NULL));
 21:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-conv", &convert, NULL));

 23:   PetscCall(MatCreate(PETSC_COMM_SELF, &A));
 24:   PetscCall(MatSetSizes(A, n * bs, n * bs, PETSC_DETERMINE, PETSC_DETERMINE));
 25:   PetscCall(MatSetBlockSize(A, bs));
 26:   PetscCall(MatSetType(A, MATSEQSBAIJ));
 27:   PetscCall(MatSetFromOptions(A));
 28:   PetscCall(MatSeqSBAIJSetPreallocation(A, bs, n, NULL));
 29:   PetscCall(MatSeqBAIJSetPreallocation(A, bs, n, NULL));
 30:   PetscCall(MatSeqAIJSetPreallocation(A, n * bs, NULL));
 31:   PetscCall(MatMPIAIJSetPreallocation(A, n * bs, NULL, n * bs, NULL));

 33:   PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &rctx));
 34:   for (i = 0; i < n; i++) {
 35:     for (j = i; j < n; j++) {
 36:       PetscCall(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:             PetscCall(PetscRandomGetValue(rctx, &v));
 43:             if (Ii == J) v = PetscRealPart(v + 3 * n * bs);
 44:             PetscCall(MatSetValue(A, Ii, J, v, INSERT_VALUES));
 45:             if (test_hermitian) {
 46:               PetscCall(MatSetValue(A, J, Ii, PetscConj(v), INSERT_VALUES));
 47:             } else {
 48:               PetscCall(MatSetValue(A, J, Ii, v, INSERT_VALUES));
 49:             }
 50:           }
 51:         }
 52:       }
 53:     }
 54:   }
 55:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 56:   PetscCall(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) PetscCall(MatSetOption(A, MAT_HERMITIAN, PETSC_TRUE));

 65:   {
 66:     Mat M;
 67:     PetscCall(MatComputeOperator(A, MATAIJ, &M));
 68:     PetscCall(MatViewFromOptions(M, NULL, "-expl_view"));
 69:     PetscCall(MatDestroy(&M));
 70:   }

 72:   A2 = NULL;
 73:   if (convert) PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &A2));

 75:   PetscCall(VecCreate(PETSC_COMM_SELF, &u));
 76:   PetscCall(VecSetSizes(u, PETSC_DECIDE, n * bs));
 77:   PetscCall(VecSetFromOptions(u));
 78:   PetscCall(VecDuplicate(u, &b));
 79:   PetscCall(VecDuplicate(b, &x));

 81:   PetscCall(VecSet(u, 1.0));
 82:   PetscCall(MatMult(A, u, b));

 84:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 85:                 Create the linear solver and set various options
 86:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 88:   /*
 89:      Create linear solver context
 90:   */
 91:   PetscCall(KSPCreate(PETSC_COMM_SELF, &ksp));

 93:   /*
 94:      Set operators.
 95:   */
 96:   PetscCall(KSPSetOperators(ksp, A2 ? A2 : A, A));

 98:   PetscCall(KSPSetFromOptions(ksp));

100:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101:                       Solve the linear system
102:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

104:   PetscCall(KSPSolve(ksp, b, x));

106:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107:                       Check solution and clean up
108:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

110:   /*
111:      Check the error
112:   */
113:   PetscCall(VecAXPY(x, -1.0, u));
114:   PetscCall(VecNorm(x, NORM_2, &norm));
115:   PetscCall(KSPGetIterationNumber(ksp, &its));

117:   /*
118:      Print convergence information.  PetscPrintf() produces a single
119:      print statement from all processes that share a communicator.
120:      An alternative is PetscFPrintf(), which prints to a file.
121:   */
122:   if (norm > 100 * PETSC_SMALL) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Norm of residual %g iterations %" PetscInt_FMT " bs %" PetscInt_FMT "\n", (double)norm, its, bs));

124:   /*
125:      Free work space.  All PETSc objects should be destroyed when they
126:      are no longer needed.
127:   */
128:   PetscCall(KSPDestroy(&ksp));
129:   PetscCall(VecDestroy(&u));
130:   PetscCall(VecDestroy(&x));
131:   PetscCall(VecDestroy(&b));
132:   PetscCall(MatDestroy(&A));
133:   PetscCall(MatDestroy(&A2));
134:   PetscCall(PetscRandomDestroy(&rctx));

136:   /*
137:      Always call PetscFinalize() before exiting a program.  This routine
138:        - finalizes the PETSc libraries as well as MPI
139:        - provides summary and diagnostic information if certain runtime
140:          options are chosen (e.g., -log_view).
141:   */
142:   PetscCall(PetscFinalize());
143:   return 0;
144: }

146: /*TEST

148:    test:
149:       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}}

151:    test:
152:       nsize: {{1 4}}
153:       suffix: cholmod
154:       requires: suitesparse
155:       args: -mat_type {{aij sbaij}} -bs 1 -pc_type cholesky -pc_factor_mat_solver_type cholmod -herm -conv {{0 1}}

157:    test:
158:       nsize: {{1 4}}
159:       suffix: superlu_dist
160:       requires: superlu_dist
161:       output_file: output/ex49_cholmod.out
162:       args: -mat_type mpiaij -bs 3 -pc_type cholesky -pc_factor_mat_solver_type superlu_dist -herm {{0 1}} -conv {{0 1}}

164:    test:
165:       suffix: mkl_pardiso
166:       requires: mkl_pardiso
167:       output_file: output/ex49_1.out
168:       args: -bs {{1 3}} -pc_type cholesky -pc_factor_mat_solver_type mkl_pardiso

170:    test:
171:       suffix: cg
172:       requires: complex
173:       output_file: output/ex49_cg.out
174:       args: -herm -ksp_cg_type hermitian -mat_type aij -ksp_type cg -pc_type jacobi -ksp_rtol 4e-07

176:    test:
177:       suffix: pipecg2
178:       requires: complex
179:       output_file: output/ex49_pipecg2.out
180:       args: -herm -mat_type aij -ksp_type pipecg2 -pc_type jacobi -ksp_rtol 4e-07 -ksp_norm_type {{preconditioned unpreconditioned natural}}

182: TEST*/