1! 2! Description: Solves a linear system with a block of right-hand sides using KSPHPDDM. 3! 4 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 21 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 KSPGetMatSolveBlockSize(ksp,M,ierr);CHKERRA(ierr) 47 if (M .ne. PETSC_DECIDE) then 48 call KSPSetMatSolveBlockSize(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 70 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_block_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 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_block_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*/ 122