Skip to content

Commit bcdfa31

Browse files
perazzperazz
committed
scale generalized eigenvalues
1 parent 4888e21 commit bcdfa31

File tree

2 files changed

+135
-21
lines changed

2 files changed

+135
-21
lines changed

‎src/stdlib_linalg.fypp

Lines changed: 13 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -836,7 +836,8 @@ module stdlib_linalg
836836
!!
837837
#:for rk,rt,ri in RC_KINDS_TYPES
838838
#:for ep,ei in EIG_PROBLEM_LIST
839-
module subroutine stdlib_linalg_eig_${ep}$_${ri}$(a#{if ei=='ggev'}#,b#{endif}#,lambda,right,left,overwrite_a,err)
839+
module subroutine stdlib_linalg_eig_${ep}$_${ri}$(a#{if ei=='ggev'}#,b#{endif}#,lambda,right,left, &
840+
overwrite_a#{if ei=='ggev'}#,overwrite_b#{endif}#,err)
840841
!! Eigendecomposition of matrix A returning an array `lambda` of eigenvalues,
841842
!! and optionally right or left eigenvectors.
842843
!> Input matrix A[m,n]
@@ -853,16 +854,21 @@ module stdlib_linalg
853854
complex(${rk}$), optional, intent(out), target :: left(:,:)
854855
!> [optional] Can A data be overwritten and destroyed?
855856
logical(lk), optional, intent(in) :: overwrite_a
857+
#:if ei=='ggev'
858+
!> [optional] Can B data be overwritten and destroyed? (default: no)
859+
logical(lk), optional, intent(in) :: overwrite_b
860+
#:endif
856861
!> [optional] state return flag. On error if not requested, the code will stop
857862
type(linalg_state_type), optional, intent(out) :: err
858863
end subroutine stdlib_linalg_eig_${ep}$_${ri}$
859864

860-
module subroutine stdlib_linalg_real_eig_${ep}$_${ri}$(a#{if ei=='ggev'}#,b#{endif}#,lambda,right,left,overwrite_a,err)
865+
module subroutine stdlib_linalg_real_eig_${ep}$_${ri}$(a#{if ei=='ggev'}#,b#{endif}#,lambda,right,left, &
866+
overwrite_a#{if ei=='ggev'}#,overwrite_b#{endif}#,err)
861867
!! Eigendecomposition of matrix A returning an array `lambda` of real eigenvalues,
862868
!! and optionally right or left eigenvectors. Returns an error if the eigenvalues had
863869
!! non-trivial imaginary parts.
864870
!> Input matrix A[m,n]
865-
${rt}$, intent(in), target :: a(:,:)
871+
${rt}$, intent(inout), target :: a(:,:)
866872
#:if ei=='ggev'
867873
!> Generalized problem matrix B[n,n]
868874
${rt}$, intent(inout), target :: b(:,:)
@@ -875,6 +881,10 @@ module stdlib_linalg
875881
complex(${rk}$), optional, intent(out), target :: left(:,:)
876882
!> [optional] Can A data be overwritten and destroyed?
877883
logical(lk), optional, intent(in) :: overwrite_a
884+
#:if ei=='ggev'
885+
!> [optional] Can B data be overwritten and destroyed? (default: no)
886+
logical(lk), optional, intent(in) :: overwrite_b
887+
#:endif
878888
!> [optional] state return flag. On error if not requested, the code will stop
879889
type(linalg_state_type), optional, intent(out) :: err
880890
end subroutine stdlib_linalg_real_eig_${ep}$_${ri}$

‎src/stdlib_linalg_eigenvalues.fypp

Lines changed: 122 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -9,9 +9,17 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues
99
use stdlib_linalg_lapack, only: geev, ggev, heev, syev
1010
use stdlib_linalg_state, only: linalg_state_type, linalg_error_handling, LINALG_ERROR, &
1111
LINALG_INTERNAL_ERROR, LINALG_VALUE_ERROR, LINALG_SUCCESS
12+
use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_positive_inf, ieee_quiet_nan
1213
implicit none(type,external)
1314

1415
character(*), parameter :: this = 'eigenvalues'
16+
17+
!> Utility function: Scale generalized eigenvalue
18+
interface scale_general_eig
19+
#:for rk,rt,ri in RC_KINDS_TYPES
20+
module procedure scale_general_eig_${ri}$
21+
#:endfor
22+
end interface scale_general_eig
1523

1624
contains
1725

@@ -157,7 +165,7 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues
157165
allocate(lambda(k))
158166

159167
!> Compute eigenvalues only
160-
call stdlib_linalg_eig_${ep}$_${ri}$(amat#{if ei=='ggev'}#,bmat#{endif}#,lambda,overwrite_a=.false.,err=err)
168+
call stdlib_linalg_eig_${ep}$_${ri}$(amat#{if ei=='ggev'}#,bmat#{endif}#,lambda,err=err)
161169

162170
end function stdlib_linalg_eigvals_${ep}$_${ri}$
163171

@@ -192,7 +200,8 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues
192200

193201
end function stdlib_linalg_eigvals_noerr_${ep}$_${ri}$
194202

195-
module subroutine stdlib_linalg_eig_${ep}$_${ri}$(a#{if ei=='ggev'}#,b#{endif}#,lambda,right,left,overwrite_a,err)
203+
module subroutine stdlib_linalg_eig_${ep}$_${ri}$(a#{if ei=='ggev'}#,b#{endif}#,lambda,right,left,&
204+
overwrite_a#{if ei=='ggev'}#,overwrite_b#{endif}#,err)
196205
!! Eigendecomposition of matrix A returning an array `lambda` of eigenvalues,
197206
!! and optionally right or left eigenvectors.
198207
!> Input matrix A[m,n]
@@ -207,15 +216,19 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues
207216
complex(${rk}$), optional, intent(out), target :: right(:,:)
208217
!> [optional] LEFT eigenvectors of A (as columns)
209218
complex(${rk}$), optional, intent(out), target :: left(:,:)
210-
!> [optional] Can A data be overwritten and destroyed?
219+
!> [optional] Can A data be overwritten and destroyed? (default: no)
211220
logical(lk), optional, intent(in) :: overwrite_a
221+
#:if ei=='ggev'
222+
!> [optional] Can B data be overwritten and destroyed? (default: no)
223+
logical(lk), optional, intent(in) :: overwrite_b
224+
#:endif
212225
!> [optional] state return flag. On error if not requested, the code will stop
213226
type(linalg_state_type), optional, intent(out) :: err
214227

215228
!> Local variables
216229
type(linalg_state_type) :: err0
217-
integer(ilp) :: m,n,lda,ldu,ldv,info,k,lwork,neig
218-
logical(lk) :: copy_a
230+
integer(ilp) :: m,n,lda,ldu,ldv,info,k,lwork,neig#{if ei=='ggev'}#,ldb,nb#{endif}#
231+
logical(lk) :: copy_a#{if ei=='ggev'}#,copy_b#{endif}#
219232
character :: task_u,task_v
220233
${rt}$, target :: work_dummy(1),u_dummy(1,1),v_dummy(1,1)
221234
${rt}$, allocatable :: work(:)
@@ -225,6 +238,10 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues
225238
#:else
226239
${rt}$, pointer :: lreal(:),limag(:)
227240
#:endif
241+
#:if ei=='ggev'
242+
${rt}$, allocatable :: beta(:)
243+
#:endif
244+
228245
!> Matrix size
229246
m = size(a,1,kind=ilp)
230247
n = size(a,2,kind=ilp)
@@ -244,18 +261,43 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues
244261
call linalg_error_handling(err0,err)
245262
return
246263
endif
264+
265+
#:if ep=='generalized'
266+
ldb = size(b,1,kind=ilp)
267+
nb = size(b,2,kind=ilp)
268+
if (ldb/=n .or. nb/=n) then
269+
err0 = linalg_state_type(this,LINALG_VALUE_ERROR,&
270+
'invalid or matrix size b=',[ldb,nb],', must be same as a=',[m,n])
271+
call linalg_error_handling(err0,err)
272+
return
273+
end if
274+
#:endif
247275

248276
! Can A be overwritten? By default, do not overwrite
249277
copy_a = .true._lk
250278
if (present(overwrite_a)) copy_a = .not.overwrite_a
251-
279+
252280
! Initialize a matrix temporary
253281
if (copy_a) then
254282
allocate(amat(m,n),source=a)
255283
else
256284
amat => a
257285
endif
258286

287+
#:if ep=='generalized'
288+
! Can B be overwritten? By default, do not overwrite
289+
copy_b = .true._lk
290+
if (present(overwrite_b)) copy_b = .not.overwrite_b
291+
292+
! Initialize a matrix temporary
293+
if (copy_b) then
294+
allocate(bmat,source=b)
295+
else
296+
bmat => b
297+
endif
298+
allocate(beta(n))
299+
#:endif
300+
259301
! Decide if U, V eigenvectors
260302
task_u = eigenvectors_task(present(left))
261303
task_v = eigenvectors_task(present(right))
@@ -304,7 +346,7 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues
304346
umat => u_dummy
305347
endif
306348

307-
get_geev: if (err0%ok()) then
349+
get_${ei}$: if (err0%ok()) then
308350

309351
ldu = size(umat,1,kind=ilp)
310352
ldv = size(vmat,1,kind=ilp)
@@ -318,11 +360,17 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues
318360

319361
lwork = -1_ilp
320362

321-
call geev(task_u,task_v,n,amat,lda,&
363+
call ${ei}$(task_u,task_v,n,amat,lda,&
364+
#:if ep=='generalized'
365+
bmat,ldb, &
366+
#:endif
322367
#{if rt.startswith('complex')}#lambda,#{else}#lreal,limag,#{endif}# &
368+
#:if ep=='generalized'
369+
beta, &
370+
#:endif
323371
umat,ldu,vmat,ldv,&
324372
work_dummy,lwork,#{if rt.startswith('complex')}#rwork,#{endif}#info)
325-
call handle_${ei}$_info(err0,info,[m,n]#{if ei=='ggev'}#,shape(bmat)#{endif}#)
373+
call handle_${ei}$_info(err0,info,shape(amat)#{if ei=='ggev'}#,shape(bmat)#{endif}#)
326374

327375
! Compute eigenvalues
328376
if (info==0) then
@@ -332,30 +380,46 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues
332380
allocate(work(lwork))
333381

334382
!> Compute eigensystem
335-
call geev(task_u,task_v,n,amat,lda,&
383+
call ${ei}$(task_u,task_v,n,amat,lda,&
384+
#:if ep=='generalized'
385+
bmat,ldb, &
386+
#:endif
336387
#{if rt.startswith('complex')}#lambda,#{else}#lreal,limag,#{endif}# &
388+
#:if ep=='generalized'
389+
beta, &
390+
#:endif
337391
umat,ldu,vmat,ldv,&
338392
work,lwork,#{if rt.startswith('complex')}#rwork,#{endif}#info)
339-
call handle_${ei}$_info(err0,info,[m,n]#{if ei=='ggev'}#,shape(bmat)#{endif}#)
393+
call handle_${ei}$_info(err0,info,shape(amat)#{if ei=='ggev'}#,shape(bmat)#{endif}#)
340394

341395
endif
342396

343397
! Finalize storage and process output flag
344398
#:if not rt.startswith('complex')
345399
lambda(:n) = cmplx(lreal(:n),limag(:n),kind=${rk}$)
400+
#:endif
401+
402+
#:if ep=='generalized'
403+
! Scale generalized eigenvalues
404+
lambda(:n) = scale_general_eig(lambda(:n),beta)
405+
#:endif
346406

407+
#:if not rt.startswith('complex')
347408
! If the j-th and (j+1)-st eigenvalues form a complex conjugate pair,
348-
! geev returns reals as:
409+
! ${ei}$ returns reals as:
349410
! u(j) = VL(:,j) + i*VL(:,j+1) and
350411
! u(j+1) = VL(:,j) - i*VL(:,j+1).
351412
! Convert these to complex numbers here.
352413
if (present(right)) call assign_real_eigenvectors_${rk}$(n,lambda,vmat,right)
353414
if (present(left)) call assign_real_eigenvectors_${rk}$(n,lambda,umat,left)
354415
#:endif
355416

356-
endif get_geev
417+
endif get_${ei}$
357418

358419
if (copy_a) deallocate(amat)
420+
#:if ep=='generalized'
421+
if (copy_b) deallocate(bmat)
422+
#:endif
359423
#:if not rt.startswith('complex')
360424
if (present(right)) deallocate(vmat)
361425
if (present(left)) deallocate(umat)
@@ -365,7 +429,7 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues
365429
end subroutine stdlib_linalg_eig_${ep}$_${ri}$
366430

367431
#:endfor
368-
432+
369433
module function stdlib_linalg_eigvalsh_${ri}$(a,upper_a,err) result(lambda)
370434
!! Return an array of eigenvalues of real symmetric / complex hermitian A
371435
!> Input matrix A[m,n]
@@ -391,7 +455,7 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues
391455
allocate(lambda(k))
392456

393457
!> Compute eigenvalues only
394-
call stdlib_linalg_eigh_${ri}$(amat,lambda,upper_a=upper_a,overwrite_a=.false.,err=err)
458+
call stdlib_linalg_eigh_${ri}$(amat,lambda,upper_a=upper_a,err=err)
395459

396460
end function stdlib_linalg_eigvalsh_${ri}$
397461

@@ -582,12 +646,18 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues
582646

583647
end subroutine assign_real_eigenvectors_${rk}$
584648

585-
module subroutine stdlib_linalg_real_eig_${ri}$(a,lambda,right,left,overwrite_a,err)
649+
#:for ep,ei in EIG_PROBLEM_LIST
650+
module subroutine stdlib_linalg_real_eig_${ep}$_${ri}$(a#{if ei=='ggev'}#,b#{endif}#,lambda,right,left, &
651+
overwrite_a#{if ei=='ggev'}#,overwrite_b#{endif}#,err)
586652
!! Eigendecomposition of matrix A returning an array `lambda` of real eigenvalues,
587653
!! and optionally right or left eigenvectors. Returns an error if the eigenvalues had
588654
!! non-trivial imaginary parts.
589655
!> Input matrix A[m,n]
590656
${rt}$, intent(inout), target :: a(:,:)
657+
#:if ei=='ggev'
658+
!> Generalized problem matrix B[n,n]
659+
${rt}$, intent(inout), target :: b(:,:)
660+
#:endif
591661
!> Array of real eigenvalues
592662
real(${rk}$), intent(out) :: lambda(:)
593663
!> The columns of RIGHT contain the right eigenvectors of A
@@ -596,6 +666,10 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues
596666
complex(${rk}$), optional, intent(out), target :: left(:,:)
597667
!> [optional] Can A data be overwritten and destroyed?
598668
logical(lk), optional, intent(in) :: overwrite_a
669+
#:if ei=='ggev'
670+
!> [optional] Can B data be overwritten and destroyed? (default: no)
671+
logical(lk), optional, intent(in) :: overwrite_b
672+
#:endif
599673
!> [optional] state return flag. On error if not requested, the code will stop
600674
type(linalg_state_type), optional, intent(out) :: err
601675

@@ -608,7 +682,8 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues
608682
n = size(lambda,dim=1,kind=ilp)
609683
allocate(clambda(n))
610684

611-
call stdlib_linalg_eig_standard_${ri}$(a,clambda,right,left,overwrite_a,err0)
685+
call stdlib_linalg_eig_${ep}$_${ri}$(a#{if ei=='ggev'}#,b#{endif}#,clambda,right,left, &
686+
overwrite_a#{if ei=='ggev'}#,overwrite_b#{endif}#,err0)
612687

613688
! Check that no eigenvalues have meaningful imaginary part
614689
if (err0%ok() .and. any(aimag(clambda)>atol+rtol*abs(abs(clambda)))) then
@@ -621,8 +696,37 @@ submodule (stdlib_linalg) stdlib_linalg_eigenvalues
621696

622697
call linalg_error_handling(err0,err)
623698

624-
end subroutine stdlib_linalg_real_eig_${ri}$
699+
end subroutine stdlib_linalg_real_eig_${ep}$_${ri}$
625700

701+
#:endfor
626702
#:endfor
703+
704+
#:for rk,rt,ri in RC_KINDS_TYPES
705+
!> Utility function: Scale generalized eigenvalue
706+
elemental complex(${rk}$) function scale_general_eig_${ri}$(alpha,beta) result(lambda)
707+
!! A generalized eigenvalue for a pair of matrices (A,B) is a scalar lambda or a ratio
708+
!! alpha/beta = lambda, such that A - lambda*B is singular. It is usually represented as the
709+
!! pair (alpha,beta), there is a reasonable interpretation for beta=0, and even for both
710+
!! being zero.
711+
complex(${rk}$), intent(in) :: alpha
712+
${rt}$, intent(in) :: beta
713+
714+
real (${rk}$), parameter :: rzero = 0.0_${rk}$
715+
complex(${rk}$), parameter :: czero = (0.0_${rk}$,0.0_${rk}$)
716+
717+
if (beta==#{if rt.startswith('real')}#r#{else}#c#{endif}#zero) then
718+
if (alpha/=czero) then
719+
lambda = cmplx(ieee_value(1.0_${rk}$, ieee_positive_inf), &
720+
ieee_value(1.0_${rk}$, ieee_positive_inf), kind=${rk}$)
721+
else
722+
lambda = ieee_value(1.0_${rk}$, ieee_quiet_nan)
723+
end if
724+
else
725+
lambda = alpha/beta
726+
end if
727+
728+
end function scale_general_eig_${ri}$
729+
730+
#:endfor
627731

628732
end submodule stdlib_linalg_eigenvalues

0 commit comments

Comments
 (0)