1: !
  2: !   Description: Solves a linear system with a block of right-hand sides using KSPHPDDM.
  3: !

  5:       program main
  6: #include <petsc/finclude/petscksp.h>
  7:       use petscksp
  8:       implicit none
  9:       Mat                            X,B
 10:       Mat                            A
 11:       KSP                            ksp
 12:       PC                             pc
 13:       Mat                            F
 14:       PetscScalar                    alpha
 15:       PetscReal                      norm
 16:       PetscInt                       m,K
 17:       PetscViewer                    viewer
 18:       character*(PETSC_MAX_PATH_LEN) name
 19:       PetscBool                      flg
 20:       PetscErrorCode                 ierr

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

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