Actual source code: ex77f.F90

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

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