1: !
  2: !   Description: Solves a linear system with a block of right-hand sides using KSPHPDDM.
  3: !
  4: #include <petsc/finclude/petscksp.h>
  5: program main
  6:   use petscksp
  7:   implicit none
  8:   Mat X, B
  9:   Mat A
 10:   KSP ksp
 11:   PC pc
 12:   Mat F
 13:   PetscScalar, parameter :: alpha = -1.0
 14:   PetscReal norm
 15:   PetscInt m, K
 16:   PetscViewer viewer
 17:   character(PETSC_MAX_PATH_LEN) name
 18:   PetscBool flg
 19:   PetscErrorCode ierr

 21:   PetscCallA(PetscInitialize(ierr))
 22:   PetscCallA(PetscOptionsGetString(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-f', name, flg, ierr))
 23:   PetscCheckA(flg .neqv. PETSC_FALSE, PETSC_COMM_WORLD, PETSC_ERR_SUP, 'Must provide a binary file for the matrix')
 24:   K = 5
 25:   PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', K, flg, ierr))
 26:   PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr))
 27:   PetscCallA(KSPCreate(PETSC_COMM_WORLD, ksp, ierr))
 28:   PetscCallA(KSPSetOperators(ksp, A, A, ierr))
 29:   PetscCallA(PetscViewerBinaryOpen(PETSC_COMM_WORLD, name, FILE_MODE_READ, viewer, ierr))
 30:   PetscCallA(MatLoad(A, viewer, ierr))
 31:   PetscCallA(PetscViewerDestroy(viewer, ierr))
 32:   PetscCallA(MatGetLocalSize(A, m, PETSC_NULL_INTEGER, ierr))
 33:   PetscCallA(MatCreateDense(PETSC_COMM_WORLD, m, PETSC_DECIDE, PETSC_DECIDE, K, PETSC_NULL_SCALAR_ARRAY, B, ierr))
 34:   PetscCallA(MatCreateDense(PETSC_COMM_WORLD, m, PETSC_DECIDE, PETSC_DECIDE, K, PETSC_NULL_SCALAR_ARRAY, X, ierr))
 35:   PetscCallA(MatSetRandom(B, PETSC_NULL_RANDOM, ierr))
 36:   PetscCallA(KSPSetFromOptions(ksp, ierr))
 37:   PetscCallA(KSPSetUp(ksp, ierr))
 38:   PetscCallA(KSPMatSolve(ksp, B, X, ierr))
 39:   PetscCallA(KSPGetMatSolveBatchSize(ksp, M, ierr))
 40:   if (M /= PETSC_DECIDE) then
 41:     PetscCallA(KSPSetMatSolveBatchSize(ksp, PETSC_DECIDE, ierr))
 42:     PetscCallA(MatZeroEntries(X, ierr))
 43:     PetscCallA(KSPMatSolve(ksp, B, X, ierr))
 44:   end if
 45:   PetscCallA(KSPGetPC(ksp, pc, ierr))
 46:   PetscCallA(PetscObjectTypeCompare(pc, PCLU, flg, ierr))
 47:   if (flg) then
 48:     PetscCallA(PCFactorGetMatrix(pc, F, ierr))
 49:     PetscCallA(MatMatSolve(F, B, B, ierr))
 50:     PetscCallA(MatAYPX(B, alpha, X, SAME_NONZERO_PATTERN, ierr))
 51:     PetscCallA(MatNorm(B, NORM_INFINITY, norm, ierr))
 52:     PetscCheckA(norm < 100*PETSC_MACHINE_EPSILON, PETSC_COMM_WORLD, PETSC_ERR_PLIB, 'KSPMatSolve() and MatMatSolve() difference has nonzero norm')
 53:   end if
 54:   PetscCallA(MatDestroy(X, ierr))
 55:   PetscCallA(MatDestroy(B, ierr))
 56:   PetscCallA(MatDestroy(A, ierr))
 57:   PetscCallA(KSPDestroy(ksp, ierr))
 58:   PetscCallA(PetscFinalize(ierr))
 59: end

 61: !/*TEST
 62: !
 63: !   testset:
 64: !      nsize: 2
 65: !      requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
 66: !      args: -ksp_converged_reason -ksp_max_it 1000 -f ${DATAFILESPATH}/matrices/hpddm/GCRODR/A_400.dat
 67: !      test:
 68: !         suffix: 1
 69: !         output_file: output/ex77_1.out
 70: !         args:
 71: !      test:
 72: !         suffix: 2a
 73: !         requires: hpddm
 74: !         output_file: output/ex77_2_ksp_hpddm_type-gmres.out
 75: !         args: -ksp_type hpddm -pc_type asm -ksp_hpddm_type gmres
 76: !      test:
 77: !         suffix: 2b
 78: !         requires: hpddm
 79: !         output_file: output/ex77_2_ksp_hpddm_type-bgmres.out
 80: !         args: -ksp_type hpddm -pc_type asm -ksp_hpddm_type bgmres
 81: !      test:
 82: !         suffix: 3a
 83: !         requires: hpddm
 84: !         output_file: output/ex77_3_ksp_hpddm_type-gcrodr.out
 85: !         args: -ksp_type hpddm -ksp_hpddm_recycle 5 -ksp_hpddm_type gcrodr
 86: !      test:
 87: !         suffix: 3b
 88: !         requires: hpddm
 89: !         output_file: output/ex77_3_ksp_hpddm_type-bgcrodr.out
 90: !         args: -ksp_type hpddm -ksp_hpddm_recycle 5 -ksp_hpddm_type bgcrodr
 91: !      test:
 92: !         nsize: 4
 93: !         suffix: 4
 94: !         requires: hpddm
 95: !         output_file: output/ex77_4.out
 96: !         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
 97: !   test:
 98: !      nsize: 1
 99: !      suffix: preonly
100: !      requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
101: !      output_file: output/empty.out
102: !      args: -N 6 -f ${DATAFILESPATH}/matrices/hpddm/GCRODR/A_400.dat -pc_type lu -ksp_type hpddm -ksp_hpddm_type preonly
103: !   test:
104: !      nsize: 4
105: !      suffix: 4_slepc
106: !      requires: hpddm datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) slepc defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
107: !      output_file: output/ex77_4.out
108: !      filter: sed "/^ksp_hpddm_recycle_ Linear eigensolve converged/d"
109: !      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 dense -ksp_hpddm_recycle_eps_converged_reason -ksp_hpddm_recycle_st_pc_type redundant
110: !
111: !TEST*/