rep_matrix.F

Go to the documentation of this file.

#include "cppdefs.h"
#if defined WEAK_CONSTRAINT && defined OBS_IMPACT
      SUBROUTINE rep_matrix (ng, model, outLoop, NinnLoop)
!
!git $Id$
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2025 The ROMS Group            Andrew M. Moore   !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.md                                               !
!=======================================================================
!                                                                      !
!  This routine estimates the TRANSPOSE of the inverse of the          !
!  stabilized representer matrix using the Lanczos vectors from        !
!  the inner-loops.                                                    !
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
      USE mod_fourdvar
      USE mod_iounits
      USE mod_scalars
 
# ifdef DISTRIBUTE
!
      USE distribute_mod, ONLY : mp_bcastf
# endif
!
      implicit none
!
!  Imported variable declarations
!
      integer, intent(in) :: ng, model, outLoop, NinnLoop
!
!  Local variable declarations.
!
      integer :: i, j, iobs, ivec, innLoop
# ifdef IMPACT_INNER
      integer :: mloop
# endif
# ifdef MINRES
      integer :: info
      integer, dimension(NinnLoop) :: ipiv
# endif
 
      real(r8) :: zbet
 
      real(r8), dimension(NinnLoop) :: zu, zgam, zt
# ifdef MINRES
      real(r8) :: zsum, zck, zgk
      real(r8), dimension(NinnLoop,NinnLoop) :: ztriT, zLT, zLTt
      real(r8), dimension(NinnLoop) :: tau, zwork1, zeref
# endif
# ifdef RPCG
      real(r8), dimension(NinnLoop) :: zfact
# endif
# ifdef IMPACT_INNER
!
!-----------------------------------------------------------------------
!  Compute observation impact on the weak constraint 4DVAR data
!  assimilation systems for each inner-loop.
!-----------------------------------------------------------------------
!
      master_thread : IF (master) THEN
 
        DO mloop=1,ninnloop
 
          DO iobs=1,ndatum(ng)
            ad_obsval(iobs,mloop)=0.0_r8
          END DO
          DO innloop=1,ninnloop
            zt(innloop)=0.0_r8
          END DO
#  ifdef RPCG
          DO innloop=1,ninnloop
            IF (innloop.eq.1) THEN
              zfact(innloop)=1.0_r8/cg_qg(1,outloop)
            ELSE
              zfact(innloop)=1.0_r8/cg_beta(innloop,outloop)
            END IF
          END DO
#  endif
!
!  Multiply TLmodVal by the tranpose matrix of Lanczos vectors.
!  Note that the factor of 1/SQRT(ObsErr) is added to convert to
!  x-space.
!
          DO innloop=1,mloop
            DO iobs=1,ndatum(ng)
              IF (obserr(iobs).NE.0.0_r8) THEN
#  ifdef RPCG
                zt(innloop)=zt(innloop)+                                &
     &                      obsscale(iobs)*                             &
     &                      zcglwk(iobs,innloop,outloop)*               &
     &                      tlmodval(iobs)
#  else
                zt(innloop)=zt(innloop)+                                &
     &                      obsscale(iobs)*                             &
     &                      zcglwk(iobs,innloop,outloop)*               &
     &                      tlmodval(iobs)/                             &
     &                      sqrt(obserr(iobs))
#  endif
              END IF
            END DO
          END DO
 
#  ifdef MINRES
!
!  Use the minimum residual method as described by Paige and Saunders
!  ("Sparse Indefinite Systems of Linear Equations", 1975, SIAM Journal
!  on Numerical Analysis, 617-619). Specifically we refer to equations
!  6.10 and 6.11 of this paper.
!
!  Perform a LQ factorization of the tridiagonal matrix.
!
          ztrit=0.0_r8
          DO i=1,mloop
            ztrit(i,i)=cg_delta(i,outloop)
          END DO
          DO i=1,mloop-1
            ztrit(i,i+1)=cg_beta(i+1,outloop)
          END DO
          DO i=2,mloop
            ztrit(i,i-1)=cg_beta(i,outloop)
          END DO
          CALL sqlq (mloop, ztrit, tau, zwork1)
!
!   Isolate L=zLT and its tranpose.
!
          zlt=0.0_r8
          zltt=0.0_r8
          DO i=1,mloop
            DO j=1,i
              zlt(i,j)=ztrit(i,j)
            END DO
          END DO
          DO j=1,mloop
            DO i=1,mloop
              zltt(i,j)=zlt(j,i)
            END DO
          END DO
!
!  Compute zck.
!
          zgk=sqrt(zlt(mloop,mloop)*zlt(mloop,mloop)+                   &
     &        cg_beta(mloop+1,outloop)*cg_beta(mloop+1,outloop))
          zck=zlt(mloop,mloop)/zgk
!
!  Now form inv(L)*zt - NOTE: we use L not L', so we use the
!  adjoint of the original solver.
!
          DO j=1,mloop
            DO i=j-1,1,-1
              zt(j)=zt(j)-zt(i)*zltt(i,j)
            END DO
            zt(j)=zt(j)/zltt(j,j)
          END DO
!
!   Now form zt=D*zt.
!
          zt(mloop)=zck*zt(mloop)
!
!   Now form zt=Q'*zt.
!
          DO i=mloop,1,-1
            DO j=1,mloop
              zeref(j)=0.0_r8
            END DO
            zeref(i)=1.0_r8
            DO j=i+1,mloop
              zeref(j)=ztrit(i,j)
            END DO
            zsum=0.0_r8
            DO j=1,mloop
              zsum=zsum+zt(j)*zeref(j)
            END DO
            DO j=1,mloop
              zt(j)=zt(j)-tau(i)*zsum*zeref(j)
            END DO
          END DO
!
!   Copy the solution zt into zu.
!
          DO i=1,mloop
            zu(i)=zt(i)
          END DO
#  else
!
!  Now multiply the result by the inverse tridiagonal matrix.
!
          zbet=cg_delta(1,outloop)
          zu(1)=zt(1)/zbet
          DO ivec=2,mloop
            zgam(ivec)=cg_beta(ivec,outloop)/zbet
            zbet=cg_delta(ivec,outloop)-cg_beta(ivec,outloop)*zgam(ivec)
            zu(ivec)=(zt(ivec)-cg_beta(ivec,outloop)*                   &
     &                zu(ivec-1))/zbet
          END DO
 
          DO ivec=mloop-1,1,-1
            zu(ivec)=zu(ivec)-zgam(ivec+1)*zu(ivec+1)
          END DO
#  endif
!
!  Finally multiply by the matrix of Lanczos vactors.
#  ifndef RPCG
!  Note that the factor of 1/SQRT(ObsErr) is added to covert to
!  x-space.
#  endif
!
          DO iobs=1,ndatum(ng)
            DO innloop=1,mloop
              IF (obserr(iobs).NE.0.0_r8) THEN
#  ifdef RPCG
                ad_obsval(iobs,mloop)=ad_obsval(iobs,mloop)+            &
     &                                obsscale(iobs)*                   &
     &                                tlmodval_s(iobs,innloop,outloop)* &
     &                                zu(innloop)*zfact(innloop)/       &
     &                                obserr(iobs)
#  else
                ad_obsval(iobs,mloop)=ad_obsval(iobs,mloop)+            &
     &                                obsscale(iobs)*                   &
     &                                zcglwk(iobs,innloop,outloop)*     &
     &                                zu(innloop)/                      &
     &                                sqrt(obserr(iobs))
#  endif
              END IF
            END DO
          END DO
 
        END DO
 
      END IF master_thread
 
# else
!
!-----------------------------------------------------------------------
!  Compute observation impact on the weak constraint 4DVAR data
!  assimilation systems.
!-----------------------------------------------------------------------
!
      master_thread : IF (master) THEN
 
        DO iobs=1,ndatum(ng)
          ad_obsval(iobs)=0.0_r8
        END DO
        DO innloop=1,ninnloop
          zt(innloop)=0.0_r8
        END DO
#  ifdef RPCG
        DO innloop=1,ninnloop
          IF (innloop.eq.1) THEN
            zfact(innloop)=1.0_r8/cg_qg(1,outloop)
          ELSE
            zfact(innloop)=1.0_r8/cg_beta(innloop,outloop)
          ENDIF
        END DO
#  endif
!
!  Multiply TLmodVal by the tranpose matrix of Lanczos vectors.
!  Note that the factor of 1/SQRT(ObsErr) is added to convert to
!  x-space.
!
        DO innloop=1,ninnloop
          DO iobs=1,ndatum(ng)
            IF (obserr(iobs).NE.0.0_r8) THEN
#  ifdef RPCG
              zt(innloop)=zt(innloop)+obsscale(iobs)*                   &
     &                    zcglwk(iobs,innloop,outloop)*tlmodval(iobs)
#  else
              zt(innloop)=zt(innloop)+obsscale(iobs)*                   &
     &                    zcglwk(iobs,innloop,outloop)*tlmodval(iobs)/  &
     &                    sqrt(obserr(iobs))
#  endif
            END IF
          END DO
        END DO
 
#  ifdef MINRES
!
!  Use the minimum residual method as described by Paige and Saunders
!  ("Sparse Indefinite Systems of Linear Equations", 1975, SIAM Journal
!  on Numerical Analysis, 617-619). Specifically we refer to equations
!  6.10 and 6.11 of this paper.
!
!  Perform a LQ factorization of the tridiagonal matrix.
!
        ztrit=0.0_r8
        DO i=1,ninnloop
          ztrit(i,i)=cg_delta(i,outloop)
        END DO
        DO i=1,ninnloop-1
          ztrit(i,i+1)=cg_beta(i+1,outloop)
        END DO
        DO i=2,ninnloop
          ztrit(i,i-1)=cg_beta(i,outloop)
        END DO
        CALL sqlq (ninnloop, ztrit, tau, zwork1)
!
!   Isolate L=zLT and its tranpose.
!
        zlt=0.0_r8
        zltt=0.0_r8
        DO i=1,ninnloop
          DO j=1,i
            zlt(i,j)=ztrit(i,j)
          END DO
        END DO
        DO j=1,ninnloop
          DO i=1,ninnloop
            zltt(i,j)=zlt(j,i)
          END DO
        END DO
!
!  Compute zck.
!
        zgk=sqrt(zlt(ninnloop,ninnloop)*zlt(ninnloop,ninnloop)+         &
     &      cg_beta(ninnloop+1,outloop)*cg_beta(ninnloop+1,outloop))
        zck=zlt(ninnloop,ninnloop)/zgk
!
!  Now form inv(L)*zt - NOTE: we use L not L', so we use the
!  adjoint of the original solver.
!
        DO j=1,ninnloop
          DO i=j-1,1,-1
            zt(j)=zt(j)-zt(i)*zltt(i,j)
          END DO
          zt(j)=zt(j)/zltt(j,j)
        END DO
!
!   Now form zt=D*zt.
!
        zt(ninnloop)=zck*zt(ninnloop)
!
!   Now form zt=Q'*zt.
!
        DO i=ninnloop,1,-1
          DO j=1,ninnloop
            zeref(j)=0.0_r8
          END DO
          zeref(i)=1.0_r8
          DO j=i+1,ninnloop
            zeref(j)=ztrit(i,j)
          END DO
          zsum=0.0_r8
          DO j=1,ninnloop
            zsum=zsum+zt(j)*zeref(j)
          END DO
          DO j=1,ninnloop
            zt(j)=zt(j)-tau(i)*zsum*zeref(j)
          END DO
        END DO
!
!   Copy the solution zt into zu.
!
        DO i=1,ninnloop
          zu(i)=zt(i)
        END DO
#  else
!
!  Now multiply the result by the inverse tridiagonal matrix.
!
        zbet=cg_delta(1,outloop)
        zu(1)=zt(1)/zbet
        DO ivec=2,ninnloop
          zgam(ivec)=cg_beta(ivec,outloop)/zbet
          zbet=cg_delta(ivec,outloop)-cg_beta(ivec,outloop)*zgam(ivec)
          zu(ivec)=(zt(ivec)-cg_beta(ivec,outloop)*                     &
     &              zu(ivec-1))/zbet
        END DO
 
        DO ivec=ninnloop-1,1,-1
          zu(ivec)=zu(ivec)-zgam(ivec+1)*zu(ivec+1)
        END DO
#  endif
!
!  Finally multiply by the matrix of Lanczos vactors.
#  ifndef RPCG
!  Note that the factor of 1/SQRT(ObsErr) is added to covert to
!  x-space.
#  endif
!
        DO iobs=1,ndatum(ng)
          DO innloop=1,ninnloop
            IF (obserr(iobs).NE.0.0_r8) THEN
#  ifdef RPCG
              ad_obsval(iobs)=ad_obsval(iobs)+                          &
     &                        obsscale(iobs)*                           &
     &                        tlmodval_s(iobs,innloop,outloop)*         &
     &                        zu(innloop)*zfact(innloop)/               &
     &                        obserr(iobs)
#  else
              ad_obsval(iobs)=ad_obsval(iobs)+                          &
     &                        obsscale(iobs)*                           &
     &                        zcglwk(iobs,innloop,outloop)*             &
     &                        zu(innloop)/                              &
     &                        sqrt(obserr(iobs))
#  endif
            END IF
          END DO
        END DO
 
      END IF master_thread
 
# endif
# ifdef DISTRIBUTE
!
!-----------------------------------------------------------------------
!  Broadcast new solution to other nodes.
!-----------------------------------------------------------------------
!
      CALL mp_bcastf (ng, model, ad_obsval)
# endif
!
      RETURN
      END SUBROUTINE rep_matrix
#else
      SUBROUTINE rep_matrix
      RETURN
      END SUBROUTINE rep_matrix
#endif