petsc_binding.f90 4.97 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

module LinearAlgebra

  implicit none

  private

#include <finclude/petsc.h90>

  type, public :: sparse_mat
     private

     !> solution 
     Vec :: x
     !> righ hand side
     Vec :: b
     !> linear system matrix
     Mat :: A
     ! linear solver context
     KSP :: ksp
     ! preconditioner context
     PC  :: pc
 
     logical :: sym = .false.

     integer(kind=4), dimension(:),  pointer :: i_indices => null()
     integer(kind=4), dimension(:),  pointer :: j_indices => null()

  end type sparse_mat

Rémy Mozul's avatar
Rémy Mozul committed
31 32 33 34 35 36
  public init   , &
         declare, &
         build  , &
         solve  , &
         erase  , &
         finalize
37 38 39

  contains

Rémy Mozul's avatar
Rémy Mozul committed
40
  subroutine init()
41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177
    implicit none
    !
    integer(kind=4)   :: ierr
    character(len=32) :: IAM
    !      12345678901234567890123456789012
    IAM = 'SparseLinearAlgebra::sparse_init'

    call PetscInitialize(PETSC_NULL_CHARACTER,ierr)

    if( ierr > 0 ) then
      print *, '[', IAM, '] PetscInitalize failed with value : ', ierr
      stop
    end if

  end subroutine

  subroutine declare(matrix,nb_dofs, nb_non_zero, i_indices, j_indices, is_sym, info)
    implicit none
    type(sparse_mat) :: matrix
    integer(kind=4)                        :: nb_dofs, nb_non_zero, info
    integer(kind=4), dimension(:), pointer :: i_indices
    integer(kind=4), dimension(:), pointer :: j_indices
    logical                                :: is_sym
    !
    PetscBool :: is_spd
    real(kind=8) :: tol
    character(len=35) :: IAM
    !      12345678901234567890123456789012345
    IAM = 'SparseLinearAlgebra::sparse_declare'

    is_spd = PETSC_FALSE
    tol = 1.e-5

    ! create right hand side
    call VecCreate(PETSC_COMM_WORLD,matrix%b,info)
    call VecSetSizes(matrix%b,PETSC_DECIDE,nb_dofs,info)
    call VecSetFromOptions(matrix%b,info)

    ! create solution
    call VecDuplicate(matrix%b,matrix%x,info)

    ! create left hand side
    call MatCreate(PETSC_COMM_WORLD,matrix%A,info)
    call MatSetSizes(matrix%A,PETSC_DECIDE,PETSC_DECIDE,nb_dofs,nb_dofs,info)
    call MatSetFromOptions(matrix%A,info)
    !if( is_sym ) then
    !  is_spd = PETSC_TRUE
    !  matrix%sym = .true.
    !end if
    !call MatSetOption(matrix%A,MAT_SPD,is_spd)
    call MatSetUp(matrix%A,info)

    ! create linear solver
    call KSPCreate(PETSC_COMM_WORLD,matrix%ksp,info)

    call KSPGetPC(matrix%ksp,matrix%pc,info)
    call PCSetType(matrix%pc,PCJACOBI,info)
    !call KSPSetTolerances(matrix%ksp,tol,PETSC_DEFAULT_REAL,PETSC_DEFAULT_REAL,PETSC_DEFAULT_REAL,info)

    matrix%i_indices => i_indices
    matrix%j_indices => j_indices

  end subroutine

  subroutine build(matrix, val, info)
    implicit none
    type(sparse_mat) :: matrix
    real(kind=8)   , dimension(:), pointer :: val
    integer(kind=4) :: info
    !
    integer(kind=4)   :: i
    character(len=33) :: IAM
    !      123456789012345678901234567890123
    IAM = 'SparseLinearAlgebra::sparse_build'

    do i = 1, size(val)
      call MatSetValues(matrix%A,1,matrix%i_indices(i)-1,1,matrix%j_indices(i)-1, &
                        val(i),INSERT_VALUES, info)
    end do
    if( .not. matrix%sym ) then
      do i = 1, size(val)
        call MatSetValues(matrix%A,1,matrix%j_indices(i)-1,1,matrix%i_indices(i)-1, &
                          val(i),INSERT_VALUES, info)
      end do
    end if

    call MatAssemblyBegin(matrix%A,MAT_FINAL_ASSEMBLY,info)
    call MatAssemblyEnd(matrix%A,MAT_FINAL_ASSEMBLY,info)

    call KSPSetOperators(matrix%ksp,matrix%A,matrix%A,info)
    call KSPSetFromOptions(matrix%ksp,info)

  end subroutine

  subroutine solve(matrix, rhs, info)
    implicit none
    type(sparse_mat) :: matrix
    real(kind=8)   , dimension(:), pointer :: rhs
    integer(kind=4) :: info
    !
    integer(kind=4) :: i
    character(len=33) :: IAM
    !      123456789012345678901234567890123
    IAM = 'SparseLinearAlgebra::sparse_solve'

    do i = 1, size(rhs)
      call VecSetValues(matrix%b,1,i-1,rhs(i),INSERT_VALUES,info)
    end do
    call VecAssemblyBegin(matrix%b,info)
    call VecAssemblyEnd(matrix%b,info)

    call KSPSolve(matrix%ksp,matrix%b,matrix%x,info)
    !activate only if verbose...
    !call KSPView(matrix%ksp,PETSC_VIEWER_STDOUT_WORLD,info)

    do i = 1, size(rhs)
      call VecGetValues(matrix%x,1,i-1,rhs(i),info)
    end do

  end subroutine

  subroutine erase(matrix)
    implicit none
    type(sparse_mat) :: matrix
    !
    integer(kind=4)   :: info
    character(len=33) :: IAM
    !      123456789012345678901234567890123
    IAM = 'SparseLinearAlgebra::sparse_erase'

    call VecDestroy(matrix%x,info)
    call VecDestroy(matrix%b,info)
    call MatDestroy(matrix%A,info)
    call KSPDestroy(matrix%ksp,info)

  end subroutine

Rémy Mozul's avatar
Rémy Mozul committed
178
  subroutine finalize()
179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195
    implicit none
    !
    integer(kind=4)   :: ierr
    character(len=36) :: IAM
    !      123456789012345678901234567890123456
    IAM = 'SparseLinearAlgebra::sparse_finalize'

    call PetsCFinalize(ierr)

    if( ierr > 0 ) then
      print *, '[', IAM, '] PetscFinalize failed with value : ', ierr
      stop
    end if

  end subroutine

end module