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;

     23:   PetscFunctionBeginUser;
     24:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
     25:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
     26:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
     27:   PetscCall(PetscOptionsGetString(NULL, NULL, "-f", name, sizeof(name), &flg));
     28:   PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_SUP, "Must provide a binary file for the matrix");
     29:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL));
     30:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-breakdown", &breakdown, NULL));
     31:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-ksp_hpddm_deflation_tol", &deflation, NULL));
     32:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
     33:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
     34:   PetscCall(KSPSetOperators(ksp, A, A));
     35:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, name, FILE_MODE_READ, &viewer));
     36:   PetscCall(MatLoad(A, viewer));
     37:   PetscCall(PetscViewerDestroy(&viewer));
     38:   PetscOptionsBegin(PETSC_COMM_WORLD, "", "", "");
     39:   PetscCall(PetscOptionsFList("-mat_type", "Matrix type", "MatSetType", MatList, deft, type, 256, &flg));
     40:   PetscOptionsEnd();
     41:   if (flg) {
     42:     PetscCall(PetscStrcmp(type, MATKAIJ, &flg));
     43:     if (!flg) {
     44:       PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));
     45:       PetscCall(MatConvert(A, type, MAT_INPLACE_MATRIX, &A));
     46:     } else {
     47:       if (size > 2) {
     48:         PetscCall(MatGetSize(A, &M, NULL));
     49:         PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
     50:         if (rank > 1) PetscCall(MatSetSizes(B, 0, 0, M, M));
     51:         else PetscCall(MatSetSizes(B, rank ? M - M / 2 : M / 2, rank ? M - M / 2 : M / 2, M, M));
     52:         PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, name, FILE_MODE_READ, &viewer));
     53:         PetscCall(MatLoad(B, viewer));
     54:         PetscCall(PetscViewerDestroy(&viewer));
     55:         PetscCall(MatHeaderReplace(A, &B));
     56:       }
     57:       PetscCall(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:       PetscCall(MatCreateKAIJ(A, N, N, S, T, &KA));
     63:     }
     64:   }
     65:   if (!flg) {
     66:     if (size > 4) {
     67:       Mat B;
     68:       PetscCall(MatGetSize(A, &M, NULL));
     69:       PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
     70:       if (rank > 3) PetscCall(MatSetSizes(B, 0, 0, M, M));
     71:       else PetscCall(MatSetSizes(B, rank == 0 ? M - 3 * (M / 4) : M / 4, rank == 0 ? M - 3 * (M / 4) : M / 4, M, M));
     72:       PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, name, FILE_MODE_READ, &viewer));
     73:       PetscCall(MatLoad(B, viewer));
     74:       PetscCall(PetscViewerDestroy(&viewer));
     75:       PetscCall(MatHeaderReplace(A, &B));
     76:     }
     77:   }
     78:   PetscCall(MatGetLocalSize(A, &m, NULL));
     79:   PetscCall(MatCreateDense(PETSC_COMM_WORLD, m, PETSC_DECIDE, PETSC_DECIDE, N, NULL, &B));
     80:   PetscCall(MatCreateDense(PETSC_COMM_WORLD, m, PETSC_DECIDE, PETSC_DECIDE, N, NULL, &X));
     81:   if (!breakdown) PetscCall(MatSetRandom(B, NULL));
     82:   PetscCall(KSPSetFromOptions(ksp));
     83:   if (!flg) {
     84:     if (!breakdown) {
     85:       PetscCall(KSPMatSolve(ksp, B, X));
     86:       PetscCall(KSPGetMatSolveBatchSize(ksp, &M));
     87:       if (M != PETSC_DECIDE) {
     88:         PetscCall(KSPSetMatSolveBatchSize(ksp, PETSC_DECIDE));
     89:         PetscCall(MatZeroEntries(X));
     90:         PetscCall(KSPMatSolve(ksp, B, X));
     91:       }
     92:       PetscCall(KSPGetPC(ksp, &pc));
     93:       PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLU, &flg));
     94:       if (flg) {
     95:         PetscCall(PCFactorGetMatrix(pc, &F));
     96:         PetscCall(MatMatSolve(F, B, B));
     97:         PetscCall(MatAYPX(B, -1.0, X, SAME_NONZERO_PATTERN));
     98:         PetscCall(MatNorm(B, NORM_INFINITY, &norm));
     99:         PetscCheck(norm < 100 * PETSC_MACHINE_EPSILON, PetscObjectComm((PetscObject)ksp), PETSC_ERR_PLIB, "KSPMatSolve() and MatMatSolve() difference has nonzero norm %g", (double)norm);
    100:       }
    101:     } else {
    102:       PetscCall(MatZeroEntries(B));
    103:       PetscCall(KSPMatSolve(ksp, B, X));
    104:       PetscCall(KSPGetConvergedReason(ksp, &reason));
    105:       PetscCheck(reason == KSP_CONVERGED_HAPPY_BREAKDOWN, PetscObjectComm((PetscObject)ksp), PETSC_ERR_PLIB, "KSPConvergedReason() %s != KSP_CONVERGED_HAPPY_BREAKDOWN", KSPConvergedReasons[reason]);
    106:       PetscCall(MatDenseGetArrayWrite(B, &x));
    107:       for (i = 0; i < m * N; ++i) x[i] = 1.0;
    108:       PetscCall(MatDenseRestoreArrayWrite(B, &x));
    109:       PetscCall(KSPMatSolve(ksp, B, X));
    110:       PetscCall(KSPGetConvergedReason(ksp, &reason));
    111:       PetscCheck(reason == KSP_DIVERGED_BREAKDOWN || deflation >= 0.0, PetscObjectComm((PetscObject)ksp), PETSC_ERR_PLIB, "KSPConvergedReason() %s != KSP_DIVERGED_BREAKDOWN", KSPConvergedReasons[reason]);
    112:     }
    113:   } else {
    114:     PetscCall(KSPSetOperators(ksp, KA, KA));
    115:     PetscCall(MatGetSize(KA, &M, NULL));
    116:     PetscCall(VecCreateFromOptions(PETSC_COMM_WORLD, NULL, 1, m * N, M, &cb));
    117:     PetscCall(VecCreateFromOptions(PETSC_COMM_WORLD, NULL, 1, m * N, M, &cx));
    118:     PetscCall(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:     PetscCall(KSPSolve(ksp, cb, cx));
    122:     PetscCall(VecDestroy(&cx));
    123:     PetscCall(VecDestroy(&cb));
    124:   }
    125:   PetscCall(MatDestroy(&X));
    126:   PetscCall(MatDestroy(&B));
    127:   PetscCall(PetscFree2(S, T));
    128:   PetscCall(MatDestroy(&KA));
    129:   PetscCall(MatDestroy(&A));
    130:   PetscCall(KSPDestroy(&ksp));
    131:   PetscCall(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*/