Actual source code: ex76f.F90

      1: !
      2: !   Description: Solves a linear systems using PCHPDDM.
      3: !

      5:       program main
      6: #include <petsc/finclude/petscksp.h>
      7:       use petscksp
      8:       use petscisdef
      9:       implicit none
     10:       Vec                            x,b
     11:       Mat                            A,aux,Y,C
     12:       KSP                            ksp
     13:       PC                             pc
     14:       IS                             is,sizes
     15:       PetscScalar                    one
     16:       PetscInt, pointer ::           idx(:)
     17:       PetscMPIInt                    rank,size
     18:       PetscInt                       m,N
     19:       PetscViewer                    viewer
     20:       character*(PETSC_MAX_PATH_LEN) dir,name
     21:       character*(8)                  fmt
     22:       character(1)                   crank,csize
     23:       PetscBool                      flg
     24:       PetscErrorCode                 ierr

     26:       PetscCallA(PetscInitialize(ierr))

     28:       PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,size,ierr))
     29:       N = 1
     30:       PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-rhs',N,flg,ierr))
     31:       PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))
     32:       PetscCallA(MatCreate(PETSC_COMM_WORLD,A,ierr))
     33:       PetscCallA(MatCreate(PETSC_COMM_SELF,aux,ierr))
     34:       PetscCallA(ISCreate(PETSC_COMM_SELF,is,ierr))
     35:       dir = '.'
     36:       PetscCallA(PetscOptionsGetString(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-load_dir',dir,flg,ierr))
     37:       fmt = '(I1)'
     38:       write (crank,fmt) rank
     39:       write (csize,fmt) size
     40:       write (name,'(a)')trim(dir)//'/sizes_'//crank//'_'//csize//'.dat'
     41:       PetscCallA(PetscViewerBinaryOpen(PETSC_COMM_SELF,name,FILE_MODE_READ, viewer,ierr))
     42:       PetscCallA(ISCreate(PETSC_COMM_SELF,sizes,ierr))
     43:       PetscCallA(ISLoad(sizes,viewer,ierr))
     44:       PetscCallA(ISGetIndicesF90(sizes,idx,ierr))
     45:       PetscCallA(MatSetSizes(A,idx(1),idx(2),idx(3),idx(4),ierr))
     46:       PetscCallA(ISRestoreIndicesF90(sizes,idx,ierr))
     47:       PetscCallA(ISDestroy(sizes,ierr))
     48:       PetscCallA(PetscViewerDestroy(viewer,ierr))
     49:       write (name,'(a)')trim(dir)//'/A.dat'
     50:       PetscCallA(PetscViewerBinaryOpen(PETSC_COMM_WORLD,name,FILE_MODE_READ,viewer,ierr))
     51:       PetscCallA(MatLoad(A,viewer,ierr))
     52:       PetscCallA(PetscViewerDestroy(viewer,ierr))
     53:       write (name,'(a)')trim(dir)//'/is_'//crank//'_'//csize//'.dat'
     54:       PetscCallA(PetscViewerBinaryOpen(PETSC_COMM_SELF,name,FILE_MODE_READ,viewer,ierr))
     55:       PetscCallA(ISLoad(is,viewer,ierr))
     56:       PetscCallA(PetscViewerDestroy(viewer,ierr))
     57:       write (name,'(a)')trim(dir)//'/Neumann_'//crank//'_'//csize//'.dat'
     58:       PetscCallA(PetscViewerBinaryOpen(PETSC_COMM_SELF,name,FILE_MODE_READ,viewer,ierr))
     59:       PetscCallA(MatSetBlockSizesFromMats(aux,A,A,ierr))
     60:       PetscCallA(MatLoad(aux,viewer,ierr))
     61:       PetscCallA(PetscViewerDestroy(viewer,ierr))
     62:       PetscCallA(MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE,ierr))
     63:       PetscCallA(MatSetOption(aux,MAT_SYMMETRIC,PETSC_TRUE,ierr))
     64:       PetscCallA(KSPCreate(PETSC_COMM_WORLD,ksp,ierr))
     65:       PetscCallA(KSPSetOperators(ksp,A,A,ierr))
     66:       PetscCallA(KSPGetPC(ksp,pc,ierr))
     67:       PetscCallA(PCSetType(pc,PCHPDDM,ierr))
     68: #if defined(PETSC_HAVE_HPDDM) && defined(PETSC_HAVE_DYNAMIC_LIBRARIES) && defined(PETSC_USE_SHARED_LIBRARIES)
     69:       PetscCallA(PCHPDDMSetAuxiliaryMat(pc,is,aux,PETSC_NULL_FUNCTION,PETSC_NULL_INTEGER,ierr))
     70:       PetscCallA(PCHPDDMHasNeumannMat(pc,PETSC_FALSE,ierr))
     71: #endif
     72:       PetscCallA(ISDestroy(is,ierr))
     73:       PetscCallA(MatDestroy(aux,ierr))
     74:       PetscCallA(KSPSetFromOptions(ksp,ierr))
     75:       PetscCallA(MatCreateVecs(A,x,b,ierr))
     76:       one = 1.0
     77:       PetscCallA(VecSet(b,one,ierr))
     78:       PetscCallA(KSPSolve(ksp,b,x,ierr))
     79:       PetscCallA(VecGetLocalSize(x,m,ierr))
     80:       PetscCallA(VecDestroy(x,ierr))
     81:       PetscCallA(VecDestroy(b,ierr))
     82:       if (N .gt. 1) then
     83:         PetscCallA(PetscOptionsClearValue(PETSC_NULL_OPTIONS,'-ksp_converged_reason',ierr))
     84:         PetscCallA(KSPSetFromOptions(ksp,ierr))
     85:         PetscCallA(MatCreateDense(PETSC_COMM_WORLD,m,PETSC_DECIDE,PETSC_DECIDE,N,PETSC_NULL_SCALAR,C,ierr))
     86:         PetscCallA(MatCreateDense(PETSC_COMM_WORLD,m,PETSC_DECIDE,PETSC_DECIDE,N,PETSC_NULL_SCALAR,Y,ierr))
     87:         PetscCallA(MatSetRandom(C,PETSC_NULL_RANDOM,ierr))
     88:         PetscCallA(KSPMatSolve(ksp,C,Y,ierr))
     89:         PetscCallA(MatDestroy(Y,ierr))
     90:         PetscCallA(MatDestroy(C,ierr))
     91:       endif
     92:       PetscCallA(KSPDestroy(ksp,ierr))
     93:       PetscCallA(MatDestroy(A,ierr))
     94:       PetscCallA(PetscFinalize(ierr))
     95:       end

     97: !/*TEST
     98: !
     99: !   test:
    100: !      requires: hpddm slepc datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
    101: !      output_file: output/ex76_1.out
    102: !      nsize: 4
    103: !      args: -ksp_rtol 1e-3 -ksp_converged_reason -pc_type {{bjacobi hpddm}shared output} -pc_hpddm_coarse_sub_pc_type lu -sub_pc_type lu -options_left no -load_dir ${DATAFILESPATH}/matrices/hpddm/GENEO
    104: !
    105: !   test:
    106: !      requires: hpddm slepc datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
    107: !      suffix: geneo
    108: !      output_file: output/ex76_geneo_pc_hpddm_levels_1_eps_nev-5.out
    109: !      nsize: 4
    110: !      args: -ksp_converged_reason -pc_type hpddm -pc_hpddm_levels_1_sub_pc_type cholesky -pc_hpddm_levels_1_eps_nev 5 -pc_hpddm_levels_1_st_pc_type cholesky -pc_hpddm_coarse_p {{1 2}shared output} -pc_hpddm_coarse_pc_type redundant -load_dir ${DATAFILESPATH}/matrices/hpddm/GENEO
    111: !
    112: !   test:
    113: !      requires: hpddm slepc datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
    114: !      suffix: fgmres_geneo_20_p_2
    115: !      output_file: output/ex76_fgmres_geneo_20_p_2.out
    116: !      nsize: 4
    117: !      args: -ksp_converged_reason -pc_type hpddm -pc_hpddm_levels_1_sub_pc_type lu -pc_hpddm_levels_1_eps_nev 20 -pc_hpddm_coarse_p 2 -pc_hpddm_coarse_pc_type redundant -ksp_type fgmres -pc_hpddm_coarse_mat_type {{baij sbaij}shared output} -load_dir ${DATAFILESPATH}/matrices/hpddm/GENEO
    118: !
    119: !   test:
    120: !      requires: hpddm slepc datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
    121: !      suffix: fgmres_geneo_20_p_2_geneo
    122: !      output_file: output/ex76_fgmres_geneo_20_p_2.out
    123: !      nsize: 4
    124: !      args: -ksp_converged_reason -pc_type hpddm -pc_hpddm_levels_1_sub_pc_type cholesky -pc_hpddm_levels_1_eps_nev 20 -pc_hpddm_levels_2_p 2 -pc_hpddm_levels_2_mat_type {{baij sbaij}shared output} -pc_hpddm_levels_2_eps_nev {{5 20}shared output} -pc_hpddm_levels_2_sub_pc_type cholesky -pc_hpddm_levels_2_ksp_type gmres -ksp_type fgmres -pc_hpddm_coarse_mat_type {{baij sbaij}shared output} -load_dir ${DATAFILESPATH}/matrices/hpddm/GENEO
    125: !   # PCHPDDM + KSPHPDDM test to exercise multilevel + multiple RHS in one go
    126: !   test:
    127: !      requires: hpddm slepc datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
    128: !      suffix: fgmres_geneo_20_p_2_geneo_rhs
    129: !      output_file: output/ex76_fgmres_geneo_20_p_2.out
    130: !      nsize: 4
    131: !      args: -ksp_converged_reason -pc_type hpddm -pc_hpddm_levels_1_sub_pc_type cholesky -pc_hpddm_levels_1_eps_nev 20 -pc_hpddm_levels_2_p 2 -pc_hpddm_levels_2_mat_type baij -pc_hpddm_levels_2_eps_nev 5 -pc_hpddm_levels_2_sub_pc_type cholesky -pc_hpddm_levels_2_ksp_max_it 10 -pc_hpddm_levels_2_ksp_type hpddm -pc_hpddm_levels_2_ksp_hpddm_type gmres -ksp_type hpddm -ksp_hpddm_variant flexible -pc_hpddm_coarse_mat_type baij -load_dir ${DATAFILESPATH}/matrices/hpddm/GENEO -rhs 4
    132: !
    133: !TEST*/