! Copyright (c) 2004-2024 Lars Nerger, lars.nerger@awi.de ! ! This routine is free software: you can redistribute it and/or modify ! it under the terms of the GNU Lesser General Public License ! as published by the Free Software Foundation, either version ! 3 of the License, or (at your option) any later version. ! ! PDAF is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! GNU Lesser General Public License for more details. ! ! You should have received a copy of the GNU Lesser General Public ! License along with this software. If not, see <http://www.gnu.org/licenses/>. ! !$Id$ !BOP ! ! !ROUTINE: PDAF_EOFCovar --- Perform EOF decomposition of state trajectory ! ! !INTERFACE: SUBROUTINE PDAF_eofcovar(dim, nstates, nfields, dim_fields, offsets, & remove_mstate, do_mv, states, stddev, svals, & svec, meanstate, verbose, status) USE PDAF_mod_filter, & ONLY: debug ! !DESCRIPTION: ! This routine performs an EOF analysis by singular value decomposition. It is ! used to prepare a covariance matrix for initializing an ensemble. For ! the decomposition a multivariate scaling can be performed by ! 'PDAF\_MVNormalize' to ensure that all fields in the state vectors have ! unit variance. ! ! To use this routine, one has to initialize the array 'states' holding in ! each column a perturbation vector (state - mean) from a state trajectory. ! Outputs are the arrays of singular values (svals) and left singular vectors ! (svec). The singular values are scaled by sqrt(1/(nstates-1)). With this, ! $svec * svals^2 * svec^T$ is the covariance matrix. In addition, the standard ! deviation of the field variations (stddev) is an output array. ! To use the multivariate normalization one has to define the number of ! different fields in the state (nfields), the dimension of each fields and ! the offset of field from the start of each state vector. ! ! The routine uses the LAPACK routine 'dgesvd' to compute the singular value ! decomposition. ! ! !REVISION HISTORY: ! 2012-09 - L. Nerger - Initial code for SANGOMA based on PDAF ! 2013-11 - L. Nerger - Adaption to SANGOMA data model ! 2016-05 - L. Nerger - Back-porting to PDAF ! 2019-11 - L. Nerger - Clarification that 'states' is destroyed by the SVD ! ! !USES: ! Include definitions for real type of different precision ! (Defines BLAS/LAPACK routines and MPI_REALTYPE) #include "typedefs.h" IMPLICIT NONE ! !ARGUMENTS: INTEGER, INTENT(in) :: dim ! Dimension of state vector INTEGER, INTENT(in) :: nstates ! Number of state vectors INTEGER, INTENT(in) :: nfields ! Number of fields in state vector INTEGER, INTENT(in) :: dim_fields(nfields) ! Size of each field INTEGER, INTENT(in) :: offsets(nfields) ! Start position of each field INTEGER, INTENT(in) :: do_mv ! 1: Do multivariate scaling; 0: no scaling ! nfields, dim_fields and offsets are only used if do_mv=1 INTEGER, INTENT(in) :: remove_mstate ! 1: subtract mean state from states ! before computing EOFs; 0: don't remove and don't touch meanstate REAL, INTENT(inout) :: states(dim, nstates) ! State perturbations ! the array content will be destroyed by the singular value decomposition REAL, INTENT(out) :: stddev(nfields) ! Standard deviation of field variability ! Without multivariate scaling (do_mv=0), it is stddev = 1.0 REAL, INTENT(out) :: svals(nstates) ! Singular values divided by sqrt(nstates-1) REAL, INTENT(out) :: svec(dim, nstates) ! Singular vectors REAL, INTENT(inout) :: meanstate(dim) ! Mean state (only changed if remove_mstate=1) INTEGER, INTENT(in) :: verbose ! Verbosity flag INTEGER, INTENT(out) :: status ! Status flag !EOP ! *** local variables *** INTEGER :: i ! Counter INTEGER :: stat ! internal status flag INTEGER :: ldwork ! variable for SVD routine REAL, ALLOCATABLE :: work(:) ! work array for SVD REAL :: svdV ! right singular values (not referenced) ! ********************** ! *** INITIALIZATION *** ! ********************** IF (verbose>0) THEN WRITE (*,'(a, 5x,a)') 'PDAF', '***********************************************' WRITE (*,'(a, 5x,a)') 'PDAF', '* PDAF_EOFCovar *' WRITE (*,'(a, 5x,a)') 'PDAF', '* *' WRITE (*,'(a, 5x,a)') 'PDAF', '* Compute EOF decomposition of a matrix *' WRITE (*,'(a, 5x,a)') 'PDAF', '* based on the Sangoma tool Sangoma_EOFCovar. *' WRITE (*,'(a, 5x,a)') 'PDAF', '***********************************************' END IF IF (debug>0) THEN WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_eofcovar -- START' WRITE (*,*) '++ PDAF-debug PDAF_eofcovar:', debug, ' dim', dim WRITE (*,*) '++ PDAF-debug PDAF_eofcovar:', debug, ' nstates', nstates WRITE (*,*) '++ PDAF-debug PDAF_eofcovar:', debug, ' nfields', nfields WRITE (*,*) '++ PDAF-debug PDAF_eofcovar:', debug, ' dim_fields', dim_fields WRITE (*,*) '++ PDAF-debug PDAF_eofcovar:', debug, ' offsets', offsets WRITE (*,*) '++ PDAF-debug PDAF_eofcovar:', debug, ' states(1,:)', states(1,:) WRITE (*,*) '++ PDAF-debug PDAF_eofcovar:', debug, & ' Note: If REAL values appear incorrect, please check if you provide them with the correct precision' END IF ! ************************* ! *** Remove mean state *** ! ************************* removemean: IF (remove_mstate == 1) THEN IF (verbose>0) & WRITE (*,'(/a,5x,a)') 'PDAF', 'EOFCOVAR: Compute and subtract mean state ----------------------' ! *** compute mean state *** meanstate = 0.0 DO i = 1, nstates meanstate(:) = meanstate(:) + states(:, i) / REAL(nstates) END DO ! *** get peturbation matrix *** DO i = 1, nstates states(:,i) = states(:,i) - meanstate(:) END DO END IF removemean ! ****************************************** ! *** Perform multivariate normalization *** ! ****************************************** stat = 0 multivar: IF (do_mv == 1) THEN IF (verbose>0) & WRITE (*,'(/a,5x,a)') 'PDAF', 'EOFCOVAR: Perform multivariate normalization -------------------' DO i = 1, nfields CALL PDAF_MVNormalize(1, dim, dim_fields(i), offsets(i), nstates, & states, stddev(i), status) if (verbose>0) & WRITE (*,'(a,5x,a,i5,a,es12.4)') 'PDAF', 'Field', i, ': standard deviation ', stddev(i) stat = stat + status END DO ELSE stddev = 1.0 END IF multivar ! ********************************************************* ! *** Singular value decomposition of covariance matrix *** ! *** *** ! *** The covariance matrix is given by the state *** ! *** sequences X of k states as *** ! *** -1 _ _ T T *** ! *** P = (k-1) (X-X) (X-X) = U L U (EVP) *** ! *** *** ! *** We compute the singular value decomposition *** ! *** _ T -1 2 T *** ! *** X-X = U S V ; P = (k-1) U S U *** ! *** *** ! *** -1/2 *** ! *** and we return U and (k-1) S. *** ! ********************************************************* IF (stat==0) THEN IF (verbose>0) & WRITE (*,'(/a,5x,a)') 'PDAF', 'EOFCOVAR: Compute SVD ------------------------------------------' ! Allocate work array and work space size ALLOCATE(work(MAX(3 * MIN(dim, nstates) + & MAX(dim, nstates), 5 * MIN(dim, nstates)))) ldwork = MAX(3 * MIN(dim, nstates) + & MAX(dim, nstates), 5 * MIN(dim, nstates)) ! Do decomposition CALL gesvdTYPE('s', 'n', dim, nstates, states, & dim, svals, svec, dim, svdV, & dim, work, ldwork, status) ! *** scale singular values *** DO i = 1, nstates svals(i) = svals(i) / SQRT(REAL(nstates - 1)) END DO stat = stat+status IF (debug>0) THEN WRITE (*,*) '++ PDAF-debug PDAF_eofcovar:', debug, ' svals', svals END IF END IF ! ******************************** ! *** Rescale singular vectors *** ! ******************************** do_rescale: IF (do_mv == 1 .AND. stat == 0) THEN IF (verbose>0) & WRITE (*,'(/a, 5x,a)') 'PDAF', 'EOFCOVAR: Re-scale singular vectors according to stddev --------' DO i = 1, nfields CALL PDAF_MVNormalize(2, dim, dim_fields(i), offsets(i), nstates-1, & svec, stddev(i), status) END DO stat = stat+status END IF do_rescale ! **************** ! *** Clean up *** ! **************** DEALLOCATE(work) status = stat IF (debug>0) THEN WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_eofcovar -- END' END IF END SUBROUTINE PDAF_eofcovar