From 7d390f737082b4da4465f3bf45330d297893195d Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 16 May 2024 09:26:33 +0200 Subject: [PATCH 01/32] add source --- src/CMakeLists.txt | 1 + src/stdlib_linalg_eigenvalues.fypp | 569 +++++++++++++++++++++++++++++ 2 files changed, 570 insertions(+) create mode 100644 src/stdlib_linalg_eigenvalues.fypp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index a031ab823..4ed16ae3b 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -25,6 +25,7 @@ set(fppFiles stdlib_linalg_outer_product.fypp stdlib_linalg_kronecker.fypp stdlib_linalg_cross_product.fypp + stdlib_linalg_eigenvalues.fypp stdlib_linalg_solve.fypp stdlib_linalg_determinant.fypp stdlib_linalg_state.fypp diff --git a/src/stdlib_linalg_eigenvalues.fypp b/src/stdlib_linalg_eigenvalues.fypp new file mode 100644 index 000000000..3b49ec968 --- /dev/null +++ b/src/stdlib_linalg_eigenvalues.fypp @@ -0,0 +1,569 @@ +#:include "common.fypp" +module stdlib_linalg_eigenvalues + use stdlib_linalg_constants + use stdlib_linalg_blas + use stdlib_linalg_lapack + use stdlib_linalg_state + use iso_fortran_env,only:real32,real64,real128,int8,int16,int32,int64,stderr => error_unit + implicit none(type,external) + private + + !> Eigendecomposition of a square matrix: return eigenvalues, and optionally eigenvectors + public :: eig + !> Eigenvalues of a square matrix + public :: eigvals + !> Eigendecomposition of a real symmetric or complex hermitian matrix + public :: eigh + !> Eigenvalues of a real symmetric or complex hermitian matrix + public :: eigvalsh + + ! Numpy: eigenvalues, eigenvectors = eig(a) + ! eigenvalues = eigvals(a) + ! Scipy: eig(a, b=None, left=False, right=True, overwrite_a=False, overwrite_b=False, check_finite=True, homogeneous_eigvals=False) + + ! Numpy: eigenvalues, eigenvectors = eigh(a, uplo='L') + ! eigenvalues = eigvalsh(a) + ! Scipy: eigh(a, b=None, *, lower=True, eigvals_only=False, overwrite_a=False, overwrite_b=False, turbo=, eigvals=, type=1, check_finite=True, subset_by_index=None, subset_by_value=None, driver=None) + + interface eig + #:for rk,rt,ri in ALL_KINDS_TYPES + module procedure stdlib_linalg_eig_${ri}$ + #:endfor + end interface eig + + interface eigvals + #:for rk,rt,ri in ALL_KINDS_TYPES + module procedure stdlib_linalg_eigvals_${ri}$ + module procedure stdlib_linalg_eigvals_noerr_${ri}$ + #:endfor + end interface eigvals + + interface eigh + #:for rk,rt,ri in ALL_KINDS_TYPES + module procedure stdlib_linalg_eigh_${ri}$ + #:endfor + end interface eigh + + interface eigvalsh + #:for rk,rt,ri in ALL_KINDS_TYPES + module procedure stdlib_linalg_eigvalsh_${ri}$ + module procedure stdlib_linalg_eigvalsh_noerr_${ri}$ + #:endfor + end interface eigvalsh + + + character(*), parameter :: this = 'eigenvalues' + + + contains + + !> Request for eigenvector calculation + elemental character function eigenvectors_flag(required) + logical(lk), intent(in) :: required + eigenvectors_flag = merge('V','N',required) + end function eigenvectors_flag + + !> Request for symmetry side (default: lower) + elemental character function symmetric_triangle(upper) + logical(lk), optional, intent(in) :: upper + if (present(upper)) then + symmetric_triangle = merge('U','L',upper) + else + symmetric_triangle = 'L' + endif + end function symmetric_triangle + + !> Process GEEV output flags + elemental subroutine geev_info(err,info,m,n) + !> Error handler + type(linalg_state), intent(inout) :: err + !> geev return flag + integer(ilp), intent(in) :: info + !> Input matrix size + integer(ilp), intent(in) :: m,n + + select case (info) + case (0) + ! Success! + err%state = LINALG_SUCCESS + case (-1) + err = linalg_state(this,LINALG_INTERNAL_ERROR,'Invalid task ID: left eigenvectors.') + case (-2) + err = linalg_state(this,LINALG_INTERNAL_ERROR,'Invalid task ID: right eigenvectors.') + case (-5,-3) + err = linalg_state(this,LINALG_VALUE_ERROR,'invalid matrix size: a=',[m,n]) + case (-9) + err = linalg_state(this,LINALG_VALUE_ERROR,'insufficient left vector matrix size.') + case (-11) + err = linalg_state(this,LINALG_VALUE_ERROR,'insufficient right vector matrix size.') + case (-13) + err = linalg_state(this,LINALG_INTERNAL_ERROR,'Insufficient work array size.') + case (1:) + err = linalg_state(this,LINALG_ERROR,'Eigenvalue computation did not converge.') + case default + err = linalg_state(this,LINALG_INTERNAL_ERROR,'Unknown error returned by geev.') + end select + + end subroutine geev_info + + + !> Process SYEV/HEEV output flags + elemental subroutine heev_info(err,info,m,n) + !> Error handler + type(linalg_state), intent(inout) :: err + !> geev return flag + integer(ilp), intent(in) :: info + !> Input matrix size + integer(ilp), intent(in) :: m,n + + select case (info) + case (0) + ! Success! + err%state = LINALG_SUCCESS + case (-1) + err = linalg_state(this,LINALG_INTERNAL_ERROR,'Invalid eigenvector request.') + case (-2) + err = linalg_state(this,LINALG_INTERNAL_ERROR,'Invalid triangular section request.') + case (-5,-3) + err = linalg_state(this,LINALG_VALUE_ERROR,'invalid matrix size: a=',[m,n]) + case (-8) + err = linalg_state(this,LINALG_INTERNAL_ERROR,'insufficient workspace size.') + case (1:) + err = linalg_state(this,LINALG_ERROR,'Eigenvalue computation did not converge.') + case default + err = linalg_state(this,LINALG_INTERNAL_ERROR,'Unknown error returned by syev/heev.') + end select + + end subroutine heev_info + + #:for rk,rt,ri in ALL_KINDS_TYPES + + !> Return an array of eigenvalues of matrix A. + function stdlib_linalg_eigvals_${ri}$(a,err) result(lambda) + !> Input matrix A[m,n] + ${rt}$, intent(in), target :: a(:,:) + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state), intent(out) :: err + !> Array of singular values + complex(${rk}$), allocatable :: lambda(:) + + !> Create + ${rt}$, pointer :: amat(:,:) + integer(ilp) :: m,n,k + + !> Create an internal pointer so the intent of A won't affect the next call + amat => a + + m = size(a,1,kind=ilp) + n = size(a,2,kind=ilp) + k = min(m,n) + + !> Allocate return storage + allocate(lambda(k)) + + !> Compute eigenvalues only + call stdlib_linalg_eig_${ri}$(amat,lambda,overwrite_a=.false.,err=err) + + end function stdlib_linalg_eigvals_${ri}$ + + !> Return an array of eigenvalues of matrix A. + function stdlib_linalg_eigvals_noerr_${ri}$(a) result(lambda) + !> Input matrix A[m,n] + ${rt}$, intent(in), target :: a(:,:) + !> Array of singular values + complex(${rk}$), allocatable :: lambda(:) + + !> Create + ${rt}$, pointer :: amat(:,:) + integer(ilp) :: m,n,k + + !> Create an internal pointer so the intent of A won't affect the next call + amat => a + + m = size(a,1,kind=ilp) + n = size(a,2,kind=ilp) + k = min(m,n) + + !> Allocate return storage + allocate(lambda(k)) + + !> Compute eigenvalues only + call stdlib_linalg_eig_${ri}$(amat,lambda,overwrite_a=.false.) + + end function stdlib_linalg_eigvals_noerr_${ri}$ + + !> Eigendecomposition of matrix A returning an array `lambda` of eigenvalues, + !> and optionally right or left eigenvectors. + subroutine stdlib_linalg_eig_${ri}$(a,lambda,right,left,overwrite_a,err) + !> Input matrix A[m,n] + ${rt}$, intent(inout), target :: a(:,:) + !> Array of eigenvalues + complex(${rk}$), intent(out) :: lambda(:) + !> The columns of RIGHT contain the right eigenvectors of A + complex(${rk}$), optional, intent(out), target :: right(:,:) + !> The columns of LEFT contain the left eigenvectors of A + complex(${rk}$), optional, intent(out), target :: left(:,:) + !> [optional] Can A data be overwritten and destroyed? + logical(lk), optional, intent(in) :: overwrite_a + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state), optional, intent(out) :: err + + !> Local variables + type(linalg_state) :: err0 + integer(ilp) :: m,n,lda,ldu,ldv,info,k,lwork,lrwork,neig + logical(lk) :: copy_a + character :: task_u,task_v + ${rt}$, target :: work_dummy(1),u_dummy(1,1),v_dummy(1,1) + ${rt}$, allocatable :: work(:) + real(${rk}$), allocatable :: rwork(:) + ${rt}$, pointer :: amat(:,:),lreal(:),limag(:),umat(:,:),vmat(:,:) + + !> Matrix size + m = size(a,1,kind=ilp) + n = size(a,2,kind=ilp) + k = min(m,n) + neig = size(lambda,kind=ilp) + lda = m + + if (.not.(k>0 .and. m==n)) then + err0 = linalg_state(this,LINALG_VALUE_ERROR,'invalid or matrix size a=',[m,n], & + ', must be square.') + goto 1 + end if + + if (.not.neig>=k) then + err0 = linalg_state(this,LINALG_VALUE_ERROR,'eigenvalue array has insufficient size:',& + ' lambda=',neig,', n=',n) + goto 1 + endif + + + + ! Can A be overwritten? By default, do not overwrite + if (present(overwrite_a)) then + copy_a = .not.overwrite_a + else + copy_a = .true._lk + endif + + ! Initialize a matrix temporary + if (copy_a) then + allocate(amat(m,n),source=a) + else + amat => a + endif + + ! Decide if U, V eigenvectors + task_u = eigenvectors_flag(present(left)) + task_v = eigenvectors_flag(present(right)) + + if (present(right)) then + + #:if rt.startswith('complex') + vmat => right + #:else + allocate(vmat(n,n)) + #:endif + + if (size(vmat,2,kind=ilp) left + #:else + allocate(umat(n,n)) + #:endif + + if (size(umat,2,kind=ilp) Prepare working storage + lwork = nint(real(work_dummy(1),kind=${rk}$), kind=ilp) + allocate(work(lwork)) + + !> Compute eigensystem + call geev(task_u,task_v,n,amat,lda,& + #{if rt.startswith('complex')}#lambda,#{else}#lreal,limag,#{endif}# & + umat,ldu,vmat,ldv,& + work,lwork,#{if rt.startswith('complex')}#rwork,#{endif}#info) + call geev_info(err0,info,m,n) + + endif + + ! Finalize storage and process output flag + #:if not rt.startswith('complex') + lambda(:n) = cmplx(lreal(:n),limag(:n),kind=${rk}$) + + ! If the j-th and (j+1)-st eigenvalues form a complex conjugate pair, + ! geev returns reals as: + ! u(j) = VL(:,j) + i*VL(:,j+1) and + ! u(j+1) = VL(:,j) - i*VL(:,j+1). + ! Convert these to complex numbers here. + if (present(right)) call assign_real_eigenvectors_${rk}$(n,lambda,vmat,right) + if (present(left)) call assign_real_eigenvectors_${rk}$(n,lambda,umat,left) + #:endif + + 2 if (copy_a) deallocate(amat) + #:if not rt.startswith('complex') + if (present(right)) deallocate(vmat) + if (present(left)) deallocate(umat) + #:endif + 1 call linalg_error_handling(err0,err) + + end subroutine stdlib_linalg_eig_${ri}$ + + !> Return an array of eigenvalues of real symmetric / complex hermitian A + function stdlib_linalg_eigvalsh_${ri}$(a,upper_a,err) result(lambda) + !> Input matrix A[m,n] + ${rt}$, intent(in), target :: a(:,:) + !> [optional] Should the upper/lower half of A be used? Default: lower + logical(lk), optional, intent(in) :: upper_a + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state), intent(out) :: err + !> Array of singular values + real(${rk}$), allocatable :: lambda(:) + + ${rt}$, pointer :: amat(:,:) + integer(ilp) :: m,n,k + + !> Create an internal pointer so the intent of A won't affect the next call + amat => a + + m = size(a,1,kind=ilp) + n = size(a,2,kind=ilp) + k = min(m,n) + + !> Allocate return storage + allocate(lambda(k)) + + !> Compute eigenvalues only + call stdlib_linalg_eigh_${ri}$(amat,lambda,upper_a=upper_a,overwrite_a=.false.,err=err) + + end function stdlib_linalg_eigvalsh_${ri}$ + + !> Return an array of eigenvalues of real symmetric / complex hermitian A + function stdlib_linalg_eigvalsh_noerr_${ri}$(a,upper_a) result(lambda) + !> Input matrix A[m,n] + ${rt}$, intent(in), target :: a(:,:) + !> [optional] Should the upper/lower half of A be used? Default: lower + logical(lk), optional, intent(in) :: upper_a + !> Array of singular values + real(${rk}$), allocatable :: lambda(:) + + ${rt}$, pointer :: amat(:,:) + integer(ilp) :: m,n,k + + !> Create an internal pointer so the intent of A won't affect the next call + amat => a + + m = size(a,1,kind=ilp) + n = size(a,2,kind=ilp) + k = min(m,n) + + !> Allocate return storage + allocate(lambda(k)) + + !> Compute eigenvalues only + call stdlib_linalg_eigh_${ri}$(amat,lambda,upper_a=upper_a,overwrite_a=.false.) + + end function stdlib_linalg_eigvalsh_noerr_${ri}$ + + !> Eigendecomposition of a real symmetric or complex Hermitian matrix A returning an array `lambda` + !> of eigenvalues, and optionally right or left eigenvectors. + subroutine stdlib_linalg_eigh_${ri}$(a,lambda,vectors,upper_a,overwrite_a,err) + !> Input matrix A[m,n] + ${rt}$, intent(inout), target :: a(:,:) + !> Array of eigenvalues + real(${rk}$), intent(out) :: lambda(:) + !> The columns of vectors contain the orthonormal eigenvectors of A + ${rt}$, optional, intent(out), target :: vectors(:,:) + !> [optional] Can A data be overwritten and destroyed? + logical(lk), optional, intent(in) :: overwrite_a + !> [optional] Should the upper/lower half of A be used? Default: lower + logical(lk), optional, intent(in) :: upper_a + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state), optional, intent(out) :: err + + !> Local variables + type(linalg_state) :: err0 + integer(ilp) :: m,n,lda,info,k,lwork,neig + logical(lk) :: copy_a + character :: triangle,task + ${rt}$, target :: work_dummy(1) + ${rt}$, allocatable :: work(:) + real(${rk}$), allocatable :: rwork(:) + ${rt}$, pointer :: amat(:,:) + + !> Matrix size + m = size(a,1,kind=ilp) + n = size(a,2,kind=ilp) + k = min(m,n) + neig = size(lambda,kind=ilp) + + if (.not.(k>0 .and. m==n)) then + err0 = linalg_state(this,LINALG_VALUE_ERROR,'invalid or matrix size a=',[m,n], & + ', must be square.') + goto 1 + end if + + if (.not.neig>=k) then + err0 = linalg_state(this,LINALG_VALUE_ERROR,'eigenvalue array has insufficient size:',& + ' lambda=',neig,' must be >= n=',n) + goto 1 + endif + + ! Check if input A can be overwritten + if (present(vectors)) then + ! No need to copy A anyways + copy_a = .false. + elseif (present(overwrite_a)) then + copy_a = .not.overwrite_a + else + copy_a = .true._lk + endif + + ! Should we use the upper or lower half of the matrix? + triangle = symmetric_triangle(upper_a) + + ! Request for eigenvectors + task = eigenvectors_flag(present(vectors)) + + if (present(vectors)) then + + ! Check size + if (any(shape(vectors,kind=ilp) Prepare working storage + lwork = nint(real(work_dummy(1),kind=${rk}$), kind=ilp) + allocate(work(lwork)) + + !> Compute eigensystem + #:if rt.startswith('complex') + call heev(task,triangle,n,amat,lda,lambda,work,lwork,rwork,info) + #:else + call syev(task,triangle,n,amat,lda,lambda,work,lwork,info) + #:endif + call heev_info(err0,info,m,n) + + endif + + ! Finalize storage and process output flag + if (copy_a) deallocate(amat) + 1 call linalg_error_handling(err0,err) + + end subroutine stdlib_linalg_eigh_${ri}$ + + #:endfor + + #:for rk,rt,ri in REAL_KINDS_TYPES + !> GEEV for real matrices returns complex eigenvalues in real arrays. + !> Convert them to complex here, following the GEEV logic. + pure subroutine assign_real_eigenvectors_${rk}$(n,lambda,lmat,out_mat) + !> Problem size + integer(ilp), intent(in) :: n + !> Array of eigenvalues + complex(${rk}$), intent(in) :: lambda(:) + !> Real matrix as returned by geev + ${rt}$, intent(in) :: lmat(:,:) + !> Complex matrix as returned by eig + complex(${rk}$), intent(out) :: out_mat(:,:) + + integer(ilp) :: i,j + + ! Copy matrix + do concurrent(i=1:n,j=1:n) + out_mat(i,j) = lmat(i,j) + end do + + ! If the j-th and (j+1)-st eigenvalues form a complex conjugate pair, + ! geev returns them as reals as: + ! u(j) = VL(:,j) + i*VL(:,j+1) and + ! u(j+1) = VL(:,j) - i*VL(:,j+1). + ! Convert these to complex numbers here. + do j=1,n-1 + if (lambda(j)==conjg(lambda(j+1))) then + out_mat(:, j) = cmplx(lmat(:,j), lmat(:,j+1),kind=${rk}$) + out_mat(:,j+1) = cmplx(lmat(:,j),-lmat(:,j+1),kind=${rk}$) + endif + end do + + end subroutine assign_real_eigenvectors_${rk}$ + #:endfor + + +end module stdlib_linalg_eigenvalues From 430c89c07839e036cdf9c23f6d86a00de77a16c3 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 16 May 2024 09:42:20 +0200 Subject: [PATCH 02/32] remove `goto`s --- src/stdlib_linalg_eigenvalues.fypp | 271 +++++++++++++++-------------- 1 file changed, 141 insertions(+), 130 deletions(-) diff --git a/src/stdlib_linalg_eigenvalues.fypp b/src/stdlib_linalg_eigenvalues.fypp index 3b49ec968..76116c84e 100644 --- a/src/stdlib_linalg_eigenvalues.fypp +++ b/src/stdlib_linalg_eigenvalues.fypp @@ -1,5 +1,7 @@ #:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES module stdlib_linalg_eigenvalues +!! Compute eigenvalues and eigenvectors use stdlib_linalg_constants use stdlib_linalg_blas use stdlib_linalg_lapack @@ -26,26 +28,26 @@ module stdlib_linalg_eigenvalues ! Scipy: eigh(a, b=None, *, lower=True, eigvals_only=False, overwrite_a=False, overwrite_b=False, turbo=, eigvals=, type=1, check_finite=True, subset_by_index=None, subset_by_value=None, driver=None) interface eig - #:for rk,rt,ri in ALL_KINDS_TYPES + #:for rk,rt,ri in RC_KINDS_TYPES module procedure stdlib_linalg_eig_${ri}$ #:endfor end interface eig interface eigvals - #:for rk,rt,ri in ALL_KINDS_TYPES + #:for rk,rt,ri in RC_KINDS_TYPES module procedure stdlib_linalg_eigvals_${ri}$ module procedure stdlib_linalg_eigvals_noerr_${ri}$ #:endfor end interface eigvals interface eigh - #:for rk,rt,ri in ALL_KINDS_TYPES + #:for rk,rt,ri in RC_KINDS_TYPES module procedure stdlib_linalg_eigh_${ri}$ #:endfor end interface eigh interface eigvalsh - #:for rk,rt,ri in ALL_KINDS_TYPES + #:for rk,rt,ri in RC_KINDS_TYPES module procedure stdlib_linalg_eigvalsh_${ri}$ module procedure stdlib_linalg_eigvalsh_noerr_${ri}$ #:endfor @@ -58,26 +60,26 @@ module stdlib_linalg_eigenvalues contains !> Request for eigenvector calculation - elemental character function eigenvectors_flag(required) + elemental character function eigenvectors_task(required) logical(lk), intent(in) :: required - eigenvectors_flag = merge('V','N',required) - end function eigenvectors_flag + eigenvectors_task = merge('V','N',required) + end function eigenvectors_task !> Request for symmetry side (default: lower) - elemental character function symmetric_triangle(upper) + elemental character function symmetric_triangle_task(upper) logical(lk), optional, intent(in) :: upper if (present(upper)) then - symmetric_triangle = merge('U','L',upper) + symmetric_triangle_task = merge('U','L',upper) else - symmetric_triangle = 'L' + symmetric_triangle_task = 'L' endif - end function symmetric_triangle + end function symmetric_triangle_task !> Process GEEV output flags - elemental subroutine geev_info(err,info,m,n) + elemental subroutine handle_geev_info(err,info,m,n) !> Error handler - type(linalg_state), intent(inout) :: err - !> geev return flag + type(linalg_state_type), intent(inout) :: err + !> GEEV return flag integer(ilp), intent(in) :: info !> Input matrix size integer(ilp), intent(in) :: m,n @@ -87,31 +89,30 @@ module stdlib_linalg_eigenvalues ! Success! err%state = LINALG_SUCCESS case (-1) - err = linalg_state(this,LINALG_INTERNAL_ERROR,'Invalid task ID: left eigenvectors.') + err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid task ID: left eigenvectors.') case (-2) - err = linalg_state(this,LINALG_INTERNAL_ERROR,'Invalid task ID: right eigenvectors.') + err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid task ID: right eigenvectors.') case (-5,-3) - err = linalg_state(this,LINALG_VALUE_ERROR,'invalid matrix size: a=',[m,n]) + err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size: a=',[m,n]) case (-9) - err = linalg_state(this,LINALG_VALUE_ERROR,'insufficient left vector matrix size.') + err = linalg_state_type(this,LINALG_VALUE_ERROR,'insufficient left vector matrix size.') case (-11) - err = linalg_state(this,LINALG_VALUE_ERROR,'insufficient right vector matrix size.') + err = linalg_state_type(this,LINALG_VALUE_ERROR,'insufficient right vector matrix size.') case (-13) - err = linalg_state(this,LINALG_INTERNAL_ERROR,'Insufficient work array size.') + err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Insufficient work array size.') case (1:) - err = linalg_state(this,LINALG_ERROR,'Eigenvalue computation did not converge.') + err = linalg_state_type(this,LINALG_ERROR,'Eigenvalue computation did not converge.') case default - err = linalg_state(this,LINALG_INTERNAL_ERROR,'Unknown error returned by geev.') + err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Unknown error returned by geev.') end select - end subroutine geev_info - + end subroutine handle_geev_info !> Process SYEV/HEEV output flags - elemental subroutine heev_info(err,info,m,n) + elemental subroutine handle_heev_info(err,info,m,n) !> Error handler - type(linalg_state), intent(inout) :: err - !> geev return flag + type(linalg_state_type), intent(inout) :: err + !> SYEV/HEEV return flag integer(ilp), intent(in) :: info !> Input matrix size integer(ilp), intent(in) :: m,n @@ -121,29 +122,29 @@ module stdlib_linalg_eigenvalues ! Success! err%state = LINALG_SUCCESS case (-1) - err = linalg_state(this,LINALG_INTERNAL_ERROR,'Invalid eigenvector request.') + err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid eigenvector request.') case (-2) - err = linalg_state(this,LINALG_INTERNAL_ERROR,'Invalid triangular section request.') + err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Invalid triangular section request.') case (-5,-3) - err = linalg_state(this,LINALG_VALUE_ERROR,'invalid matrix size: a=',[m,n]) + err = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid matrix size: a=',[m,n]) case (-8) - err = linalg_state(this,LINALG_INTERNAL_ERROR,'insufficient workspace size.') + err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'insufficient workspace size.') case (1:) - err = linalg_state(this,LINALG_ERROR,'Eigenvalue computation did not converge.') + err = linalg_state_type(this,LINALG_ERROR,'Eigenvalue computation did not converge.') case default - err = linalg_state(this,LINALG_INTERNAL_ERROR,'Unknown error returned by syev/heev.') + err = linalg_state_type(this,LINALG_INTERNAL_ERROR,'Unknown error returned by syev/heev.') end select - end subroutine heev_info + end subroutine handle_heev_info - #:for rk,rt,ri in ALL_KINDS_TYPES + #:for rk,rt,ri in RC_KINDS_TYPES - !> Return an array of eigenvalues of matrix A. function stdlib_linalg_eigvals_${ri}$(a,err) result(lambda) + !! Return an array of eigenvalues of matrix A. !> Input matrix A[m,n] ${rt}$, intent(in), target :: a(:,:) !> [optional] state return flag. On error if not requested, the code will stop - type(linalg_state), intent(out) :: err + type(linalg_state_type), intent(out) :: err !> Array of singular values complex(${rk}$), allocatable :: lambda(:) @@ -166,8 +167,8 @@ module stdlib_linalg_eigenvalues end function stdlib_linalg_eigvals_${ri}$ - !> Return an array of eigenvalues of matrix A. function stdlib_linalg_eigvals_noerr_${ri}$(a) result(lambda) + !! Return an array of eigenvalues of matrix A. !> Input matrix A[m,n] ${rt}$, intent(in), target :: a(:,:) !> Array of singular values @@ -192,9 +193,9 @@ module stdlib_linalg_eigenvalues end function stdlib_linalg_eigvals_noerr_${ri}$ - !> Eigendecomposition of matrix A returning an array `lambda` of eigenvalues, - !> and optionally right or left eigenvectors. subroutine stdlib_linalg_eig_${ri}$(a,lambda,right,left,overwrite_a,err) + !! Eigendecomposition of matrix A returning an array `lambda` of eigenvalues, + !! and optionally right or left eigenvectors. !> Input matrix A[m,n] ${rt}$, intent(inout), target :: a(:,:) !> Array of eigenvalues @@ -206,10 +207,10 @@ module stdlib_linalg_eigenvalues !> [optional] Can A data be overwritten and destroyed? logical(lk), optional, intent(in) :: overwrite_a !> [optional] state return flag. On error if not requested, the code will stop - type(linalg_state), optional, intent(out) :: err + type(linalg_state_type), optional, intent(out) :: err !> Local variables - type(linalg_state) :: err0 + type(linalg_state_type) :: err0 integer(ilp) :: m,n,lda,ldu,ldv,info,k,lwork,lrwork,neig logical(lk) :: copy_a character :: task_u,task_v @@ -226,19 +227,17 @@ module stdlib_linalg_eigenvalues lda = m if (.not.(k>0 .and. m==n)) then - err0 = linalg_state(this,LINALG_VALUE_ERROR,'invalid or matrix size a=',[m,n], & + err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid or matrix size a=',[m,n], & ', must be square.') - goto 1 - end if - - if (.not.neig>=k) then - err0 = linalg_state(this,LINALG_VALUE_ERROR,'eigenvalue array has insufficient size:',& + call linalg_error_handling(err0,err) + return + elseif (.not.neig>=k) then + err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'eigenvalue array has insufficient size:',& ' lambda=',neig,', n=',n) - goto 1 + call linalg_error_handling(err0,err) + return endif - - ! Can A be overwritten? By default, do not overwrite if (present(overwrite_a)) then copy_a = .not.overwrite_a @@ -254,22 +253,25 @@ module stdlib_linalg_eigenvalues endif ! Decide if U, V eigenvectors - task_u = eigenvectors_flag(present(left)) - task_v = eigenvectors_flag(present(right)) + task_u = eigenvectors_task(present(left)) + task_v = eigenvectors_task(present(right)) if (present(right)) then - + #:if rt.startswith('complex') + ! For a complex matrix, GEEV returns complex arrays. + ! Point directly to output. vmat => right #:else + ! For a real matrix, GEEV returns real arrays. + ! Allocate temporary reals, will be converted to complex vectors at the end. allocate(vmat(n,n)) #:endif if (size(vmat,2,kind=ilp) left #:else + ! For a real matrix, GEEV returns real arrays. + ! Allocate temporary reals, will be converted to complex vectors at the end. allocate(umat(n,n)) - #:endif - + #:endif + if (size(umat,2,kind=ilp) Prepare working storage - lwork = nint(real(work_dummy(1),kind=${rk}$), kind=ilp) - allocate(work(lwork)) - - !> Compute eigensystem - call geev(task_u,task_v,n,amat,lda,& - #{if rt.startswith('complex')}#lambda,#{else}#lreal,limag,#{endif}# & - umat,ldu,vmat,ldv,& - work,lwork,#{if rt.startswith('complex')}#rwork,#{endif}#info) - call geev_info(err0,info,m,n) + ! Compute workspace size + #:if rt.startswith('complex') + allocate(rwork(2*n)) + #:else + allocate(lreal(n),limag(n)) + #:endif - endif - - ! Finalize storage and process output flag - #:if not rt.startswith('complex') - lambda(:n) = cmplx(lreal(:n),limag(:n),kind=${rk}$) + lwork = -1_ilp + + call geev(task_u,task_v,n,amat,lda,& + #{if rt.startswith('complex')}#lambda,#{else}#lreal,limag,#{endif}# & + umat,ldu,vmat,ldv,& + work_dummy,lwork,#{if rt.startswith('complex')}#rwork,#{endif}#info) + call handle_geev_info(err0,info,m,n) + + ! Compute eigenvalues + if (info==0) then + + !> Prepare working storage + lwork = nint(real(work_dummy(1),kind=${rk}$), kind=ilp) + allocate(work(lwork)) + + !> Compute eigensystem + call geev(task_u,task_v,n,amat,lda,& + #{if rt.startswith('complex')}#lambda,#{else}#lreal,limag,#{endif}# & + umat,ldu,vmat,ldv,& + work,lwork,#{if rt.startswith('complex')}#rwork,#{endif}#info) + call handle_geev_info(err0,info,m,n) + + endif + + ! Finalize storage and process output flag + #:if not rt.startswith('complex') + lambda(:n) = cmplx(lreal(:n),limag(:n),kind=${rk}$) + + ! If the j-th and (j+1)-st eigenvalues form a complex conjugate pair, + ! geev returns reals as: + ! u(j) = VL(:,j) + i*VL(:,j+1) and + ! u(j+1) = VL(:,j) - i*VL(:,j+1). + ! Convert these to complex numbers here. + if (present(right)) call assign_real_eigenvectors_${rk}$(n,lambda,vmat,right) + if (present(left)) call assign_real_eigenvectors_${rk}$(n,lambda,umat,left) + #:endif - ! If the j-th and (j+1)-st eigenvalues form a complex conjugate pair, - ! geev returns reals as: - ! u(j) = VL(:,j) + i*VL(:,j+1) and - ! u(j+1) = VL(:,j) - i*VL(:,j+1). - ! Convert these to complex numbers here. - if (present(right)) call assign_real_eigenvectors_${rk}$(n,lambda,vmat,right) - if (present(left)) call assign_real_eigenvectors_${rk}$(n,lambda,umat,left) - #:endif + endif get_geev - 2 if (copy_a) deallocate(amat) + if (copy_a) deallocate(amat) #:if not rt.startswith('complex') if (present(right)) deallocate(vmat) if (present(left)) deallocate(umat) #:endif - 1 call linalg_error_handling(err0,err) + call linalg_error_handling(err0,err) end subroutine stdlib_linalg_eig_${ri}$ - !> Return an array of eigenvalues of real symmetric / complex hermitian A function stdlib_linalg_eigvalsh_${ri}$(a,upper_a,err) result(lambda) + !! Return an array of eigenvalues of real symmetric / complex hermitian A !> Input matrix A[m,n] ${rt}$, intent(in), target :: a(:,:) !> [optional] Should the upper/lower half of A be used? Default: lower logical(lk), optional, intent(in) :: upper_a !> [optional] state return flag. On error if not requested, the code will stop - type(linalg_state), intent(out) :: err + type(linalg_state_type), intent(out) :: err !> Array of singular values real(${rk}$), allocatable :: lambda(:) @@ -379,9 +388,9 @@ module stdlib_linalg_eigenvalues call stdlib_linalg_eigh_${ri}$(amat,lambda,upper_a=upper_a,overwrite_a=.false.,err=err) end function stdlib_linalg_eigvalsh_${ri}$ - - !> Return an array of eigenvalues of real symmetric / complex hermitian A + function stdlib_linalg_eigvalsh_noerr_${ri}$(a,upper_a) result(lambda) + !! Return an array of eigenvalues of real symmetric / complex hermitian A !> Input matrix A[m,n] ${rt}$, intent(in), target :: a(:,:) !> [optional] Should the upper/lower half of A be used? Default: lower @@ -407,9 +416,9 @@ module stdlib_linalg_eigenvalues end function stdlib_linalg_eigvalsh_noerr_${ri}$ - !> Eigendecomposition of a real symmetric or complex Hermitian matrix A returning an array `lambda` - !> of eigenvalues, and optionally right or left eigenvectors. subroutine stdlib_linalg_eigh_${ri}$(a,lambda,vectors,upper_a,overwrite_a,err) + !! Eigendecomposition of a real symmetric or complex Hermitian matrix A returning an array `lambda` + !! of eigenvalues, and optionally right or left eigenvectors. !> Input matrix A[m,n] ${rt}$, intent(inout), target :: a(:,:) !> Array of eigenvalues @@ -421,10 +430,10 @@ module stdlib_linalg_eigenvalues !> [optional] Should the upper/lower half of A be used? Default: lower logical(lk), optional, intent(in) :: upper_a !> [optional] state return flag. On error if not requested, the code will stop - type(linalg_state), optional, intent(out) :: err + type(linalg_state_type), optional, intent(out) :: err !> Local variables - type(linalg_state) :: err0 + type(linalg_state_type) :: err0 integer(ilp) :: m,n,lda,info,k,lwork,neig logical(lk) :: copy_a character :: triangle,task @@ -440,15 +449,15 @@ module stdlib_linalg_eigenvalues neig = size(lambda,kind=ilp) if (.not.(k>0 .and. m==n)) then - err0 = linalg_state(this,LINALG_VALUE_ERROR,'invalid or matrix size a=',[m,n], & + err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid or matrix size a=',[m,n], & ', must be square.') - goto 1 - end if - - if (.not.neig>=k) then - err0 = linalg_state(this,LINALG_VALUE_ERROR,'eigenvalue array has insufficient size:',& + call linalg_error_handling(err0,err) + return + elseif (.not.neig>=k) then + err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'eigenvalue array has insufficient size:',& ' lambda=',neig,' must be >= n=',n) - goto 1 + call linalg_error_handling(err0,err) + return endif ! Check if input A can be overwritten @@ -462,19 +471,20 @@ module stdlib_linalg_eigenvalues endif ! Should we use the upper or lower half of the matrix? - triangle = symmetric_triangle(upper_a) + triangle = symmetric_triangle_task(upper_a) ! Request for eigenvectors - task = eigenvectors_flag(present(vectors)) + task = eigenvectors_task(present(vectors)) if (present(vectors)) then ! Check size if (any(shape(vectors,kind=ilp), eigvals=, type=1, check_finite=True, subset_by_index=None, subset_by_value=None, driver=None) - - interface eig - #:for rk,rt,ri in RC_KINDS_TYPES - #:if rk!="xdp" - module procedure stdlib_linalg_eig_${ri}$ - #:endif - #:endfor - end interface eig - - interface eigvals - #:for rk,rt,ri in RC_KINDS_TYPES - #:if rk!="xdp" - module procedure stdlib_linalg_eigvals_${ri}$ - module procedure stdlib_linalg_eigvals_noerr_${ri}$ - #:endif - #:endfor - end interface eigvals - interface eigh - #:for rk,rt,ri in RC_KINDS_TYPES - #:if rk!="xdp" - module procedure stdlib_linalg_eigh_${ri}$ - #:endif - #:endfor - end interface eigh - - interface eigvalsh - #:for rk,rt,ri in RC_KINDS_TYPES - #:if rk!="xdp" - module procedure stdlib_linalg_eigvalsh_${ri}$ - module procedure stdlib_linalg_eigvalsh_noerr_${ri}$ - #:endif - #:endfor - end interface eigvalsh - - character(*), parameter :: this = 'eigenvalues' - contains !> Request for eigenvector calculation @@ -236,13 +181,14 @@ module stdlib_linalg_eigenvalues lda = m if (.not.(k>0 .and. m==n)) then - err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'invalid or matrix size a=',[m,n], & - ', must be square.') + err0 = linalg_state_type(this,LINALG_VALUE_ERROR,& + 'invalid or matrix size a=',[m,n],', must be square.') call linalg_error_handling(err0,err) return elseif (.not.neig>=k) then - err0 = linalg_state_type(this,LINALG_VALUE_ERROR,'eigenvalue array has insufficient size:',& - ' lambda=',neig,', n=',n) + err0 = linalg_state_type(this,LINALG_VALUE_ERROR,& + 'eigenvalue array has insufficient size:',& + ' lambda=',neig,', n=',n) call linalg_error_handling(err0,err) return endif @@ -590,4 +536,4 @@ module stdlib_linalg_eigenvalues #:endfor -end module stdlib_linalg_eigenvalues +end submodule stdlib_linalg_eigenvalues From 2d58c682ef7fabe0044043611b1f39bf3b627906 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 16 May 2024 10:11:45 +0200 Subject: [PATCH 05/32] add tests --- src/stdlib_linalg.fypp | 4 + test/linalg/CMakeLists.txt | 2 + test/linalg/test_linalg_eigenvalues.fypp | 208 +++++++++++++++++++++++ 3 files changed, 214 insertions(+) create mode 100644 test/linalg/test_linalg_eigenvalues.fypp diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 3de02fef7..c7f1b80d6 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -19,6 +19,10 @@ module stdlib_linalg public :: det public :: operator(.det.) public :: diag + public :: eig + public :: eigh + public :: eigvals + public :: eigvalsh public :: eye public :: lstsq public :: lstsq_space diff --git a/test/linalg/CMakeLists.txt b/test/linalg/CMakeLists.txt index 9ec242213..54eb1a545 100644 --- a/test/linalg/CMakeLists.txt +++ b/test/linalg/CMakeLists.txt @@ -2,6 +2,7 @@ set( fppFiles "test_linalg.fypp" "test_blas_lapack.fypp" + "test_linalg_eigenvalues.fypp" "test_linalg_solve.fypp" "test_linalg_lstsq.fypp" "test_linalg_determinant.fypp" @@ -11,6 +12,7 @@ fypp_f90("${fyppFlags}" "${fppFiles}" outFiles) ADDTEST(linalg) ADDTEST(linalg_determinant) +ADDTEST(linalg_eigenvalues) ADDTEST(linalg_matrix_property_checks) ADDTEST(linalg_solve) ADDTEST(linalg_lstsq) diff --git a/test/linalg/test_linalg_eigenvalues.fypp b/test/linalg/test_linalg_eigenvalues.fypp new file mode 100644 index 000000000..3af893617 --- /dev/null +++ b/test/linalg/test_linalg_eigenvalues.fypp @@ -0,0 +1,208 @@ +#:include "common.fypp" +! Test eigenvalues and eigendecompositions +module test_linalg_eigenvalues + use stdlib_linalg_constants + use stdlib_linalg_state + use stdlib_linalg, only: eig, eigh, eigvals, eigvalsh, diag + use testdrive, only: error_type, check, new_unittest, unittest_type + + implicit none (type,external) + private + + public :: test_eig_eigh + + contains + + !> SVD tests + subroutine test_eig_eigh(tests) + !> Collection of tests + type(unittest_type), allocatable, intent(out) :: tests(:) + + allocate(tests(0)) + + #:for rk,rt,ri in REAL_KINDS_TYPES + #:if rk!="xdp" + tests = [tests,new_unittest("test_eig_real_${ri}$",test_eig_real_${ri}$), & + new_unittest("test_eigh_real_${ri}$",test_eigh_real_${ri}$)] + #:endif + #: endfor + + #:for ck,ct,ci in CMPLX_KINDS_TYPES + tests = [tests,new_unittest("test_eig_complex_${ci}$",test_eig_complex_${ci}$)] + #: endfor + + end subroutine test_eig_eigh + + !> Simple real matrix eigenvalues + #:for rk,rt,ri in REAL_KINDS_TYPES + #:if rk!="xdp" + subroutine test_eig_real_${ri}$(error) + type(error_type), allocatable, intent(out) :: error + + !> Reference solution + real(${rk}$), parameter :: zero = 0.0_${rk}$ + real(${rk}$), parameter :: two = 2.0_${rk}$ + real(${rk}$), parameter :: sqrt2o2 = sqrt(two)*0.5_${rk}$ + real(${rk}$), parameter :: tol = sqrt(epsilon(zero)) + + !> Local variables + type(linalg_state_type) :: state + ${rt}$ :: A(3,3),B(2,2) + complex(${rk}$) :: lambda(3),Bvec(2,2),Bres(2,2) + + !> Matrix with real eigenvalues + A = reshape([1,0,0, & + 0,2,0, & + 0,0,3],[3,3]) + + call eig(A,lambda,err=state) + + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + call check(error, all(aimag(lambda)==zero.and.real(lambda,kind=${rk}$)==[1,2,3]),'expected results') + if (allocated(error)) return + + !> Matrix with complex eigenvalues + B = transpose(reshape([1, -1, & + 1, 1],[2,2])) + + !> Expected right eigenvectors + Bres(1,1:2) = sqrt2o2 + Bres(2,1) = cmplx(zero,-sqrt2o2,kind=${rk}$) + Bres(2,2) = cmplx(zero,+sqrt2o2,kind=${rk}$) + + call eig(B,lambda,right=Bvec,err=state) + + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + call check(error, all(abs(Bres-Bvec)<=tol),'expected results') + if (allocated(error)) return + + end subroutine test_eig_real_${ri}$ + + ! Symmetric matrix eigenvalues + subroutine test_eigh_real_${ri}$(error) + type(error_type), allocatable, intent(out) :: error + + !> Reference solution + real(${rk}$), parameter :: zero = 0.0_${rk}$ + real(${rk}$), parameter :: tol = sqrt(epsilon(zero)) + real(${rk}$), parameter :: A(4,4) = reshape([6,3,1,5, & + 3,0,5,1, & + 1,5,6,2, & + 5,1,2,2],[4,4]) + + !> Local variables + real(${rk}$) :: Amat(4,4),lambda(4),vect(4,4),Av(4,4),lv(4,4) + type(linalg_state_type) :: state + + Amat = A + + call eigh(Amat,lambda,vect,err=state) + + Av = matmul(A,vect) + lv = matmul(vect,diag(lambda)) + + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + call check(error, all(abs(Av-lv)<=tol*abs(Av)),'expected results') + if (allocated(error)) return + + !> Test functional versions: no state interface + lambda = eigvalsh(Amat) + + !> State interface + lambda = eigvalsh(Amat,err=state) + + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + !> Functional version, lower A + Amat = A + lambda = eigvalsh(Amat,upper_a=.false.,err=state) + + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + end subroutine test_eigh_real_${ri}$ + + #:endif + #:endfor + + !> Simple complex matrix eigenvalues + #:for ck,ct,ci in CMPLX_KINDS_TYPES + #:if ck!="xdp" + subroutine test_eig_complex_${ci}$(error) + type(error_type), allocatable, intent(out) :: error + + !> Reference solution + real(${ck}$), parameter :: zero = 0.0_${ck}$ + real(${ck}$), parameter :: two = 2.0_${ck}$ + real(${ck}$), parameter :: sqrt2o2 = sqrt(two)*0.5_${ck}$ + real(${ck}$), parameter :: tol = sqrt(epsilon(zero)) + ${ct}$, parameter :: cone = (1.0_${ck}$,0.0_${ck}$) + ${ct}$, parameter :: cimg = (0.0_${ck}$,1.0_${ck}$) + ${ct}$, parameter :: czero = (0.0_${ck}$,0.0_${ck}$) + + !> Local vaciables + type(linalg_state_type) :: state + ${ct}$ :: A(2,2),lambda(2),Avec(2,2),Ares(2,2),lres(2) + + !> Matcix with real eigenvalues + A = transpose(reshape([ cone, cimg, & + -cimg, cone], [2,2])) + + call eig(A,lambda,right=Avec,err=state) + + !> Expected eigenvalues and eigenvectors + lres(1) = two + lres(2) = zero + + !> Eigenvectors may vary: do not use for error + Ares(1,1) = cmplx(zero,sqrt2o2,kind=${ck}$) + Ares(1,2) = cmplx(sqrt2o2,zero,kind=${ck}$) + Ares(2,1) = cmplx(sqrt2o2,zero,kind=${ck}$) + Ares(2,2) = cmplx(zero,sqrt2o2,kind=${ck}$) + + call check(error,state%ok(),state%print()) + if (allocated(error)) return + + call check(error, all(abs(lambda-lres)<=tol), 'results match expected') + if (allocated(error)) return + + end subroutine test_eig_complex_${ci}$ + + #:endif + #:endfor + + +end module test_linalg_eigenvalues + +program test_eigenvalues + use, intrinsic :: iso_fortran_env, only : error_unit + use testdrive, only : run_testsuite, new_testsuite, testsuite_type + use test_linalg_eigenvalues, only : test_eig_eigh + implicit none + integer :: stat, is + type(testsuite_type), allocatable :: testsuites(:) + character(len=*), parameter :: fmt = '("#", *(1x, a))' + + stat = 0 + + testsuites = [ & + new_testsuite("linalg_eigenvalues", test_eig_eigh) & + ] + + do is = 1, size(testsuites) + write(error_unit, fmt) "Testing:", testsuites(is)%name + call run_testsuite(testsuites(is)%collect, error_unit, stat) + end do + + if (stat > 0) then + write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!" + error stop + end if +end program test_eigenvalues From 15b73fc3090978dec8d76a304c330041e59dfe21 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 16 May 2024 10:31:44 +0200 Subject: [PATCH 06/32] test: remove `xdp` --- test/linalg/test_linalg_eigenvalues.fypp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/linalg/test_linalg_eigenvalues.fypp b/test/linalg/test_linalg_eigenvalues.fypp index 3af893617..ecec0be3b 100644 --- a/test/linalg/test_linalg_eigenvalues.fypp +++ b/test/linalg/test_linalg_eigenvalues.fypp @@ -28,7 +28,9 @@ module test_linalg_eigenvalues #: endfor #:for ck,ct,ci in CMPLX_KINDS_TYPES + #:if ck!="xdp" tests = [tests,new_unittest("test_eig_complex_${ci}$",test_eig_complex_${ci}$)] + #:endif #: endfor end subroutine test_eig_eigh From b5f7f2125a7e5e2e156c480d589f87ac00bf744c Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 16 May 2024 10:58:02 +0200 Subject: [PATCH 07/32] add examples: `eig`, `eigh`, `eigvals`, `eigvalsh` --- example/linalg/CMakeLists.txt | 4 +++ example/linalg/example_eig.f90 | 27 ++++++++++++++++++++ example/linalg/example_eigh.f90 | 39 +++++++++++++++++++++++++++++ example/linalg/example_eigvals.f90 | 25 ++++++++++++++++++ example/linalg/example_eigvalsh.f90 | 26 +++++++++++++++++++ 5 files changed, 121 insertions(+) create mode 100644 example/linalg/example_eig.f90 create mode 100644 example/linalg/example_eigh.f90 create mode 100644 example/linalg/example_eigvals.f90 create mode 100644 example/linalg/example_eigvalsh.f90 diff --git a/example/linalg/CMakeLists.txt b/example/linalg/CMakeLists.txt index 2d80f067c..ebbb0f387 100644 --- a/example/linalg/CMakeLists.txt +++ b/example/linalg/CMakeLists.txt @@ -13,6 +13,10 @@ ADD_EXAMPLE(is_square) ADD_EXAMPLE(is_symmetric) ADD_EXAMPLE(is_triangular) ADD_EXAMPLE(outer_product) +ADD_EXAMPLE(eig) +ADD_EXAMPLE(eigh) +ADD_EXAMPLE(eigvals) +ADD_EXAMPLE(eigvalsh) ADD_EXAMPLE(trace) ADD_EXAMPLE(state1) ADD_EXAMPLE(state2) diff --git a/example/linalg/example_eig.f90 b/example/linalg/example_eig.f90 new file mode 100644 index 000000000..beb1ef271 --- /dev/null +++ b/example/linalg/example_eig.f90 @@ -0,0 +1,27 @@ +! Eigendecomposition of a real square matrix +program example_eig + use stdlib_linalg_constants, only: sp + use stdlib_linalg, only: eig + implicit none + + integer :: i + real(sp), allocatable :: A(:,:) + complex(sp), allocatable :: lambda(:),vectors(:,:) + + ! Decomposition of this square matrix + ! NB Fortran is column-major -> transpose input + A = transpose(reshape( [ [2, 2, 4], & + [1, 3, 5], & + [2, 3, 4] ], [3,3] )) + + ! Get eigenvalues and right eigenvectors + allocate(lambda(3),vectors(3,3)) + + call eig(A, lambda, right=vectors) + + do i=1,3 + print *, 'eigenvalue ',i,': ',lambda(i) + print *, 'eigenvector ',i,': ',vectors(:,i) + end do + +end program example_eig diff --git a/example/linalg/example_eigh.f90 b/example/linalg/example_eigh.f90 new file mode 100644 index 000000000..1cf30138d --- /dev/null +++ b/example/linalg/example_eigh.f90 @@ -0,0 +1,39 @@ +! Eigendecomposition of a real symmetric matrix +program example_eigh + use stdlib_linalg_constants, only: sp + use stdlib_linalg, only: eigh + implicit none + + integer :: i + real(sp), allocatable :: A(:,:),lambda(:),v(:,:) + complex(sp), allocatable :: cA(:,:),cv(:,:) + + ! Decomposition of this symmetric matrix + ! NB Fortran is column-major -> transpose input + A = transpose(reshape( [ [2, 1, 4], & + [1, 3, 5], & + [4, 5, 4] ], [3,3] )) + + ! Note: real symmetric matrices have real (orthogonal) eigenvalues and eigenvectors + allocate(lambda(3),v(3,3)) + call eigh(A, lambda, vectors=v) + + print *, 'Real matrix' + do i=1,3 + print *, 'eigenvalue ',i,': ',lambda(i) + print *, 'eigenvector ',i,': ',v(:,i) + end do + + ! Complex hermitian matrices have real (orthogonal) eigenvalues and complex eigenvectors + cA = A + + allocate(cv(3,3)) + call eigh(cA, lambda, vectors=cv) + + print *, 'Complex matrix' + do i=1,3 + print *, 'eigenvalue ',i,': ',lambda(i) + print *, 'eigenvector ',i,': ',cv(:,i) + end do + +end program example_eigh diff --git a/example/linalg/example_eigvals.f90 b/example/linalg/example_eigvals.f90 new file mode 100644 index 000000000..973d50666 --- /dev/null +++ b/example/linalg/example_eigvals.f90 @@ -0,0 +1,25 @@ +! Eigenvalues of a general real / complex matrix +program example_eigvals + use stdlib_linalg_constants, only: sp + use stdlib_linalg, only: eigvals + implicit none + + integer :: i + real(sp), allocatable :: A(:,:),lambda(:) + complex(sp), allocatable :: cA(:,:),clambda(:) + + ! NB Fortran is column-major -> transpose input + A = transpose(reshape( [ [2, 8, 4], & + [1, 3, 5], & + [9, 5,-2] ], [3,3] )) + + ! Note: real symmetric matrix + lambda = eigvals(A) + print *, 'Real matrix eigenvalues: ',lambda + + ! Complex general matrix + cA = cmplx(A, -2*A, kind=sp) + clambda = eigvals(cA) + print *, 'Complex matrix eigenvalues: ',clambda + +end program example_eigvals diff --git a/example/linalg/example_eigvalsh.f90 b/example/linalg/example_eigvalsh.f90 new file mode 100644 index 000000000..6853d9051 --- /dev/null +++ b/example/linalg/example_eigvalsh.f90 @@ -0,0 +1,26 @@ +! Eigenvalues of a real symmetric / complex hermitian matrix +program example_eigvalsh + use stdlib_linalg_constants, only: sp + use stdlib_linalg, only: eigvalsh + implicit none + + integer :: i + real(sp), allocatable :: A(:,:),lambda(:) + complex(sp), allocatable :: cA(:,:) + + ! Decomposition of this symmetric matrix + ! NB Fortran is column-major -> transpose input + A = transpose(reshape( [ [2, 1, 4], & + [1, 3, 5], & + [4, 5, 4] ], [3,3] )) + + ! Note: real symmetric matrices have real (orthogonal) eigenvalues and eigenvectors + lambda = eigvalsh(A) + print *, 'Symmetric matrix eigenvalues: ',lambda + + ! Complex hermitian matrices have real (orthogonal) eigenvalues and complex eigenvectors + cA = A + lambda = eigvalsh(cA) + print *, 'Hermitian matrix eigenvalues: ',lambda + +end program example_eigvalsh From 892308ff7ea439b494fcbd859c4314966842fcc4 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 16 May 2024 10:59:51 +0200 Subject: [PATCH 08/32] fix: `module` procedures --- src/stdlib_linalg_eigenvalues.fypp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/stdlib_linalg_eigenvalues.fypp b/src/stdlib_linalg_eigenvalues.fypp index 53c9abd41..db5282904 100644 --- a/src/stdlib_linalg_eigenvalues.fypp +++ b/src/stdlib_linalg_eigenvalues.fypp @@ -93,7 +93,7 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues #:for rk,rt,ri in RC_KINDS_TYPES #:if rk!="xdp" - function stdlib_linalg_eigvals_${ri}$(a,err) result(lambda) + module function stdlib_linalg_eigvals_${ri}$(a,err) result(lambda) !! Return an array of eigenvalues of matrix A. !> Input matrix A[m,n] ${rt}$, intent(in), target :: a(:,:) @@ -121,7 +121,7 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues end function stdlib_linalg_eigvals_${ri}$ - function stdlib_linalg_eigvals_noerr_${ri}$(a) result(lambda) + module function stdlib_linalg_eigvals_noerr_${ri}$(a) result(lambda) !! Return an array of eigenvalues of matrix A. !> Input matrix A[m,n] ${rt}$, intent(in), target :: a(:,:) @@ -147,7 +147,7 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues end function stdlib_linalg_eigvals_noerr_${ri}$ - subroutine stdlib_linalg_eig_${ri}$(a,lambda,right,left,overwrite_a,err) + module subroutine stdlib_linalg_eig_${ri}$(a,lambda,right,left,overwrite_a,err) !! Eigendecomposition of matrix A returning an array `lambda` of eigenvalues, !! and optionally right or left eigenvectors. !> Input matrix A[m,n] @@ -315,7 +315,7 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues end subroutine stdlib_linalg_eig_${ri}$ - function stdlib_linalg_eigvalsh_${ri}$(a,upper_a,err) result(lambda) + module function stdlib_linalg_eigvalsh_${ri}$(a,upper_a,err) result(lambda) !! Return an array of eigenvalues of real symmetric / complex hermitian A !> Input matrix A[m,n] ${rt}$, intent(in), target :: a(:,:) @@ -344,7 +344,7 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues end function stdlib_linalg_eigvalsh_${ri}$ - function stdlib_linalg_eigvalsh_noerr_${ri}$(a,upper_a) result(lambda) + module function stdlib_linalg_eigvalsh_noerr_${ri}$(a,upper_a) result(lambda) !! Return an array of eigenvalues of real symmetric / complex hermitian A !> Input matrix A[m,n] ${rt}$, intent(in), target :: a(:,:) @@ -371,7 +371,7 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues end function stdlib_linalg_eigvalsh_noerr_${ri}$ - subroutine stdlib_linalg_eigh_${ri}$(a,lambda,vectors,upper_a,overwrite_a,err) + module subroutine stdlib_linalg_eigh_${ri}$(a,lambda,vectors,upper_a,overwrite_a,err) !! Eigendecomposition of a real symmetric or complex Hermitian matrix A returning an array `lambda` !! of eigenvalues, and optionally right or left eigenvectors. !> Input matrix A[m,n] From 1ebf374111365804a60b70fec83ce245511f1873 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 16 May 2024 11:24:57 +0200 Subject: [PATCH 09/32] add documentation --- src/stdlib_linalg.fypp | 90 ++++++++++++++++++++++++++++-- src/stdlib_linalg_eigenvalues.fypp | 16 +++--- 2 files changed, 93 insertions(+), 13 deletions(-) diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index c7f1b80d6..90b920395 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -556,8 +556,27 @@ module stdlib_linalg #:endfor end interface - interface eig - !! Eigendecomposition of a square matrix: return eigenvalues, and optionally eigenvectors + ! Eigendecomposition of a square matrix: eigenvalues, and optionally eigenvectors + interface eig + !! version: experimental + !! + !! Solves the eigendecomposition \( A \cdot \bar{v} - \lambda \cdot \bar{v} \) for square matrix \( A \). + !! ([Specification](../page/specs/stdlib_linalg.html#eig-eigenvalues-and-eigenvectors-of-a-square-matrix)) + !! + !!### Summary + !! Subroutine interface for computing eigenvalues and eigenvectors of a square matrix. + !! + !!### Description + !! + !! This interface provides methods for computing the eigenvalues, and optionally eigenvectors, + !! of a general square matrix. Supported data types include `real` and `complex`, and no assumption is + !! made on the matrix structure. The user may request either left, right, or both + !! eigenvectors to be returned. They are returned as columns of a square matrix with the same size as `A`. + !! Preallocated space for both eigenvalues `lambda` and the eigenvector matrices must be user-provided. + !! + !!@note The solution is based on LAPACK's general eigenproblem solvers `*GEEV`. + !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``). + !! #:for rk,rt,ri in RC_KINDS_TYPES #:if rk!="xdp" module subroutine stdlib_linalg_eig_${ri}$(a,lambda,right,left,overwrite_a,err) @@ -580,8 +599,26 @@ module stdlib_linalg #:endfor end interface eig + ! Eigenvalues of a square matrix interface eigvals - !! Eigenvalues of a square matrix + !! version: experimental + !! + !! Returns the eigenvalues \( lambda \), \( A \cdot \bar{v} - \lambda \cdot \bar{v} \), for square matrix \( A \). + !! ([Specification](../page/specs/stdlib_linalg.html#eigvals-eigenvalues-of-a-square-matrix)) + !! + !!### Summary + !! Function interface for computing the eigenvalues of a square matrix. + !! + !!### Description + !! + !! This interface provides functions for returning the eigenvalues of a general square matrix. + !! Supported data types include `real` and `complex`, and no assumption is made on the matrix structure. + !! An `error stop` is thrown in case of failure; otherwise, error information can be returned + !! as an optional `type(linalg_state_type)` output flag. + !! + !!@note The solution is based on LAPACK's general eigenproblem solvers `*GEEV`. + !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``). + !! #:for rk,rt,ri in RC_KINDS_TYPES #:if rk!="xdp" module function stdlib_linalg_eigvals_${ri}$(a,err) result(lambda) @@ -605,8 +642,30 @@ module stdlib_linalg #:endfor end interface eigvals + ! Eigendecomposition of a real symmetric or complex hermitian matrix interface eigh - !! Eigendecomposition of a real symmetric or complex hermitian matrix + !! version: experimental + !! + !! Solves the eigendecomposition \( A \cdot \bar{v} - \lambda \cdot \bar{v} \) for a real symmetric + !! \( A = A^T \) or complex Hermitian \( A = A^H \) square matrix. + !! ([Specification](../page/specs/stdlib_linalg.html#eigh-eigenvalues-and-eigenvectors-of-a-real-symmetric-or-complex-hermitian-square-matrix)) + !! + !!### Summary + !! Subroutine interface for computing eigenvalues and eigenvectors of a real symmetric or complex Hermitian square matrix. + !! + !!### Description + !! + !! This interface provides methods for computing the eigenvalues, and optionally eigenvectors, + !! of a real symmetric or complex Hermitian square matrix. Supported data types include `real` and `complex`. + !! The matrix must be symmetric (if `real`) or Hermitian (if `complex`). Only the lower or upper + !! half of the matrix is accessed, and the user can select which using the optional `upper_a` + !! flag (default: use lower half). The vectors are orthogonal, and may be returned as columns of an optional + !! matrix `vectors` with the same kind and size as `A`. + !! Preallocated space for both eigenvalues `lambda` and the eigenvector matrix must be user-provided. + !! + !!@note The solution is based on LAPACK's eigenproblem solvers `*SYEV`/`*HEEV`. + !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``). + !! #:for rk,rt,ri in RC_KINDS_TYPES #:if rk!="xdp" module subroutine stdlib_linalg_eigh_${ri}$(a,lambda,vectors,upper_a,overwrite_a,err) @@ -629,8 +688,29 @@ module stdlib_linalg #:endfor end interface eigh + ! Eigenvalues of a real symmetric or complex hermitian matrix interface eigvalsh - !! Eigenvalues of a real symmetric or complex hermitian matrix + !! version: experimental + !! + !! Returns the eigenvalues \( lambda \), \( A \cdot \bar{v} - \lambda \cdot \bar{v} \), for a real + !! symmetric \( A = A^T \) or complex Hermitian \( A = A^H \) square matrix. + !! ([Specification](../page/specs/stdlib_linalg.html#eigvalsh-eigenvalues-of-a-real-symmetric-or-complex-hermitian-square-matrix)) + !! + !!### Summary + !! Function interface for computing the eigenvalues of a real symmetric or complex hermitian square matrix. + !! + !!### Description + !! + !! This interface provides functions for returning the eigenvalues of a real symmetric or complex Hermitian + !! square matrix. Supported data types include `real` and `complex`. The matrix must be symmetric + !! (if `real`) or Hermitian (if `complex`). Only the lower or upper half of the matrix is accessed, + !! and the user can select which using the optional `upper_a` flag (default: use lower half). + !! An `error stop` is thrown in case of failure; otherwise, error information can be returned + !! as an optional `type(linalg_state_type)` output flag. + !! + !!@note The solution is based on LAPACK's eigenproblem solvers `*SYEV`/`*HEEV`. + !!@note BLAS/LAPACK backends do not currently support extended precision (``xdp``). + !! #:for rk,rt,ri in RC_KINDS_TYPES #:if rk!="xdp" module function stdlib_linalg_eigvalsh_${ri}$(a,upper_a,err) result(lambda) diff --git a/src/stdlib_linalg_eigenvalues.fypp b/src/stdlib_linalg_eigenvalues.fypp index db5282904..ae7ac482c 100644 --- a/src/stdlib_linalg_eigenvalues.fypp +++ b/src/stdlib_linalg_eigenvalues.fypp @@ -99,7 +99,7 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues ${rt}$, intent(in), target :: a(:,:) !> [optional] state return flag. On error if not requested, the code will stop type(linalg_state_type), intent(out) :: err - !> Array of singular values + !> Array of eigenvalues complex(${rk}$), allocatable :: lambda(:) !> Create @@ -125,7 +125,7 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues !! Return an array of eigenvalues of matrix A. !> Input matrix A[m,n] ${rt}$, intent(in), target :: a(:,:) - !> Array of singular values + !> Array of eigenvalues complex(${rk}$), allocatable :: lambda(:) !> Create @@ -154,9 +154,9 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues ${rt}$, intent(inout), target :: a(:,:) !> Array of eigenvalues complex(${rk}$), intent(out) :: lambda(:) - !> The columns of RIGHT contain the right eigenvectors of A + !> [optional] RIGHT eigenvectors of A (as columns) complex(${rk}$), optional, intent(out), target :: right(:,:) - !> The columns of LEFT contain the left eigenvectors of A + !> [optional] LEFT eigenvectors of A (as columns) complex(${rk}$), optional, intent(out), target :: left(:,:) !> [optional] Can A data be overwritten and destroyed? logical(lk), optional, intent(in) :: overwrite_a @@ -323,7 +323,7 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues logical(lk), optional, intent(in) :: upper_a !> [optional] state return flag. On error if not requested, the code will stop type(linalg_state_type), intent(out) :: err - !> Array of singular values + !> Array of eigenvalues real(${rk}$), allocatable :: lambda(:) ${rt}$, pointer :: amat(:,:) @@ -348,9 +348,9 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues !! Return an array of eigenvalues of real symmetric / complex hermitian A !> Input matrix A[m,n] ${rt}$, intent(in), target :: a(:,:) - !> [optional] Should the upper/lower half of A be used? Default: lower + !> [optional] Should the upper/lower half of A be used? Default: use lower logical(lk), optional, intent(in) :: upper_a - !> Array of singular values + !> Array of eigenvalues real(${rk}$), allocatable :: lambda(:) ${rt}$, pointer :: amat(:,:) @@ -382,7 +382,7 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues ${rt}$, optional, intent(out), target :: vectors(:,:) !> [optional] Can A data be overwritten and destroyed? logical(lk), optional, intent(in) :: overwrite_a - !> [optional] Should the upper/lower half of A be used? Default: lower + !> [optional] Should the upper/lower half of A be used? Default: use lower logical(lk), optional, intent(in) :: upper_a !> [optional] state return flag. On error if not requested, the code will stop type(linalg_state_type), optional, intent(out) :: err From f9279b39920f4603fc2ebc67a592f8d86493684a Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 16 May 2024 11:50:23 +0200 Subject: [PATCH 10/32] specs: `eig`, `eigh` --- doc/specs/stdlib_linalg.md | 96 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 96 insertions(+) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index c414c071a..8a77defda 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -687,6 +687,8 @@ Expert (`Pure`) interface: `overwrite_a` (optional): Shall be an input logical flag. if `.true.`, input matrix `a` will be used as temporary storage and overwritten. This avoids internal data allocation. This is an `intent(in)` argument. +`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. + ### Return value For a full-rank matrix, returns an array value that represents the solution to the linear system of equations. @@ -898,3 +900,97 @@ Exceptions trigger an `error stop`. ```fortran {!example/linalg/example_determinant2.f90!} ``` + +## `eig` - Eigenvalues and Eigenvectors of a Square Matrix + +### Status + +Experimental + +### Description + +This subroutine computes the solution to the eigenproblem \( A \cdot \bar{v} - \lambda \cdot \bar{v} \), where \( A \) is a square, full-rank, `real` or `complex` matrix. + +Result array `lambda` returns the eigenvalues of \( A \). The user can request eigenvectors to be returned: +on output `left` may contain the left eigenvectors, `right` the right eigenvectors of \( A \). +Each eigenvector is returned as a column of `left` of `right`. +The solver is based on LAPACK's `*GEEV` backends. + +### Syntax + +`call ` [[stdlib_linalg(module):eig(interface)]] `(a, lambda [, right] [,left] [,overwrite_a] [,err])` + +### Arguments + +`a` : `real` or `complex` square array containing the coefficient matrix. It is normally an `intent(in)` argument. If `overwrite_a=.true.`, it is an `intent(inout)` argument and is destroyed by the call. + +`lambda`: Shall be a `complex` rank-1 array of the same kind as `a`, containing the eigenvalues. It is an `intent(out)` argument. + +`right` (optional): Shall be a `complex` rank-2 array of the same size and kind as `a`, containing the right eigenvectors of `a`. It is an `intent(out)` argument. + +`left` (optional): Shall be a `complex` rank-2 array of the same size and kind as `a`, containing the left eigenvectors of `a`. It is an `intent(out)` argument. + +`overwrite_a` (optional): Shall be an input logical flag. if `.true.`, input matrix `a` will be used as temporary storage and overwritten. This avoids internal data allocation. This is an `intent(in)` argument. + +`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. + +### Return value + +Raises `LINALG_ERROR` if the matrix is singular to working precision. +Raises `LINALG_VALUE_ERROR` if any matrix or arrays have invalid/incompatible sizes. +If `err` is not present, exceptions trigger an `error stop`. + +### Example + +```fortran +{!example/linalg/example_eig.f90!} +``` + +## `eigh` - Eigenvalues and Eigenvectors of a Real symmetric or Complex Hermitian Square Matrix + +### Status + +Experimental + +### Description + +This subroutine computes the solution to the eigendecomposition \( A \cdot \bar{v} - \lambda \cdot \bar{v} \), +where \( A \) is a square, full-rank, `real` symmetric \( A = A^T \) or `complex` Hermitian \( A = A^H \) matrix. + +Result array `lambda` returns the `real` eigenvalues of \( A \). The user can request the orthogonal eigenvectors +to be returned: on output `vectors` may contain the matrix of eigenvectors, returned as a columns. + +Normally, only the lower triangular part of \( A \) is accessed. On input, `logical` flag `upper_a` +allows the user to request what triangular part of the matrix should be used. + +The solver is based on LAPACK's `*SYEV` and `*HEEV` backends. + +### Syntax + +`call ` [[stdlib_linalg(module):eigh(interface)]] `(a, lambda [, vectors] [, upper_a] [, overwrite_a] [,err])` + +### Arguments + +`a` : `real` or `complex` square array containing the coefficient matrix. It is normally an `intent(in)` argument. If `overwrite_a=.true.`, it is an `intent(inout)` argument and is destroyed by the call. + +`lambda`: Shall be a `complex` rank-1 array of the same kind as `a`, containing the eigenvalues. It is an `intent(out)` argument. + +`vectors` (optional): Shall be a rank-2 array of the same type, size and kind as `a`, containing the eigenvectors of `a`. It is an `intent(out)` argument. + +`upper_a` (optional): Shall be an input logical flag. if `.true.`, the upper triangular part of `a` will be used accessed. Otherwise, the lower triangular part will be accessed. It is an `intent(in)` argument. + +`overwrite_a` (optional): Shall be an input logical flag. if `.true.`, input matrix `a` will be used as temporary storage and overwritten. This avoids internal data allocation. This is an `intent(in)` argument. + +`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. + +### Return value + +Raises `LINALG_ERROR` if the matrix is singular to working precision. +Raises `LINALG_VALUE_ERROR` if any matrix or arrays have invalid/incompatible sizes. +If `err` is not present, exceptions trigger an `error stop`. + +### Example + +```fortran +{!example/linalg/example_eigh.f90!} +``` From 56d89a8f2fb8f55e8f23aea6c3db3b1558829b63 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Thu, 16 May 2024 11:58:09 +0200 Subject: [PATCH 11/32] specs: `eigvals`, `eigvalsh` --- doc/specs/stdlib_linalg.md | 83 +++++++++++++++++++++++++++++++++++++- 1 file changed, 81 insertions(+), 2 deletions(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 8a77defda..7d8ef1b99 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -936,7 +936,7 @@ The solver is based on LAPACK's `*GEEV` backends. ### Return value -Raises `LINALG_ERROR` if the matrix is singular to working precision. +Raises `LINALG_ERROR` if the calculation did not converge. Raises `LINALG_VALUE_ERROR` if any matrix or arrays have invalid/incompatible sizes. If `err` is not present, exceptions trigger an `error stop`. @@ -985,7 +985,7 @@ The solver is based on LAPACK's `*SYEV` and `*HEEV` backends. ### Return value -Raises `LINALG_ERROR` if the matrix is singular to working precision. +Raises `LINALG_ERROR` if the calculation did not converge. Raises `LINALG_VALUE_ERROR` if any matrix or arrays have invalid/incompatible sizes. If `err` is not present, exceptions trigger an `error stop`. @@ -994,3 +994,82 @@ If `err` is not present, exceptions trigger an `error stop`. ```fortran {!example/linalg/example_eigh.f90!} ``` + +## `eigvals` - Eigenvalues of a Square Matrix + +### Status + +Experimental + +### Description + +This function returns the eigenvalues to matrix \( A \): a square, full-rank, `real` or `complex` matrix. +The eigenvalues are solutions to the eigenproblem \( A \cdot \bar{v} - \lambda \cdot \bar{v} \). + +Result array `lambda` is `complex`, and returns the eigenvalues of \( A \). +The solver is based on LAPACK's `*GEEV` backends. + +### Syntax + +`lambda = ` [[stdlib_linalg(module):eigvals(interface)]] `(a, [,err])` + +### Arguments + +`a` : `real` or `complex` square array containing the coefficient matrix. It is normally an `intent(in)` argument. + +`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. + +### Return value + +Returns a `complex` array containing the eigenvalues of `a`. + +Raises `LINALG_ERROR` if the calculation did not converge. +Raises `LINALG_VALUE_ERROR` if any matrix or arrays have invalid/incompatible sizes. +If `err` is not present, exceptions trigger an `error stop`. + +### Example + +```fortran +{!example/linalg/example_eigvals.f90!} +``` + +## `eigvalsh` - Eigenvalues of a Real Symmetric or Complex Hermitian Square Matrix + +### Status + +Experimental + +### Description + +This function returns the eigenvalues to matrix \( A \): a where \( A \) is a square, full-rank, +`real` symmetric \( A = A^T \) or `complex` Hermitian \( A = A^H \) matrix. +The eigenvalues are solutions to the eigenproblem \( A \cdot \bar{v} - \lambda \cdot \bar{v} \). + +Result array `lambda` is `real`, and returns the eigenvalues of \( A \). +The solver is based on LAPACK's `*SYEV` and `*HEEV` backends. + +### Syntax + +`lambda = ` [[stdlib_linalg(module):eigvalsh(interface)]] `(a, [, upper_a] [,err])` + +### Arguments + +`a` : `real` or `complex` square array containing the coefficient matrix. It is normally an `intent(in)` argument. + +`upper_a` (optional): Shall be an input logical flag. if `.true.`, the upper triangular part of `a` will be used accessed. Otherwise, the lower triangular part will be accessed (default). It is an `intent(in)` argument. + +`err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. + +### Return value + +Returns a `real` array containing the eigenvalues of `a`. + +Raises `LINALG_ERROR` if the calculation did not converge. +Raises `LINALG_VALUE_ERROR` if any matrix or arrays have invalid/incompatible sizes. +If `err` is not present, exceptions trigger an `error stop`. + +### Example + +```fortran +{!example/linalg/example_eigvalsh.f90!} +``` From f13e75d76f12af4c2084efc0b962d05b1474ebc2 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 5 Jun 2024 08:31:24 +0200 Subject: [PATCH 12/32] clarify svd --- doc/specs/stdlib_linalg.md | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 7d8ef1b99..aa2b39a1e 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -911,9 +911,8 @@ Experimental This subroutine computes the solution to the eigenproblem \( A \cdot \bar{v} - \lambda \cdot \bar{v} \), where \( A \) is a square, full-rank, `real` or `complex` matrix. -Result array `lambda` returns the eigenvalues of \( A \). The user can request eigenvectors to be returned: -on output `left` may contain the left eigenvectors, `right` the right eigenvectors of \( A \). -Each eigenvector is returned as a column of `left` of `right`. +Result array `lambda` returns the eigenvalues of \( A \). The user can request eigenvectors to be returned: if provided, on output `left` will contain the left eigenvectors, `right` the right eigenvectors of \( A \). +Both `left` and `right` are rank-2 arrays, where eigenvectors are stored as columns. The solver is based on LAPACK's `*GEEV` backends. ### Syntax From 34bb7d7f4a72682cbf8c21661d0afa802f798655 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 5 Jun 2024 08:32:22 +0200 Subject: [PATCH 13/32] Update doc/specs/stdlib_linalg.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 7d8ef1b99..fd5dfbc22 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -922,7 +922,7 @@ The solver is based on LAPACK's `*GEEV` backends. ### Arguments -`a` : `real` or `complex` square array containing the coefficient matrix. It is normally an `intent(in)` argument. If `overwrite_a=.true.`, it is an `intent(inout)` argument and is destroyed by the call. +`a` : `real` or `complex` square array containing the coefficient matrix. If `overwrite_a=.false.`, it is an `intent(in)` argument. Otherwise, it is an `intent(inout)` argument and is destroyed by the call. `lambda`: Shall be a `complex` rank-1 array of the same kind as `a`, containing the eigenvalues. It is an `intent(out)` argument. From cdef45f142bc40da7599574c87fe554088298901 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 5 Jun 2024 08:34:05 +0200 Subject: [PATCH 14/32] remove `(sp)` from the examples --- example/linalg/example_eig.f90 | 5 ++--- example/linalg/example_eigh.f90 | 5 ++--- example/linalg/example_eigvals.f90 | 5 ++--- example/linalg/example_eigvalsh.f90 | 5 ++--- 4 files changed, 8 insertions(+), 12 deletions(-) diff --git a/example/linalg/example_eig.f90 b/example/linalg/example_eig.f90 index beb1ef271..ca7e6728b 100644 --- a/example/linalg/example_eig.f90 +++ b/example/linalg/example_eig.f90 @@ -1,12 +1,11 @@ ! Eigendecomposition of a real square matrix program example_eig - use stdlib_linalg_constants, only: sp use stdlib_linalg, only: eig implicit none integer :: i - real(sp), allocatable :: A(:,:) - complex(sp), allocatable :: lambda(:),vectors(:,:) + real, allocatable :: A(:,:) + complex, allocatable :: lambda(:),vectors(:,:) ! Decomposition of this square matrix ! NB Fortran is column-major -> transpose input diff --git a/example/linalg/example_eigh.f90 b/example/linalg/example_eigh.f90 index 1cf30138d..d5cf9793f 100644 --- a/example/linalg/example_eigh.f90 +++ b/example/linalg/example_eigh.f90 @@ -1,12 +1,11 @@ ! Eigendecomposition of a real symmetric matrix program example_eigh - use stdlib_linalg_constants, only: sp use stdlib_linalg, only: eigh implicit none integer :: i - real(sp), allocatable :: A(:,:),lambda(:),v(:,:) - complex(sp), allocatable :: cA(:,:),cv(:,:) + real, allocatable :: A(:,:),lambda(:),v(:,:) + complex, allocatable :: cA(:,:),cv(:,:) ! Decomposition of this symmetric matrix ! NB Fortran is column-major -> transpose input diff --git a/example/linalg/example_eigvals.f90 b/example/linalg/example_eigvals.f90 index 973d50666..eff28b9d9 100644 --- a/example/linalg/example_eigvals.f90 +++ b/example/linalg/example_eigvals.f90 @@ -1,12 +1,11 @@ ! Eigenvalues of a general real / complex matrix program example_eigvals - use stdlib_linalg_constants, only: sp use stdlib_linalg, only: eigvals implicit none integer :: i - real(sp), allocatable :: A(:,:),lambda(:) - complex(sp), allocatable :: cA(:,:),clambda(:) + real, allocatable :: A(:,:),lambda(:) + complex, allocatable :: cA(:,:),clambda(:) ! NB Fortran is column-major -> transpose input A = transpose(reshape( [ [2, 8, 4], & diff --git a/example/linalg/example_eigvalsh.f90 b/example/linalg/example_eigvalsh.f90 index 6853d9051..de672b0df 100644 --- a/example/linalg/example_eigvalsh.f90 +++ b/example/linalg/example_eigvalsh.f90 @@ -1,12 +1,11 @@ ! Eigenvalues of a real symmetric / complex hermitian matrix program example_eigvalsh - use stdlib_linalg_constants, only: sp use stdlib_linalg, only: eigvalsh implicit none integer :: i - real(sp), allocatable :: A(:,:),lambda(:) - complex(sp), allocatable :: cA(:,:) + real, allocatable :: A(:,:),lambda(:) + complex, allocatable :: cA(:,:) ! Decomposition of this symmetric matrix ! NB Fortran is column-major -> transpose input From fe1c89e0956e5c7dd57def880087d526e84bae22 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 5 Jun 2024 08:34:32 +0200 Subject: [PATCH 15/32] Update doc/specs/stdlib_linalg.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index fd5dfbc22..84ce319b5 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -971,7 +971,7 @@ The solver is based on LAPACK's `*SYEV` and `*HEEV` backends. ### Arguments -`a` : `real` or `complex` square array containing the coefficient matrix. It is normally an `intent(in)` argument. If `overwrite_a=.true.`, it is an `intent(inout)` argument and is destroyed by the call. +`a` : `real` or `complex` square array containing the coefficient matrix. It is an `intent(in)` argument. If `overwrite_a=.true.`, it is an `intent(inout)` argument and is destroyed by the call. `lambda`: Shall be a `complex` rank-1 array of the same kind as `a`, containing the eigenvalues. It is an `intent(out)` argument. From 65099e3397cb1aff13fd7264dbe633b3f8db2f41 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 5 Jun 2024 08:34:49 +0200 Subject: [PATCH 16/32] Update doc/specs/stdlib_linalg.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 84ce319b5..ba831685f 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -977,7 +977,7 @@ The solver is based on LAPACK's `*SYEV` and `*HEEV` backends. `vectors` (optional): Shall be a rank-2 array of the same type, size and kind as `a`, containing the eigenvectors of `a`. It is an `intent(out)` argument. -`upper_a` (optional): Shall be an input logical flag. if `.true.`, the upper triangular part of `a` will be used accessed. Otherwise, the lower triangular part will be accessed. It is an `intent(in)` argument. +`upper_a` (optional): Shall be an input `logical` flag. if `.true.`, the upper triangular part of `a` will be accessed. Otherwise, the lower triangular part will be accessed. It is an `intent(in)` argument. `overwrite_a` (optional): Shall be an input logical flag. if `.true.`, input matrix `a` will be used as temporary storage and overwritten. This avoids internal data allocation. This is an `intent(in)` argument. From 0e13fe59f54cfabeaedaf032800308873c395386 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 5 Jun 2024 08:36:02 +0200 Subject: [PATCH 17/32] kind -> precision --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index aa2b39a1e..bb15bd5c8 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -972,7 +972,7 @@ The solver is based on LAPACK's `*SYEV` and `*HEEV` backends. `a` : `real` or `complex` square array containing the coefficient matrix. It is normally an `intent(in)` argument. If `overwrite_a=.true.`, it is an `intent(inout)` argument and is destroyed by the call. -`lambda`: Shall be a `complex` rank-1 array of the same kind as `a`, containing the eigenvalues. It is an `intent(out)` argument. +`lambda`: Shall be a `complex` rank-1 array of the same precision as `a`, containing the eigenvalues. It is an `intent(out)` argument. `vectors` (optional): Shall be a rank-2 array of the same type, size and kind as `a`, containing the eigenvectors of `a`. It is an `intent(out)` argument. From d0f385ab7a07adef9c68a5d76d57ace5b92e6edf Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 5 Jun 2024 08:36:16 +0200 Subject: [PATCH 18/32] Update doc/specs/stdlib_linalg.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index ba831685f..7dc87f336 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -979,7 +979,7 @@ The solver is based on LAPACK's `*SYEV` and `*HEEV` backends. `upper_a` (optional): Shall be an input `logical` flag. if `.true.`, the upper triangular part of `a` will be accessed. Otherwise, the lower triangular part will be accessed. It is an `intent(in)` argument. -`overwrite_a` (optional): Shall be an input logical flag. if `.true.`, input matrix `a` will be used as temporary storage and overwritten. This avoids internal data allocation. This is an `intent(in)` argument. +`overwrite_a` (optional): Shall be an input `logical` flag. If `.true.`, input matrix `a` will be used as temporary storage and overwritten. This avoids internal data allocation. This is an `intent(in)` argument. `err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. From 36c8d5afaf4d4921695da36adccc78adac2523ff Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 5 Jun 2024 09:16:09 +0200 Subject: [PATCH 19/32] fix `sp` in example --- example/linalg/example_eigvals.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/example/linalg/example_eigvals.f90 b/example/linalg/example_eigvals.f90 index eff28b9d9..7a8727e75 100644 --- a/example/linalg/example_eigvals.f90 +++ b/example/linalg/example_eigvals.f90 @@ -17,7 +17,7 @@ program example_eigvals print *, 'Real matrix eigenvalues: ',lambda ! Complex general matrix - cA = cmplx(A, -2*A, kind=sp) + cA = cmplx(A, -2*A) clambda = eigvals(cA) print *, 'Complex matrix eigenvalues: ',clambda From e6aa0f50e95eefe96cd5e6ade2c50e42fa158cdc Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 5 Jun 2024 09:19:15 +0200 Subject: [PATCH 20/32] Option to output `real` eigenvalues only --- doc/specs/stdlib_linalg.md | 3 ++- src/stdlib_linalg.fypp | 21 +++++++++++++++ src/stdlib_linalg_eigenvalues.fypp | 42 +++++++++++++++++++++++++++++- 3 files changed, 64 insertions(+), 2 deletions(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index cc3ee2778..aa50aab51 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -923,7 +923,7 @@ The solver is based on LAPACK's `*GEEV` backends. `a` : `real` or `complex` square array containing the coefficient matrix. If `overwrite_a=.false.`, it is an `intent(in)` argument. Otherwise, it is an `intent(inout)` argument and is destroyed by the call. -`lambda`: Shall be a `complex` rank-1 array of the same kind as `a`, containing the eigenvalues. It is an `intent(out)` argument. +`lambda`: Shall be a `complex` or `real` rank-1 array of the same kind as `a`, containing the eigenvalues, or their `real` component only. It is an `intent(out)` argument. `right` (optional): Shall be a `complex` rank-2 array of the same size and kind as `a`, containing the right eigenvectors of `a`. It is an `intent(out)` argument. @@ -937,6 +937,7 @@ The solver is based on LAPACK's `*GEEV` backends. Raises `LINALG_ERROR` if the calculation did not converge. Raises `LINALG_VALUE_ERROR` if any matrix or arrays have invalid/incompatible sizes. +Raises `LINALG_VALUE_ERROR` if the `real` component is only requested, but the eigenvalues have non-trivial imaginary parts. If `err` is not present, exceptions trigger an `error stop`. ### Example diff --git a/src/stdlib_linalg.fypp b/src/stdlib_linalg.fypp index 90b920395..2e8d9ef8a 100644 --- a/src/stdlib_linalg.fypp +++ b/src/stdlib_linalg.fypp @@ -597,6 +597,27 @@ module stdlib_linalg end subroutine stdlib_linalg_eig_${ri}$ #:endif #:endfor + #:for rk,rt,ri in REAL_KINDS_TYPES + #:if rk!="xdp" + module subroutine stdlib_linalg_real_eig_${ri}$(a,lambda,right,left,overwrite_a,err) + !! Eigendecomposition of matrix A returning an array `lambda` of real eigenvalues, + !! and optionally right or left eigenvectors. Returns an error if the eigenvalues had + !! non-trivial imaginary parts. + !> Input matrix A[m,n] + ${rt}$, intent(inout), target :: a(:,:) + !> Array of real eigenvalues + real(${rk}$), intent(out) :: lambda(:) + !> The columns of RIGHT contain the right eigenvectors of A + complex(${rk}$), optional, intent(out), target :: right(:,:) + !> The columns of LEFT contain the left eigenvectors of A + complex(${rk}$), optional, intent(out), target :: left(:,:) + !> [optional] Can A data be overwritten and destroyed? + logical(lk), optional, intent(in) :: overwrite_a + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), optional, intent(out) :: err + end subroutine stdlib_linalg_real_eig_${ri}$ + #:endif + #:endfor end interface eig ! Eigenvalues of a square matrix diff --git a/src/stdlib_linalg_eigenvalues.fypp b/src/stdlib_linalg_eigenvalues.fypp index ae7ac482c..37f40ff8e 100644 --- a/src/stdlib_linalg_eigenvalues.fypp +++ b/src/stdlib_linalg_eigenvalues.fypp @@ -315,6 +315,47 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues end subroutine stdlib_linalg_eig_${ri}$ + module subroutine stdlib_linalg_real_eig_${ri}$(a,lambda,right,left,overwrite_a,err) + !! Eigendecomposition of matrix A returning an array `lambda` of real eigenvalues, + !! and optionally right or left eigenvectors. Returns an error if the eigenvalues had + !! non-trivial imaginary parts. + !> Input matrix A[m,n] + ${rt}$, intent(inout), target :: a(:,:) + !> Array of real eigenvalues + real(${rk}$), intent(out) :: lambda(:) + !> The columns of RIGHT contain the right eigenvectors of A + complex(${rk}$), optional, intent(out), target :: right(:,:) + !> The columns of LEFT contain the left eigenvectors of A + complex(${rk}$), optional, intent(out), target :: left(:,:) + !> [optional] Can A data be overwritten and destroyed? + logical(lk), optional, intent(in) :: overwrite_a + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), optional, intent(out) :: err + + type(linalg_state_type) :: err0 + integer(ilp) :: n + complex(${rk}$), allocatable :: clambda(:) + real(${rk}$), parameter :: rtol = epsilon(0.0_${rk}$) + real(${rk}$), parameter :: atol = tiny(0.0_${rk}$) + + n = size(lambda,dim=1,kind=ilp) + allocate(clambda(n)) + + call stdlib_linalg_eig_${ri}$(a,clambda,right,left,overwrite_a,err0) + + ! Check that no eigenvalues have meaningful imaginary part + if (err0%ok() .and. any(aimag(clambda)>atol+rtol*abs(abs(clambda)))) then + err0 = linalg_state_type(this,LINALG_VALUE_ERROR, & + 'complex eigenvalues detected: max(imag(lambda))=',maxval(aimag(clambda))) + endif + + ! Return real components only + lambda(:n) = real(clambda,kind=${rk}$) + + call linalg_error_handling(err0,err) + + end subroutine stdlib_linalg_real_eig_${ri}$ + module function stdlib_linalg_eigvalsh_${ri}$(a,upper_a,err) result(lambda) !! Return an array of eigenvalues of real symmetric / complex hermitian A !> Input matrix A[m,n] @@ -535,5 +576,4 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues #:endif #:endfor - end submodule stdlib_linalg_eigenvalues From 17ebed27c66df0128a6721f6f90f2f72fa8c87d2 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 21 Jun 2024 14:22:11 +0200 Subject: [PATCH 21/32] reorganize `if(present(x))` --- src/stdlib_linalg_eigenvalues.fypp | 14 ++++---------- 1 file changed, 4 insertions(+), 10 deletions(-) diff --git a/src/stdlib_linalg_eigenvalues.fypp b/src/stdlib_linalg_eigenvalues.fypp index 37f40ff8e..d5111d435 100644 --- a/src/stdlib_linalg_eigenvalues.fypp +++ b/src/stdlib_linalg_eigenvalues.fypp @@ -21,11 +21,8 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues !> Request for symmetry side (default: lower) elemental character function symmetric_triangle_task(upper) logical(lk), optional, intent(in) :: upper - if (present(upper)) then - symmetric_triangle_task = merge('U','L',upper) - else - symmetric_triangle_task = 'L' - endif + symmetric_triangle_task = 'L' + if (present(upper)) symmetric_triangle_task = merge('U','L',upper) end function symmetric_triangle_task !> Process GEEV output flags @@ -194,11 +191,8 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues endif ! Can A be overwritten? By default, do not overwrite - if (present(overwrite_a)) then - copy_a = .not.overwrite_a - else - copy_a = .true._lk - endif + copy_a = .true._lk + if (present(overwrite_a)) copy_a = .not.overwrite_a ! Initialize a matrix temporary if (copy_a) then From bac318798d2096f59f659ef0c0b7ecc1e0ecbcc0 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 21 Jun 2024 14:23:34 +0200 Subject: [PATCH 22/32] no negated checks --- src/stdlib_linalg_eigenvalues.fypp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/stdlib_linalg_eigenvalues.fypp b/src/stdlib_linalg_eigenvalues.fypp index d5111d435..07f2c7fef 100644 --- a/src/stdlib_linalg_eigenvalues.fypp +++ b/src/stdlib_linalg_eigenvalues.fypp @@ -177,12 +177,12 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues neig = size(lambda,kind=ilp) lda = m - if (.not.(k>0 .and. m==n)) then + if (k<=0 .or. m/=n)) then err0 = linalg_state_type(this,LINALG_VALUE_ERROR,& - 'invalid or matrix size a=',[m,n],', must be square.') + 'invalid or matrix size a=',[m,n],', must be nonempty square.') call linalg_error_handling(err0,err) return - elseif (.not.neig>=k) then + elseif (neig=k) then + elseif (neig= n=',n) call linalg_error_handling(err0,err) From 4f503f74e3c751f8b79855484b4de08c1b6e2e4a Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 21 Jun 2024 14:24:44 +0200 Subject: [PATCH 23/32] copy_a: simplify --- src/stdlib_linalg_eigenvalues.fypp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/stdlib_linalg_eigenvalues.fypp b/src/stdlib_linalg_eigenvalues.fypp index 07f2c7fef..d33af9902 100644 --- a/src/stdlib_linalg_eigenvalues.fypp +++ b/src/stdlib_linalg_eigenvalues.fypp @@ -451,13 +451,12 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues endif ! Check if input A can be overwritten + copy_a = .true._lk if (present(vectors)) then ! No need to copy A anyways copy_a = .false. elseif (present(overwrite_a)) then copy_a = .not.overwrite_a - else - copy_a = .true._lk endif ! Should we use the upper or lower half of the matrix? From f232a763a131517b244a79bad2101741c583defa Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 21 Jun 2024 15:37:31 +0200 Subject: [PATCH 24/32] typo --- src/stdlib_linalg_eigenvalues.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stdlib_linalg_eigenvalues.fypp b/src/stdlib_linalg_eigenvalues.fypp index d33af9902..ac2b0b170 100644 --- a/src/stdlib_linalg_eigenvalues.fypp +++ b/src/stdlib_linalg_eigenvalues.fypp @@ -177,7 +177,7 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues neig = size(lambda,kind=ilp) lda = m - if (k<=0 .or. m/=n)) then + if (k<=0 .or. m/=n) then err0 = linalg_state_type(this,LINALG_VALUE_ERROR,& 'invalid or matrix size a=',[m,n],', must be nonempty square.') call linalg_error_handling(err0,err) From d19425ab99d8f9e912bf6707f08ff893b54297d7 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Fri, 21 Jun 2024 15:55:38 +0200 Subject: [PATCH 25/32] fix intel `module subroutine` error --- src/stdlib_linalg_eigenvalues.fypp | 82 +++++++++++++++--------------- 1 file changed, 41 insertions(+), 41 deletions(-) diff --git a/src/stdlib_linalg_eigenvalues.fypp b/src/stdlib_linalg_eigenvalues.fypp index ac2b0b170..2c8ee8783 100644 --- a/src/stdlib_linalg_eigenvalues.fypp +++ b/src/stdlib_linalg_eigenvalues.fypp @@ -309,47 +309,6 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues end subroutine stdlib_linalg_eig_${ri}$ - module subroutine stdlib_linalg_real_eig_${ri}$(a,lambda,right,left,overwrite_a,err) - !! Eigendecomposition of matrix A returning an array `lambda` of real eigenvalues, - !! and optionally right or left eigenvectors. Returns an error if the eigenvalues had - !! non-trivial imaginary parts. - !> Input matrix A[m,n] - ${rt}$, intent(inout), target :: a(:,:) - !> Array of real eigenvalues - real(${rk}$), intent(out) :: lambda(:) - !> The columns of RIGHT contain the right eigenvectors of A - complex(${rk}$), optional, intent(out), target :: right(:,:) - !> The columns of LEFT contain the left eigenvectors of A - complex(${rk}$), optional, intent(out), target :: left(:,:) - !> [optional] Can A data be overwritten and destroyed? - logical(lk), optional, intent(in) :: overwrite_a - !> [optional] state return flag. On error if not requested, the code will stop - type(linalg_state_type), optional, intent(out) :: err - - type(linalg_state_type) :: err0 - integer(ilp) :: n - complex(${rk}$), allocatable :: clambda(:) - real(${rk}$), parameter :: rtol = epsilon(0.0_${rk}$) - real(${rk}$), parameter :: atol = tiny(0.0_${rk}$) - - n = size(lambda,dim=1,kind=ilp) - allocate(clambda(n)) - - call stdlib_linalg_eig_${ri}$(a,clambda,right,left,overwrite_a,err0) - - ! Check that no eigenvalues have meaningful imaginary part - if (err0%ok() .and. any(aimag(clambda)>atol+rtol*abs(abs(clambda)))) then - err0 = linalg_state_type(this,LINALG_VALUE_ERROR, & - 'complex eigenvalues detected: max(imag(lambda))=',maxval(aimag(clambda))) - endif - - ! Return real components only - lambda(:n) = real(clambda,kind=${rk}$) - - call linalg_error_handling(err0,err) - - end subroutine stdlib_linalg_real_eig_${ri}$ - module function stdlib_linalg_eigvalsh_${ri}$(a,upper_a,err) result(lambda) !! Return an array of eigenvalues of real symmetric / complex hermitian A !> Input matrix A[m,n] @@ -566,6 +525,47 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues end subroutine assign_real_eigenvectors_${rk}$ + module subroutine stdlib_linalg_real_eig_${ri}$(a,lambda,right,left,overwrite_a,err) + !! Eigendecomposition of matrix A returning an array `lambda` of real eigenvalues, + !! and optionally right or left eigenvectors. Returns an error if the eigenvalues had + !! non-trivial imaginary parts. + !> Input matrix A[m,n] + ${rt}$, intent(inout), target :: a(:,:) + !> Array of real eigenvalues + real(${rk}$), intent(out) :: lambda(:) + !> The columns of RIGHT contain the right eigenvectors of A + complex(${rk}$), optional, intent(out), target :: right(:,:) + !> The columns of LEFT contain the left eigenvectors of A + complex(${rk}$), optional, intent(out), target :: left(:,:) + !> [optional] Can A data be overwritten and destroyed? + logical(lk), optional, intent(in) :: overwrite_a + !> [optional] state return flag. On error if not requested, the code will stop + type(linalg_state_type), optional, intent(out) :: err + + type(linalg_state_type) :: err0 + integer(ilp) :: n + complex(${rk}$), allocatable :: clambda(:) + real(${rk}$), parameter :: rtol = epsilon(0.0_${rk}$) + real(${rk}$), parameter :: atol = tiny(0.0_${rk}$) + + n = size(lambda,dim=1,kind=ilp) + allocate(clambda(n)) + + call stdlib_linalg_eig_${ri}$(a,clambda,right,left,overwrite_a,err0) + + ! Check that no eigenvalues have meaningful imaginary part + if (err0%ok() .and. any(aimag(clambda)>atol+rtol*abs(abs(clambda)))) then + err0 = linalg_state_type(this,LINALG_VALUE_ERROR, & + 'complex eigenvalues detected: max(imag(lambda))=',maxval(aimag(clambda))) + endif + + ! Return real components only + lambda(:n) = real(clambda,kind=${rk}$) + + call linalg_error_handling(err0,err) + + end subroutine stdlib_linalg_real_eig_${ri}$ + #:endif #:endfor From cc95fc425449b3e59d0984310514b90996737869 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 26 Jun 2024 08:36:43 +0200 Subject: [PATCH 26/32] Update doc/specs/stdlib_linalg.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index c6d2521a1..af6c68ba0 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -958,7 +958,7 @@ This subroutine computes the solution to the eigendecomposition \( A \cdot \bar{ where \( A \) is a square, full-rank, `real` symmetric \( A = A^T \) or `complex` Hermitian \( A = A^H \) matrix. Result array `lambda` returns the `real` eigenvalues of \( A \). The user can request the orthogonal eigenvectors -to be returned: on output `vectors` may contain the matrix of eigenvectors, returned as a columns. +to be returned: on output `vectors` may contain the matrix of eigenvectors, returned as a column. Normally, only the lower triangular part of \( A \) is accessed. On input, `logical` flag `upper_a` allows the user to request what triangular part of the matrix should be used. From 18b95afc1e9c6d8baafe6eaf1a1791ac3642036d Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 26 Jun 2024 08:36:52 +0200 Subject: [PATCH 27/32] Update doc/specs/stdlib_linalg.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index af6c68ba0..414b9e1a3 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -977,7 +977,7 @@ The solver is based on LAPACK's `*SYEV` and `*HEEV` backends. `vectors` (optional): Shall be a rank-2 array of the same type, size and kind as `a`, containing the eigenvectors of `a`. It is an `intent(out)` argument. -`upper_a` (optional): Shall be an input `logical` flag. if `.true.`, the upper triangular part of `a` will be accessed. Otherwise, the lower triangular part will be accessed. It is an `intent(in)` argument. +`upper_a` (optional): Shall be an input `logical` flag. If `.true.`, the upper triangular part of `a` will be accessed. Otherwise, the lower triangular part will be accessed. It is an `intent(in)` argument. `overwrite_a` (optional): Shall be an input `logical` flag. If `.true.`, input matrix `a` will be used as temporary storage and overwritten. This avoids internal data allocation. This is an `intent(in)` argument. From 1fb2ba90d24255d5ff773551eb4134a9bd8883c0 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 26 Jun 2024 08:37:01 +0200 Subject: [PATCH 28/32] Update doc/specs/stdlib_linalg.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index 414b9e1a3..c87c21ffb 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -1015,7 +1015,7 @@ The solver is based on LAPACK's `*GEEV` backends. ### Arguments -`a` : `real` or `complex` square array containing the coefficient matrix. It is normally an `intent(in)` argument. +`a` : `real` or `complex` square array containing the coefficient matrix. It is an `intent(in)` argument. `err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. From f5303bf84996da935d4c136406add9203fc8bb7f Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 26 Jun 2024 08:37:12 +0200 Subject: [PATCH 29/32] Update doc/specs/stdlib_linalg.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index c87c21ffb..c2d1bbfc4 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -1055,7 +1055,7 @@ The solver is based on LAPACK's `*SYEV` and `*HEEV` backends. ### Arguments -`a` : `real` or `complex` square array containing the coefficient matrix. It is normally an `intent(in)` argument. +`a` : `real` or `complex` square array containing the coefficient matrix. It is an `intent(in)` argument. `upper_a` (optional): Shall be an input logical flag. if `.true.`, the upper triangular part of `a` will be used accessed. Otherwise, the lower triangular part will be accessed (default). It is an `intent(in)` argument. From 66f1f17763b722324d4ac42ca5ba8ba0993a88cb Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 26 Jun 2024 08:37:23 +0200 Subject: [PATCH 30/32] Update doc/specs/stdlib_linalg.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_linalg.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_linalg.md b/doc/specs/stdlib_linalg.md index c2d1bbfc4..4e8fe3b19 100644 --- a/doc/specs/stdlib_linalg.md +++ b/doc/specs/stdlib_linalg.md @@ -1057,7 +1057,7 @@ The solver is based on LAPACK's `*SYEV` and `*HEEV` backends. `a` : `real` or `complex` square array containing the coefficient matrix. It is an `intent(in)` argument. -`upper_a` (optional): Shall be an input logical flag. if `.true.`, the upper triangular part of `a` will be used accessed. Otherwise, the lower triangular part will be accessed (default). It is an `intent(in)` argument. +`upper_a` (optional): Shall be an input logical flag. If `.true.`, the upper triangular part of `a` will be used accessed. Otherwise, the lower triangular part will be accessed (default). It is an `intent(in)` argument. `err` (optional): Shall be a `type(linalg_state_type)` value. This is an `intent(out)` argument. From 5501b83f92706d99e6235821cf3aeb76de6792cb Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 26 Jun 2024 09:18:17 +0200 Subject: [PATCH 31/32] example_eigh: v->vectors --- example/linalg/example_eigh.f90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/example/linalg/example_eigh.f90 b/example/linalg/example_eigh.f90 index d5cf9793f..417e73cff 100644 --- a/example/linalg/example_eigh.f90 +++ b/example/linalg/example_eigh.f90 @@ -4,8 +4,8 @@ program example_eigh implicit none integer :: i - real, allocatable :: A(:,:),lambda(:),v(:,:) - complex, allocatable :: cA(:,:),cv(:,:) + real, allocatable :: A(:,:),lambda(:),vectors(:,:) + complex, allocatable :: cA(:,:),cvectors(:,:) ! Decomposition of this symmetric matrix ! NB Fortran is column-major -> transpose input @@ -15,24 +15,24 @@ program example_eigh ! Note: real symmetric matrices have real (orthogonal) eigenvalues and eigenvectors allocate(lambda(3),v(3,3)) - call eigh(A, lambda, vectors=v) + call eigh(A, lambda, vectors=vectors) print *, 'Real matrix' do i=1,3 print *, 'eigenvalue ',i,': ',lambda(i) - print *, 'eigenvector ',i,': ',v(:,i) + print *, 'eigenvector ',i,': ',vectors(:,i) end do ! Complex hermitian matrices have real (orthogonal) eigenvalues and complex eigenvectors cA = A allocate(cv(3,3)) - call eigh(cA, lambda, vectors=cv) + call eigh(cA, lambda, vectors=cvectors) print *, 'Complex matrix' do i=1,3 print *, 'eigenvalue ',i,': ',lambda(i) - print *, 'eigenvector ',i,': ',cv(:,i) + print *, 'eigenvector ',i,': ',cvectors(:,i) end do end program example_eigh From d433869b29af59183d639ab5c839063d7386a327 Mon Sep 17 00:00:00 2001 From: Federico Perini Date: Wed, 26 Jun 2024 09:31:30 +0200 Subject: [PATCH 32/32] fix vectors --- example/linalg/example_eigh.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/example/linalg/example_eigh.f90 b/example/linalg/example_eigh.f90 index 417e73cff..7f19946ae 100644 --- a/example/linalg/example_eigh.f90 +++ b/example/linalg/example_eigh.f90 @@ -14,7 +14,7 @@ program example_eigh [4, 5, 4] ], [3,3] )) ! Note: real symmetric matrices have real (orthogonal) eigenvalues and eigenvectors - allocate(lambda(3),v(3,3)) + allocate(lambda(3),vectors(3,3)) call eigh(A, lambda, vectors=vectors) print *, 'Real matrix' @@ -26,7 +26,7 @@ program example_eigh ! Complex hermitian matrices have real (orthogonal) eigenvalues and complex eigenvectors cA = A - allocate(cv(3,3)) + allocate(cvectors(3,3)) call eigh(cA, lambda, vectors=cvectors) print *, 'Complex matrix'