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

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