Actual source code: ex77.c

  1: #include <petsc.h>

  3: static char help[] = "Solves a linear system with a block of right-hand sides using KSPHPDDM.\n\n";

  5: int main(int argc, char **args)
  6: {
  7:   Mat                X, B;         /* computed solutions and RHS */
  8:   Vec                cx, cb;       /* columns of X and B */
  9:   Mat                A, KA = NULL; /* linear system matrix */
 10:   KSP                ksp;          /* linear solver context */
 11:   PC                 pc;           /* preconditioner context */
 12:   Mat                F;            /* factored matrix from the preconditioner context */
 13:   PetscScalar       *x, *S = NULL, *T = NULL;
 14:   PetscReal          norm, deflation = -1.0;
 15:   PetscInt           m, M, N = 5, i;
 16:   PetscMPIInt        rank, size;
 17:   const char        *deft = MATAIJ;
 18:   PetscViewer        viewer;
 19:   char               name[PETSC_MAX_PATH_LEN], type[256];
 20:   PetscBool          breakdown = PETSC_FALSE, flg;
 21:   KSPConvergedReason reason;

 24:   PetscInitialize(&argc, &args, NULL, help);
 25:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
 26:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 27:   PetscOptionsGetString(NULL, NULL, "-f", name, sizeof(name), &flg);
 29:   PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL);
 30:   PetscOptionsGetBool(NULL, NULL, "-breakdown", &breakdown, NULL);
 31:   PetscOptionsGetReal(NULL, NULL, "-ksp_hpddm_deflation_tol", &deflation, NULL);
 32:   MatCreate(PETSC_COMM_WORLD, &A);
 33:   KSPCreate(PETSC_COMM_WORLD, &ksp);
 34:   KSPSetOperators(ksp, A, A);
 35:   PetscViewerBinaryOpen(PETSC_COMM_WORLD, name, FILE_MODE_READ, &viewer);
 36:   MatLoad(A, viewer);
 37:   PetscViewerDestroy(&viewer);
 38:   PetscOptionsBegin(PETSC_COMM_WORLD, "", "", "");
 39:   PetscOptionsFList("-mat_type", "Matrix type", "MatSetType", MatList, deft, type, 256, &flg);
 40:   PetscOptionsEnd();
 41:   if (flg) {
 42:     PetscStrcmp(type, MATKAIJ, &flg);
 43:     if (!flg) {
 44:       MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE);
 45:       MatConvert(A, type, MAT_INPLACE_MATRIX, &A);
 46:     } else {
 47:       if (size > 2) {
 48:         MatGetSize(A, &M, NULL);
 49:         MatCreate(PETSC_COMM_WORLD, &B);
 50:         if (rank > 1) MatSetSizes(B, 0, 0, M, M);
 51:         else MatSetSizes(B, rank ? M - M / 2 : M / 2, rank ? M - M / 2 : M / 2, M, M);
 52:         PetscViewerBinaryOpen(PETSC_COMM_WORLD, name, FILE_MODE_READ, &viewer);
 53:         MatLoad(B, viewer);
 54:         PetscViewerDestroy(&viewer);
 55:         MatHeaderReplace(A, &B);
 56:       }
 57:       PetscCalloc2(N * N, &S, N * N, &T);
 58:       for (i = 0; i < N; i++) { /* really easy problem used for testing */
 59:         S[i * (N + 1)] = 1e+6;
 60:         T[i * (N + 1)] = 1e-2;
 61:       }
 62:       MatCreateKAIJ(A, N, N, S, T, &KA);
 63:     }
 64:   }
 65:   if (!flg) {
 66:     if (size > 4) {
 67:       Mat B;
 68:       MatGetSize(A, &M, NULL);
 69:       MatCreate(PETSC_COMM_WORLD, &B);
 70:       if (rank > 3) MatSetSizes(B, 0, 0, M, M);
 71:       else MatSetSizes(B, rank == 0 ? M - 3 * (M / 4) : M / 4, rank == 0 ? M - 3 * (M / 4) : M / 4, M, M);
 72:       PetscViewerBinaryOpen(PETSC_COMM_WORLD, name, FILE_MODE_READ, &viewer);
 73:       MatLoad(B, viewer);
 74:       PetscViewerDestroy(&viewer);
 75:       MatHeaderReplace(A, &B);
 76:     }
 77:   }
 78:   MatGetLocalSize(A, &m, NULL);
 79:   MatCreateDense(PETSC_COMM_WORLD, m, PETSC_DECIDE, PETSC_DECIDE, N, NULL, &B);
 80:   MatCreateDense(PETSC_COMM_WORLD, m, PETSC_DECIDE, PETSC_DECIDE, N, NULL, &X);
 81:   if (!breakdown) MatSetRandom(B, NULL);
 82:   KSPSetFromOptions(ksp);
 83:   if (!flg) {
 84:     if (!breakdown) {
 85:       KSPMatSolve(ksp, B, X);
 86:       KSPGetMatSolveBatchSize(ksp, &M);
 87:       if (M != PETSC_DECIDE) {
 88:         KSPSetMatSolveBatchSize(ksp, PETSC_DECIDE);
 89:         MatZeroEntries(X);
 90:         KSPMatSolve(ksp, B, X);
 91:       }
 92:       KSPGetPC(ksp, &pc);
 93:       PetscObjectTypeCompare((PetscObject)pc, PCLU, &flg);
 94:       if (flg) {
 95:         PCFactorGetMatrix(pc, &F);
 96:         MatMatSolve(F, B, B);
 97:         MatAYPX(B, -1.0, X, SAME_NONZERO_PATTERN);
 98:         MatNorm(B, NORM_INFINITY, &norm);
100:       }
101:     } else {
102:       MatZeroEntries(B);
103:       KSPMatSolve(ksp, B, X);
104:       KSPGetConvergedReason(ksp, &reason);
106:       MatDenseGetArrayWrite(B, &x);
107:       for (i = 0; i < m * N; ++i) x[i] = 1.0;
108:       MatDenseRestoreArrayWrite(B, &x);
109:       KSPMatSolve(ksp, B, X);
110:       KSPGetConvergedReason(ksp, &reason);
112:     }
113:   } else {
114:     KSPSetOperators(ksp, KA, KA);
115:     MatGetSize(KA, &M, NULL);
116:     VecCreateMPI(PETSC_COMM_WORLD, m * N, M, &cb);
117:     VecCreateMPI(PETSC_COMM_WORLD, m * N, M, &cx);
118:     VecSetRandom(cb, NULL);
119:     /* solving with MatKAIJ is equivalent to block solving with row-major RHS and solutions */
120:     /* only applies if MatKAIJGetScaledIdentity() returns true                              */
121:     KSPSolve(ksp, cb, cx);
122:     VecDestroy(&cx);
123:     VecDestroy(&cb);
124:   }
125:   MatDestroy(&X);
126:   MatDestroy(&B);
127:   PetscFree2(S, T);
128:   MatDestroy(&KA);
129:   MatDestroy(&A);
130:   KSPDestroy(&ksp);
131:   PetscFinalize();
132:   return 0;
133: }

135: /*TEST

137:    testset:
138:       nsize: 2
139:       requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
140:       args: -ksp_converged_reason -ksp_max_it 500 -f ${DATAFILESPATH}/matrices/hpddm/GCRODR/A_400.dat -mat_type {{aij sbaij}shared output}
141:       test:
142:          suffix: 1
143:          args:
144:       test:
145:          suffix: 2
146:          requires: hpddm
147:          args: -ksp_type hpddm -pc_type asm -ksp_hpddm_type {{gmres bgmres}separate output}
148:       test:
149:          suffix: 3
150:          requires: hpddm
151:          args: -ksp_type hpddm -ksp_hpddm_recycle 5 -ksp_hpddm_type {{gcrodr bgcrodr}separate output}
152:       test:
153:          nsize: 4
154:          suffix: 4
155:          requires: hpddm
156:          args: -ksp_rtol 1e-4 -ksp_type hpddm -ksp_hpddm_recycle 5 -ksp_hpddm_type bgcrodr -ksp_view_final_residual -N 12 -ksp_matsolve_batch_size 5

158:    test:
159:       nsize: 1
160:       requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
161:       suffix: preonly
162:       args: -N 6 -f ${DATAFILESPATH}/matrices/hpddm/GCRODR/A_400.dat -pc_type lu -ksp_type hpddm -ksp_hpddm_type preonly

164:    testset:
165:       nsize: 1
166:       requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
167:       args: -N 3 -f ${DATAFILESPATH}/matrices/hpddm/GCRODR/A_400.dat -ksp_type hpddm -breakdown
168:       test:
169:          suffix: breakdown_wo_deflation
170:          output_file: output/ex77_preonly.out
171:          args: -pc_type none -ksp_hpddm_type {{bcg bgmres bgcrodr bfbcg}shared output}
172:       test:
173:          suffix: breakdown_w_deflation
174:          output_file: output/ex77_deflation.out
175:          filter: sed "s/BGCRODR/BGMRES/g"
176:          args: -pc_type lu -ksp_hpddm_type {{bgmres bgcrodr}shared output} -ksp_hpddm_deflation_tol 1e-1 -info :ksp

178:    test:
179:       nsize: 2
180:       requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
181:       args: -N 12 -ksp_converged_reason -ksp_max_it 500 -f ${DATAFILESPATH}/matrices/hpddm/GCRODR/A_400.dat -mat_type kaij -pc_type pbjacobi -ksp_type hpddm -ksp_hpddm_type {{gmres bgmres}separate output}

183:    test:
184:       nsize: 3
185:       requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
186:       suffix: kaij_zero
187:       output_file: output/ex77_ksp_hpddm_type-bgmres.out
188:       args: -N 12 -ksp_converged_reason -ksp_max_it 500 -f ${DATAFILESPATH}/matrices/hpddm/GCRODR/A_400.dat -mat_type kaij -pc_type pbjacobi -ksp_type hpddm -ksp_hpddm_type bgmres

190:    test:
191:       nsize: 4
192:       requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) slepc defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
193:       suffix: 4_slepc
194:       output_file: output/ex77_4.out
195:       filter: sed "/^ksp_hpddm_recycle_ Linear eigensolve converged/d"
196:       args: -ksp_converged_reason -ksp_max_it 500 -f ${DATAFILESPATH}/matrices/hpddm/GCRODR/A_400.dat -ksp_rtol 1e-4 -ksp_type hpddm -ksp_hpddm_recycle 5 -ksp_hpddm_type bgcrodr -ksp_view_final_residual -N 12 -ksp_matsolve_batch_size 5 -ksp_hpddm_recycle_redistribute 2 -ksp_hpddm_recycle_mat_type {{aij dense}shared output} -ksp_hpddm_recycle_eps_converged_reason -ksp_hpddm_recycle_st_pc_type redundant

198:    testset:
199:       nsize: 4
200:       requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) slepc defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
201:       filter: sed "/^ksp_hpddm_recycle_ Linear eigensolve converged/d"
202:       args: -ksp_converged_reason -ksp_max_it 500 -f ${DATAFILESPATH}/matrices/hpddm/GCRODR/A_400.dat -ksp_rtol 1e-4 -ksp_type hpddm -ksp_hpddm_recycle 5 -ksp_hpddm_type bgcrodr -ksp_view_final_residual -N 12 -ksp_matsolve_batch_size 5 -ksp_hpddm_recycle_redistribute 2 -ksp_hpddm_recycle_eps_converged_reason
203:       test:
204:          requires: elemental
205:          suffix: 4_elemental
206:          output_file: output/ex77_4.out
207:          args: -ksp_hpddm_recycle_mat_type elemental
208:       test:
209:          requires: elemental
210:          suffix: 4_elemental_matmat
211:          output_file: output/ex77_4.out
212:          args: -ksp_hpddm_recycle_mat_type elemental -ksp_hpddm_recycle_eps_type subspace -ksp_hpddm_recycle_eps_target 0 -ksp_hpddm_recycle_eps_target_magnitude -ksp_hpddm_recycle_st_type sinvert -ksp_hpddm_recycle_bv_type mat -ksp_hpddm_recycle_bv_orthog_block svqb
213:       test:
214:          requires: scalapack
215:          suffix: 4_scalapack
216:          output_file: output/ex77_4.out
217:          args: -ksp_hpddm_recycle_mat_type scalapack

219:    test:
220:       nsize: 5
221:       requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
222:       suffix: 4_zero
223:       output_file: output/ex77_4.out
224:       args: -ksp_converged_reason -ksp_max_it 500 -f ${DATAFILESPATH}/matrices/hpddm/GCRODR/A_400.dat -ksp_rtol 1e-4 -ksp_type hpddm -ksp_hpddm_recycle 5 -ksp_hpddm_type bgcrodr -ksp_view_final_residual -N 12 -ksp_matsolve_batch_size 5

226: TEST*/