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;

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

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

 34:   PetscRandomCreate(PETSC_COMM_SELF, &rctx);
 35:   for (i = 0; i < n; i++) {
 36:     for (j = i; j < n; j++) {
 37:       PetscRandomGetValue(rctx, &v);
 38:       if (PetscRealPart(v) < .1 || i == j) {
 39:         for (k = 0; k < bs; k++) {
 40:           for (l = 0; l < bs; l++) {
 41:             Ii = i * bs + k;
 42:             J  = j * bs + l;
 43:             PetscRandomGetValue(rctx, &v);
 44:             if (Ii == J) v = PetscRealPart(v + 3 * n * bs);
 45:             MatSetValue(A, Ii, J, v, INSERT_VALUES);
 46:             if (test_hermitian) {
 47:               MatSetValue(A, J, Ii, PetscConj(v), INSERT_VALUES);
 48:             } else {
 49:               MatSetValue(A, J, Ii, v, INSERT_VALUES);
 50:             }
 51:           }
 52:         }
 53:       }
 54:     }
 55:   }
 56:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 57:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

 59:   /* With complex numbers:
 60:      - PETSc cholesky does not support hermitian matrices
 61:      - CHOLMOD only supports hermitian matrices
 62:      - SUPERLU_DIST seems supporting both
 63:   */
 64:   if (test_hermitian) MatSetOption(A, MAT_HERMITIAN, PETSC_TRUE);

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

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

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

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

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

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

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

 99:   KSPSetFromOptions(ksp);

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

105:   KSPSolve(ksp, b, x);

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

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

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

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

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

147: /*TEST

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

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

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

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

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

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

183: TEST*/