From 5f204a1822593a6525915d61b746f1612dfc60a4 Mon Sep 17 00:00:00 2001 From: "lin.gan" Date: Mon, 1 Mar 2021 19:27:23 +0000 Subject: [PATCH 1/4] Part of Doxygen updates for sorc/grid_tools.fd #366 Doxygen update for regional_esg_grid.fd/psym2.f90 --- .../regional_esg_grid.fd/psym2.f90 | 227 +++++++++--------- 1 file changed, 110 insertions(+), 117 deletions(-) diff --git a/sorc/grid_tools.fd/regional_esg_grid.fd/psym2.f90 b/sorc/grid_tools.fd/regional_esg_grid.fd/psym2.f90 index 3fd573cb9..7f637e1f6 100644 --- a/sorc/grid_tools.fd/regional_esg_grid.fd/psym2.f90 +++ b/sorc/grid_tools.fd/regional_esg_grid.fd/psym2.f90 @@ -1,5 +1,6 @@ !> @file -!! @author R. J. Purser @date September 2018 +!! @author R. J. Purser +!! @date September 2018 !! !! A suite of routines to perform the eigen-decomposition of symmetric 2*2 !! matrices and to deliver basic analytic functions, and the derivatives @@ -35,18 +36,17 @@ module psym2 contains -!============================================================================= +!> Get the orthogonal eigenvectors, vv, and diagonal matrix of eigenvalues, oo, +!! of the symmetric 2*2 matrix, em. +!! @author R. J. Purser +!! @param em symmetric 2*2 matrix +!! @param vv orthogonal eigenvectors +!! @param oo diagonal matrix of eigenvalues subroutine eigensym2(em,vv,oo)! [eigensym2] -!============================================================================= -! Get the orthogonal eigenvectors, vv, and diagonal matrix of eigenvalues, oo, -! of the symmetric 2*2 matrix, em. -!============================================================================= implicit none real(dp),dimension(2,2),intent(in ):: em real(dp),dimension(2,2),intent(out):: vv,oo -!----------------------------------------------------------------------------- real(dp):: a,b,c,d,e,f,g,h -!============================================================================= a=em(1,1); b=em(1,2); c=em(2,2) d=a*c-b*b! <- det(em) e=(a+c)*o2; f=(a-c)*o2 @@ -60,26 +60,29 @@ subroutine eigensym2(em,vv,oo)! [eigensym2] oo=matmul(transpose(vv),matmul(em,vv)) oo(1,2)=u0; oo(2,1)=u0 end subroutine eigensym2 -!============================================================================= + +!> For a symmetric 2*2 matrix, em, return the normalized eigenvectors, vv, and +!! the diagonal matrix of eigenvalues, oo. If the two eigenvalues are equal, +!! proceed no further and raise the logical failure flag, ff, to .true.; +!! otherwise, return with vvd=d(vv)/d(em) and ood=d(oo)/d(em) and ff=.false., +!! and maintain the symmetries between the last two of the indices of +!! these derivatives. +!! @author R. J. Purser +!! @param em symmetric 2*2 matrix +!! @param vv normalized eigenvectors +!! @param oo diagonal matrix of eigenvalues +!! @param ff logical failure flag +!! @param vvd vvd=d(vv)/d(em) +!! @param ood ood=d(oo)/d(em) subroutine eigensym2d(em,vv,oo,vvd,ood,ff)! [eigensym2] -!============================================================================= -! For a symmetric 2*2 matrix, em, return the normalized eigenvectors, vv, and -! the diagonal matrix of eigenvalues, oo. If the two eigenvalues are equal, -! proceed no further and raise the logical failure flag, ff, to .true.; -! otherwise, return with vvd=d(vv)/d(em) and ood=d(oo)/d(em) and ff=.false., -! and maintain the symmetries between the last two of the indices of -! these derivatives. -!============================================================================= implicit none real(dp),dimension(2,2), intent(in ):: em real(dp),dimension(2,2), intent(out):: vv,oo real(dp),dimension(2,2,2,2),intent(out):: vvd,ood logical, intent(out):: ff -!----------------------------------------------------------------------------- real(dp),dimension(2,2):: emd,tt,vvr real(dp) :: oodif,dtheta integer(spi) :: i,j -!============================================================================= call eigensym2(em,vv,oo); vvr(1,:)=-vv(2,:); vvr(2,:)=vv(1,:) oodif=oo(1,1)-oo(2,2); ff=oodif==u0; if(ff)return ood=0 @@ -101,37 +104,35 @@ subroutine eigensym2d(em,vv,oo,vvd,ood,ff)! [eigensym2] enddo end subroutine eigensym2d -!============================================================================= +!> Get the inverse of a 2*2 matrix (need not be symmetric in this case). +!! @author R. J. Purser +!! @param em 2*2 matrix +!! @param z inverse of a 2*2 matrix subroutine invsym2(em,z)! [invsym2] -!============================================================================= -! Get the inverse of a 2*2 matrix (need not be symmetric in this case). -!============================================================================= implicit none real(dp),dimension(2,2),intent(in ):: em real(dp),dimension(2,2),intent(out):: z -!----------------------------------------------------------------------------- real(dp):: detem -!============================================================================= z(1,1)=em(2,2); z(2,1)=-em(2,1); z(1,2)=-em(1,2); z(2,2)=em(1,1) detem=em(1,1)*em(2,2)-em(2,1)*em(1,2) z=z/detem end subroutine invsym2 -!============================================================================= + +!> Get the inverse, z,of a 2*2 symmetric matrix, em, and its derivative, zd, +!! with respect to symmetric variations of its components. I.e., for a +!! symmetric infinitesimal change, delta_em, in em, the resulting +!1 infinitesimal change in z would be: +!! delta_z(i,j) = matmul(zd(i,j,:,:),delta_em) +!! @author R. J. Purser +!! @param z inverse of a 2*2 symmetric matrix +!! @param em 2*2 symmetric matrix +!! @param zd derivative of the 2*2 symmetric matrix subroutine invsym2d(em,z,zd)! [invsym2] -!============================================================================= -! Get the inverse, z,of a 2*2 symmetric matrix, em, and its derivative, zd, -! with respect to symmetric variations of its components. I.e., for a -! symmetric infinitesimal change, delta_em, in em, the resulting -! infinitesimal change in z would be: -! delta_z(i,j) = matmul(zd(i,j,:,:),delta_em) -!============================================================================= implicit none real(dp),dimension(2,2) ,intent(in ):: em real(dp),dimension(2,2) ,intent(out):: z real(dp),dimension(2,2,2,2),intent(out):: zd -!----------------------------------------------------------------------------- integer(spi):: k,l -!============================================================================= call invsym2(em,z) call id2222(zd) do l=1,2; do k=1,2 @@ -139,18 +140,16 @@ subroutine invsym2d(em,z,zd)! [invsym2] enddo; enddo end subroutine invsym2d -!============================================================================= +!> Get the sqrt of a symmetric positive-definite 2*2 matrix +!! @author R. J. Purser +!! @param em 2*2 symmetric matrix +!! @param z sqrt of a symmetric positive-definite 2*2 matrix subroutine sqrtsym2(em,z)! [sqrtsym2] -!============================================================================= -! Get the sqrt of a symmetric positive-definite 2*2 matrix -!============================================================================= implicit none real(dp),dimension(2,2),intent(in ):: em real(dp),dimension(2,2),intent(out):: z -!----------------------------------------------------------------------------- real(dp),dimension(2,2):: vv,oo integer(spi) :: i -!============================================================================= call eigensym2(em,vv,oo) do i=1,2 if(oo(i,i)<0)stop 'In sqrtsym2; matrix em is not non-negative' @@ -158,52 +157,52 @@ subroutine sqrtsym2(em,z)! [sqrtsym2] z=matmul(vv,matmul(oo,transpose(vv))) end subroutine sqrtsym2 -!============================================================================= +!> General routine to evaluate z=sqrt(x), and the symmetric +!! derivative, zd = dz/dx, where x is a symmetric 2*2 positive-definite +!! matrix. If the eigenvalues are very close together, extract their +!! geometric mean for "preconditioning" a scaled version, px, of x, whose +!! sqrt, and hence its derivative, can be easily obtained by the series +!! expansion method. Otherwise, use the eigen-method (which entails dividing +!! by the difference in the eignevalues to get zd, and which therefore +!! fails when the eigenvalues become too similar). +!! @author R. J. Purser +!! @param z sqrt(x) result +!! @param zd symmetric derivative +!! @param x symmetric 2*2 positive-definite matrix subroutine sqrtsym2d(x,z,zd)! [sqrtsym2] -!============================================================================= -! General routine to evaluate z=sqrt(x), and the symmetric -! derivative, zd = dz/dx, where x is a symmetric 2*2 positive-definite -! matrix. If the eigenvalues are very close together, extract their -! geometric mean for "preconditioning" a scaled version, px, of x, whose -! sqrt, and hence its derivative, can be easily obtained by the series -! expansion method. Otherwise, use the eigen-method (which entails dividing -! by the difference in the eignevalues to get zd, and which therefore -! fails when the eigenvalues become too similar). -!============================================================================= implicit none real(dp),dimension(2,2), intent(in ):: x real(dp),dimension(2,2), intent(out):: z real(dp),dimension(2,2,2,2),intent(out):: zd -!----------------------------------------------------------------------------- real(dp),dimension(2,2):: px real(dp) :: rdetx,lrdetx,htrpxs,q -!============================================================================= rdetx=sqrt(x(1,1)*x(2,2)-x(1,2)*x(2,1)) ! <- sqrt(determinant of x) lrdetx=sqrt(rdetx) -px=x/rdetx ! <- preconditioned x (has unit determinant) -htrpxs= ((px(1,1)+px(2,2))/2)**2 ! <- {half-trace-px}-squared +px=x/rdetx ! - preconditioned x (has unit determinant) +htrpxs= ((px(1,1)+px(2,2))/2)**2 ! - {half-trace-px}-squared q=htrpxs-u1 -if(q<.05_dp)then ! <- Taylor expansion method +if(q<.05_dp)then ! - Taylor expansion method call sqrtsym2d_t(px,z,zd) z=z*lrdetx; zd=zd/lrdetx else - call sqrtsym2d_e(x,z,zd) ! <- Eigen-method + call sqrtsym2d_e(x,z,zd) ! - Eigen-method endif end subroutine sqrtsym2d -!============================================================================= +!> Eigen-method +!! @author R. J. Purser +!! @param z sqrt(x) result +!! @param zd symmetric derivative +!! @param x symmetric 2*2 positive-definite matrix subroutine sqrtsym2d_e(x,z,zd)! [sqrtsym2d_e] -!============================================================================= implicit none real(dp),dimension(2,2), intent(in ):: x real(dp),dimension(2,2), intent(out):: z real(dp),dimension(2,2,2,2),intent(out):: zd -!----------------------------------------------------------------------------- real(dp),dimension(2,2,2,2):: vvd,ood real(dp),dimension(2,2) :: vv,oo,oori,tt integer(spi) :: i,j logical :: ff -!============================================================================= call eigensym2(x,vv,oo,vvd,ood,ff) z=u0; z(1,1)=sqrt(oo(1,1)); z(2,2)=sqrt(oo(2,2)) z=matmul(matmul(vv,z),transpose(vv)) @@ -217,26 +216,25 @@ subroutine sqrtsym2d_e(x,z,zd)! [sqrtsym2d_e] enddo end subroutine sqrtsym2d_e -!============================================================================= +!> Use the Taylor-series method (eigenvalues both fairly close to unity). +!! For a 2*2 positive definite symmetric matrix x, try to get both the z=sqrt(x) +!! and dz/dx using the binomial-expansion method applied to the intermediate +!! matrix, r = (x-1). ie z=sqrt(x) = (1+r)^{1/2} = I + (1/2)*r -(1/8)*r^2 ... +!! + [(-)^n *(2n)!/{(n+1)! * n! *2^{2*n-1}} ]*r^{n+1} +!! @author R. J. Purser +!! @param z sqrt(x) result +!! @param zd symmetric derivative +!! @param x symmetric 2*2 positive-definite matrix subroutine sqrtsym2d_t(x,z,zd)! [sqrtsym2d_t] -!============================================================================= -! Use the Taylor-series method (eigenvalues both fairly close to unity). -! For a 2*2 positive definite symmetric matrix x, try to get both the z=sqrt(x) -! and dz/dx using the binomial-expansion method applied to the intermediate -! matrix, r = (x-1). ie z=sqrt(x) = (1+r)^{1/2} = I + (1/2)*r -(1/8)*r^2 ... -! + [(-)^n *(2n)!/{(n+1)! * n! *2^{2*n-1}} ]*r^{n+1} -!============================================================================= implicit none real(dp),dimension(2,2), intent(in ):: x real(dp),dimension(2,2), intent(out):: z real(dp),dimension(2,2,2,2),intent(out):: zd -!----------------------------------------------------------------------------- integer(spi),parameter :: nit=300 ! number of iterative increments allowed real(dp),parameter :: crit=1.e-17 real(dp),dimension(2,2) :: r,rp,rd,rpd real(dp) :: c integer(spi) :: i,j,n -!============================================================================= r=x; r(1,1)=x(1,1)-1; r(2,2)=x(2,2)-1 z=u0; z(1,1)=u1; z(2,2)=u1 rp=r @@ -263,33 +261,32 @@ subroutine sqrtsym2d_t(x,z,zd)! [sqrtsym2d_t] enddo; enddo end subroutine sqrtsym2d_t -!============================================================================= +!> Get the exp of a symmetric 2*2 matrix +!! @author R. J. Purser +!! @param em symmetric 2*2 matrix +!! @param expem exp of a symmetric 2*2 matrix subroutine expsym2(em,expem)! [expsym2] -!============================================================================= -! Get the exp of a symmetric 2*2 matrix -!============================================================================= implicit none real(dp),dimension(2,2),intent(in ):: em real(dp),dimension(2,2),intent(out):: expem -!----------------------------------------------------------------------------- real(dp),dimension(2,2):: vv,oo integer(spi) :: i -!============================================================================= call eigensym2(em,vv,oo) do i=1,2; oo(i,i)=exp(oo(i,i)); enddo expem=matmul(vv,matmul(oo,transpose(vv))) end subroutine expsym2 -!============================================================================= + +!> Sub Process for process symmetric 2*2 matrix +!! @author R. J. Purser +!! @param zd symmetric derivative +!! @param x symmetric 2*2 positive-definite matrix subroutine expsym2d(x,z,zd)! [expsym2] -!============================================================================= implicit none real(dp),dimension(2,2), intent(in ):: x real(dp),dimension(2,2), intent(out):: z real(dp),dimension(2,2,2,2),intent(out):: zd -!----------------------------------------------------------------------------- real(dp),dimension(2,2):: px real(dp) :: trxh,detpx -!============================================================================= trxh=(x(1,1)+x(2,2))*o2 px=x;px(1,1)=x(1,1)-trxh;px(2,2)=x(2,2)-trxh detpx=abs(px(1,1)*px(2,2)-px(1,2)*px(2,1)) @@ -299,19 +296,19 @@ subroutine expsym2d(x,z,zd)! [expsym2] z=z*exp(trxh) end subroutine expsym2d -!============================================================================= +!> Sub Process for process symmetric 2*2 matrix +!! @author R. J. Purser +!! @param zd symmetric derivative +!! @param x symmetric 2*2 positive-definite matrix subroutine expsym2d_e(x,z,zd)! [expsym2d_e] -!============================================================================= implicit none real(dp),dimension(2,2), intent(in ):: x real(dp),dimension(2,2), intent(out):: z real(dp),dimension(2,2,2,2),intent(out):: zd -!----------------------------------------------------------------------------- real(dp),dimension(2,2,2,2):: vvd,ood real(dp),dimension(2,2) :: vv,oo,ooe,tt integer(spi) :: i,j logical :: ff -!============================================================================= call eigensym2(x,vv,oo,vvd,ood,ff) z=u0; z(1,1)=exp(oo(1,1)); z(2,2)=exp(oo(2,2)) z=matmul(matmul(vv,z),transpose(vv)) @@ -325,24 +322,23 @@ subroutine expsym2d_e(x,z,zd)! [expsym2d_e] enddo end subroutine expsym2d_e -!============================================================================= +!> Use the Taylor-series method (eigenvalues both fairly close to zero). +!! For a 2*2 symmetric matrix x, try to get both the z=exp(x) +!! and dz/dx using the Taylor series expansion method. +!! @author R. J. Purser +!! @param z Taylor series expansion method exp(x) +!! @param zd symmetric derivative +!! @param x symmetric 2*2 positive-definite matrix subroutine expsym2d_t(x,z,zd)! [expsym2d_t] -!============================================================================= -! Use the Taylor-series method (eigenvalues both fairly close to zero). -! For a 2*2 symmetric matrix x, try to get both the z=exp(x) -! and dz/dx using the Taylor series expansion method. -!============================================================================= implicit none real(dp),dimension(2,2), intent(in ):: x real(dp),dimension(2,2), intent(out):: z real(dp),dimension(2,2,2,2),intent(out):: zd -!----------------------------------------------------------------------------- integer(spi),parameter :: nit=100 ! number of iterative increments allowed real(dp),parameter :: crit=1.e-17_dp real(dp),dimension(2,2) :: xp,xd,xpd real(dp) :: c integer(spi) :: i,j,n -!============================================================================= z=0; z(1,1)=u1; z(2,2)=u1 xp=x c=u1 @@ -368,18 +364,16 @@ subroutine expsym2d_t(x,z,zd)! [expsym2d_t] enddo; enddo end subroutine expsym2d_t -!============================================================================= +!> Get the log of a symmetric positive-definite 2*2 matrix +!! @author R. J. Purser +!! @param em symmetric 2*2 matrix +!! @param logem log of a symmetric positive-definite 2*2 matrix subroutine logsym2(em,logem)! [logsym2] -!============================================================================= -! Get the log of a symmetric positive-definite 2*2 matrix -!============================================================================= implicit none real(dp),dimension(2,2),intent(in ):: em real(dp),dimension(2,2),intent(out):: logem -!----------------------------------------------------------------------------- real(dp),dimension(2,2):: vv,oo integer(spi) :: i -!============================================================================= call eigensym2(em,vv,oo) do i=1,2 if(oo(i,i)<=u0)stop 'In logsym2; matrix em is not positive definite' @@ -387,23 +381,23 @@ subroutine logsym2(em,logem)! [logsym2] enddo logem=matmul(vv,matmul(oo,transpose(vv))) end subroutine logsym2 -!============================================================================= + +!> General routine to evaluate the logarithm, z=log(x), and the symmetric +!! derivative, zd = dz/dx, where x is a symmetric 2*2 positive-definite +!! matrix. +!! @author R. J. Purser +!! @param zd the symmetric derivative +!! @param x a symmetric 2*2 positive-definite matrix +!! @param z evaluate the logarithm log(x) subroutine logsym2d(x,z,zd)! [logsym2] -!============================================================================= -! General routine to evaluate the logarithm, z=log(x), and the symmetric -! derivative, zd = dz/dx, where x is a symmetric 2*2 positive-definite -! matrix. -!============================================================================= use pfun, only: sinhox implicit none real(dp),dimension(2,2), intent(in ):: x real(dp),dimension(2,2), intent(out):: z real(dp),dimension(2,2,2,2),intent(out):: zd -!----------------------------------------------------------------------------- real(dp),dimension(2,2):: vv,oo,d11,d12,d22,pqr real(dp) :: c,s,cc,cs,ss,c2h,p,q,r,lp,lq,L integer(spi) :: i -!============================================================================= call eigensym2(x,vv,oo) if(oo(1,1)<=u0 .or. oo(2,2)<=u0)stop 'In logsym2; matrix x is not positive definite' c=vv(1,1); s=vv(1,2); cc=c*c; cs=c*s; ss=s*s; c2h=(cc-ss)*o2 @@ -420,9 +414,9 @@ subroutine logsym2d(x,z,zd)! [logsym2] zd(:,:,2,1)=zd(:,:,1,2) end subroutine logsym2d -!============================================================================= +!> General routine for Effective identity +!! @author R. J. Purser subroutine id2222(em)! [id2222] -!============================================================================= implicit none real(dp),dimension(2,2,2,2),intent(out):: em real(dp),dimension(2,2,2,2) :: id @@ -430,20 +424,19 @@ subroutine id2222(em)! [id2222] em=id end subroutine id2222 -!=========================================================================== +!> Return the cholesky lower triangular factor, C, of the 2X2 symmetric +!! matrix, S, or raise the failure flag, FF, if S is not positive-definite. +!! @author R. J. Purser +!! @param c cholesky lower triangular factor +!! @param s 2X2 symmetric matrix +!! @param ff raise the failure flag subroutine chol2(s,c,ff)! [chol2] -!=========================================================================== -! Return the cholesky lower triangular factor, C, of the 2X2 symmetric -! matrix, S, or raise the failure flag, FF, if S is not positive-definite. -!=========================================================================== use pietc, only: u0 implicit none real(dp),dimension(2,2),intent(in ):: s real(dp),dimension(2,2),intent(out):: c logical ,intent(out):: ff -!--------------------------------------------------------------------------- real(dp):: r -!=========================================================================== ff=s(1,1)<=u0; if(ff)return c(1,1)=sqrt(s(1,1)) c(1,2)=u0 From 2f5d83e4b58db8a71f9c999e322c7dce2589cf3b Mon Sep 17 00:00:00 2001 From: "lin.gan" Date: Mon, 1 Mar 2021 20:24:59 +0000 Subject: [PATCH 2/4] Part of Doxygen updates for sorc/grid_tools.fd #366 modify grid_tools.fd/regional_esg_grid.fd/psym2.f90 with doxygen improvement --- sorc/grid_tools.fd/regional_esg_grid.fd/psym2.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sorc/grid_tools.fd/regional_esg_grid.fd/psym2.f90 b/sorc/grid_tools.fd/regional_esg_grid.fd/psym2.f90 index 7f637e1f6..68d74d59a 100644 --- a/sorc/grid_tools.fd/regional_esg_grid.fd/psym2.f90 +++ b/sorc/grid_tools.fd/regional_esg_grid.fd/psym2.f90 @@ -121,7 +121,7 @@ end subroutine invsym2 !> Get the inverse, z,of a 2*2 symmetric matrix, em, and its derivative, zd, !! with respect to symmetric variations of its components. I.e., for a !! symmetric infinitesimal change, delta_em, in em, the resulting -!1 infinitesimal change in z would be: +!! infinitesimal change in z would be: !! delta_z(i,j) = matmul(zd(i,j,:,:),delta_em) !! @author R. J. Purser !! @param z inverse of a 2*2 symmetric matrix From 28071dafd282e9b1f284627ae989f85e7347e56e Mon Sep 17 00:00:00 2001 From: "Lin.Gan" Date: Tue, 2 Mar 2021 14:41:04 +0000 Subject: [PATCH 3/4] part of Doxygen updates for sorc/grid_tools.fd #366 sorc/grid_tools.fd/regional_esg_grid.fd/psym2.f90 --- .../regional_esg_grid.fd/psym2.f90 | 40 +++++++++---------- 1 file changed, 19 insertions(+), 21 deletions(-) diff --git a/sorc/grid_tools.fd/regional_esg_grid.fd/psym2.f90 b/sorc/grid_tools.fd/regional_esg_grid.fd/psym2.f90 index 68d74d59a..50cac9c04 100644 --- a/sorc/grid_tools.fd/regional_esg_grid.fd/psym2.f90 +++ b/sorc/grid_tools.fd/regional_esg_grid.fd/psym2.f90 @@ -1,18 +1,16 @@ !> @file -!! @author R. J. Purser -!! @date September 2018 -!! !! A suite of routines to perform the eigen-decomposition of symmetric 2*2 !! matrices and to deliver basic analytic functions, and the derivatives !! of these functions, of such matrices. !! In addition, we include a simple cholesky routine -!! !! DIRECT DEPENDENCIES !! Library: pfun !! Module: pkind, pietc, pfun -!! +!! @author R. J. Purser @date September 2018 + +!> Module of 11 subroutines +!! @author R. J. Purser module psym2 -!============================================================================= use pkind, only: spi,dp use pietc, only: u0,u1,o2 implicit none @@ -38,10 +36,10 @@ module psym2 !> Get the orthogonal eigenvectors, vv, and diagonal matrix of eigenvalues, oo, !! of the symmetric 2*2 matrix, em. -!! @author R. J. Purser !! @param em symmetric 2*2 matrix !! @param vv orthogonal eigenvectors !! @param oo diagonal matrix of eigenvalues +!! @author R. J. Purser subroutine eigensym2(em,vv,oo)! [eigensym2] implicit none real(dp),dimension(2,2),intent(in ):: em @@ -67,13 +65,13 @@ end subroutine eigensym2 !! otherwise, return with vvd=d(vv)/d(em) and ood=d(oo)/d(em) and ff=.false., !! and maintain the symmetries between the last two of the indices of !! these derivatives. -!! @author R. J. Purser !! @param em symmetric 2*2 matrix !! @param vv normalized eigenvectors !! @param oo diagonal matrix of eigenvalues !! @param ff logical failure flag !! @param vvd vvd=d(vv)/d(em) !! @param ood ood=d(oo)/d(em) +!! @author R. J. Purser subroutine eigensym2d(em,vv,oo,vvd,ood,ff)! [eigensym2] implicit none real(dp),dimension(2,2), intent(in ):: em @@ -105,9 +103,9 @@ subroutine eigensym2d(em,vv,oo,vvd,ood,ff)! [eigensym2] end subroutine eigensym2d !> Get the inverse of a 2*2 matrix (need not be symmetric in this case). -!! @author R. J. Purser !! @param em 2*2 matrix !! @param z inverse of a 2*2 matrix +!! @author R. J. Purser subroutine invsym2(em,z)! [invsym2] implicit none real(dp),dimension(2,2),intent(in ):: em @@ -123,10 +121,10 @@ end subroutine invsym2 !! symmetric infinitesimal change, delta_em, in em, the resulting !! infinitesimal change in z would be: !! delta_z(i,j) = matmul(zd(i,j,:,:),delta_em) -!! @author R. J. Purser !! @param z inverse of a 2*2 symmetric matrix !! @param em 2*2 symmetric matrix !! @param zd derivative of the 2*2 symmetric matrix +!! @author R. J. Purser subroutine invsym2d(em,z,zd)! [invsym2] implicit none real(dp),dimension(2,2) ,intent(in ):: em @@ -141,9 +139,9 @@ subroutine invsym2d(em,z,zd)! [invsym2] end subroutine invsym2d !> Get the sqrt of a symmetric positive-definite 2*2 matrix -!! @author R. J. Purser !! @param em 2*2 symmetric matrix !! @param z sqrt of a symmetric positive-definite 2*2 matrix +!! @author R. J. Purser subroutine sqrtsym2(em,z)! [sqrtsym2] implicit none real(dp),dimension(2,2),intent(in ):: em @@ -165,10 +163,10 @@ end subroutine sqrtsym2 !! expansion method. Otherwise, use the eigen-method (which entails dividing !! by the difference in the eignevalues to get zd, and which therefore !! fails when the eigenvalues become too similar). -!! @author R. J. Purser !! @param z sqrt(x) result !! @param zd symmetric derivative !! @param x symmetric 2*2 positive-definite matrix +!! @author R. J. Purser subroutine sqrtsym2d(x,z,zd)! [sqrtsym2] implicit none real(dp),dimension(2,2), intent(in ):: x @@ -190,10 +188,10 @@ subroutine sqrtsym2d(x,z,zd)! [sqrtsym2] end subroutine sqrtsym2d !> Eigen-method -!! @author R. J. Purser !! @param z sqrt(x) result !! @param zd symmetric derivative !! @param x symmetric 2*2 positive-definite matrix +!! @author R. J. Purser subroutine sqrtsym2d_e(x,z,zd)! [sqrtsym2d_e] implicit none real(dp),dimension(2,2), intent(in ):: x @@ -221,10 +219,10 @@ end subroutine sqrtsym2d_e !! and dz/dx using the binomial-expansion method applied to the intermediate !! matrix, r = (x-1). ie z=sqrt(x) = (1+r)^{1/2} = I + (1/2)*r -(1/8)*r^2 ... !! + [(-)^n *(2n)!/{(n+1)! * n! *2^{2*n-1}} ]*r^{n+1} -!! @author R. J. Purser !! @param z sqrt(x) result !! @param zd symmetric derivative !! @param x symmetric 2*2 positive-definite matrix +!! @author R. J. Purser subroutine sqrtsym2d_t(x,z,zd)! [sqrtsym2d_t] implicit none real(dp),dimension(2,2), intent(in ):: x @@ -262,9 +260,9 @@ subroutine sqrtsym2d_t(x,z,zd)! [sqrtsym2d_t] end subroutine sqrtsym2d_t !> Get the exp of a symmetric 2*2 matrix -!! @author R. J. Purser !! @param em symmetric 2*2 matrix !! @param expem exp of a symmetric 2*2 matrix +!! @author R. J. Purser subroutine expsym2(em,expem)! [expsym2] implicit none real(dp),dimension(2,2),intent(in ):: em @@ -277,9 +275,9 @@ subroutine expsym2(em,expem)! [expsym2] end subroutine expsym2 !> Sub Process for process symmetric 2*2 matrix -!! @author R. J. Purser !! @param zd symmetric derivative !! @param x symmetric 2*2 positive-definite matrix +!! @author R. J. Purser subroutine expsym2d(x,z,zd)! [expsym2] implicit none real(dp),dimension(2,2), intent(in ):: x @@ -297,9 +295,9 @@ subroutine expsym2d(x,z,zd)! [expsym2] end subroutine expsym2d !> Sub Process for process symmetric 2*2 matrix -!! @author R. J. Purser !! @param zd symmetric derivative !! @param x symmetric 2*2 positive-definite matrix +!! @author R. J. Purser subroutine expsym2d_e(x,z,zd)! [expsym2d_e] implicit none real(dp),dimension(2,2), intent(in ):: x @@ -325,10 +323,10 @@ end subroutine expsym2d_e !> Use the Taylor-series method (eigenvalues both fairly close to zero). !! For a 2*2 symmetric matrix x, try to get both the z=exp(x) !! and dz/dx using the Taylor series expansion method. -!! @author R. J. Purser !! @param z Taylor series expansion method exp(x) !! @param zd symmetric derivative !! @param x symmetric 2*2 positive-definite matrix +!! @author R. J. Purser subroutine expsym2d_t(x,z,zd)! [expsym2d_t] implicit none real(dp),dimension(2,2), intent(in ):: x @@ -365,9 +363,9 @@ subroutine expsym2d_t(x,z,zd)! [expsym2d_t] end subroutine expsym2d_t !> Get the log of a symmetric positive-definite 2*2 matrix -!! @author R. J. Purser !! @param em symmetric 2*2 matrix !! @param logem log of a symmetric positive-definite 2*2 matrix +!! @author R. J. Purser subroutine logsym2(em,logem)! [logsym2] implicit none real(dp),dimension(2,2),intent(in ):: em @@ -385,10 +383,10 @@ end subroutine logsym2 !> General routine to evaluate the logarithm, z=log(x), and the symmetric !! derivative, zd = dz/dx, where x is a symmetric 2*2 positive-definite !! matrix. -!! @author R. J. Purser !! @param zd the symmetric derivative !! @param x a symmetric 2*2 positive-definite matrix !! @param z evaluate the logarithm log(x) +!! @author R. J. Purser subroutine logsym2d(x,z,zd)! [logsym2] use pfun, only: sinhox implicit none @@ -426,10 +424,10 @@ end subroutine id2222 !> Return the cholesky lower triangular factor, C, of the 2X2 symmetric !! matrix, S, or raise the failure flag, FF, if S is not positive-definite. -!! @author R. J. Purser !! @param c cholesky lower triangular factor !! @param s 2X2 symmetric matrix !! @param ff raise the failure flag +!! @author R. J. Purser subroutine chol2(s,c,ff)! [chol2] use pietc, only: u0 implicit none From df3ea0a36c5bef68e5867ccadc494f489845f203 Mon Sep 17 00:00:00 2001 From: Edward Hartnett Date: Tue, 2 Mar 2021 12:06:52 -0700 Subject: [PATCH 4/4] doxygen fixes --- .../regional_esg_grid.fd/psym2.f90 | 203 ++++++++++-------- 1 file changed, 112 insertions(+), 91 deletions(-) diff --git a/sorc/grid_tools.fd/regional_esg_grid.fd/psym2.f90 b/sorc/grid_tools.fd/regional_esg_grid.fd/psym2.f90 index 50cac9c04..f8ed4293a 100644 --- a/sorc/grid_tools.fd/regional_esg_grid.fd/psym2.f90 +++ b/sorc/grid_tools.fd/regional_esg_grid.fd/psym2.f90 @@ -1,14 +1,12 @@ !> @file -!! A suite of routines to perform the eigen-decomposition of symmetric 2*2 -!! matrices and to deliver basic analytic functions, and the derivatives -!! of these functions, of such matrices. -!! In addition, we include a simple cholesky routine -!! DIRECT DEPENDENCIES -!! Library: pfun -!! Module: pkind, pietc, pfun +!! @brief Matrix routines. !! @author R. J. Purser @date September 2018 -!> Module of 11 subroutines +!> A suite of routines to perform the eigen-decomposition of symmetric 2*2 +!! matrices and to deliver basic analytic functions, and the +!! derivatives of these functions, of such matrices. In addition, we +!! include a simple cholesky routine. +!! !! @author R. J. Purser module psym2 use pkind, only: spi,dp @@ -17,7 +15,7 @@ module psym2 private public:: eigensym2,invsym2,sqrtsym2,expsym2,logsym2,id2222,chol2 -real(dp),dimension(2,2,2,2):: id +real(dp),dimension(2,2,2,2):: id !< ID. data id/u1,u0,u0,u0, u0,o2,o2,u0, u0,o2,o2,u0, u0,u0,u0,u1/! Effective identity interface eigensym2; module procedure eigensym2,eigensym2d; end interface @@ -34,11 +32,12 @@ module psym2 contains -!> Get the orthogonal eigenvectors, vv, and diagonal matrix of eigenvalues, oo, -!! of the symmetric 2*2 matrix, em. -!! @param em symmetric 2*2 matrix -!! @param vv orthogonal eigenvectors -!! @param oo diagonal matrix of eigenvalues +!> Get the orthogonal eigenvectors, vv, and diagonal matrix of +!! eigenvalues, oo, of the symmetric 2*2 matrix, em. +!! +!! @param[in] em symmetric 2*2 matrix +!! @param[out] vv orthogonal eigenvectors +!! @param[out] oo diagonal matrix of eigenvalues !! @author R. J. Purser subroutine eigensym2(em,vv,oo)! [eigensym2] implicit none @@ -59,18 +58,19 @@ subroutine eigensym2(em,vv,oo)! [eigensym2] oo(1,2)=u0; oo(2,1)=u0 end subroutine eigensym2 -!> For a symmetric 2*2 matrix, em, return the normalized eigenvectors, vv, and -!! the diagonal matrix of eigenvalues, oo. If the two eigenvalues are equal, -!! proceed no further and raise the logical failure flag, ff, to .true.; -!! otherwise, return with vvd=d(vv)/d(em) and ood=d(oo)/d(em) and ff=.false., -!! and maintain the symmetries between the last two of the indices of -!! these derivatives. -!! @param em symmetric 2*2 matrix -!! @param vv normalized eigenvectors -!! @param oo diagonal matrix of eigenvalues -!! @param ff logical failure flag -!! @param vvd vvd=d(vv)/d(em) -!! @param ood ood=d(oo)/d(em) +!> For a symmetric 2*2 matrix, em, return the normalized eigenvectors, +!! vv, and the diagonal matrix of eigenvalues, oo. If the two +!! eigenvalues are equal, proceed no further and raise the logical +!! failure flag, ff, to .true.; otherwise, return with vvd=d(vv)/d(em) +!! and ood=d(oo)/d(em) and ff=.false., and maintain the symmetries +!! between the last two of the indices of these derivatives. +!! +!! @param[in] em symmetric 2*2 matrix +!! @param[out] vv normalized eigenvectors +!! @param[out] oo diagonal matrix of eigenvalues +!! @param[out] vvd vvd=d(vv)/d(em) +!! @param[out] ood ood=d(oo)/d(em) +!! @param[out] ff logical failure flag !! @author R. J. Purser subroutine eigensym2d(em,vv,oo,vvd,ood,ff)! [eigensym2] implicit none @@ -102,9 +102,11 @@ subroutine eigensym2d(em,vv,oo,vvd,ood,ff)! [eigensym2] enddo end subroutine eigensym2d -!> Get the inverse of a 2*2 matrix (need not be symmetric in this case). -!! @param em 2*2 matrix -!! @param z inverse of a 2*2 matrix +!> Get the inverse of a 2*2 matrix (need not be symmetric in this +!! case). +!! +!! @param[in] em 2*2 matrix +!! @param[out] z inverse of a 2*2 matrix !! @author R. J. Purser subroutine invsym2(em,z)! [invsym2] implicit none @@ -116,14 +118,15 @@ subroutine invsym2(em,z)! [invsym2] z=z/detem end subroutine invsym2 -!> Get the inverse, z,of a 2*2 symmetric matrix, em, and its derivative, zd, -!! with respect to symmetric variations of its components. I.e., for a -!! symmetric infinitesimal change, delta_em, in em, the resulting -!! infinitesimal change in z would be: -!! delta_z(i,j) = matmul(zd(i,j,:,:),delta_em) -!! @param z inverse of a 2*2 symmetric matrix -!! @param em 2*2 symmetric matrix -!! @param zd derivative of the 2*2 symmetric matrix +!> Get the inverse, z,of a 2*2 symmetric matrix, em, and its +!! derivative, zd, with respect to symmetric variations of its +!! components. I.e., for a symmetric infinitesimal change, delta_em, +!! in em, the resulting infinitesimal change in z would be: +!! 
delta_z(i,j) = matmul(zd(i,j,:,:),delta_em)
+!! +!! @param[in] em 2*2 symmetric matrix +!! @param[out] z inverse of a 2*2 symmetric matrix +!! @param[out] zd derivative of the 2*2 symmetric matrix !! @author R. J. Purser subroutine invsym2d(em,z,zd)! [invsym2] implicit none @@ -138,9 +141,10 @@ subroutine invsym2d(em,z,zd)! [invsym2] enddo; enddo end subroutine invsym2d -!> Get the sqrt of a symmetric positive-definite 2*2 matrix -!! @param em 2*2 symmetric matrix -!! @param z sqrt of a symmetric positive-definite 2*2 matrix +!> Get the sqrt of a symmetric positive-definite 2*2 matrix. +!! +!! @param[in] em 2*2 symmetric matrix +!! @param[out] z sqrt of a symmetric positive-definite 2*2 matrix !! @author R. J. Purser subroutine sqrtsym2(em,z)! [sqrtsym2] implicit none @@ -155,17 +159,19 @@ subroutine sqrtsym2(em,z)! [sqrtsym2] z=matmul(vv,matmul(oo,transpose(vv))) end subroutine sqrtsym2 -!> General routine to evaluate z=sqrt(x), and the symmetric -!! derivative, zd = dz/dx, where x is a symmetric 2*2 positive-definite -!! matrix. If the eigenvalues are very close together, extract their -!! geometric mean for "preconditioning" a scaled version, px, of x, whose -!! sqrt, and hence its derivative, can be easily obtained by the series -!! expansion method. Otherwise, use the eigen-method (which entails dividing -!! by the difference in the eignevalues to get zd, and which therefore -!! fails when the eigenvalues become too similar). -!! @param z sqrt(x) result -!! @param zd symmetric derivative -!! @param x symmetric 2*2 positive-definite matrix +!> General routine to evaluate z=sqrt(x), and the symmetric +!! derivative, zd = dz/dx, where x is a symmetric 2*2 +!! positive-definite matrix. If the eigenvalues are very close +!! together, extract their geometric mean for "preconditioning" a +!! scaled version, px, of x, whose sqrt, and hence its derivative, can +!! be easily obtained by the series expansion method. Otherwise, use +!! the eigen-method (which entails dividing by the difference in the +!! eignevalues to get zd, and which therefore fails when the +!! eigenvalues become too similar). +!! +!! @param[in] x symmetric 2*2 positive-definite matrix +!! @param[out] z sqrt(x) result +!! @param[out] zd symmetric derivative !! @author R. J. Purser subroutine sqrtsym2d(x,z,zd)! [sqrtsym2] implicit none @@ -187,10 +193,11 @@ subroutine sqrtsym2d(x,z,zd)! [sqrtsym2] endif end subroutine sqrtsym2d -!> Eigen-method -!! @param z sqrt(x) result -!! @param zd symmetric derivative -!! @param x symmetric 2*2 positive-definite matrix +!> Eigen-method. +!! +!! @param[in] x symmetric 2*2 positive-definite matrix +!! @param[out] z sqrt(x) result +!! @param[out] zd symmetric derivative !! @author R. J. Purser subroutine sqrtsym2d_e(x,z,zd)! [sqrtsym2d_e] implicit none @@ -215,13 +222,15 @@ subroutine sqrtsym2d_e(x,z,zd)! [sqrtsym2d_e] end subroutine sqrtsym2d_e !> Use the Taylor-series method (eigenvalues both fairly close to unity). -!! For a 2*2 positive definite symmetric matrix x, try to get both the z=sqrt(x) -!! and dz/dx using the binomial-expansion method applied to the intermediate -!! matrix, r = (x-1). ie z=sqrt(x) = (1+r)^{1/2} = I + (1/2)*r -(1/8)*r^2 ... -!! + [(-)^n *(2n)!/{(n+1)! * n! *2^{2*n-1}} ]*r^{n+1} -!! @param z sqrt(x) result -!! @param zd symmetric derivative -!! @param x symmetric 2*2 positive-definite matrix +!! For a 2*2 positive definite symmetric matrix x, try to get both the +!! z=sqrt(x) and dz/dx using the binomial-expansion method applied to +!! the intermediate matrix, +!! 
r = (x-1). ie z=sqrt(x) = (1+r)^{1/2} = I + (1/2)*r -(1/8)*r^2 ...
+!!  + [(-)^n *(2n)!/{(n+1)! * n! *2^{2*n-1}} ]*r^{n+1}
+!! +!! @param[in] x symmetric 2*2 positive-definite matrix +!! @param[out] z sqrt(x) result +!! @param[out] zd symmetric derivative !! @author R. J. Purser subroutine sqrtsym2d_t(x,z,zd)! [sqrtsym2d_t] implicit none @@ -259,9 +268,10 @@ subroutine sqrtsym2d_t(x,z,zd)! [sqrtsym2d_t] enddo; enddo end subroutine sqrtsym2d_t -!> Get the exp of a symmetric 2*2 matrix -!! @param em symmetric 2*2 matrix -!! @param expem exp of a symmetric 2*2 matrix +!> Get the exp of a symmetric 2*2 matrix. +!! +!! @param[in] em symmetric 2*2 matrix +!! @param[out] expem exp of a symmetric 2*2 matrix !! @author R. J. Purser subroutine expsym2(em,expem)! [expsym2] implicit none @@ -274,9 +284,11 @@ subroutine expsym2(em,expem)! [expsym2] expem=matmul(vv,matmul(oo,transpose(vv))) end subroutine expsym2 -!> Sub Process for process symmetric 2*2 matrix -!! @param zd symmetric derivative -!! @param x symmetric 2*2 positive-definite matrix +!> Sub Process for process symmetric 2*2 matrix. +!! +!! @param[in] x symmetric 2*2 positive-definite matrix +!! @param[out] z ??? +!! @param[out] zd symmetric derivative !! @author R. J. Purser subroutine expsym2d(x,z,zd)! [expsym2] implicit none @@ -294,9 +306,11 @@ subroutine expsym2d(x,z,zd)! [expsym2] z=z*exp(trxh) end subroutine expsym2d -!> Sub Process for process symmetric 2*2 matrix -!! @param zd symmetric derivative -!! @param x symmetric 2*2 positive-definite matrix +!> Sub Process for process symmetric 2*2 matrix. +!! +!! @param[in] x symmetric 2*2 positive-definite matrix +!! @param[out] z ??? +!! @param[out] zd symmetric derivative !! @author R. J. Purser subroutine expsym2d_e(x,z,zd)! [expsym2d_e] implicit none @@ -320,12 +334,13 @@ subroutine expsym2d_e(x,z,zd)! [expsym2d_e] enddo end subroutine expsym2d_e -!> Use the Taylor-series method (eigenvalues both fairly close to zero). -!! For a 2*2 symmetric matrix x, try to get both the z=exp(x) +!> Use the Taylor-series method (eigenvalues both fairly close to +!! zero). For a 2*2 symmetric matrix x, try to get both the z=exp(x) !! and dz/dx using the Taylor series expansion method. -!! @param z Taylor series expansion method exp(x) -!! @param zd symmetric derivative -!! @param x symmetric 2*2 positive-definite matrix +!! +!! @param[in] x symmetric 2*2 positive-definite matrix +!! @param[out] z Taylor series expansion method exp(x) +!! @param[out] zd symmetric derivative !! @author R. J. Purser subroutine expsym2d_t(x,z,zd)! [expsym2d_t] implicit none @@ -362,9 +377,10 @@ subroutine expsym2d_t(x,z,zd)! [expsym2d_t] enddo; enddo end subroutine expsym2d_t -!> Get the log of a symmetric positive-definite 2*2 matrix -!! @param em symmetric 2*2 matrix -!! @param logem log of a symmetric positive-definite 2*2 matrix +!> Get the log of a symmetric positive-definite 2*2 matrix. +!! +!! @param[in] em symmetric 2*2 matrix +!! @param[out] logem log of a symmetric positive-definite 2*2 matrix !! @author R. J. Purser subroutine logsym2(em,logem)! [logsym2] implicit none @@ -380,12 +396,13 @@ subroutine logsym2(em,logem)! [logsym2] logem=matmul(vv,matmul(oo,transpose(vv))) end subroutine logsym2 -!> General routine to evaluate the logarithm, z=log(x), and the symmetric -!! derivative, zd = dz/dx, where x is a symmetric 2*2 positive-definite -!! matrix. -!! @param zd the symmetric derivative -!! @param x a symmetric 2*2 positive-definite matrix -!! @param z evaluate the logarithm log(x) +!> General routine to evaluate the logarithm, z=log(x), and the +!! symmetric derivative, zd = dz/dx, where x is a symmetric 2*2 +!! positive-definite matrix. +!! +!! @param[in] zd the symmetric derivative +!! @param[out] x a symmetric 2*2 positive-definite matrix +!! @param[out] z evaluate the logarithm log(x) !! @author R. J. Purser subroutine logsym2d(x,z,zd)! [logsym2] use pfun, only: sinhox @@ -412,7 +429,9 @@ subroutine logsym2d(x,z,zd)! [logsym2] zd(:,:,2,1)=zd(:,:,1,2) end subroutine logsym2d -!> General routine for Effective identity +!> General routine for Effective identity. +!! +!! @param[out] em ??? !! @author R. J. Purser subroutine id2222(em)! [id2222] implicit none @@ -422,11 +441,13 @@ subroutine id2222(em)! [id2222] em=id end subroutine id2222 -!> Return the cholesky lower triangular factor, C, of the 2X2 symmetric -!! matrix, S, or raise the failure flag, FF, if S is not positive-definite. -!! @param c cholesky lower triangular factor -!! @param s 2X2 symmetric matrix -!! @param ff raise the failure flag +!> Return the cholesky lower triangular factor, C, of the 2X2 +!! symmetric matrix, S, or raise the failure flag, FF, if S is not +!! positive-definite. +!! +!! @param[in] s 2X2 symmetric matrix +!! @param[out] c cholesky lower triangular factor +!! @param[out] ff raise the failure flag !! @author R. J. Purser subroutine chol2(s,c,ff)! [chol2] use pietc, only: u0