! (C) Copyright 2015- ECMWF. ! (C) Copyright 2015- Meteo-France. ! ! This software is licensed under the terms of the Apache Licence Version 2.0 ! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. ! In applying this licence, ECMWF does not waive the privileges and immunities ! granted to it by virtue of its status as an intergovernmental organisation ! nor does it submit to any jurisdiction. ! MODULE INTERPOL_DECOMP_MOD ! Compute Interpolative Decomposions (ID) ! See Cheng,H., Gimbutas,Z., Martinsson,P.G. and Rokhlin,V. (2005) ! "On the compression of low rank matrices", SIAM.J.Sci.Comput., ! Vol. 26, No. 4, pp1389-1404 ! Also lecture notes "Mulilevel compression of Linear Operators: ! Descendents of Fast Multiple Methods and Calderon-Zygmund Theory" ! P.G.Martinsson and Mark Tygert, 2011. Chapter 7. ! Author: Mats Hamrud USE EC_PARKIND, ONLY : JPIM, JPRD, JPIB IMPLICIT NONE CONTAINS !=========================================================================== SUBROUTINE COMPUTE_ID(PEPS,KM,KN,PMAT,KRANK,KBCLIST,PNONIM) IMPLICIT NONE ! Compute ID REAL(KIND=JPRD),INTENT(IN) :: PEPS ! Precision for computation ! of numerical rank INTEGER(KIND=JPIM),INTENT(IN) :: KM ! Number of rows in matrix pmat INTEGER(KIND=JPIM),INTENT(IN) :: KN ! Number of columns in matrix pmat REAL(KIND=JPRD) ,INTENT(IN) :: PMAT(:,:) ! Original matrix INTEGER(KIND=JPIM),INTENT(OUT) :: KRANK ! Numerical rank INTEGER(KIND=JPIM),INTENT(OUT) :: KBCLIST(:) ! List of columns REAL(KIND=JPRD) ,INTENT(OUT) :: PNONIM(:,:) ! Non-identity part of projection ! matrix INTEGER(KIND=JPIM) :: JM,JN REAL(KIND=JPRD) :: ZR(KM,KN) REAL(KIND=JPRD),ALLOCATABLE :: ZS(:,:),ZT(:,:) !---------------------------------------------------------------------------- !Avoid destroying input matrix ZR(:,:) = PMAT(1:KM,1:KN) ! Householder QR CALL ALG541(PEPS,KM,KN,ZR,KRANK,KBCLIST) DO JN=1,KN DO JM=JN+1,KM ZR(JM,JN) = 0.0_JPRD ENDDO ENDDO ! S leftmost kxk block of R ALLOCATE(ZS(KRANK,KRANK)) DO JN=1,KRANK DO JM=1,KRANK IF(JM <= KM ) THEN ZS(JM,JN) = ZR(JM,JN) ELSE ZS(JM,JN) = 0.0_JPRD ENDIF ENDDO ENDDO ! T Rightmost kx(k-n) block of R ALLOCATE(ZT(KRANK,KN-KRANK)) DO JN=1,KN-KRANK DO JM=1,KRANK IF(JM <= KM ) THEN ZT(JM,JN) = ZR(JM,JN+KRANK) ELSE ZT(JM,JN) = 0.0_JPRD ENDIF ENDDO ENDDO !Solve linear equation (BLAS level 3 routine) IF( KRANK <= 0 ) THEN write(0,*) 'warning: KRANK DTRSM ', KRANK, KM, KN CALL ABOR1('DTRSM : KRANK <=0 not allowed') ENDIF ! IF (JPRB == JPRD) THEN CALL DTRSM('Left','Upper','No transpose','Non-unit',KRANK,KN-KRANK,1.0_JPRD, & & ZS,KRANK,ZT,KRANK) ! ELSE ! CALL STRSM('Left','Upper','No transpose','Non-unit',KRANK,KN-KRANK,1.0_JPRD, & ! & ZS,KRANK,ZT,KRANK) ! ENDIF DO JM=1,KRANK DO JN=1,KN-KRANK PNONIM(JM,JN) = ZT(JM,JN) ENDDO ENDDO DEALLOCATE(ZS,ZT) !!$IF(KRANK < KN) THEN !!$ PRINT *,'MAXVAL PNONIM ',KM,KM,KRANK,MAXVAL( PNONIM(1:KRANK,1:KN-KRANK)) !!$ENDIF END SUBROUTINE COMPUTE_ID !============================================================================== SUBROUTINE ALG541(PEPS,KM,KN,PA,KRANK,KLIST) IMPLICIT NONE ! Householder QR with Column Pivoting ! Algorithm 5.4.1 from Matrix Computations, G.H.Golub & C.F van Loen, third ed. ! Algorithm modified to terminate at numerical precision "peps" REAL(KIND=JPRD),INTENT(IN) :: PEPS ! Precision INTEGER(KIND=JPIM),INTENT(IN) :: KM ! Number of rows in matrix pa INTEGER(KIND=JPIM),INTENT(IN) :: KN ! Number of columns in matrix pa REAL(KIND=JPRD),INTENT(INOUT) :: PA(:,:) ! On input : original matrix ! on output : R in upper triangle etc ! see Golub&Van Loen INTEGER(KIND=JPIM),INTENT(OUT) :: KRANK ! Numerical rank of matrix INTEGER(KIND=JPIM),INTENT(OUT) :: KLIST(:) ! List of columns (pivots) INTEGER(KIND=JPIM) :: JN,ISWAP,IK,IM,IN,IMIN,ILIST(KN) REAL(KIND=JPRD) :: ZC(KN),ZTAU,ZSWAPA(KM),ZSWAP,ZV(KM),ZBETA,ZWORK(KN),ZTAU_IN REAL(KIND=JPRD) :: ZTAU_REC,ZEPS !------------------------------------------------------------------------------- ZEPS = 10000.0_JPRD*EPSILON(ZEPS) IMIN=MIN(KM,KN) ! Compute initial column norms,its max and the first column where c=tau IK = 0 ZTAU = 0._JPRD DO JN=1,KN ZC(JN) = DOT_PRODUCT(PA(1:KM,JN),PA(1:KM,JN)) IF(ZC(JN) > ZTAU) THEN IK = JN ZTAU = ZC(JN) ENDIF ENDDO ZTAU_IN = ZTAU ZTAU_REC= ZTAU KRANK = 0 DO WHILE (ZTAU > PEPS**2*ZTAU_IN) KRANK = KRANK+1 IF( KRANK <= IMIN ) THEN ILIST(KRANK) = IK ! Column swap KRANK with IK ZSWAPA(:) = PA(:,KRANK) PA(:,KRANK) = PA(:,IK) PA(:,IK) = ZSWAPA(:) ZSWAP = ZC(KRANK) ZC(KRANK) = ZC(IK) ZC(IK) = ZSWAP ! Compute Householder vector ZBETA=0.0_JPRD IF( KM-KRANK >= 0 ) THEN CALL ALG511(ZEPS,KM-KRANK+1,PA(KRANK:KM,KRANK),ZV,ZBETA) ENDIF ! Apply Householder matrix IM = KM-KRANK+1 IN = KN-KRANK+1 ! LAPACK CALL DLARF('Left',IM,IN,ZV,1,ZBETA,PA(KRANK,KRANK),KM,ZWORK) ENDIF ! Update column norms ZTAU = 0.0_JPRD IF(KRANK < IMIN) THEN PA(KRANK+1:KM,KRANK) = ZV(2:IM) DO JN=KRANK+1,KN ZC(JN) = ZC(JN)-PA(KRANK,JN)**2 IF(ZC(JN) > ZTAU) THEN IK = JN ZTAU = ZC(JN) ENDIF ENDDO ! Re-compute column norms due to round-off error IF(ZTAU < ZEPS*ZTAU_REC .OR. ZTAU < 0._JPRD .or. (KN-KRANK) > 100 ) THEN DO JN=KRANK+1,KN ZC(JN) = DOT_PRODUCT(PA(KRANK+1:,JN),PA(KRANK+1:,JN)) IF(ZC(JN) > ZTAU) THEN IK = JN ZTAU = ZC(JN) ENDIF ENDDO !write(0,*) 'RECOMPUTE TAU ',KRANK,ZTAU_REC,ZTAU ZTAU_REC = ZTAU ENDIF ENDIF ENDDO ! Make sure klist is filled also beyond krank DO JN=1,KN KLIST(JN) = JN ENDDO DO JN=1,KRANK ISWAP = KLIST(JN) KLIST(JN) = KLIST(ILIST(JN)) KLIST(ILIST(JN)) = ISWAP ENDDO END SUBROUTINE ALG541 !============================================================================== SUBROUTINE ALG511(PEPS,KSIZE,PX,PV,PBETA) IMPLICIT NONE ! Compute Householder vector ! Algorithm 5.1.1 from Matrix Computations, G.H.Golub & C.F van Loen, third ed. REAL(KIND=JPRD),INTENT(IN) :: PEPS ! Precision REAL(KIND=JPRD),INTENT(IN) :: PX(:) INTEGER(KIND=JPIM), INTENT(IN) :: KSIZE REAL(KIND=JPRD),INTENT(OUT) :: PV(:) REAL(KIND=JPRD),INTENT(OUT) :: PBETA INTEGER(KIND=JPIM) :: IL REAL(KIND=JPRD) :: ZSIGMA,ZMU, ZNORM REAL(KIND=JPRD) :: ZX(KSIZE) !------------------------------------------------------------------------------- ! normalize ZNORM=0._JPRD DO IL=1,KSIZE ZNORM = ZNORM + PX(IL)*PX(IL) ENDDO ZNORM=SQRT(ZNORM) ZX(:)=PX(1:KSIZE) IF( ZNORM > PEPS ) ZX(:)=PX(1:KSIZE)/ZNORM ZSIGMA=0._JPRD IF( KSIZE > 1 ) ZSIGMA = DOT_PRODUCT(ZX(2:KSIZE),ZX(2:KSIZE)) PV(1) = 1.0_JPRD IF( KSIZE > 1 ) PV(2:KSIZE) = ZX(2:KSIZE) IF(ABS(ZSIGMA) < PEPS**2) THEN PBETA = 0.0_JPRD ELSE ZMU = SQRT(ZX(1)**2+ZSIGMA) IF(ZX(1) <= 0.0_JPRD) THEN PV(1) = ZX(1)-ZMU ELSE PV(1) = -ZSIGMA/(ZX(1)+ZMU) ENDIF PBETA = 2.0_JPRD*PV(1)**2/(ZSIGMA+PV(1)**2) PV(:) = PV(:)/(PV(1)) ENDIF END SUBROUTINE ALG511 !================================================================================ END MODULE INTERPOL_DECOMP_MOD