! Copyright (c) 2004-2024 Lars Nerger ! ! This file is part of PDAF. ! ! PDAF 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 PDAF. If not, see <http://www.gnu.org/licenses/>. ! !$Id$ !BOP ! ! !ROUTINE: PDAF_seek_rediag --- Perform rediagonalization of P in SEEK ! ! !INTERFACE: SUBROUTINE PDAF_seek_rediag(dim_p, dim_eof, Uinv, eofV_p, subtype, & screen, flag) ! !DESCRIPTION: ! Re-orthogonalization of the modes V of the ! low-rank approximated covariance matrix in ! its decomposed form P = V U V$^T$. ! ! Compute eigenmodes of the matrix B = L$^T$ L = C D C$^T$ ! where L = V U$^{1/2}$ (from Cholesky decomposition) ! and get new modes V as V = L C D$^{-1/2}$, ! $D = diag(\lambda_1,...,\lambda_r),\ \lambda_i > \lambda_i+1$ ! and C matrix of corresponding eigenvectors. ! The new U is given by the matrix D. ! ! Variant for domain decomposed states. ! ! New version to compute matrix B. More efficient for ! dim $>>$ dim\_eof ! ! ! This is a core routine of PDAF and ! should not be changed by the user ! ! ! !REVISION HISTORY: ! 2003-10 - Lars Nerger - Initial code ! Later revisions - see svn log ! ! !USES: ! Include definitions for real type of different precision ! (Defines BLAS/LAPACK routines and MPI_REALTYPE) #include "typedefs.h" USE mpi USE PDAF_timer, & ONLY: PDAF_timeit USE PDAF_memcounting, & ONLY: PDAF_memcount USE PDAF_mod_filtermpi, & ONLY: mype, MPIerr, COMM_filter IMPLICIT NONE ! !ARGUMENTS: INTEGER, INTENT(in) :: dim_p ! PE-Local state dimension INTEGER, INTENT(in) :: dim_eof ! Number of EOFs REAL, INTENT(inout) :: Uinv(dim_eof,dim_eof) ! Inverse of matrix U REAL, INTENT(inout) :: eofV_p(dim_p,dim_eof) ! PE-local matrix V INTEGER, INTENT(in) :: subtype ! Filter subtype INTEGER, INTENT(in) :: screen ! Verbosity flag INTEGER, INTENT(inout) :: flag ! Status Flag ! !CALLING SEQUENCE: ! Called by: PDAF_seek_update ! Calls: PDAF_timeit ! Calls: PDAF_memcount ! Calls: gemmTYPE (BLAS; dgemm or sgemm dependent on precision) ! Calls: potrfTYPE (LAPACK; dpotrf ot spotrf dependent on precision) ! Calls: gesvTYPE (LAPACK; dgesv or sgesv dependent on precision) ! Calls: MPI_allreduce (MPI) !EOP ! *** local variables *** INTEGER :: row, col ! counters INTEGER, SAVE :: allocflag = 0 ! Flag whether first time allocation is done INTEGER, ALLOCATABLE :: ipiv(:) ! vector of pivot indices for SGESV REAL, ALLOCATABLE :: ev(:) ! vector of eigenvalues REAL, ALLOCATABLE :: rwork(:) ! workarray for eigenproblem REAL, ALLOCATABLE :: U(:,:) ! temporary for covar matrix REAL, ALLOCATABLE :: L(:,:) ! covariance matrix REAL, ALLOCATABLE :: Temp1(:,:) ! temporary for covar matrix REAL, ALLOCATABLE :: B(:,:) ! temporary for covar matrix INTEGER :: syev_info, gesv_info, potrf_info ! Info flags for LAPACK calls ! ********************** ! *** INITIALIZATION *** ! ********************** IF (mype == 0 .AND. screen > 0) THEN WRITE (*, '(a, 5x, a)') 'PDAF', 'Re-orthogonalize covariance matrix modes' END IF ! ************************************** ! *** Compute matrix B = A^T V^T V A *** ! ************************************** CALL PDAF_timeit(20, 'new') ! *** Get U by inversion of Uinv *** ALLOCATE(U(dim_eof, dim_eof)) ALLOCATE(ipiv(dim_eof)) IF (allocflag == 0) CALL PDAF_memcount(4, 'r', dim_eof * dim_eof) IF (allocflag == 0) CALL PDAF_memcount(4, 'i', dim_eof) ! Initialize U U = 0.0 DO row = 1, dim_eof U(row, row) = 1.0 END DO ! call solver CALL PDAF_timeit(32, 'new') CALL gesvTYPE(dim_eof, dim_eof, Uinv, dim_eof, ipiv, & U, dim_eof, gesv_info) CALL PDAF_timeit(32, 'old') DEALLOCATE(ipiv) ! Check if inversion was successful INVok: IF (gesv_info == 0) THEN ! *** Cholesky decomposition of U: U = A A^T *** CALL PDAF_timeit(33, 'new') CALL potrfTYPE('l', dim_eof, U, dim_eof, potrf_info) CALL PDAF_timeit(33, 'old') ! *** set upper elements to zero *** DO col = 2, dim_eof DO row = 1, col - 1 U(row, col) = 0.0 END DO END DO !*** Compute B = A^T V^T V A *** ALLOCATE(Temp1(dim_eof, dim_eof)) ALLOCATE(B(dim_eof, dim_eof)) IF (allocflag == 0) CALL PDAF_memcount(4, 'r', 2 * dim_eof * dim_eof) ! local V^T V CALL gemmTYPE('t', 'n', dim_eof, dim_eof, dim_p, & 1.0, eofV_p, dim_p, eofV_p, dim_p, & 0.0, Temp1, dim_eof) CALL PDAF_timeit(20, 'old') CALL PDAF_timeit(21, 'new') CALL MPI_allreduce(Temp1, B, dim_eof * dim_eof, MPI_REALTYPE, & MPI_SUM, COMM_filter, MPIerr) CALL PDAF_timeit(21, 'old') CALL PDAF_timeit(20, 'new') ! (V^T V) A (A stored in U) CALL gemmTYPE('n', 'n', dim_eof, dim_eof, dim_eof, & 1.0, B, dim_eof, U, dim_eof, & 0.0, Temp1, dim_eof) ! B = A^T (V^T V A) (A stored in U) CALL gemmTYPE('t', 'n', dim_eof, dim_eof, dim_eof, & 1.0, U, dim_eof, Temp1, dim_eof, & 0.0, B, dim_eof) CALL PDAF_timeit(31, 'old') ! ******************************* ! *** Eigendecomposition of B *** ! *** B = C D C^T *** ! ******************************* ALLOCATE(ev(dim_eof)) ALLOCATE(rwork(3 * dim_eof)) IF (allocflag == 0) CALL PDAF_memcount(4, 'r', 4 * dim_eof) CALL syevTYPE('v', 'u', dim_eof, B, dim_eof, & ev, rwork, 3 * dim_eof, syev_info) CALL PDAF_timeit(20, 'old') DEALLOCATE(rwork) ! check if eigendecomposition was successful EVPok: IF (syev_info == 0) THEN ! Eigendecomposition OK, continue ! *** Reorder matrix of eigenvectors *** ! eof_mid = floor(real(dim_eof)/2.0) ! do col=1,eof_mid ! rwork(1:dim_eof) = U(:,col) ! U(:,col) = U(:,dim_eof-col+1) ! U(:,dim_eof-col+1) = rwork(1:dim_eof) ! end do ! ***************************************************** ! *** Initialize covar with re-orthogonalized modes *** ! ***************************************************** CALL PDAF_timeit(22, 'new') ! AC = A C (AC stored in Temp1, A stored in U, C stored in B) CALL gemmTYPE('n', 'n', dim_eof, dim_eof, dim_eof, & 1.0, U, dim_eof, B, dim_eof, & 0.0, Temp1, dim_eof) ! initialize L from V ALLOCATE(L(dim_p, dim_eof)) IF (allocflag == 0) CALL PDAF_memcount(4, 'r', dim_eof * dim_p) L = eofV_p ! V = L AC (AC stored in Temp1) CALL gemmTYPE('n', 'n', dim_p, dim_eof, dim_eof, & 1.0, L, dim_p, Temp1, dim_eof, & 0.0, eofV_p, dim_p) DEALLOCATE(L, U, Temp1, B) unitmodes: IF (subtype /= 1) THEN ! use eigenvalues in U with unit modes in V ! U = diag(lambda_1,...,lambda_r) IF (mype == 0 .AND. screen > 0) THEN WRITE (*, '(a, 5x, a)') 'PDAF', '--- Use normalized modes' END IF Uinv = 0.0 DO row = 1, dim_eof Uinv(row, row) = 1.0 / ev(row) ev(row) = SQRT(1.0 / ev(row)) END DO ! Rescale modes V DO col = 1, dim_eof DO row = 1, dim_p eofV_p(row, col) = eofV_p(row, col) * ev(col) END DO END DO ELSE unitmodes ! use unit U with modes in V scaled by eigenvalues IF (mype==0 .AND. screen > 0) THEN WRITE (*, '(a, 5x, a)') 'PDAF', '--- Use unit U' END IF Uinv = 0.0 DO row = 1, dim_eof Uinv(row, row) = 1.0 END DO END IF unitmodes CALL PDAF_timeit(22, 'old') ELSE ! eigendecomposition failed IF (mype == 0) WRITE (*, '(/5x, a/)') & 'PDAF-ERROR(1): Problem with EOF decomposition of matrix Lt L !!!' flag = 1 ENDIF EVPok DEALLOCATE(ev) ELSE INVok ! Inversion failed IF (mype == 0) WRITE (*, '(/5x, a/)') & 'PDAF-ERROR(2): Problem with inversion of Uinv !!!' flag = 2 ENDIF INVok ! ******************** ! *** Finishing up *** ! ******************** IF (allocflag == 0) allocflag = 1 END SUBROUTINE PDAF_seek_rediag