diff --git a/fpm.toml b/fpm.toml index 66ea2c3..7cc1a76 100644 --- a/fpm.toml +++ b/fpm.toml @@ -12,13 +12,14 @@ auto-tests = true auto-examples = true module-naming = false link = ["blas", "lapack"] +#link = ["openblas"] [install] -library = false +library = true [fortran] implicit-typing = false -implicit-external = false +implicit-external = true source-form = "free" [dependencies] diff --git a/src/AbstractVector.f90 b/src/AbstractVector.f90 index 94c81b5..e06a2dd 100644 --- a/src/AbstractVector.f90 +++ b/src/AbstractVector.f90 @@ -13,41 +13,41 @@ module AbstractVector contains private - !> Basic operations. + !> Basic operations. procedure, pass(from), public :: copy generic , public :: assignment(=) => copy procedure(abstract_zero), deferred, public :: zero - !> Scalar-vector product. + !> Scalar-vector product. procedure(abstract_scal), deferred, public :: scal - !> Vector-vector operations. + !> Vector-vector operations. procedure(axpby_interface), deferred, pass(self), public :: axpby procedure, pass(self), public :: add procedure, pass(self), public :: sub - !> Reduction operations. + !> Reduction operations. procedure(abstract_dot), pass(self), deferred, public :: dot procedure, pass(self), public :: norm end type abstract_vector - ! --> Abstract interfaces for the type-bound procedures. + !> Abstract interfaces for the type-bound procedures. abstract interface - !> Basic operations. + !> Basic operations. subroutine abstract_zero(self) import abstract_vector class(abstract_vector), intent(inout) :: self end subroutine abstract_zero - !> Scalar-vector product. + !> Scalar-vector product. subroutine abstract_scal(self, alpha) import abstract_vector, wp class(abstract_vector), intent(inout) :: self real(kind=wp), intent(in) :: alpha end subroutine abstract_scal - !> Vector-vector operations. + !> Vector-vector operations. subroutine axpby_interface(self, alpha, vec, beta) import abstract_vector, wp class(abstract_vector), intent(inout) :: self @@ -55,7 +55,7 @@ subroutine axpby_interface(self, alpha, vec, beta) real(kind=wp) , intent(in) :: alpha, beta end subroutine axpby_interface - !> Reduction operations. + !> Reduction operations. real(kind=wp) function abstract_dot(self, vec) result(alpha) import abstract_vector, wp class(abstract_vector), intent(in) :: self, vec @@ -67,6 +67,7 @@ end function abstract_dot subroutine copy(out, from) class(abstract_vector), intent(in) :: from class(abstract_vector), allocatable, intent(out) :: out + if (allocated(out)) deallocate(out) allocate(out, source=from) return end subroutine copy @@ -74,7 +75,7 @@ end subroutine copy subroutine add(self, y) class(abstract_vector), intent(inout) :: self class(abstract_vector), intent(in) :: y - ! --> Vector addition. + !> Vector addition. call self%axpby(1.0_wp, y, 1.0_wp) return end subroutine add @@ -82,7 +83,7 @@ end subroutine add subroutine sub(self, y) class(abstract_vector), intent(inout) :: self class(abstract_vector), intent(in) :: y - ! --> Vector subtraction. + !> Vector subtraction. call self%axpby(1.0_wp, y, -1.0_wp) return end subroutine sub @@ -93,29 +94,32 @@ real(kind=wp) function norm(self) result(alpha) end function norm subroutine get_vec(y, X, v) - !> Output vector. + !> Output vector. class(abstract_vector), allocatable, intent(out) :: y - !> Krylov basis. + !> Krylov basis. class(abstract_vector), intent(in) :: X(:) - !> Coordinates of the output vector y in the Krylov basis X. + !> Coordinates of the output vector y in the Krylov basis X. real(kind=wp), intent(in) :: v(:) - !> Temporary Krylov vector. + !> Temporary Krylov vector. class(abstract_vector), allocatable :: wrk - !> Miscellaneous. + !> Miscellaneous. integer :: i - ! --> Check sizes. + !> Check sizes. if (size(X) .ne. size(v)) then write(*, *) "INFO : Krylov basis X and low-dimension vector v have different sizes." - call exit() + return endif - ! --> Initialize output vector. + !> Initialize output vector. if (allocated(y) .eqv. .false.) allocate(y, source=X(1)) ; call y%zero() - ! --> Compute output vector. + !> Compute output vector. + if (.not. allocated(wrk)) allocate(wrk, source=X(1)) do i = 1, size(X) - wrk = X(i) ; call wrk%scal(v(i)) + wrk = X(i) + call wrk%scal(v(i)) call y%add(wrk) enddo + if (allocated(wrk)) deallocate(wrk) return end subroutine get_vec -end module AbstractVector +end module AbstractVector \ No newline at end of file