From 5df2c4be22f79bfccdd3ae9b0adb5daf2fb362c0 Mon Sep 17 00:00:00 2001
From: ronesy <rlt@nilu.no>
Date: Mon, 8 Jul 2024 15:55:31 +0200
Subject: [PATCH] Change to calculation of eigendecomposition of prior error
covariance to only calculate eigenvectors for the kept eigenvalues
---
source/congrad.f90 | 3 +-
source/error_cov.f90 | 69 ++++++++++++++++++++++++++++++++++++++++----
2 files changed, 65 insertions(+), 7 deletions(-)
diff --git a/source/congrad.f90 b/source/congrad.f90
index e1ef950..8c1373f 100644
--- a/source/congrad.f90
+++ b/source/congrad.f90
@@ -14,7 +14,8 @@
! You should have received a copy of the GNU General Public License
! along with FLEXINVERT. If not, see <http://www.gnu.org/licenses/>.
!
-! Adapted from the code of Mike Fisher (ECMWF) by Rona Thompson, 2017
+! Based on the code from Mike Fisher (ECMWF, 2002) and the adaptation by
+! Frederic Chevallier (LSCE, 2004) written in Fortran90 by Rona Thompson, 2017
!---------------------------------------------------------------------------------------
!
!> congrad
diff --git a/source/error_cov.f90 b/source/error_cov.f90
index 62c2014..7d9105d 100644
--- a/source/error_cov.f90
+++ b/source/error_cov.f90
@@ -74,6 +74,7 @@ subroutine error_cov(config, files, states, covar, corr, xerr)
real(kind=8), dimension(:,:), allocatable :: evecs
integer :: ierr, info, neval, ind
real, parameter :: smallnum=1.e-15
+ integer :: timei, timef, countrate
! calculate error covariance
! --------------------------
@@ -106,6 +107,13 @@ subroutine error_cov(config, files, states, covar, corr, xerr)
! calculate eigendecomposition of cov
! -----------------------------------
+ ! call to dsyevr to find only ILth to IUth eigenvalues and vectors
+ ! thus reducing time needed for eigen decomposition for large B
+ ! IL = index of smallest eigenvalue to be found
+ ! IU = index of largest eigenvalue to be found
+ ! where 1 <= IL <= IU <=N
+ ! IU = N and IL = max(1, N - f*N) where f is fraction of eigenvalues to keep
+
! B = U*D*U^(T)
lwork = 26*ndvar
liwork = 10*ndvar
@@ -121,24 +129,73 @@ subroutine error_cov(config, files, states, covar, corr, xerr)
if ( ierr.ne.0 ) stop 'ERROR error_cov: not enough memory'
allocate( zwork(ndvar,ndvar), stat=ierr )
if ( ierr.ne.0 ) stop 'ERROR error_cov: not enough memory'
+ zwork = cov
+
+ ! arguments to dsyevr
+ ! jobsz = N : N only eigenvalues, V both
+ ! range = A : A all eigenvalues, V values in interval (VL,VU], I from ILth to IUth values
+ ! uplo = L : lower triangle of matrix A stored
+ ! n = ndvar : order of matrix A
+ ! matrix = zwork : matrix A
+ ! lda = ndvar : leading dimension of A
+ ! vl = used if range is V : lower bound of interval to search for eigenvalues
+ ! vu = used if range is V : upper bound of interval to search for eigenvalues
+ ! il = used if range is I : index of smallest eigenvalue to be found
+ ! iu = used if range is I : index of largest eigenvalue to be found
+ ! absol = 0 : absolute tolerance for eigenvalues
+ ! m = neval : number of eigenvalues found
+ ! w = evals : eigenvalues in ascending order
+ ! z = evecs : found eigenvectors as columns
+ ! ldz = ndvar : leading dimension of z
+ ! isuppz = isuppz : matrix of indices to non-zero values in z
+ ! work = temporary working array
+ ! lwork = dimension of work
+ ! iwork = temporary working array
+ ! liwork = dimension of iwork
+ ! info : 0 success, < 0 if info = -i ith argument had an illegal value, > 0 internal error
+
+ ! first calculate only eigenvalues
+ call system_clock(timei, countrate)
+ call dsyevr('N','A','L',ndvar,zwork,ndvar,1.,1.,1,1,0.,neval,evals,evecs,ndvar,isuppz,work,lwork,iwork,liwork,info)
+ if ( info.ne.0 ) then
+ write(logid,*) 'ERROR error_cov: eigenvalue decomposition'
+ stop
+ endif
+ call system_clock(timef, countrate)
+ print*, 'error_cov: time to calculate eigenvalues (s) = ',(timef-timei)/countrate
+ ! now calculate eigenvalues in range IU = N and IL = max(1, N - f*N)
+ ! where f is fraction of eigenvalues to keep
+ ! and their corresponding eigenvectors
+ ind = minloc(evals, dim=1, mask=evals.ge.(config%trunc*maxval(evals)))
+ print*, 'error_cov: ind, evals(ind) = ',ind, evals(ind)
+ call system_clock(timei, countrate)
+ ! because on exit the lower triangle (if UPLO='L') or the upper triangle
+ ! (if UPLO='U') of A, including the diagonal, is destroyed
+ ! need to reassign matrix A
zwork = cov
- call dsyevr('V','A','L',ndvar,zwork,ndvar,1.,1.,1,1,0.,neval,evals,evecs,ndvar,isuppz,work,lwork,iwork,liwork,info)
+ evals(:) = 0.
+ evecs(:,:) = 0.
+ call dsyevr('V','I','L',ndvar,zwork,ndvar,1.,1.,ind,ndvar,0.,neval,evals,evecs,ndvar,isuppz,work,lwork,iwork,liwork,info)
if ( info.ne.0 ) then
write(logid,*) 'ERROR error_cov: eigenvalue decomposition'
stop
endif
+ call system_clock(timef, countrate)
+ print*, 'error_cov: time to calculate eigenvalues and eigenvectors (s) = ',(timef-timei)/countrate
+ print*, 'error_cov: neval = ',neval
! truncate eigenvalues
! eigenvalues are sorted in ascending order
- ind = minloc(evals, dim=1, mask=evals.ge.(config%trunc*maxval(evals)))
- write(logid,*) 'Truncating eigenvalues at: ',config%trunc*maxval(evals)
- neig = ndvar - ind + 1
+ neig = neval
write(logid,*) 'Number of retained eigenvalues: ',neig
call alloc_covar(covar)
- covar%evals = evals(ind:ndvar)
- covar%evecs = evecs(:,ind:ndvar)
+ ! eigenvalues in ascending order
+! covar%evals = evals(ind:ndvar)
+ covar%evals = evals(1:neig)
+ ! calculated eigenvectors in columns 1:neig
+ covar%evecs = evecs(:,1:neig)
deallocate(work)
deallocate(zwork)
--
GitLab