MODULE MODULE_pmat1
!$$$ module documentation block
! . . . .
! module: module_pmat1
!
! abstract: Routines for basic algebraic operations on general matrices
! and vectors
!
! additional notes:
! These routines, perform basic algebraic operations on real vectors and
! matrices. The task performed by each routine is, as far as possible,
! encoded in each routine's name; three letters describe the
! operation, the remainder defining the type of operand and, if needed to
! resolve an ambiguity, the type of result.
!
! program history log:
! 2011-07-04 todling - set to double precision to allow running GSI in
! in either single or double precision
!
! remarks:
! 1. routines here must work under REAL*8 (double precision)
!
! OPERATIONS:
! DET evaluate log-determinant
! DIF differentiate
! INT integrate
! INV invert the matrix, or linear system involving the matrix operand
! L1L Cholesky LU decomposition, where U is just L-transpose
! L1U L-U decomposition of first arg, with 1's along diagonal of L and U
! LDL Cholesky LDU decomposition, where U is just L-transpose and D diag.
! LDU LDU decomposition
! NOR evaluate norm of operand
! POL polynomial (first argument) of second argument
! POW raise operand to some integer power
! SWP swap first two operands
! TRC evaluate trace of operand
! U1L back substitution with matrix decomposed into LU form, 1's on diag.
! UDL back substitution with matrix decomposed into LDU form
! WRT write out
! ZER set operand to zero
!
! OPERAND TYPES:
! B banded matrix
! C circulant matrix
! D diagonal matrix
! H symmetric or hermitian matrix
! L lower triangular matrix
! M matrix (rectangular, in general)
! P polynomial or power-series coefficient vector
! Q sQuare matrix with Fortran dimension same as logical dimension
! R row of a matrix
! S scalar
! T transpose of the matrix
! U upper triangular matrix
! V vector, or column of a matrix
! X field of parallel X-vectors (aligned like "columns" of a matrix)
! Y field of parallel Y-vectors (aligned like "rows" of a matrix)
!
! program history:
! 1994- - R.J.Purser - initial coding
! 2008-04-25 safford - add standard documentation blocks
!
! subroutines included:
! pro333
! dpro333
! cro33
! dcro33
! norv
! dnorv
! norq
! dnorq
! swpvv
! dswpvv
! mulmd
! dmulmd
! multd
! dmultd
! muldm
! dmuldm
! muldt
! dmuldt
! mulpp
! dmulpp
! madpp
! dmadpp
! msbpp
! dmsbpp
! difp
! ddifp
! intp
! dintp
! invp
! dinvp
! prgv
! dprgv
! mulcc
! dmulcc
! madcc
! dmadcc
! msbcc
! dmsbcc
! zerl
! dzerl
! zeru
! dzeru
! ldum
! dldum
! udlmm, udlmv
! dudlmm,dudlmv
! linvan
! dlinvan
! copdm
! dcopdm
! condm
! dcondm
! copsm
! dcopsm
! consm
! dconsm
! addmd
! daddmd
! submd
! dsubmd
! addms
! daddms
! subms
! dsubms
! l1lm
! dl1lm
! ldlm
! dldlm
! invh
! dinvh
! invl
! dinvl
! linlv
! dlinlv
! linuv
! dlinuv
! powp
! dpowp
! polps
! dpolps
! polpp
! dpolpp
! trcm
! dtrcm
! invmt, linmmt, linmvt
! dinvmt,dlinmmt,dlinmvt
!
! variable definitions:
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: i_kind
use kinds, only: r_kind => r_double
use constants,only: zero,one
IMPLICIT NONE
! set default to private
private
! set subroutines/interfaces to public
public :: pro333
public :: pro333_d
public :: cro33
public :: cro33_d
public :: norv
public :: norv_d
public :: norq
public :: norq_d
public :: swpvv
public :: swpvv_d
public :: mulmd
public :: mulmd_d
public :: multd
public :: multd_d
public :: muldm
public :: muldm_d
public :: muldt
public :: muldt_d
public :: mulpp
public :: mulpp_d
public :: madpp
public :: madpp_d
public :: msbpp
public :: msbpp_d
public :: difp
public :: difp_d
public :: intp
public :: intp_d
public :: invp
public :: invp_d
public :: prgv
public :: prgv_d
public :: mulcc
public :: mulcc_d
public :: madcc
public :: madcc_d
public :: msbcc
public :: msbcc_d
public :: zerl
public :: zerl_d
public :: zeru
public :: zeru_d
public :: ldum
public :: ldum_d
public :: udlmm
public :: udlmm_d
public :: linvan
public :: linvan_d
public :: copdm
public :: copdm_d
public :: condm
public :: condm_d
public :: copsm
public :: copsm_d
public :: consm
public :: consm_d
public :: addmd
public :: addmd_d
public :: submd
public :: submd_d
public :: addms
public :: addms_d
public :: subms
public :: subms_d
public :: l1lm
public :: l1lm_d
public :: ldlm
public :: ldlm_d
public :: invh
public :: invh_d
public :: invl
public :: invl_d
public :: linlv
public :: linlv_d
public :: linuv
public :: linuv_d
public :: powp
public :: powp_d
public :: polps
public :: polps_d
public :: polpp
public :: polpp_d
public :: trcm
public :: trcm_d
public :: inv
public :: inv_d
INTERFACE pro333 ; MODULE PROCEDURE pro333; END INTERFACE
INTERFACE pro333_d; MODULE PROCEDURE dpro333; END INTERFACE
INTERFACE cro33 ; MODULE PROCEDURE cro33; END INTERFACE
INTERFACE cro33_d; MODULE PROCEDURE dcro33; END INTERFACE
INTERFACE norv; MODULE PROCEDURE norv; END INTERFACE
INTERFACE norv_d; MODULE PROCEDURE dnorv; END INTERFACE
INTERFACE norq; MODULE PROCEDURE norq; END INTERFACE
INTERFACE norq_d; MODULE PROCEDURE dnorq; END INTERFACE
INTERFACE swpvv; MODULE PROCEDURE swpvv; END INTERFACE
INTERFACE swpvv_d; MODULE PROCEDURE dswpvv; END INTERFACE
INTERFACE mulmd; MODULE PROCEDURE mulmd; END INTERFACE
INTERFACE mulmd_d; MODULE PROCEDURE dmulmd; END INTERFACE
INTERFACE multd; MODULE PROCEDURE multd; END INTERFACE
INTERFACE multd_d; MODULE PROCEDURE dmultd; END INTERFACE
INTERFACE muldm; MODULE PROCEDURE muldm; END INTERFACE
INTERFACE muldm_d; MODULE PROCEDURE dmuldm; END INTERFACE
INTERFACE muldt; MODULE PROCEDURE muldt; END INTERFACE
INTERFACE muldt_d; MODULE PROCEDURE dmuldt; END INTERFACE
INTERFACE mulpp; MODULE PROCEDURE mulpp; END INTERFACE
INTERFACE mulpp_d; MODULE PROCEDURE dmulpp; END INTERFACE
INTERFACE madpp; MODULE PROCEDURE madpp; END INTERFACE
INTERFACE madpp_d; MODULE PROCEDURE dmadpp; END INTERFACE
INTERFACE msbpp; MODULE PROCEDURE msbpp; END INTERFACE
INTERFACE msbpp_d; MODULE PROCEDURE dmsbpp; END INTERFACE
INTERFACE difp; MODULE PROCEDURE difp; END INTERFACE
INTERFACE difp_d; MODULE PROCEDURE ddifp; END INTERFACE
INTERFACE intp; MODULE PROCEDURE intp; END INTERFACE
INTERFACE intp_d; MODULE PROCEDURE dintp; END INTERFACE
INTERFACE invp; MODULE PROCEDURE invp; END INTERFACE
INTERFACE invp_d; MODULE PROCEDURE dinvp; END INTERFACE
INTERFACE prgv; MODULE PROCEDURE prgv; END INTERFACE
INTERFACE prgv_d; MODULE PROCEDURE dprgv; END INTERFACE
INTERFACE mulcc; MODULE PROCEDURE mulcc; END INTERFACE
INTERFACE mulcc_d; MODULE PROCEDURE dmulcc; END INTERFACE
INTERFACE madcc; MODULE PROCEDURE madcc; END INTERFACE
INTERFACE madcc_d; MODULE PROCEDURE dmadcc; END INTERFACE
INTERFACE msbcc; MODULE PROCEDURE msbcc; END INTERFACE
INTERFACE msbcc_d; MODULE PROCEDURE dmsbcc; END INTERFACE
INTERFACE zerl; MODULE PROCEDURE zerl; END INTERFACE
INTERFACE zerl_d; MODULE PROCEDURE dzerl; END INTERFACE
INTERFACE zeru; MODULE PROCEDURE zeru; END INTERFACE
INTERFACE zeru_d; MODULE PROCEDURE dzeru; END INTERFACE
INTERFACE ldum; MODULE PROCEDURE ldum; END INTERFACE
INTERFACE ldum_d; MODULE PROCEDURE dldum; END INTERFACE
INTERFACE udlmm; MODULE PROCEDURE udlmm, udlmv; END INTERFACE
INTERFACE udlmm_d; MODULE PROCEDURE dudlmm,dudlmv; END INTERFACE
INTERFACE linvan; MODULE PROCEDURE linvan; END INTERFACE
INTERFACE linvan_d; MODULE PROCEDURE dlinvan; END INTERFACE
INTERFACE copdm; MODULE PROCEDURE copdm; END INTERFACE
INTERFACE copdm_d; MODULE PROCEDURE dcopdm; END INTERFACE
INTERFACE condm; MODULE PROCEDURE condm; END INTERFACE
INTERFACE condm_d; MODULE PROCEDURE dcondm; END INTERFACE
INTERFACE copsm; MODULE PROCEDURE copsm; END INTERFACE
INTERFACE copsm_d; MODULE PROCEDURE dcopsm; END INTERFACE
INTERFACE consm; MODULE PROCEDURE consm; END INTERFACE
INTERFACE consm_d; MODULE PROCEDURE dconsm; END INTERFACE
INTERFACE addmd; MODULE PROCEDURE addmd; END INTERFACE
INTERFACE addmd_d; MODULE PROCEDURE daddmd; END INTERFACE
INTERFACE submd; MODULE PROCEDURE submd; END INTERFACE
INTERFACE submd_d; MODULE PROCEDURE dsubmd; END INTERFACE
INTERFACE addms; MODULE PROCEDURE addms; END INTERFACE
INTERFACE addms_d; MODULE PROCEDURE daddms; END INTERFACE
INTERFACE subms; MODULE PROCEDURE subms; END INTERFACE
INTERFACE subms_d; MODULE PROCEDURE dsubms; END INTERFACE
INTERFACE l1lm; MODULE PROCEDURE l1lm; END INTERFACE
INTERFACE l1lm_d; MODULE PROCEDURE dl1lm; END INTERFACE
INTERFACE ldlm; MODULE PROCEDURE ldlm; END INTERFACE
INTERFACE ldlm_d; MODULE PROCEDURE dldlm; END INTERFACE
INTERFACE invh; MODULE PROCEDURE invh; END INTERFACE
INTERFACE invh_d; MODULE PROCEDURE dinvh; END INTERFACE
INTERFACE invl; MODULE PROCEDURE invl; END INTERFACE
INTERFACE invl_d; MODULE PROCEDURE dinvl; END INTERFACE
INTERFACE linlv; MODULE PROCEDURE linlv; END INTERFACE
INTERFACE linlv_d; MODULE PROCEDURE dlinlv; END INTERFACE
INTERFACE linuv; MODULE PROCEDURE linuv; END INTERFACE
INTERFACE linuv_d; MODULE PROCEDURE dlinuv; END INTERFACE
INTERFACE powp; MODULE PROCEDURE powp; END INTERFACE
INTERFACE powp_d; MODULE PROCEDURE dpowp; END INTERFACE
INTERFACE polps; MODULE PROCEDURE polps; END INTERFACE
INTERFACE polps_d; MODULE PROCEDURE dpolps; END INTERFACE
INTERFACE polpp; MODULE PROCEDURE polpp; END INTERFACE
INTERFACE polpp_d; MODULE PROCEDURE dpolpp; END INTERFACE
INTERFACE trcm; MODULE PROCEDURE trcm; END INTERFACE
INTERFACE trcm_d; MODULE PROCEDURE dtrcm; END INTERFACE
INTERFACE inv; MODULE PROCEDURE invmt, linmmt, linmvt; END INTERFACE
INTERFACE inv_d; MODULE PROCEDURE dinvmt,dlinmmt,dlinmvt; END INTERFACE
CONTAINS
FUNCTION pro333(d,e,f) RESULT(pro_res)
!$$$ subprogram documentation block
! . . .
! subprogram: pro333
!
! prgrmmr:
!
! abstract: triple product of 3 3-vectors
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! d(3), e(3), f(3) - input vectors
!
! output argument list:
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(IN ) :: d(3), e(3), f(3)
REAL(r_kind) :: pro_res
REAL(r_kind) :: g(3)
CALL CRO33(E,F,G)
pro_res=DOT_PRODUCT(d,g)
END FUNCTION pro333
FUNCTION dpro333(d,e,f) RESULT(pro_res)
!$$$ subprogram documentation block
! . . .
! subprogram: dpro333
!
! prgrmmr:
!
! abstract: triple product of 3 3-vectors
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! d(3), e(3), f(3) - input vectors
!
! output argument list:
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(IN ) :: d(3), e(3), f(3)
REAL(r_kind) :: pro_res
REAL(r_kind) :: g(3)
CALL CRO33_d(E,F,G)
pro_res=DOT_PRODUCT(d,g)
END FUNCTION dpro333
SUBROUTINE cro33(a,b,c)
!$$$ subprogram documentation block
! . . .
! subprogram: cro33
!
! prgrmmr:
!
! abstract: special case of 3-dimensions: cross-product
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a(3), b(3) - input vectors
!
! output argument list:
! c(3) - resulting cross-product
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block$
implicit none
REAL(r_kind), INTENT(IN ) :: a(3), b(3)
REAL(r_kind), INTENT( OUT) :: c(3)
c(1)=a(2)*b(3)-a(3)*b(2)
c(2)=a(3)*b(1)-a(1)*b(3)
c(3)=a(1)*b(2)-a(2)*b(1)
END SUBROUTINE cro33
SUBROUTINE dcro33(a,b,c)
!$$$ subprogram documentation block
! . . .
! subprogram: dcro33
!
! prgrmmr:
!
! abstract: special case of 3-dimensions: cross-product
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a(3), b(3) - input vectors
!
! output argument list:
! c(3) - resulting cross-product
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(IN ) :: a(3), b(3)
REAL(r_kind), INTENT( OUT) :: c(3)
c(1)=a(2)*b(3)-a(3)*b(2)
c(2)=a(3)*b(1)-a(1)*b(3)
c(3)=a(1)*b(2)-a(2)*b(1)
END SUBROUTINE dcro33
FUNCTION norv(d) RESULT(norv_res)
!$$$ subprogram documentation block
! . . .
! subprogram: norv
!
! prgrmmr:
!
! abstract: norm of vector
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! d - input vector
!
! output argument list:
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(IN ) :: d(:)
REAL(r_kind) :: norv_res
norv_res=SQRT(DOT_PRODUCT(D,D))
END FUNCTION norv
FUNCTION dnorv(d)
!$$$ subprogram documentation block
! . . .
! subprogram: dnorv
!
! prgrmmr:
!
! abstract: norm of vector
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! d - input vector
!
! output argument list:
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(IN ) :: d(:)
REAL(r_kind):: dnorv
dnorv=SQRT(DOT_PRODUCT(d,d))
END FUNCTION dnorv
FUNCTION norq(d)
!$$$ subprogram documentation block
! . . .
! subprogram: norq
!
! prgrmmr:
!
! abstract: norm of a matrix
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! d - input matrix
!
! output argument list:
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind),INTENT(IN ) :: d(:,:)
REAL(r_kind):: norq
INTEGER(i_kind) m2,i2
m2=SIZE(d,2)
norq=zero; DO i2=1,m2; norq=norq+dot_PRODUCT(d(:,i2),d(:,i2)); ENDDO
norq=SQRT(norq)
END FUNCTION norq
FUNCTION dnorq(d) ! norm of a matrix
!$$$ subprogram documentation block
! . . .
! subprogram: dnorq
!
! prgrmmr:
!
! abstract: norm of a matrix
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! d - input matrix
!
! output argument list:
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind),INTENT(IN ) :: d(:,:)
REAL(r_kind):: dnorq
INTEGER(i_kind) m2,i2
m2=SIZE(d,2)
dnorq=zero; DO i2=1,m2; dnorq=dnorq+dot_PRODUCT(d(:,i2),d(:,i2)); ENDDO
dnorq=SQRT(dnorq)
END FUNCTION dnorq
SUBROUTINE swpvv(d,e)
!$$$ subprogram documentation block
! . . .
! subprogram: swpvv
!
! prgrmmr:
!
! abstract: swap first two operands of input vectors
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! d, e -
!
! output argument list:
! d, e -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(INOUT) :: d(:), e(:)
REAL(r_kind) :: t(SIZE(d))
t = d; d = e; e = t
END SUBROUTINE swpvv
SUBROUTINE dswpvv(d,e)
!$$$ subprogram documentation block
! . . .
! subprogram: dswpvv
!
! prgrmmr:
!
! abstract: swap first two operads of input vectors
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! d, e -
!
! output argument list:
! d, e -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(INOUT) :: d(:), e(:)
REAL(r_kind) :: t(SIZE(d))
t = d; d = e; e = t
END SUBROUTINE dswpvv
SUBROUTINE mulmd(a,d,b)
!$$$ subprogram documentation block
! . . .
! subprogram: mulmd
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a, b -
! d -
! output argument list:
! a, b -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(INOUT) :: a(:,:),b(:,:)
REAL(r_kind), INTENT(IN ) :: d(*)
INTEGER(i_kind):: m2,j
m2=SIZE(a,2)
DO j=1,m2; b(:,j)=a(:,j)*d(j); ENDDO
END SUBROUTINE mulmd
SUBROUTINE dmulmd(a,d,b)
!$$$ subprogram documentation block
! . . .
! subprogram: dmulmd
!
! prgrmmr:
!
! abstract: special case of 3-dimensions: cross-product
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a, b -
! d -
!
! output argument list:
! a, b -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(INOUT) :: a(:,:),b(:,:)
REAL(r_kind), INTENT(IN ) :: d(*)
INTEGER(i_kind):: m2,j
m2=SIZE(a,2)
DO j=1,m2; b(:,j)=a(:,j)*d(j); ENDDO
END SUBROUTINE dmulmd
SUBROUTINE multd(a,d,b)
!$$$ subprogram documentation block
! . . .
! subprogram: multd
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a, b -
! d -
!
! output argument list:
! a, b -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(INOUT) :: a(:,:),b(:,:)
REAL(r_kind), INTENT(IN ) :: d(*)
INTEGER(i_kind):: m2,j
m2=SIZE(a,1)
DO j=1,m2; b(:,j) = a(j,:) * d(j); ENDDO
END SUBROUTINE multd
SUBROUTINE dmultd(a,d,b)
!$$$ subprogram documentation block
! . . .
! subprogram: dmultd
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a, b -
! d -
!
! output argument list:
! a, b -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(INOUT) :: a(:,:),b(:,:)
REAL(r_kind), INTENT(IN ) :: d(*)
INTEGER(i_kind):: m2,j
m2=SIZE(a,1)
DO j=1,m2; b(:,j) = a(j,:) * d(j); ENDDO
END SUBROUTINE dmultd
SUBROUTINE muldm(d,a,b)
!$$$ subprogram documentation block
! . . .
! subprogram: muldm
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a, b -
! d -
!
! output argument list:
! a, b -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(INOUT) :: a(:,:),b(:,:)
REAL(r_kind), INTENT(IN ) :: d(*)
INTEGER(i_kind) :: m1,i
m1=SIZE(a,1)
DO i=1,m1; b(i,:) = d(i)*a(i,:); ENDDO
END SUBROUTINE muldm
SUBROUTINE dmuldm(d,a,b)
!$$$ subprogram documentation block
! . . .
! subprogram: dmuldm
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a, b -
! d -
!
! output argument list:
! a, b -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(INOUT) :: a(:,:),b(:,:)
REAL(r_kind), INTENT(IN ) :: d(*)
INTEGER(i_kind) :: m1,i
m1=SIZE(a,1)
DO i=1,m1; b(i,:) = d(i)*a(i,:); ENDDO
END SUBROUTINE dmuldm
SUBROUTINE muldt(d,a,b)
!$$$ subprogram documentation block
! . . .
! subprogram: muldt
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a, b -
! d -
!
! output argument list:
! a, b -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(INOUT) :: a(:,:),b(:,:)
REAL(r_kind), INTENT(IN ) :: d(*)
INTEGER(i_kind) :: m1,i
m1=SIZE(a,2)
DO i=1,m1; b(i,:) = d(i)*a(:,i); ENDDO
END SUBROUTINE muldt
SUBROUTINE dmuldt(d,a,b)
!$$$ subprogram documentation block
! . . .
! subprogram: dmuldt
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a, b -
! d -
!
! output argument list:
! a, b -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
REAL(r_kind), INTENT(INOUT) :: a(:,:),b(:,:)
REAL(r_kind), INTENT(IN ) :: d(*)
INTEGER(i_kind):: m1,i
m1=SIZE(a,2)
DO i=1,m1; b(i,:) = d(i)*a(:,i); ENDDO
END SUBROUTINE dmuldt
SUBROUTINE mulpp(a,b,c)
!$$$ subprogram documentation block
! . . .
! subprogram: mulpp
!
! prgrmmr:
!
! abstract: multiply polynomials, possibly in place
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block, rm unused vars
!
! input argument list:
! a, b -
! c -
!
! output argument list:
! c -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(IN ) :: a(0:), b(0:)
REAL(r_kind), INTENT(INOUT) :: c(0:)
INTEGER(i_kind) :: m,mcp, j
REAL(r_kind) :: s
m=SIZE(a)-1
mcp=mcmax(a,b,m)
c(mcp:m) = zero
DO j=mcp,1,-1
s = SUM(a(j-1:0:-1)*b(0:j-1))
c(j-1)=s
ENDDO
RETURN
ENTRY madpp(a,b,c)
m=SIZE(a)-1
mcp=mcmax(a,b,m)
DO j=mcp,1,-1
s = SUM(a(j-1:0:-1)*b(0:j-1))
c(j-1)=c(j-1)+s
ENDDO
RETURN
ENTRY msbpp(a,b,c)
m=SIZE(a)-1
mcp=mcmax(a,b,m)
DO j=mcp,1,-1
s = SUM(a(j-1:0:-1)*b(0:j-1))
c(j-1)=c(j-1)-s
ENDDO
RETURN
CONTAINS
FUNCTION mcmax(a,b,m) RESULT(mmx_res) ! This fn can be contained in mulpp().
!$$$ subprogram documentation block
! . . . .
! subprogram: mcmax
! prgmmr:
!
! abstract:
!
! program history log:
! 2009-08-26 lueken - added subprogram doc block
!
! input argument list:
! m
! a,b
!
! output argument list:
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m
REAL(r_kind), INTENT(IN ) :: a(0:m), b(0:m)
INTEGER(i_kind) :: mmx_res
INTEGER(i_kind) :: ma, mb
mmx_res=0 ! default for when ALL elements of c are zero
DO ma=m,0,-1 ! seek last nonzero coefficient of polynomial a
IF(a(ma) /= zero)THEN
DO mb=m,0,-1 ! seek last nonzero coefficient of polynomial b
IF(b(mb) /= zero)THEN
mmx_res=MIN(m,ma+mb)+1 ! hence, 1+last non-0 element of their product
RETURN
ENDIF
ENDDO
RETURN
ENDIF
ENDDO
END FUNCTION mcmax
END SUBROUTINE mulpp
SUBROUTINE difp(a,b) ! Symbolically differentiate polynomial
!$$$ subprogram documentation block
! . . .
! subprogram: difp
!
! prgrmmr:
!
! abstract: Symbolically differentiate polynomial
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block, rm unused vars
!
! input argument list:
! a -
!
! output argument list:
! b -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(IN ) :: a(0:)
REAL(r_kind), INTENT( OUT) :: b(0:)
INTEGER(i_kind) :: m, i
REAL(r_kind) :: s, b0
m=SIZE(a)-1
DO i=1,m ! possibly with coincident storage for a and b
b(i-1)=i*a(i)
ENDDO
b(m)=zero
RETURN
ENTRY intp(a,b) ! Symbolically integrate polynomial
m=SIZE(a)-1
DO i=m,1,-1 ! possibly with coincident storage for a and b
b(i)=a(i-1)/i
ENDDO
b(0)=zero
RETURN
ENTRY invp(a,b) ! Invert polynomial or power-series
m=SIZE(a)-1
b0=one/a(0) ! storage of a and b must not be the same
b(0)=b0
DO i=1,m
s = SUM(b(i-1:0:-1)*a(1:i))
b(i)=-b0*s
ENDDO
END SUBROUTINE difp
SUBROUTINE dmulpp(a,b,c) ! multiply polynomials, possibly in place
!$$$ subprogram documentation block
! . . .
! subprogram: dmulpp
!
! prgrmmr:
!
! abstract: multiply polynomials, possibly in place
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block, rm unused vars
!
! input argument list:
! a, b -
! c -
!
! output argument list:
! c -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(IN ) :: a(0:), b(0:)
REAL(r_kind), INTENT(INOUT) :: c(0:)
INTEGER(i_kind) :: m,mcp, j
REAL(r_kind) :: s
m=SIZE(a)-1
mcp=mcmax(a,b,m)
c(mcp:m) = zero
DO j=mcp,1,-1
s = SUM(a(j-1:0:-1)*b(0:j-1))
c(j-1)=s
ENDDO
RETURN
ENTRY dmadpp(a,b,c)
m=SIZE(a)-1
mcp=mcmax(a,b,m)
DO j=mcp,1,-1
s = SUM(a(j-1:0:-1)*b(0:j-1))
c(j-1)=c(j-1)+s
ENDDO
RETURN
ENTRY dmsbpp(a,b,c)
m=SIZE(a)-1
mcp=mcmax(a,b,m)
DO j=mcp,1,-1
s = SUM(a(j-1:0:-1)*b(0:j-1))
c(j-1)=c(j-1)-s
ENDDO
RETURN
CONTAINS
FUNCTION mcmax(a,b,m) RESULT(mmx_res)
!$$$ subprogram documentation block
! . . . .
! subprogram: mcmax
! prgmmr:
!
! abstract:
!
! program history log:
! 2009-08-26 lueken - added subprogram doc block
!
! input argument list:
! m
! a,b
!
! output argument list:
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m
REAL(r_kind) , INTENT(IN ) :: a(0:m), b(0:m)
INTEGER(i_kind) :: mmx_res
INTEGER(i_kind) :: ma, mb
mmx_res=0 ! default for when all elements of c are zero
DO ma=m,0,-1 ! seek last nonzero coefficient of polynomial a
IF(a(ma) /= zero)THEN
DO mb=m,0,-1 ! seek last nonzero coefficient of polynomial b
IF(b(mb) /= zero)THEN
mmx_res=MIN(m,ma+mb)+1 ! hence, 1+last non-0 element of their product
RETURN
ENDIF
ENDDO
RETURN
ENDIF
ENDDO
RETURN
END FUNCTION mcmax
END SUBROUTINE dmulpp
SUBROUTINE ddifp(a,b) ! Symbolically differentiate polynomial
!$$$ subprogram documentation block
! . . .
! subprogram: ddifp
!
! prgrmmr:
!
! abstract: Symbolically differentiate polynomial
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block, rm unused vars
!
! input argument list:
! a -
! b -
!
! output argument list:
! b -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(IN ) :: a(0:)
REAL(r_kind), INTENT(INOUT) :: b(0:)
INTEGER(i_kind) :: m, i
REAL(r_kind) :: s, b0
m=SIZE(a)-1
DO i=1,m ! possibly with coincident storage for a and b
b(i-1)=i*a(i)
ENDDO
b(m)=zero
RETURN
ENTRY dintp(a,b) ! Symbolically integrate polynomial
m=SIZE(a)-1
DO i=m,1,-1 ! possibly with coincident storage for a and b
b(i)=a(i-1)/i
ENDDO
b(0)=zero
RETURN
ENTRY dinvp(a,b) ! Invert polynomial or power-series
m=SIZE(a)-1
b0=one/a(0) ! storage of a and b must not be the same
b(0)=b0
DO i=1,m
s = SUM(b(i-1:0:-1)*a(1:i))
b(i)=-b0*s
ENDDO
END SUBROUTINE ddifp
SUBROUTINE prgv(d)
!$$$ subprogram documentation block
! . . .
! subprogram: prgv
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! d -
!
! output argument list:
! d -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(INOUT) :: d(:)
REAL(r_kind), PARAMETER :: crit=1.E-30_r_kind
INTEGER(i_kind) :: i,m
m=SIZE(d)
DO i=1,m; IF(ABS(d(i)) <= crit)d(i)=zero; ENDDO
END SUBROUTINE prgv
SUBROUTINE dprgv(d)
!$$$ subprogram documentation block
! . . .
! subprogram: dprgv
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! d -
!
! output argument list:
! d -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(INOUT) :: d(:)
REAL(r_kind), PARAMETER :: crit=1.E-30_r_kind
INTEGER(i_kind) :: i,m
m=SIZE(d)
DO i=1,m; IF(ABS(d(i)) <= crit)d(i)=zero; ENDDO
END SUBROUTINE dprgv
SUBROUTINE mulcc(a,b,c,m)
!$$$ subprogram documentation block
! . . .
! subprogram: mulcc
!
! prgrmmr:
!
! abstract: Multiply circulant matrices of period M
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a, b, c -
! m -
!
! output argument list:
! a, b, c -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m
REAL(r_kind) , INTENT(INOUT) :: a(0:m-1), b(0:m-1), c(0:m-1)
INTEGER(i_kind) :: mm, j
c(0:m-1) = zero
ENTRY madcc(a,b,c,m)
mm=m-1
DO j=0,mm
c(j:m-1) = c(j:m-1) + a(0:m-j-1)*b(j)
c(0:j-1) = c(0:j-1) + a(m-j:m-1)*b(j)
ENDDO
RETURN
ENTRY msbcc(a,b,c,m)
mm=m-1
DO j=0,mm
c(j:m-1) = c(j:m-1) - a(0:m-j-1)*b(j)
c(0:j-1) = c(0:j-1) - a(m-j:m-1)*b(j)
ENDDO
END SUBROUTINE mulcc
SUBROUTINE dmulcc(a,b,c,m) ! Multiply circulant matrices of period M
!$$$ subprogram documentation block
! . . .
! subprogram: dmulcc
!
! prgrmmr:
!
! abstract: Multiply circulant matrices of period M
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a, b, c -
! m -
!
! output argument list:
! a, b, c -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m
REAL(r_kind) , INTENT(INOUT) :: a(0:m-1), b(0:m-1), c(0:m-1)
INTEGER(i_kind) :: mm, j
c(0:m-1) = zero
ENTRY dmadcc(a,b,c,m)
mm=m-1
DO j=0,mm
c(j:m-1) = c(j:m-1) + a(0:m-j-1)*b(j)
c(0:j-1) = c(0:j-1) + a(m-j:m-1)*b(j)
ENDDO
RETURN
ENTRY dmsbcc(a,b,c,m)
mm=m-1
DO j=0,mm
c(j:m-1) = c(j:m-1) - a(0:m-j-1)*b(j)
c(0:j-1) = c(0:j-1) - a(m-j:m-1)*b(j)
ENDDO
END SUBROUTINE dmulcc
SUBROUTINE zerl(a)
!$$$ subprogram documentation block
! . . .
! subprogram: zerl
!
! prgrmmr:
!
! abstract: Zero out the strictly lower triangle of elements
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a -
!
! output argument list:
! a -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind),INTENT(INOUT) :: a(:,:)
INTEGER(i_kind) :: m,j
m=SIZE(a,1); DO j=1,m; a(j+1:m,j) = zero; ENDDO; RETURN
ENTRY zeru(a) ! Zero out the strictly upper triangle of elements
m=SIZE(a,1); DO j=1,m; a(1:j-1,j) = zero; ENDDO
END SUBROUTINE zerl
SUBROUTINE dzerl(a)
!$$$ subprogram documentation block
! . . .
! subprogram: dmuldm
!
! prgrmmr:
!
! abstract: Zero out the strictly lower triangle of elements
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a -
!
! output argument list:
! a -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind),INTENT(INOUT) :: a(:,:)
INTEGER(i_kind) :: m,j
m=SIZE(a,1); DO j=1,m; a(j+1:m,j) = zero; ENDDO; RETURN
ENTRY dzeru(a) ! Zero out the strictly upper triangle of elements
m=SIZE(a,1); DO j=1,m; a(1:j-1,j) = zero; ENDDO
END SUBROUTINE dzerl
SUBROUTINE ldum(a,ipiv,d)
!$$$ subprogram documentation block
! . . .
! subprogram: ldum
!
! prgrmmr: R.J.Purser, NCEP, Washington D.C. 1996
!
! abstract: perform l-d-u decomposition of square matrix a in place with
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a - square matrix to be factorized
!
! output argument list:
! a - square matrix to be factorized
! ipiv - ipiv array encoding the pivoting sequence
! d - indicator for possible sign change of determinant
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(INOUT) :: a(:,:)
REAL(r_kind), INTENT(OUT ) :: d
INTEGER(i_kind), INTENT(OUT ) :: ipiv(:)
INTEGER(i_kind) :: m,i, j, jp, ibig, jm
REAL(r_kind) :: s(SIZE(a,1)), aam, aa, abig, ajj, ajji, aij
m=SIZE(a,1)
DO i=1,m
aam=zero
DO j=1,m
aa=ABS(a(i,j))
IF(aa > aam)aam=aa
ENDDO
IF(aam == zero)THEN
PRINT '(" row ",i3," of matrix in ldum vanishes")',i
STOP
ENDIF
s(i)=one/aam
ENDDO
d=one
ipiv(m)=m
DO j=1,m-1
jp=j+1
abig=s(j)*ABS(a(j,j))
ibig=j
DO i=jp,m
aa=s(i)*ABS(a(i,j))
IF(aa > abig)THEN
ibig=i
abig=aa
ENDIF
ENDDO
! swap rows, recording changed sign of determinant
ipiv(j)=ibig
IF(ibig /= j)THEN
d=-d
CALL swpvv(a(j,:),a(ibig,:))
s(ibig)=s(j)
ENDIF
ajj=a(j,j)
IF(ajj == zero)THEN
jm=j-1
PRINT '(" failure in ldum:"/" matrix singular, rank=",i3)',jm
STOP
ENDIF
ajji=one/ajj
DO i=jp,m
aij=ajji*a(i,j)
a(i,j)=aij
a(i,jp:m) = a(i,jp:m) - aij*a(j,jp:m)
ENDDO
ENDDO
END SUBROUTINE ldum
SUBROUTINE DLDUM(A,IPIV,D)
!$$$ subprogram documentation block
! . . .
! subprogram: dldum
!
! prgrmmr: R.J.Purser, NCEP, Washington D.C. 1996
!
! abstract: perform l-d-u decomposition of square matrix a in place with
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a - square matrix to be factorized
!
! output argument list:
! a - square matrix to be factorized
! ipiv - ipiv array encoding the pivoting sequence
! d - indicator for possible sign change of determinant
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind) , INTENT(INOUT) :: a(:,:)
REAL(r_kind) , INTENT( OUT) :: d
INTEGER(i_kind), INTENT( OUT) :: ipiv(:)
INTEGER(i_kind) :: m,i, j, jp, ibig, jm
REAL(r_kind) :: s(SIZE(a,1)), aam, aa, abig, ajj, ajji, aij
m=SIZE(a,1)
DO i=1,m
aam=zero
DO j=1,m
aa=ABS(a(i,j))
IF(aa > aam)aam=aa
ENDDO
IF(aam == zero)THEN
PRINT '(" row ",i3," of matrix in dldum vanishes")',i
STOP
ENDIF
s(i)=one/aam
ENDDO
d=one
ipiv(m)=m
DO j=1,m-1
jp=j+1
abig=s(j)*ABS(a(j,j))
ibig=j
DO i=jp,m
aa=s(i)*ABS(a(i,j))
IF(aa > abig)THEN
ibig=i
abig=aa
ENDIF
ENDDO
! swap rows, recording changed sign of determinant
ipiv(j)=ibig
IF(ibig /= j)THEN
d=-d
CALL swpvv_d(a(j,:),a(ibig,:))
s(ibig)=s(j)
ENDIF
ajj=a(j,j)
IF(ajj == zero)THEN
jm=j-1
PRINT '(" Failure in dldum:"/" matrix singular, rank=",i3)',jm
STOP
ENDIF
ajji=one/ajj
DO i=jp,m
aij=ajji*a(i,j)
a(i,j)=aij
a(i,jp:m) = a(i,jp:m) - aij*a(j,jp:m)
ENDDO
ENDDO
END SUBROUTINE dldum
SUBROUTINE udlmm(a,b,ipiv)
!$$$ subprogram documentation block
! . . .
! subprogram: udlmm
!
! prgrmmr: R.J.Purser, NCEP, Washington D.C. 1993
!
! abstract: use l-u factors in A to back-substitute for mm rhs in B,
! using ipiv to define the pivoting permutation used in the l-u
! decomposition.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a - L-D-U factorization of linear system matrux
! b - right-hand-sides on entry, corresponding matrix of solution
! vectors on return
! ipiv - ipiv array encoding the pivoting sequence
!
! output argument list:
! b - right-hand-sides on entry, corresponding matrix of solution
! vectors on return
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: ipiv(:)
REAL(r_kind), INTENT(IN ) :: a(:,:)
REAL(r_kind), INTENT(INOUT) :: b(:,:)
INTEGER(i_kind) :: m,mm,i, k, l
REAL(r_kind) :: s,aiii
m=SIZE(a,1); mm=SIZE(b,2)
DO k=1,mm !loop over columns of b
DO i=1,m
l=ipiv(i)
s=b(l,k)
b(l,k)=b(i,k)
s = s - SUM(b(1:i-1,k)*a(i,1:i-1))
b(i,k)=s
ENDDO
b(m,k)=b(m,k)/a(m,m)
DO i=m-1,1,-1
aiii=one/a(i,i)
b(i,k) = b(i,k) - SUM(b(i+1:m,k)*a(i,i+1:m))
b(i,k)=b(i,k)*aiii
ENDDO
ENDDO
END SUBROUTINE udlmm
SUBROUTINE dudlmm(a,b,ipiv)
!$$$ subprogram documentation block
! . . .
! subprogram: dudlmm
!
! prgrmmr: R.J.Purser, NCEP, Washington D.C. 1993
!
! abstract: use l-u factors in A to back-substitute for mm rhs in B,
! using ipiv to define the pivoting permutation used in the l-u
! decomposition.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a - L-D-U factorization of linear system matrux
! b - right-hand-sides on entry, corresponding matrix of solution
! vectors on return
! ipiv - ipiv array encoding the pivoting sequence
!
! output argument list:
! b - right-hand-sides on entry, corresponding matrix of solution
! vectors on return
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: ipiv(:)
REAL(r_kind) , INTENT(IN ) :: a(:,:)
REAL(r_kind) , INTENT(INOUT) :: b(:,:)
INTEGER(i_kind) :: m,mm,i, k, l
REAL(r_kind) :: s,aiii
m=SIZE(a,1); mm=SIZE(b,2)
DO k=1,mm !loop over columns of b
DO i=1,m
l=ipiv(i)
s=b(l,k)
b(l,k)=b(i,k)
s = s - SUM(b(1:i-1,k)*a(i,1:i-1))
b(i,k)=s
ENDDO
b(m,k)=b(m,k)/a(m,m)
DO i=m-1,1,-1
aiii=one/a(i,i)
b(i,k) = b(i,k) - SUM(b(i+1:m,k)*a(i,i+1:m))
b(i,k)=b(i,k)*aiii
ENDDO
ENDDO
END SUBROUTINE dudlmm
SUBROUTINE udlmv(a,b,ipiv)
!$$$ subprogram documentation block
! . . .
! subprogram: udlmv
!
! prgrmmr: R.J.Purser, NCEP, Washington D.C. 1993
!
! abstract: use l-u factors in A to back-substitute for mm rhs in B, using
! ipiv to define the pivoting permutation used in the l-u
! decomposition.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a - L-D-U factorization of linear system matrux
! b - right-hand-side on entry, corresponding vector solution
! on return
! ipiv - ipiv array encoding the pivoting sequence
!
! output argument list:
! b - right-hand-side on entry, corresponding vector solution
! on return
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: ipiv(:)
REAL(r_kind), INTENT(IN ) :: a(:,:)
REAL(r_kind), INTENT(INOUT) :: b(:)
INTEGER(i_kind) :: m,i, l
REAL(r_kind) :: s,aiii
m=SIZE(a,1)
DO i=1,m
l=ipiv(i)
s=b(l)
b(l)=b(i)
s = s - SUM(b(1:i-1)*a(i,1:i-1))
b(i)=s
ENDDO
b(m)=b(m)/a(m,m)
DO i=m-1,1,-1
aiii=one/a(i,i)
b(i) = b(i) - SUM(b(i+1:m)*a(i,i+1:m))
b(i)=b(i)*aiii
ENDDO
END SUBROUTINE udlmv
SUBROUTINE dudlmv(a,b,ipiv)
!$$$ subprogram documentation block
! . . .
! subprogram: dudlmv
!
! prgrmmr: R.J.Purser, NCEP, Washington D.C. 1993
!
! abstract: use l-u factors in A to back-substitute for mm rhs in B, using
! ipiv to define the pivoting permutation used in the l-u
! decomposition.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a - L-D-U factorization of linear system matrux
! b - right-hand-side on entry, corresponding vector solution
! on return
! ipiv - ipiv array encoding the pivoting sequence
!
! output argument list:
! b - right-hand-side on entry, corresponding vector solution
! on return
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: ipiv(:)
REAL(r_kind) , INTENT(IN ) :: a(:,:)
REAL(r_kind) , INTENT(INOUT) :: b(:)
INTEGER(i_kind) :: m,i, l
REAL(r_kind) :: s,aiii
m=SIZE(a,1)
DO i=1,m
l=ipiv(i)
s=b(l)
b(l)=b(i)
s = s - SUM(b(1:i-1)*a(i,1:i-1))
b(i)=s
ENDDO
b(m)=b(m)/a(m,m)
DO i=m-1,1,-1
aiii=one/a(i,i)
b(i) = b(i) - SUM(b(i+1:m)*a(i,i+1:m))
b(i)=b(i)*aiii
ENDDO
END SUBROUTINE dudlmv
SUBROUTINE linvan(w,ab)
!$$$ subprogram documentation block
! . . .
! subprogram: linvan
!
! prgrmmr: R.J.Purser, NCEP, Washington D.C. 1993
!
! abstract:
! Take square matrix W and seek row and column scalings to produce non-
! vanishing elements of rescaled W having magnitudes as close to unity
! as possible. The approach is make the geometric mean of the nonvanishing
! elements of each row and of each column +1 or -1. Having rescaled the
! matrix and the r.h.s. vector AB, compute the product P of row-vector
! norms, then compute the determinant D and solve the linear system.
! Rescale the solution vector (now AB) and put the conditioning indicator
! formed by the ratio D/P into the first element of W.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block, rm unused vars
!
! input argument list:
! W - Generalized Vandermonde matrix in, conditioning indicator out.
! AB - R.h.s. vector in, solution vector of numerical coefficients out.
!
! output argument list:
! W - Generalized Vandermonde matrix in, conditioning indicator out.
! AB - R.h.s. vector in, solution vector of numerical coefficients out.
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(INOUT) :: w(:,:), ab(:)
INTEGER(i_kind), PARAMETER :: nit=20
REAL(r_kind) :: d1(SIZE(w,1)), d2(SIZE(w,1)), &
w2(SIZE(w,1),SIZE(w,1)),v(SIZE(w,1))
INTEGER(i_kind) :: i, j, it, jt, ipiv(SIZE(w,1)), nc
REAL(r_kind) :: p, e, dw, c, d, d2j
REAL(r_kind),ALLOCATABLE :: wv(:,:) ! work variable for ab(nc) and v(nn)
nc = SIZE(w,DIM=1)
ALLOCATE(wv(nc,1))
w2=w ! Preserve original W and AB for use
v = ab(1:nc) ! in later "clean-up" operation.
d1 = one ! Row scaling factors set to default
d2 = one ! Column scaling factors set to default
C=1.E-16_r_kind ! Set initial criterion for "negligible" elements of W
! In first attempt to estimate row and column scalings, use logarithms
! to avoid the risk of under- or over-flows of the line products of W:
DO i=1,nc
p=zero
e=zero
DO j=1,nc
dw=ABS(w(i,j))
IF(dw > c)THEN
e=e+one
p=p+LOG(dw)
ENDIF
ENDDO
IF(E == zero)STOP 'W effectively singular in LINVAN'
d1(i)=EXP(-p/e)
ENDDO
CALL muldm(d1,w2,w)
DO j=1,nc
p=zero
e=zero
DO i=1,nc
dw=ABS(w(i,j))
IF(dw > c)THEN
e=e+one
p=p+LOG(dw)
ENDIF
ENDDO
IF(E == zero)STOP 'W effectively singular in LINVAN'
d2(j)=EXP(-p/e)
ENDDO
CALL mulmd(w,d2,w)
c=1.e-8_r_kind ! reset the criterion for "negligible" elements
! revert to iterations of the more efficient method without logarithms:
DO jt=1,2
DO it=1,nit ! perform nit relaxation iterations
DO i=1,nc ! do rows:
p=one
e=zero
DO j=1,nc
dw=ABS(w(i,j))
IF(dw > c)THEN
e=e+one
p=p*dw
ENDIF
ENDDO
p=one/(p**(one/e))
w(i,:) = w(i,:) * p ! rescale this row of w..
d1(i)=d1(i)*p ! ..and update d1 consistently
ENDDO
DO j=1,nc ! do columns:
p=one
e=zero
d2j=d2(j)
DO i=1,nc
dw=ABS(w(i,j))
IF(dw > c)THEN
e=e+one
p=p*dw
ENDIF
ENDDO
p=one/(p**(one/e))
w(:,j) = w(:,j) * p ! rescale this column of w..
d2(j)=d2(j)*p ! ..and update d2 consistently
ENDDO
ENDDO
c=1.e-3_r_kind ! final setting for criterion for "negligible" elements
ENDDO
ab(1:nc) = d1(1:nc) * ab(1:nc) ! rescale r.h.s vector by d1
p=one ! p becomes product of row-lengths:
DO i=1,nc
p=p*SQRT(dot_PRODUCT(w(i,:),w(i,:)))
ENDDO
CALL ldum(w,ipiv,d)
DO i=1,nc
d=d*w(i,i) ! d becomes the determinant of w
ENDDO
wv(:,1) = ab ! convert shape of array
CALL udlmm(w,wv(:,1:1),ipiv)
ab = d2 * wv(:,1) ! rescale solution vector by d2
! ab(1:nc) = d2(1:nc) * ab(1:nc) ! rescale solution vector by d2
! note: it is very likely that round-off errors have accumulated during
! the iterative rescaling of w. we invoke original matrix elements w2 and
! substitute the tentative solution vector into the original (unscaled)
! equation in order to estimate the residual components of roundoff error.
! begin "clean-up" process. substitute solution vector in original
! equation and leave the residual difference in v
v=v-MATMUL(w2,ab)
v = d1 * v ! rescale the residual vector by d1
wv(:,1) = v ! convert shape of array
CALL udlmm(w,wv(:,1:1),ipiv) ! solve linear system with this rhs.
ab=ab+wv(:,1)*d2 ! add residual solution vector,
! scaled, to ab
DEALLOCATE(wv)
w(1,1)=d/p ! this ratio is an indicator of the overall conditioning
! when d/p is very small, treat the results with suspicion!
END SUBROUTINE linvan
SUBROUTINE dlinvan(w,ab)
!$$$ subprogram documentation block
! . . .
! subprogram: dlinvan
! . . .
! prgrmmr: R.J.Purser, NCEP, Washington D.C. 1996
!
! abstract:
! Take square matrix W and seek row and column scalings to produce non-
! vanishing elements of rescaled W having magnitudes as close to unity
! as possible. The approach is make the geometric mean of the nonvanishing
! elements of each row and of each column +1 or -1. Having rescaled the
! matrix and the r.h.s. vector AB, compute the product P of row-vector
! norms, then compute the determinant D and solve the linear system.
! Rescale the solution vector (now AB) and put the conditioning indicator
! formed by the ratio D/P into the first element of W.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block, rm unused vars
!
! input argument list:
! W - Generalized Vandermonde matrix in, conditioning indicator out.
! AB - R.h.s. vector in, solution vector of numerical coefficients out.
!
! output argument list:
! W - Generalized Vandermonde matrix in, conditioning indicator out.
! AB - R.h.s. vector in, solution vector of numerical coefficients out.
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(INOUT) :: w(:,:), ab(:)
INTEGER(i_kind), PARAMETER :: nit=20
REAL(r_kind) :: d1(SIZE(w,1)), d2(SIZE(w,1)), &
w2(SIZE(w,1),SIZE(w,1)),v(SIZE(w,1))
INTEGER(i_kind) :: i, j, it, jt, ipiv(SIZE(w,1)), nc
REAL(r_kind) :: p, e, dw, c, d, d2j
REAL(r_kind),ALLOCATABLE :: wv(:,:) ! work variable for ab(nc) and v(nn)
nc = SIZE(w,DIM=1)
ALLOCATE(wv(nc,1))
w2=w ! Preserve original W and AB for use
v = ab(1:nc) ! in later "clean-up" operation.
d1 = one ! Row scaling factors set to default
d2 = one ! Column scaling factors set to default
C=1.E-16_r_kind ! Set initial criterion for "negligible" elements of W
! In first attempt to estimate row and column scalings, use logarithms
! to avoid the risk of under- or over-flows of the line products of W:
DO i=1,nc
p=zero
e=zero
DO j=1,nc
dw=ABS(w(i,j))
IF(dw > c)THEN
e=e+one
p=p+LOG(dw)
ENDIF
ENDDO
IF(e == zero)STOP 'w effectively singular in linvan'
d1(i)=EXP(-p/e)
ENDDO
CALL muldm_d(d1,w2,w)
DO j=1,nc
p=zero
e=zero
DO i=1,nc
dw=ABS(w(i,j))
IF(dw > c)THEN
e=e+one
p=p+LOG(dw)
ENDIF
ENDDO
IF(e == zero)STOP 'w effectively singular in linvan'
d2(j)=EXP(-p/e)
ENDDO
CALL mulmd_d(w,d2,w)
c=1.e-8_r_kind ! reset the criterion for "negligible" elements
! revert to iterations of the more efficient method without logarithms:
DO jt=1,2
DO it=1,nit ! perform nit relaxation iterations
DO i=1,nc ! do rows:
p=one
e=zero
DO j=1,nc
dw=ABS(w(i,j))
IF(dw > c)THEN
e=e+one
p=p*dw
ENDIF
ENDDO
p=one/(p**(one/e))
w(i,:) = w(i,:) * p ! rescale this row of w..
d1(i)=d1(i)*p ! ..and update d1 consistently
ENDDO
DO j=1,nc ! do columns:
p=one
e=zero
d2j=d2(j)
DO i=1,nc
dw=ABS(w(i,j))
IF(dw > c)THEN
e=e+one
p=p*dw
ENDIF
ENDDO
p=one/(p**(one/e))
w(:,j) = w(:,j) * p ! rescale this column of w..
d2(j)=d2(j)*p ! ..and update d2 consistently
ENDDO
ENDDO
c=1.e-3_r_kind ! final setting for criterion for "negligible" elements
ENDDO
ab(1:nc) = d1(1:nc) * ab(1:nc) ! rescale r.h.s vector by d1
p=one ! p becomes product of row-lengths:
DO i=1,nc
p=p*SQRT(dot_PRODUCT(w(i,:),w(i,:)))
ENDDO
CALL ldum_d(w,ipiv,d)
DO i=1,nc
d=d*w(i,i) ! d becomes the determinant of w
ENDDO
wv(:,1) = ab ! convert shape of array
CALL udlmm_d(w,wv(:,1:1),ipiv)
ab = d2 * wv(:,1) ! rescale solution vector by d2
! ab(1:nc) = d2(1:nc) * ab(1:nc) ! Rescale solution vector by D2
! Note: it is very likely that round-off errors have accumulated during
! the iterative rescaling of W. We invoke original matrix elements W2 and
! substitute the tentative solution vector into the original (unscaled)
! equation in order to estimate the residual components of roundoff error.
! Begin "clean-up" process. Substitute solution vector in original
! equation and leave the residual difference in V
v=v-MATMUL(w2,ab)
v = d1 * v ! Rescale the residual vector by D1
wv(:,1) = v ! Convert shape of array
CALL UDLMM_d(w,wv(:,1:1),ipiv) ! Solve linear system with THIS rhs.
ab=ab+wv(:,1)*d2 ! Add residual solution vector,
! scaled, to AB
DEALLOCATE(wv)
w(1,1)=d/p ! this ratio is an indicator of the overall conditioning
! When D/P is very small, treat the results with suspicion!
END SUBROUTINE dlinvan
SUBROUTINE copdm(d,a)
!$$$ subprogram documentation block
! . . .
! subprogram: copdm
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! d -
!
! output argument list:
! a -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind),DIMENSION(:),INTENT(IN)::d; REAL(r_kind),DIMENSION(:,:),INTENT(OUT)::a; INTEGER(i_kind) i
a=zero; DO i=1,SIZE(a,1); a(i,i)= d(i); ENDDO; RETURN
ENTRY condm(d,a); a=zero; DO i=1,SIZE(a,1); a(i,i)=-d(i); ENDDO
END SUBROUTINE copdm
SUBROUTINE dcopdm(d,a)
!$$$ subprogram documentation block
! . . .
! subprogram: dcopdm
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! d -
!
! output argument list:
! a -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind),DIMENSION(:),INTENT(IN)::d; REAL(r_kind),DIMENSION(:,:),INTENT(OUT)::a
INTEGER(i_kind) i
a=zero; DO i=1,SIZE(a,1); a(i,i)= d(i); ENDDO; RETURN
ENTRY dcondm(d,a); a=zero; DO i=1,SIZE(a,1); a(i,i)=-d(i); ENDDO
END SUBROUTINE dcopdm
SUBROUTINE copsm(s,a)
!$$$ subprogram documentation block
! . . .
! subprogram: copsm
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! s -
!
! output argument list:
! a -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind),INTENT(IN) :: s; REAL(r_kind),DIMENSION(:,:),INTENT(OUT):: a; INTEGER(i_kind) i
a=zero; DO i=1,SIZE(a,1); a(i,i)= s; ENDDO; RETURN
ENTRY consm(s,a); a=zero; DO i=1,SIZE(a,1); a(i,i)=-s; ENDDO
END SUBROUTINE copsm
SUBROUTINE dcopsm(s,a)
!$$$ subprogram documentation block
! . . .
! subprogram: dcopsm
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! s -
!
! output argument list:
! a -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind),INTENT(IN) :: s; REAL(r_kind),DIMENSION(:,:),INTENT(OUT):: a; INTEGER(i_kind) i
a=zero; DO i=1,SIZE(a,1); a(i,i)= s; ENDDO; RETURN
ENTRY dconsm(s,a); a=zero; DO i=1,SIZE(a,1); a(i,i)=-s; ENDDO
END SUBROUTINE dcopsm
SUBROUTINE addmd(a,b,d)
!$$$ subprogram documentation block
! . . .
! subprogram: addmd
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a, b -
! d -
!
! output argument list:
! a, b -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind),DIMENSION(:,:),INTENT(INOUT):: a,b; REAL(r_kind),DIMENSION(:),INTENT(IN):: d
REAL(r_kind) s; INTEGER(i_kind) i
b=a; DO i=1,SIZE(a,1); b(i,i)=b(i,i)+d(i); ENDDO; RETURN
ENTRY submd(a,b,d);b=a; DO i=1,SIZE(a,1); b(i,i)=b(i,i)-d(i); ENDDO; RETURN
ENTRY addms(a,b,s);b=a; DO I=1,SIZE(a,1); b(i,i)=b(i,i)+s; ENDDO; RETURN
ENTRY SUBMS(A,B,S);b=a; DO I=1,SIZE(a,1); B(I,I)=B(I,I)-S; ENDDO;
END SUBROUTINE addmd
SUBROUTINE daddmd(a,b,d)
!$$$ subprogram documentation block
! . . .
! subprogram: daddmd
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a, b -
! d -
!
! output argument list:
! a, b -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
REAL(r_kind),DIMENSION(:,:),INTENT(INOUT)::A,B;REAL(r_kind),DIMENSION(:),INTENT(IN)::D
REAL(r_kind) s; INTEGER(i_kind) i
b=a; DO i=1,SIZE(a,1); b(i,i)=b(i,i)+d(i); ENDDO; RETURN
ENTRY DSUBMD(A,B,D); b=a; DO i=1,SIZE(a,1); b(i,i)=b(i,i)-d(i); ENDDO; RETURN
ENTRY DADDMS(A,B,S); b=a; DO i=1,SIZE(a,1); b(i,i)=b(i,i)+s; ENDDO; RETURN
ENTRY DSUBMS(A,B,S); b=a; DO i=1,SIZE(a,1); b(i,i)=b(i,i)-s; ENDDO;
END SUBROUTINE daddmd
SUBROUTINE l1lm(a,b)
!$$$ subprogram documentation block
! . . .
! subprogram: l1lm
!
! prgrmmr:
!
! abstract: Cholesky, M -> L*U, U(i,j)=L(j,i)
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a -
! b -
!
! output argument list:
! b -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(IN ) :: a(:,:)
REAL(r_kind), INTENT(INOUT) :: b(:,:)
INTEGER(i_kind) :: m,j, jm, jp, i
REAL(r_kind) :: s, bjji
m=SIZE(a,1)
DO j=1,m
jm=j-1
jp=j+1
s = a(j,j) - SUM(b(j,1:jm)*b(j,1:jm))
IF(S <= zero)THEN
PRINT '(" L1LM detects non-positivity at diagonal index",i2)',J
STOP
ENDIF
b(j,j)=SQRT(s)
bjji=one/b(j,j)
DO i=jp,m
s = a(i,j) - SUM(b(i,1:jm)*b(j,1:jm))
b(i,j)=s*bjji
ENDDO
b(1:jm,j) = zero
ENDDO
END SUBROUTINE l1lm
SUBROUTINE DL1LM(A,B)
!$$$ subprogram documentation block
! . . .
! subprogram: dl1lm
!
! prgrmmr:
!
! abstract: Cholesky, M -> L*U, U(i,j)=L(j,i)
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a -
! b -
!
! output argument list:
! b -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(IN ) :: a(:,:)
REAL(r_kind), INTENT(INOUT) :: b(:,:)
INTEGER(i_kind) :: m,j, jm, jp, i
REAL(r_kind) :: s, bjji
m=SIZE(a,1)
DO j=1,m
jm=j-1
jp=j+1
s = a(j,j) - SUM(b(j,1:jm)*b(j,1:jm))
IF(s <= zero)THEN
PRINT '(" L1LM detects non-positivity at diagonal index",i2)',J
STOP
ENDIF
b(j,j)=SQRT(s)
bjji=one/b(j,j)
DO i=jp,m
s = a(i,j) - SUM(b(i,1:jm)*b(j,1:jm))
b(i,j)=s*bjji
ENDDO
b(1:jm,j) = zero
ENDDO
RETURN
END SUBROUTINE dl1lm
SUBROUTINE ldlm(a,b,d) ! Modified Cholesky decompose Q --> L*D*U, U(i,j)=L(j,i)
!$$$ subprogram documentation block
! . . .
! subprogram: ldlm
!
! prgrmmr:
!
! abstract: Modified Cholesky decompose Q --> L*D*U, U(i,j)=L(j,i)
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block, rm unused vars
!
! input argument list:
! a -
! b -
!
! output argument list:
! b -
! d -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(IN ) :: a(:,:)
REAL(r_kind), INTENT(INOUT) :: b(:,:)
REAL(r_kind), INTENT( OUT) :: d(:)
INTEGER(i_kind) :: m,j, jm, jp, i
REAL(r_kind) :: bjji
m=SIZE(a,1)
DO j=1,m
jm=j-1
jp=j+1
d(j)=a(j,j) - SUM(b(1:jm,j)*b(j,1:jm))
b(j,j) = one
IF(d(j) == zero)THEN
PRINT '(" LDLM detects singularity at diagonal index",i2)',J
STOP
ENDIF
bjji=one/d(j)
DO i=jp,m
b(j,i)= a(i,j) - dot_PRODUCT(b(1:jm,j),b(i,1:jm))
b(i,j)=b(j,i)*bjji
ENDDO
ENDDO
CALL zeru(b)
RETURN
END SUBROUTINE ldlm
SUBROUTINE dldlm(a,b,d) ! Modified Cholesky Q --> L*D*U, U(i,j)=L(j,i)
!$$$ subprogram documentation block
! . . .
! subprogram: dldlm
!
! prgrmmr:
!
! abstract: Modified Cholesky Q --> L*D*U, U(i,j)=L(j,i)
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block, rm unused vars
!
! input argument list:
! a -
! b -
!
! output argument list:
! b -
! d -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(IN ) :: a(:,:)
REAL(r_kind), INTENT(INOUT) :: b(:,:)
REAL(r_kind), INTENT( OUT) :: d(:)
INTEGER(i_kind) :: m,j, jm, jp, i
REAL(r_kind) :: bjji
m=SIZE(a,1)
DO j=1,m; jm=j-1; jp=j+1
d(j)=a(j,j) - SUM(b(1:jm,j)*b(j,1:jm))
b(j,j) = one
IF(d(j) == zero)THEN
PRINT '(" DLDLM detects singularity at diagonal index",i2)',J
STOP
ENDIF
bjji=one/d(j)
DO i=jp,m
b(j,i)= a(i,j) - dot_PRODUCT(b(1:jm,j),b(i,1:jm))
b(i,j)=b(j,i)*bjji
ENDDO
ENDDO
CALL zeru_d(b)
RETURN
END SUBROUTINE dldlm
SUBROUTINE invh(a)
!$$$ subprogram documentation block
! . . .
! subprogram: invh
!
! prgrmmr: R.J.Purser, National Meteorological Center, Washington D.C. 1993
!
! abstract: Inver,t in place, a symmetric matrix
!
! limitation: This routine incorporates no pivoting - it is intended for matrices
! that are already diagonally dominant
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! A - symmetric square matrix, output as inverse of input
!
! output argument list:
! A - symmetric square matrix, output as inverse of input
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(INOUT) :: a(:,:)
INTEGER(i_kind) :: m,k, kp, i, ip, j
REAL(r_kind),DIMENSION(SIZE(a,1)):: d
m=SIZE(a,1)
! PERFORM L.D.U DECOMPOSITION OF THE SYMMETRIC MATRIX:
CALL ldlm(a,a,d)
! INVERT (IN PLACE) THE LOWER TRIANGULAR PART OF A, (ASSUMING UNIT
! DIAGONAL ELEMENTS), AND INVERT THE DIAGONAL PART OF A (ASSUMING
! ZERO OFF-DIAGONAL ELEMENTS). PUT TRANSPOSE OF LOWER, TIMES DIAGONAL,
! INTO UPPER PART OF A.
DO k=1,m; kp=k+1
a(k,k)=one/d(k)
DO i=kp,m
a(i,k) = a(i,k) + SUM(a(kp:i-1,k)*a(i,kp:i-1)) ! really??
a(i,k) =-a(i,k)
ENDDO
ENDDO
! MULTIPLY: THE TRANSPOSE OF THE LOWER PART OF A (ASSUMING UNIT DIAGS),
! TIMES THE DIAGONAL PART (ASSUMING ZERO OFF-DIAGS), TIMES THE LOWER
! PART. THIS PRODUCT IS THE SYMMETRIC INVERSE OF THE ORIGINAL B.
DO i=2,m
a(1:i-1,i) = a(i,1:i-1) * a(i,i) ! Really?
ENDDO
DO i=1,m
ip=i+1
DO j=1,i-1
a(j,i) = a(j,i) + SUM(a(ip:ip+m-i-1,i)*a(j,ip:ip+m-i-1))
a(i,j) = a(j,i)
ENDDO
a(i,i) = a(i,i) + SUM(a(ip:ip+m-i-1,i)*a(i,ip:ip+m-i-1))
ENDDO
END SUBROUTINE invh
SUBROUTINE dinvh(a)
!$$$ subprogram documentation block
! . . .
! subprogram: dinvh
!
! prgrmmr: R.J.Purser, National Meteorological Center, Washington D.C. 1993
!
! abstract: Inver,t in place, a symmetric matrix
!
! limitation: This routine incorporates no pivoting - it is intended for matrices
! that are already diagonally dominant
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! A - symmetric square matrix, output as inverse of input
!
! output argument list:
! A - symmetric square matrix, output as inverse of input
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(INOUT) :: a(:,:)
INTEGER(i_kind) :: m,k, kp, i, ip, j
REAL(r_kind),DIMENSION(SIZE(a,1)):: d
m=SIZE(a,1)
! PERFORM L.D.U DECOMPOSITION OF THE SYMMETRIC MATRIX:
CALL ldlm_d(a,a,d)
! INVERT (IN PLACE) THE LOWER TRIANGULAR PART OF A, (ASSUMING UNIT
! DIAGONAL ELEMENTS), AND INVERT THE DIAGONAL PART OF A (ASSUMING
! ZERO OFF-DIAGONAL ELEMENTS). PUT TRANSPOSE OF LOWER, TIMES DIAGONAL,
! INTO UPPER PART OF A.
DO k=1,m
kp=k+1
a(k,k)=one/d(k)
DO i=kp,m
a(i,k) = a(i,k) + SUM(a(kp:i-1,k)*a(i,kp:i-1)) ! really??
a(i,k) =-a(i,k)
ENDDO
ENDDO
! MULTIPLY: THE TRANSPOSE OF THE LOWER PART OF A (ASSUMING UNIT DIAGS),
! TIMES THE DIAGONAL PART (ASSUMING ZERO OFF-DIAGS), TIMES THE LOWER
! PART. THIS PRODUCT IS THE SYMMETRIC INVERSE OF THE ORIGINAL B.
DO i=2,m
a(1:i-1,i) = a(i,1:i-1) * a(i,i) ! really?
ENDDO
DO i=1,m
ip=i+1
DO j=1,i-1
a(j,i) = a(j,i) + SUM(a(ip:ip+m-i-1,i)*a(j,ip:ip+m-i-1))
a(i,j) = a(j,i)
ENDDO
a(i,i) = a(i,i) + SUM(a(ip:ip+m-i-1,i)*a(i,ip:ip+m-i-1))
ENDDO
END SUBROUTINE dinvh
SUBROUTINE invl(a)
!$$$ subprogram documentation block
! . . .
! subprogram: invl
!
! prgrmmr: R.J.Purser, National Meteorological Center, Washington D.C. 1994
!
! abstract: Invert lower triangular matrix in place if A are same
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a -
!
! output argument list:
! a -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(INOUT) :: a(:,:)
INTEGER(i_kind) :: m,j, i
REAL(r_kind) :: s
m=SIZE(a,1)
DO j=m,1,-1
a(1:j-1,j) = zero
a(j,j)=one/a(j,j)
DO i=j+1,m
s = SUM(a(j:i-1,j)*a(i,j:i-1))
a(i,j)=-a(i,i)*s
ENDDO
ENDDO
END SUBROUTINE invl
SUBROUTINE dinvl(a)
!$$$ subprogram documentation block
! . . .
! subprogram: dinvl
!
! prgrmmr: R.J.Purser, National Meteorological Center, Washington D.C. 1994
!
! abstract: Invert lower triangular matrix in place if A are same
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a -
!
! output argument list:
! a -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(INOUT) :: a(:,:)
INTEGER(i_kind) :: m,j, i
REAL(r_kind) :: s
m=SIZE(a,1)
DO j=m,1,-1
a(1:j-1,j) = zero
a(j,j)=one/a(j,j)
DO i=j+1,m
s = SUM(a(j:i-1,j)*a(i,j:i-1))
a(i,j)=-a(i,i)*s
ENDDO
ENDDO
END SUBROUTINE dinvl
SUBROUTINE linlv(a,u)
!$$$ subprogram documentation block
! . . .
! subprogram: linvlv
!
! prgrmmr: R.J.Purser, National Meteorological Center, Washington D.C. 1994
!
! abstract: Solve linear system involving lower triangular (LINLV) or upper
! triangular (LINUV) matrix. u is input as right-hand-side, output
! as the solution vector.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a -
! u -
!
! output argument list:
! u -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(IN ):: a(:,:)
REAL(r_kind), INTENT(INOUT):: u(:)
INTEGER(i_kind) :: m,i, j, jp
DO i=1,SIZE(a,1); u(i)=(u(i) - SUM(u(1:i-1)*a(i,1:i-1)))/a(i,i); ENDDO
RETURN
ENTRY linuv(a,u); m=SIZE(a,1)
DO j=m,1,-1; jp=j+1; u(j)=(u(j) - SUM(a(jp:m,j)*u(jp:m))) /a(j,j); ENDDO
END SUBROUTINE linlv
SUBROUTINE dlinlv(a,u)
!$$$ subprogram documentation block
! . . .
! subprogram: dlinvlv
!
! prgrmmr: R.J.Purser, National Meteorological Center, Washington D.C. 1994
!
! abstract: Invert lower triangular matrix in place if A are same
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a -
! u -
!
! output argument list:
! u -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(IN ) :: a(:,:)
REAL(r_kind), INTENT(INOUT) :: u(:)
INTEGER(i_kind) :: m,i, j, jp
DO i=1,SIZE(a,1); u(i)= (u(i) - SUM(u(1:i-1)*a(i,1:i-1)))/a(i,i); ENDDO
RETURN
ENTRY dlinuv(a,u); m=SIZE(a,1)
DO j=m,1,-1; jp=j+1; u(j) = (u(j) - SUM(a(jp:m,j)*u(jp:m)))/a(j,j); ENDDO
END SUBROUTINE dlinlv
SUBROUTINE powp(a,b,n)
!$$$ subprogram documentation block
! . . .
! subprogram: powp
!
! prgrmmr:
!
! abstract: Raise power series A to the power
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a - power series
! n - power to raise to
!
! output argument list:
! b - output power series
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: n ! of N and output as B
REAL(r_kind), INTENT(IN ) :: a(0:)
REAL(r_kind), INTENT( OUT) :: b(0:)
REAL(r_kind),DIMENSION(0:SIZE(a)-1):: t; INTEGER(i_kind) :: k
b(0)=one; b(1:) = zero; DO k=1,n; CALL mulpp(a,b,t); b=t; ENDDO
END SUBROUTINE powp
SUBROUTINE DPOWP(A,B,N) ! Raise power series A to the power
!$$$ subprogram documentation block
! . . .
! subprogram: dpowp
!
! prgrmmr:
!
! abstract: Raise power series A to the power
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a - power series
! n - power to raise to
!
! output argument list:
! b - output power series
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: n ! of N and output as B
REAL(r_kind) , INTENT(IN ) :: a(0:)
REAL(r_kind) , INTENT( OUT) :: b(0:)
REAL(r_kind),DIMENSION(0:SIZE(a)-1):: t; INTEGER(i_kind) :: k
B(0)=one; b(1:) = zero; DO k=1,n; CALL mulpp_d(a,b,t); b=t; ENDDO
END SUBROUTINE dpowp
SUBROUTINE polps(a,s1,s2) ! Apply series A to scalar S1 to obtain S2
!$$$ subprogram documentation block
! . . .
! subprogram: polps
!
! prgrmmr:
!
! abstract: Apply series A to scalar S1 to obtain S2
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a -
! s1 -
!
! output argument list:
! s2 -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind),INTENT(IN ) :: a(0:)
REAL(r_kind),INTENT(IN ) :: s1
REAL(r_kind),INTENT( OUT) :: s2
INTEGER(i_kind) m,k
m=SIZE(a)-1; s2=a(m); DO k=m-1,0,-1; s2=s2*s1+a(k); ENDDO
END SUBROUTINE polps
SUBROUTINE dpolps(a,s1,s2) ! Apply series A to scalar S1 to obtain S2
!$$$ subprogram documentation block
! . . .
! subprogram: dpolps
!
! prgrmmr:
!
! abstract: Apply series A to scalar S1 to obtain S2
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a -
! s1 -
!
! output argument list:
! s2 -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind),INTENT(IN ) :: a(0:)
REAL(r_kind),INTENT(IN ) :: s1
REAL(r_kind),INTENT( OUT) :: s2
INTEGER(i_kind) m,k
m=SIZE(a)-1; s2=a(m); DO k=m-1,0,-1; s2=s2*s1+a(k); ENDDO
END SUBROUTINE dpolps
SUBROUTINE polpp(a,b,c)
!$$$ subprogram documentation block
! . . .
! subprogram: polpp
!
! prgrmmr:
!
! abstract: Apply power series A to power series B and put
! the result out as power-series C.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a,b,c -
!
! output argument list:
! a,b,c -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind),INTENT(INOUT) :: a(0:),b(0:),c(0:)
REAL(r_kind),DIMENSION(0:SIZE(a)-1):: t
INTEGER(i_kind) m,k
m=SIZE(a)-1; c(0)=a(m); c(1:m) = zero
DO k=m-1,0,-1; CALL mulpp(b,c,t); c=t; c(0)=c(0)+a(k); ENDDO
END SUBROUTINE polpp
SUBROUTINE dpolpp(a,b,c)
!$$$ subprogram documentation block
! . . .
! subprogram: dpolpp
!
! prgrmmr:
!
! abstract: Apply power series A to power series B and put
! the result out as power-series C.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a,b,c -
!
! output argument list:
! a,b,c -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind),INTENT(INOUT) :: a(0:),b(0:),c(0:)
REAL(r_kind),DIMENSION(0:SIZE(a)-1):: t
INTEGER(i_kind) m,k
m=SIZE(a)-1
c(0)=a(m); c(1:m) = zero
DO k=m-1,0,-1; CALL mulpp_d(b,c,t); c=t; c(0)=c(0)+a(k); ENDDO
END SUBROUTINE dpolpp
FUNCTION trcm(a) RESULT(trc_res)
!$$$ subprogram documentation block
! . . .
! subprogram: trcm
!
! prgrmmr:
!
! abstract: Trace of square matrix A
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block, rm unused vars
!
! input argument list:
! a -
!
! output argument list:
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(IN ) :: a(:,:)
REAL(r_kind) :: trc_res
INTEGER(i_kind) :: i
trc_res=zero; DO i=1,SIZE(a,1); trc_res=trc_res+a(i,i); ENDDO
END FUNCTION trcm
FUNCTION dtrcm(a) RESULT(trc_res) ! Trace of square matrix A
!$$$ subprogram documentation block
! . . . .
! subprogram: dtrcm
! prgmmr:
!
! abstract:
!
! program history log:
! 2009-08-26 lueken - added subprogram doc block
!
! input argument list:
! a
!
! output argument list:
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
REAL(r_kind), INTENT(IN ) :: a(:,:)
REAL(r_kind) :: trc_res
INTEGER(i_kind) :: i
trc_res=zero; DO i=1,SIZE(a,1); trc_res=trc_res+a(i,i); ENDDO
END FUNCTION dtrcm
SUBROUTINE invmt(a)
!$$$ subprogram documentation block
! . . .
! subprogram: invmt
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a -
!
! output argument list:
! a -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind),DIMENSION(:,:),INTENT(INOUT) :: a
INTEGER(i_kind) m,i,j,jp,l
REAL(r_kind) d
INTEGER(i_kind),DIMENSION(SIZE(a,1)):: ipiv
m=SIZE(a,1)
IF(m /= SIZE(a,2))STOP 'matrix passed to invmt is not square'
! Perform a pivoted L-D-U decomposition on matrix a:
CALL ldum(a,ipiv,d)
! Invert upper triangular portion U in place:
DO i=1,m; a(i,i)=one/a(i,i); ENDDO
DO i=1,m-1
DO j=i+1,m; a(i,j)=-a(j,j)*DOT_PRODUCT(a(i:j-1,j),a(i,i:j-1)); ENDDO
ENDDO
! Invert lower triangular portion L in place:
DO j=1,m-1; jp=j+1
DO i=jp,m; a(i,j)=-a(i,j)-DOT_PRODUCT(a(jp:i-1,j),a(i,jp:i-1)); ENDDO
ENDDO
! Form the product of U**-1 and L**-1 in place
DO j=1,m-1; jp=j+1
DO i=1,j; a(i,j)=a(i,j)+DOT_PRODUCT(a(jp:m,j),a(i,jp:m)); ENDDO
DO i=jp,m; a(i,j)=DOT_PRODUCT(a(i:m,j),a(i,i:m)); ENDDO
ENDDO
! Permute columns according to ipiv
DO j=m-1,1,-1; l=ipiv(j); CALL swpvv(a(:,j),a(:,l)); ENDDO
END SUBROUTINE invmt
SUBROUTINE dinvmt(a)
!$$$ subprogram documentation block
! . . .
! subprogram: dinvmt
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a -
!
! output argument list:
! a -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind),DIMENSION(:,:),INTENT(INOUT) :: a
INTEGER(i_kind) :: m,i,j,jp,l
REAL(r_kind) :: d
INTEGER(i_kind),DIMENSION(SIZE(a,1)) :: ipiv
m=SIZE(a,1)
IF(m /= SIZE(a,2))STOP 'matrix passed to dinvmt is not square'
! Perform a pivoted L-D-U decomposition on matrix a:
CALL ldum_d(a,ipiv,d)
! Invert upper triangular portion U in place:
DO i=1,m; a(i,i)=one/a(i,i); ENDDO
DO i=1,m-1
DO j=i+1,m; a(i,j)=-a(j,j)*DOT_PRODUCT(a(i:j-1,j),a(i,i:j-1)); ENDDO
ENDDO
! Invert lower triangular portion L in place:
DO j=1,m-1; jp=j+1
DO i=jp,m; a(i,j)=-a(i,j)-DOT_PRODUCT(a(jp:i-1,j),a(i,jp:i-1)); ENDDO
ENDDO
! Form the product of U**-1 and L**-1 in place
DO j=1,m-1; jp=j+1
DO i=1,j; a(i,j)=a(i,j)+DOT_PRODUCT(a(jp:m,j),a(i,jp:m)); ENDDO
DO i=jp,m; a(i,j)=DOT_PRODUCT(a(i:m,j),a(i,i:m)); ENDDO
ENDDO
! Permute columns according to ipiv
DO j=m-1,1,-1; l=ipiv(j); CALL swpvv_d(a(:,j),a(:,l)); ENDDO
END SUBROUTINE dinvmt
SUBROUTINE linmmt(a,b)
!$$$ subprogram documentation block
! . . .
! subprogram: linmmt
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a,b -
!
! output argument list:
! a,b -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind),DIMENSION(:,:),INTENT(INOUT) :: a,b
INTEGER(i_kind),DIMENSION(SIZE(a,1)) :: ipiv
INTEGER(i_kind) :: m
REAL(r_kind) :: d
m=SIZE(a,1)
IF(m /= SIZE(a,2))STOP 'matrix passed to linmmt is not square'
IF(m /= SIZE(b,1))STOP 'matrix and vectors in linmmt have unmatched sizes'
CALL ldum(a,ipiv,d); CALL udlmm(a,b,ipiv)
END SUBROUTINE linmmt
SUBROUTINE dlinmmt(a,b)
!$$$ subprogram documentation block
! . . .
! subprogram: dlinmmt
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a,b -
!
! output argument list:
! a,b -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind),DIMENSION(:,:),INTENT(INOUT) :: a,b
INTEGER(i_kind),DIMENSION(SIZE(a,1)) :: ipiv
INTEGER(i_kind) :: m
REAL(r_kind) :: d
m=SIZE(a,1)
IF(m /= SIZE(a,2))STOP 'matrix passed to linmmt_d is not square'
IF(m /= SIZE(b,1))STOP 'matrix and vectors in linmmt_d have unmatched sizes'
CALL ldum_d(a,ipiv,d); CALL udlmm_d(a,b,ipiv)
END SUBROUTINE dlinmmt
SUBROUTINE linmvt(a,b)
!$$$ subprogram documentation block
! . . .
! subprogram: linmvt
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a,b -
!
! output argument list:
! a,b -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind),DIMENSION(:,:),INTENT(INOUT) :: a
REAL(r_kind),DIMENSION(:), INTENT(INOUT) :: b
INTEGER(i_kind),DIMENSION(SIZE(a,1)) :: ipiv
INTEGER(i_kind) :: m
REAL(r_kind) :: d
m=SIZE(a,1)
IF(m /= SIZE(a,2))STOP 'matrix passed to linmvt is not square'
IF(m /= SIZE(b))STOP 'matrix and vectors in linmvt have unmatched sizes'
CALL ldum(a,ipiv,d); CALL udlmm(a,b,ipiv)
END SUBROUTINE linmvt
SUBROUTINE dlinmvt(a,b)
!$$$ subprogram documentation block
! . . .
! subprogram: dlinmvt
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a,b -
!
! output argument list:
! a,b -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
REAL(r_kind),DIMENSION(:,:),INTENT(INOUT) :: a
REAL(r_kind),DIMENSION(:), INTENT(INOUT) :: b
INTEGER(i_kind),DIMENSION(SIZE(a,1)) :: ipiv
INTEGER(i_kind) m; REAL(r_kind) d
m=SIZE(a,1)
IF(m /= SIZE(a,2))STOP 'matrix passed to linmvt_d is not square'
IF(m /= SIZE(b))STOP 'matrix and vectors in linmvt_d have unmatched sizes'
CALL ldum_d(a,ipiv,d); CALL udlmm_d(a,b,ipiv)
END SUBROUTINE dlinmvt
end module module_pmat1
MODULE MODULE_pmat2
!$$$ module documentation block
! . . . .
! module: module_pmat2
! prgmmr: purser
!
! abstract:
!
! program history log:
! 1994- - purser
! 2008-04-25 safford - add documentation block
!
! subroutines included:
! avco
! davco
! dfco
! ddfco
! dfco2
! ddfco2
! clib
! dclib
! cad1b
! csb1b
! cad2b
! csb2b
! copbt
! conbt
! copmb
! conmb
! copbm
! conbm
! mulbb
! madbb
! msbbb
! ldub
! dldub
! l1ubb
! dl1ubb
! l1ueb
! dl1ueb
! l1lb
! ldlb
! dldlb
! udub
! dudub
! mulbv
! madbv
! msbbv
! mulbx
! madbx
! msbbx
! mulby
! madby
! msbby
! mulvb
! madvb
! msbvb
! mulxb
! madxb
! msbxb
! mulyb
! madyb
! msbyb
! mulbd
! madbd
! msbbd
! muldb
! maddb
! msbdb
! udlbv
! udlbx
! udlby
! udlvb
! udlxb
! udlyb
! u1lbv
! u1lbx
! u1lby
! u1lvb
! u1lxb
! u1lyb
! linbv
! wrtb
!
! variable definitions:
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
USE MODULE_pmat1
use kinds, only: i_kind
use kinds, only: r_kind => r_double
use constants, only: zero,one,two
IMPLICIT NONE
! set default to private
private
! set subroutines/interfaces to public
public :: avco
public :: avco_d
public :: dfco
public :: dfco_d
public :: dfco2
public :: dfco2_d
public :: clib
public :: clib_d
public :: cad1b
public :: csb1b
public :: cad2b
public :: csb2b
public :: copbt
public :: conbt
public :: copmb
public :: conmb
public :: mulbb
public :: madbb
public :: msbbb
public :: ldub
public :: ldub_d
public :: l1ubb
public :: l1ubb_d
public :: l1ueb
public :: l1ueb_d
public :: l1lb
public :: ldlb
public :: ldlb_d
public :: udub
public :: udub_d
public :: mulbv
public :: madbv
public :: msbbv
public :: mulbx
public :: madbx
public :: msbbx
public :: mulby
public :: madby
public :: msbby
public :: mulvb
public :: madvb
public :: msbvb
public :: mulxb
public :: madxb
public :: msbxb
public :: mulyb
public :: madyb
public :: msbyb
public :: mulbd
public :: madbd
public :: msbbd
public :: muldb
public :: maddb
public :: msbdb
public :: udlbv
public :: udlbx
public :: udlby
public :: udlvb
public :: udlxb
public :: udlyb
public :: u1lbv
public :: u1lbx
public :: u1lby
public :: u1lvb
public :: u1lxb
public :: u1lyb
public :: linbv
public :: wrtb
INTERFACE avco; MODULE PROCEDURE avco; END INTERFACE
INTERFACE avco_d; MODULE PROCEDURE davco; END INTERFACE
INTERFACE dfco; MODULE PROCEDURE dfco; END INTERFACE
INTERFACE dfco_d; MODULE PROCEDURE ddfco; END INTERFACE
INTERFACE dfco2; MODULE PROCEDURE dfco2; END INTERFACE
INTERFACE dfco2_d;MODULE PROCEDURE ddfco2; END INTERFACE
INTERFACE clib; MODULE PROCEDURE clib; END INTERFACE
INTERFACE clib_d; MODULE PROCEDURE dclib; END INTERFACE
INTERFACE cad1b; MODULE PROCEDURE cad1b; END INTERFACE
INTERFACE csb1b; MODULE PROCEDURE csb1b; END INTERFACE
INTERFACE cad2b; MODULE PROCEDURE cad2b; END INTERFACE
INTERFACE csb2b; MODULE PROCEDURE csb2b; END INTERFACE
INTERFACE copbt; MODULE PROCEDURE copbt; END INTERFACE
INTERFACE conbt; MODULE PROCEDURE conbt; END INTERFACE
INTERFACE copmb; MODULE PROCEDURE copmb; END INTERFACE
INTERFACE conmb; MODULE PROCEDURE conmb; END INTERFACE
INTERFACE copbm; MODULE PROCEDURE copbm; END INTERFACE
INTERFACE conbm; MODULE PROCEDURE conbm; END INTERFACE
INTERFACE mulbb; MODULE PROCEDURE mulbb; END INTERFACE
INTERFACE madbb; MODULE PROCEDURE madbb; END INTERFACE
INTERFACE msbbb; MODULE PROCEDURE msbbb; END INTERFACE
INTERFACE ldub; MODULE PROCEDURE ldub; END INTERFACE
INTERFACE ldub_d; MODULE PROCEDURE dldub; END INTERFACE
INTERFACE l1ubb; MODULE PROCEDURE l1ubb; END INTERFACE
INTERFACE l1ubb_d;MODULE PROCEDURE dl1ubb; END INTERFACE
INTERFACE l1ueb; MODULE PROCEDURE l1ueb; END INTERFACE
INTERFACE l1ueb_d;MODULE PROCEDURE dl1ueb; END INTERFACE
INTERFACE l1lb; MODULE PROCEDURE l1lb; END INTERFACE
INTERFACE ldlb; MODULE PROCEDURE ldlb; END INTERFACE
INTERFACE ldlb_d; MODULE PROCEDURE dldlb; END INTERFACE
INTERFACE udub; MODULE PROCEDURE udub; END INTERFACE
INTERFACE udub_d; MODULE PROCEDURE dudub; END INTERFACE
INTERFACE mulbv; MODULE PROCEDURE mulbv; END INTERFACE
INTERFACE madbv; MODULE PROCEDURE madbv; END INTERFACE
INTERFACE msbbv; MODULE PROCEDURE msbbv; END INTERFACE
INTERFACE mulbx; MODULE PROCEDURE mulbx; END INTERFACE
INTERFACE madbx; MODULE PROCEDURE madbx; END INTERFACE
INTERFACE msbbx; MODULE PROCEDURE msbbx; END INTERFACE
INTERFACE mulby; MODULE PROCEDURE mulby; END INTERFACE
INTERFACE madby; MODULE PROCEDURE madby; END INTERFACE
INTERFACE msbby; MODULE PROCEDURE msbby; END INTERFACE
INTERFACE mulvb; MODULE PROCEDURE mulvb; END INTERFACE
INTERFACE madvb; MODULE PROCEDURE madvb; END INTERFACE
INTERFACE msbvb; MODULE PROCEDURE msbvb; END INTERFACE
INTERFACE mulxb; MODULE PROCEDURE mulxb; END INTERFACE
INTERFACE madxb; MODULE PROCEDURE madxb; END INTERFACE
INTERFACE msbxb; MODULE PROCEDURE msbxb; END INTERFACE
INTERFACE mulyb; MODULE PROCEDURE mulyb; END INTERFACE
INTERFACE madyb; MODULE PROCEDURE madyb; END INTERFACE
INTERFACE msbyb; MODULE PROCEDURE msbyb; END INTERFACE
INTERFACE mulbd; MODULE PROCEDURE mulbd; END INTERFACE
INTERFACE madbd; MODULE PROCEDURE madbd; END INTERFACE
INTERFACE msbbd; MODULE PROCEDURE msbbd; END INTERFACE
INTERFACE muldb; MODULE PROCEDURE muldb; END INTERFACE
INTERFACE maddb; MODULE PROCEDURE maddb; END INTERFACE
INTERFACE msbdb; MODULE PROCEDURE msbdb; END INTERFACE
INTERFACE udlbv; MODULE PROCEDURE udlbv; END INTERFACE
INTERFACE udlbx; MODULE PROCEDURE udlbx; END INTERFACE
INTERFACE udlby; MODULE PROCEDURE udlby; END INTERFACE
INTERFACE udlvb; MODULE PROCEDURE udlvb; END INTERFACE
INTERFACE udlxb; MODULE PROCEDURE udlxb; END INTERFACE
INTERFACE udlyb; MODULE PROCEDURE udlyb; END INTERFACE
INTERFACE u1lbv; MODULE PROCEDURE u1lbv; END INTERFACE
INTERFACE u1lbx; MODULE PROCEDURE u1lbx; END INTERFACE
INTERFACE u1lby; MODULE PROCEDURE u1lby; END INTERFACE
INTERFACE u1lvb; MODULE PROCEDURE u1lvb; END INTERFACE
INTERFACE u1lxb; MODULE PROCEDURE u1lxb; END INTERFACE
INTERFACE u1lyb; MODULE PROCEDURE u1lyb; END INTERFACE
INTERFACE linbv; MODULE PROCEDURE linbv; END INTERFACE
INTERFACE wrtb; MODULE PROCEDURE wrtb; END INTERFACE
CONTAINS
!=============================================================================
SUBROUTINE davco(na,nb,za,zb,z0,a,b)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: davco
!
! prgrmmr: purser 1999
!
! abstract: Compute one row of the coefficients for the compact mid-interval
! interpolation scheme characterized by matrix equation of the form,
! A.t = B.s (*)
! Where s is the vector of "source" values, t the staggered "target"
! values.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! NA - number of t-points operated on by this row of the A of (*)
! NB - number of s-points operated on by this row of the B of (*)
! ZA - coordinates of t-points used in this row of (*)
! ZB - coordinates of s-points used in this row of (*)
! Z0 - nominal point of application of this row of (*)
!
! output argument list:
! A - the NA coefficients A for this scheme
! B - the NB coefficients B for this scheme
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: na,nb
REAL(r_kind) , INTENT(IN ) :: za(na),zb(nb),z0
REAL(r_kind) , INTENT( OUT) :: a(na),b(nb)
!-----------------------------------------------------------------------------
INTEGER(i_kind) :: na1,nab,i
REAL(r_kind),DIMENSION(na+nb,na+nb):: w
REAL(r_kind),DIMENSION(na) :: za0,pa
REAL(r_kind),DIMENSION(nb) :: zb0,pb
REAL(r_kind),DIMENSION(na+nb) :: ab
!=============================================================================
na1=na+1; nab=na+nb
za0=za-z0 ; zb0=zb-z0
pa=one; pb=-one
w=zero; ab=zero
w(1,1:na)=one; ab(1)=one
DO i=2,nab; w(i,1:na)=pa; pa=pa*za0; w(i,na1:nab)=pb; pb=pb*zb0; ENDDO
CALL inv_d(w,ab)
a=ab(1:na); b=ab(na1:nab)
END SUBROUTINE davco
!=============================================================================
SUBROUTINE avco(na,nb,za,zb,z0,a,b)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: avco
!
! prgrmmr: purser 1999
!
! abstract: Compute one row of the coefficients for the compact mid-interval
! interpolation scheme characterized by matrix equation of the form,
! A.t = B.s (*)
! Where s is the vector of "source" values, t the staggered "target"
! values.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! NA - number of t-points operated on by this row of the A of (*)
! NB - number of s-points operated on by this row of the B of (*)
! ZA - coordinates of t-points used in this row of (*)
! ZB - coordinates of s-points used in this row of (*)
! Z0 - nominal point of application of this row of (*)
!
! output argument list:
! A - the NA coefficients A for this scheme
! B - the NB coefficients B for this scheme
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: na,nb
REAL(r_kind), INTENT(IN ) :: za(na),zb(nb),z0
REAL(r_kind), INTENT( OUT) :: a(na),b(nb)
!-----------------------------------------------------------------------------
INTEGER(i_kind) :: na1,nab,i
REAL(r_kind), DIMENSION(na+nb,na+nb):: w
REAL(r_kind), DIMENSION(na) :: za0,pa
REAL(r_kind), DIMENSION(nb) :: zb0,pb
REAL(r_kind), DIMENSION(na+nb) :: ab
!=============================================================================
na1=na+1; nab=na+nb
za0=za-z0 ; zb0=zb-z0
pa=one; pb=-one
w=zero; ab=zero
w(1,1:na)=one; ab(1)=one
DO i=2,nab; w(i,1:na)=pa; pa=pa*za0; w(i,na1:nab)=pb; pb=pb*zb0; ENDDO
CALL inv(w,ab)
a=ab(1:na); b=ab(na1:nab)
END SUBROUTINE avco
!=============================================================================
SUBROUTINE ddfco(na,nb,za,zb,z0,a,b)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: ddfco
!
! prgrmmr: purser 1999
!
! abstract: Compute one row of the coefficients for either the compact
! differencing or quadrature scheme characterized by matrix
! equation of the form,
! A.d = B.c (*)
! In either case, d is the derivative of c.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! NA - number of d-points operated on by this row of the A of (*)
! NB - number of c-points operated on by this row of the B of (*)
! ZA - coordinates of d-points used in this row of (*)
! ZB - coordinates of c-points used in this row of (*)
! Z0 - nominal point of application of this row of (*)
!
! output argument list:
! A - the A-coefficients for this scheme
! B - the B-coefficients for this scheme
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: na,nb
REAL(r_kind) , INTENT(IN ) :: za(na),zb(nb),z0
REAL(r_kind) , INTENT( OUT) :: a(na),b(nb)
!-----------------------------------------------------------------------------
INTEGER(i_kind) :: na1,nab,i
REAL(r_kind), DIMENSION(na+nb,na+nb):: w
REAL(r_kind), DIMENSION(na) :: za0,pa
REAL(r_kind), DIMENSION(nb) :: zb0,pb
REAL(r_kind), DIMENSION(na+nb) :: ab
!=============================================================================
na1=na+1; nab=na+nb
za0=za-z0 ; zb0=zb-z0
pa=one; pb=-one
w=zero; ab=zero
w(1,1:na)=one; ab(1)=one
DO i=3,nab; w(i,1:na) =pa*(i-2); pa=pa*za0; ENDDO
DO i=2,nab; w(i,na1:nab)=pb; pb=pb*zb0; ENDDO
CALL inv_d(w,ab)
a=ab(1:na); b=ab(na1:nab)
END SUBROUTINE ddfco
!=============================================================================
SUBROUTINE dfco(na,nb,za,zb,z0,a,b)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: dfco
!
! prgrmmr: purser 1999
!
! abstract: Compute one row of the coefficients for either the compact
! differencing or quadrature scheme characterized by matrix
! equation of the form,
! A.d = B.c (*)
! In either case, d is the derivative of c.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! NA - number of d-points operated on by this row of the A of (*)
! NB - number of c-points operated on by this row of the B of (*)
! ZA - coordinates of d-points used in this row of (*)
! ZB - coordinates of c-points used in this row of (*)
! Z0 - nominal point of application of this row of (*)
!
! output argument list:
! A - the A-coefficients for this scheme
! B - the B-coefficients for this scheme
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: na,nb
REAL(r_kind), INTENT(IN ) :: za(na),zb(nb),z0
REAL(r_kind), INTENT( OUT) :: a(na),b(nb)
!-----------------------------------------------------------------------------
INTEGER(i_kind) :: na1,nab,i
REAL(r_kind), DIMENSION(na+nb,na+nb):: w
REAL(r_kind), DIMENSION(na) :: za0,pa
REAL(r_kind), DIMENSION(nb) :: zb0,pb
REAL(r_kind), DIMENSION(na+nb) :: ab
!=============================================================================
na1=na+1; nab=na+nb
za0=za-z0 ; zb0=zb-z0
pa=one; pb=-one
w=zero; ab=zero
w(1,1:na)=one; ab(1)=one
DO i=3,nab; w(i,1:na) =pa*(i-2); pa=pa*za0; ENDDO
DO i=2,nab; w(i,na1:nab)=pb; pb=pb*zb0; ENDDO
CALL inv(w,ab)
a=ab(1:na); b=ab(na1:nab)
END SUBROUTINE dfco
!=============================================================================
SUBROUTINE ddfco2(na,nb,za,zb,z0,a,b)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: ddfco2
!
! prgrmmr: purser 1999
!
! abstract: Compute one row of the coefficients for either the compact second-
! differencing scheme characterized by matrix equation of the form,
! A.d = B.c (*)
! Where d is the second-derivative of c.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! NA - number of d-points operated on by this row of the A of (*)
! NB - number of c-points operated on by this row of the B of (*)
! ZA - coordinates of d-points used in this row of (*)
! ZB - coordinates of c-points used in this row of (*)
! Z0 - nominal point of application of this row of (*)
!
! output argument list:
! A - the NA coefficients A for this scheme
! B - the NB coefficients B for this scheme
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: na,nb
REAL(r_kind) , INTENT(IN ) :: za(na),zb(nb),z0
REAL(r_kind) , INTENT( OUT) :: a(na),b(nb)
!-----------------------------------------------------------------------------
INTEGER(i_kind) :: na1,nab,i
REAL(r_kind), DIMENSION(na+nb,na+nb):: w
REAL(r_kind), DIMENSION(na) :: za0,pa
REAL(r_kind), DIMENSION(nb) :: zb0,pb
REAL(r_kind), DIMENSION(na+nb) :: ab
!=============================================================================
na1=na+1; nab=na+nb
za0=za-z0 ; zb0=zb-z0
pa=one; pb=-one
w=zero; ab=zero
w(1,1:na)=one; ab(1)=one
DO i=4,nab; w(i,1:na) =pa*(i-2)*(i-3); pa=pa*za0; ENDDO
DO i=2,nab; w(i,na1:nab)=pb; pb=pb*zb0; ENDDO
CALL inv_d(w,ab)
a=ab(1:na); b=ab(na1:nab)
END SUBROUTINE ddfco2
!=============================================================================
SUBROUTINE dfco2(na,nb,za,zb,z0,a,b)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: dfco2
!
! prgrmmr: purser 1999
!
! abstract: Compute one row of the coefficients for either the compact second-
! differencing scheme characterized by matrix equation of the form,
! A.d = B.c (*)
! Where d is the second-derivative of c.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! NA - number of d-points operated on by this row of the A of (*)
! NB - number of c-points operated on by this row of the B of (*)
! ZA - coordinates of d-points used in this row of (*)
! ZB - coordinates of c-points used in this row of (*)
! Z0 - nominal point of application of this row of (*)
!
! output argument list:
! A - the NA coefficients A for this scheme
! B - the NB coefficients B for this scheme
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: na,nb
REAL(r_kind), INTENT(IN ) :: za(na),zb(nb),z0
REAL(r_kind), INTENT( OUT) :: a(na),b(nb)
!-----------------------------------------------------------------------------
INTEGER(i_kind) :: na1,nab,i
REAL(r_kind), DIMENSION(na+nb,na+nb):: w
REAL(r_kind), DIMENSION(na) :: za0,pa
REAL(r_kind), DIMENSION(nb) :: zb0,pb
REAL(r_kind), DIMENSION(na+nb) :: ab
!=============================================================================
na1=na+1; nab=na+nb
za0=za-z0 ; zb0=zb-z0
pa=one; pb=-one
w=zero; ab=zero
w(1,1:na)=one; ab(1)=one
DO i=4,nab; w(i,1:na) =pa*(i-2)*(i-3); pa=pa*za0; ENDDO
DO i=2,nab; w(i,na1:nab)=pb; pb=pb*zb0; ENDDO
CALL inv(w,ab)
a=ab(1:na); b=ab(na1:nab)
END SUBROUTINE dfco2
!=============================================================================
SUBROUTINE clib(a,m1,m2,mah1,mah2) ! Clip the dead space of the band matrix, a
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: clib
!
! prgrmmr:
!
! abstract: Clip the dead space of the band matrix
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! m1, m2, mah1, mah2 -
!
! output argument list:
! a -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m1, m2, mah1, mah2
REAL(r_kind) , INTENT(INOUT) :: a(m1,-mah1:mah2)
INTEGER(i_kind) :: j
IF(m2-m1+mah1 < 0)STOP 'In CLIB, form of band matrix implies redundant rows'
DO j=1,mah1; a(1:min(m1,j),-j)=zero; ENDDO; DO j=m2-m1+1,mah2; a(max(1,m2-j+1):m1,j)=zero; ENDDO
END SUBROUTINE clib
!=============================================================================
SUBROUTINE dclib(a,m1,m2,mah1,mah2)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: dclib
!
! prgrmmr:
!
! abstract: Clip the dead space of the band matrix
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! m1, m2, mah1, mah2 -
!
! output argument list:
! a -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m1, m2, mah1, mah2
REAL(r_kind) , INTENT(INOUT) :: a(m1,-mah1:mah2)
INTEGER(i_kind) :: j
IF(m2-m1+mah1 < 0)STOP 'In CLIB_d, form of band matrix implies redundant rows'
DO j=1,mah1; a(1:min(m1,j),-j)=zero; ENDDO; DO j=m2-m1+1,mah2; a(max(1,m2-j+1):m1,j)=zero; ENDDO
END SUBROUTINE dclib
!=============================================================================
SUBROUTINE cad1b(a,m1,mah1,mah2,mirror2)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: cad1b
!
! prgrmmr:
!
! abstract: Incorporate operand symmetry near end-1 of a band matrix operator
!
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! A - Input as unclipped operator, output as symmetrized and clipped.
! m1 - Sizes of implied full matrix
! mah1, mah2 - Left and right semi-bandwidths of A.
! mirror2 - 2*location of symmetry axis relative to end-1 operand element.
!
! output argument list:
! A - Input as unclipped operator, output as symmetrized and clipped.
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m1,mah1,mah2,mirror2
REAL(r_kind), INTENT(INOUT) :: a(0:m1-1,-mah1:mah2)
INTEGER(i_kind) :: i,i2,jm,jp,jpmax
IF(mirror2+mah1 > mah2)STOP 'In cad1b, mah2 insufficient'
DO i=0,m1-1; i2=i*2; jpmax=mirror2+mah1-i2; IF(jpmax <= -mah1)EXIT
DO jm=-mah1,mah2; jp=mirror2-jm-i2; IF(jp <= jm)EXIT
a(i,jp)=a(i,jp)+a(i,jm) ! Reflect and add
a(i,jm)=zero ! zero the exterior part
ENDDO
ENDDO
RETURN
!=============================================================================
ENTRY csb1b(a,m1,mah1,mah2,mirror2)
!=============================================================================
! Like cad1b, but for antisymmetric operand
IF(mirror2+mah1 > mah2)STOP 'In csb1b, mah2 insufficient'
DO i=0,m1-1; i2=i*2; jpmax=mirror2+mah1-i2; IF(jpmax < -mah1)EXIT
DO jm=-mah1,mah2; jp=mirror2-jm-i2; IF(jp < jm)EXIT
a(i,jp)=a(i,jp)-a(i,jm) ! Reflect and subtract
a(i,jm)=zero ! zero the exterior part
ENDDO
ENDDO
END SUBROUTINE cad1b
!=============================================================================
SUBROUTINE cad2b(a,m1,m2,mah1,mah2,mirror2)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: cad2b
!
! prgrmmr:
!
! abstract: Incorporate operand symmetry near end-2 of a band matrix operator
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! A - Input as unclipped operator, output as symmetrized and clipped.
! m1, m2 - Sizes of implied full matrix
! mah1, mah2 - Left and right semi-bandwidths of A.
! mirror2 - 2*location of symmetry axis relative to end-2 operand element.
!
! output argument list:
! A - Input as unclipped operator, output as symmetrized and clipped.
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m1,m2,mah1,mah2,mirror2
REAL(r_kind), INTENT(INOUT) :: a(1-m1:0,m1-m2-mah1:m1-m2+mah2)
INTEGER(i_kind) :: i,i2,jm,jp,jmmin,nah1,nah2,mirror,j0
nah1=mah1+m2-m1; nah2=mah2+m1-m2 ! Effective 2nd-index bounds of A
IF(mirror2-nah1 > -nah2)STOP 'In cad2b, mah1 insufficient'
DO i=0,1-m1,-1; i2=i*2; jmmin=mirror2-nah2-i2; IF(jmmin >= nah2)EXIT
DO jp=nah2,nah1,-1; jm=mirror2-jp-i2; IF(jm >= jp)EXIT
a(i,jm)=a(i,jm)+a(i,jp) ! Reflect and add
a(i,jp)=zero ! zero the exterior part
ENDDO
ENDDO
RETURN
!=============================================================================
ENTRY csb2b(a,m1,m2,mah1,mah2,mirror2)
!=============================================================================
nah1=mah1+m2-m1; nah2=mah2+m1-m2 ! Effective 2nd-index bounds of A
IF(mirror2-nah1 > -nah2)STOP 'In csb2b, mah1 insufficient'
DO i=0,1-m1,-1; i2=i*2; jmmin=mirror2-nah2-i2; IF(jmmin > nah2)EXIT
DO jp=nah2,nah1,-1; jm=mirror2-jp-i2; IF(jm > jp)EXIT
a(i,jm)=a(i,jm)-a(i,jp) ! Reflect and subtract
a(i,jp)=zero ! zero the exterior part
ENDDO
ENDDO
!=============================================================================
ENTRY cex2b(a,m1,m2,mah1,mah2,mirror2)
!=============================================================================
nah1=mah1+m2-m1; nah2=mah2+m1-m2 ! Effective 2nd-index bounds of A
IF(mirror2-nah1 > -nah2)STOP 'In cex2b, mah1 insufficient'
mirror=mirror2/2
IF(mirror*2 /= mirror2)STOP 'In cex2b, mirror2 is not even'
DO i=0,1-m1,-1; i2=i*2; jmmin=mirror2-nah2-i2; IF(jmmin >= nah2)EXIT
j0=mirror-i
DO jp=nah2,nah1,-1; jm=mirror2-jp-i2; IF(jm >= jp)EXIT
a(i,jm)=a(i,jm)-a(i,jp) ! Reflect and subtract
a(i,j0)=a(i,j0)+two*a(i,jp)! Apply double the coefficient to end
a(i,jp)=zero ! zero the exterior part
ENDDO
ENDDO
END SUBROUTINE cad2b
!=============================================================================
SUBROUTINE copbt(a,b,m1,m2,mah1,mah2)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: copbt
!
! prgrmmr: R.J.Purser, National Meteorological Center, Washington D.C. 1994
!
! abstract: Copy transpose of rectangular banded matrix A to B
!
! Note: This routine expects A and B always to occupy separate storage.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! A - input matrix in banded format
! M1 - number of rows of A, columns of B
! M2 - number of columns of A, rows of B
! MAH1 - left-half-bandwidth of A, right-half-bandwidth of B
! MAH2 - right-half-bandwidth of A, left-half-bandwidth of B
!
! output argument list:
! B - output matrix in banded format
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m1, m2, mah1, mah2
REAL(r_kind), INTENT(IN ) :: a(m1,-mah1:mah2)
REAL(r_kind), INTENT( OUT) :: b(m2,-mah2:mah1)
INTEGER(i_kind) :: j, i
CALL clib(b,mah2,mah1,m2,m1)
DO j=-mah1,mah2
DO i=MAX(1,1-j),MIN(m1,m2-j); b(j+i,-j)=a(i,j); ENDDO
ENDDO
RETURN
ENTRY conbt(a,b,m1,m2,mah1,mah2)
CALL clib(b,mah2,mah1,m2,m1)
DO j=-mah1,mah2
DO i=MAX(1,1-j),MIN(m1,m2-j); b(j+i,-j)=-a(i,j); ENDDO
ENDDO
END SUBROUTINE copbt
!=============================================================================
SUBROUTINE copmb(afull,aband,m1,m2,mah1,mah2)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: copmb
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! m1, m2, mah1, mah2 -
! afull -
!
! output argument list:
! aband -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m1, m2, mah1, mah2
REAL(r_kind), DIMENSION(m1,m2), INTENT(IN ) :: afull
REAL(r_kind), DIMENSION(m1,-mah1:mah2),INTENT( OUT) :: aband
INTEGER(i_kind) :: i1,i2, i, j
CALL clib(aband,m1,m2,mah1,mah2)
DO j=1,m1; i1=MAX(1,1-j); i2=MIN(m1,m2-j)
DO i=i1,i2; aband(i,j)= afull(i,j+i); ENDDO
ENDDO
RETURN
!=============================================================================
ENTRY conmb(afull,aband,m1,m2,mah1,mah2)
!=============================================================================
CALL clib(aband,m1,m2,mah1,mah2)
DO j=1,m1; i1=MAX(1,1-j); i2=MIN(m1,m2-j)
DO i=i1,i2; aband(i,j)=-afull(i,j+i); ENDDO
ENDDO
END SUBROUTINE copmb
!=============================================================================
SUBROUTINE copbm(aband,afull,m1,m2,mah1,mah2)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: copbm
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! m1, m2, mah1, mah2 -
! aband -
!
! output argument list:
! afull -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m1, m2, mah1, mah2
REAL(r_kind), DIMENSION(m1,-mah1:mah2),INTENT(IN ) :: aband
REAL(r_kind), DIMENSION(m1,m2), INTENT( OUT) :: afull
INTEGER(i_kind) :: i1,i2, i, j
afull=zero
DO j=1,m1; i1=MAX(1,1-j); i2=MIN(m1,m2-j)
DO i=i1,i2; afull(i,j+i)= aband(i,j); ENDDO
ENDDO
RETURN
!=============================================================================
ENTRY conbm(aband,afull,m1,m2,mah1,mah2)
!=============================================================================
afull=zero
DO j=1,m1; i1=MAX(1,1-j); i2=MIN(m1,m2-j)
DO i=i1,i2; afull(i,j+i)=-aband(i,j); ENDDO
ENDDO
END SUBROUTINE copbm
!=============================================================================
SUBROUTINE mulbb(a,b,c,m1,m2,mah1,mah2,mbh1,mbh2,mch1,mch2)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: mulbb
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! m1, m2, mah1, mah2 -
! mbh1, mbh2, mch1, mch2 -
! a, b -
! c -
!
! output argument list:
! c -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m1, m2, mah1, mah2, mbh1, mbh2, mch1, mch2
REAL(r_kind), INTENT(IN ) :: a(m1,-mah1:mah2), b(m2,-mbh1:mbh2)
REAL(r_kind), INTENT(INOUT) :: c(m1,-mch1:mch2)
INTEGER(i_kind) :: nch1, nch2, j, k, jpk, i1,i2
c=zero
ENTRY madbb(a,b,c,m1,m2,mah1,mah2,mbh1,mbh2,mch1,mch2)
nch1=mah1+mbh1; nch2=mah2+mbh2
IF(nch1 /= mch1 .OR. nch2 /= mch2)STOP 'In MULBB, dimensions inconsistent'
DO j=-mah1,mah2
DO k=-mbh1,mbh2; jpk=j+k; i1=MAX(1,1-j); i2=MIN(m1,m2-j)
c(i1:i2,jpk)=c(i1:i2,jpk)+a(i1:i2,j)*b(j+i1:j+i2,k)
ENDDO
ENDDO
END SUBROUTINE mulbb
!=============================================================================
SUBROUTINE msbbb(a,b,c,m1,m2,mah1,mah2,mbh1,mbh2,mch1,mch2)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: msbbb
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! m1, m2, mah1, mah2 -
! mbh1, mbh2, mch1, mch2 -
! a, b -
!
! output argument list:
! c -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m1, m2, mah1, mah2, mbh1, mbh2, mch1, mch2
REAL(r_kind), INTENT(IN ) :: a(m1,-mah1:mah2), b(m2,-mbh1:mbh2)
REAL(r_kind), INTENT( OUT) :: c(m1,-mch1:mch2)
INTEGER(i_kind) :: nch1, nch2, j, k, jpk, i1,i2
nch1=mah1+mbh1; nch2=mah2+mbh2
IF(nch1 /= mch1 .OR. nch2 /= mch2)STOP 'In MSBBB, dimensions inconsistent'
DO j=-mah1,mah2
DO k=-mbh1,mbh2; jpk=j+k; i1=MAX(1,1-j); i2=MIN(m1,m2-j)
c(i1:i2,jpk)=c(i1:i2,jpk)-a(i1:i2,j)*b(j+i1:j+i2,k)
ENDDO
ENDDO
END SUBROUTINE msbbb
!=============================================================================
SUBROUTINE LDUB(a,m,mah1,mah2)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: ldub
!
! prgrmmr: R.J.Purser, 1994
!
! abstract: Compute [L]*[D**-1]*[U] decomposition of asymmetric band-matrix
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! A - input as the asymmetric band matrix. On output, it contains
! the [L]*[D**-1]*[U] factorization of the input matrix, where
! [L] is lower triangular with unit main diagonal
! [D] is a diagonal matrix
! [U] is upper triangular with unit main diagonal
! M - The number of rows of array A
! MAH1 - The left half-bandwidth of fortran array A
! MAH2 - The right half-bandwidth of fortran array A
!
! output argument list:
! A - input as the asymmetric band matrix. On output, it contains
! the [L]*[D**-1]*[U] factorization of the input matrix, where
! [L] is lower triangular with unit main diagonal
! [D] is a diagonal matrix
! [U] is upper triangular with unit main diagonal
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m,mah1, mah2
REAL(r_kind), INTENT(INOUT) :: a(m,-mah1:mah2)
INTEGER(i_kind) :: j, imost, jmost, jp, i
REAL(r_kind) :: ajj, ajji, aij
DO j=1,m
imost=MIN(m,j+mah1)
jmost=MIN(m,j+mah2)
jp=j+1
ajj=a(j,0)
IF(ajj == zero)THEN
PRINT '(" Failure in LDUB:"/" Matrix requires pivoting or is singular")'
STOP
ENDIF
ajji=one/ajj
a(j,0)=ajji
DO i=jp,imost
aij=ajji*a(i,j-i)
a(i,j-i)=aij
a(i,jp-i:jmost-i)=a(i,jp-i:jmost-i)-aij*a(j,jp-j:jmost-j)
ENDDO
a(j,jp-j:jmost-j)=ajji*a(j,jp-j:jmost-j)
ENDDO
END SUBROUTINE LDUB
!=============================================================================
SUBROUTINE DLDUB(a,m,mah1,mah2)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: dldub
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! m, mah1, mah2 -
! a -
!
! output argument list:
! a -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
INTEGER(i_kind), INTENT(IN ) :: m,mah1, mah2
REAL(r_kind) , INTENT(INOUT) :: a(m,-mah1:mah2)
INTEGER(i_kind) :: j, imost, jmost, jp, i
REAL(r_kind) :: ajj, ajji, aij
DO j=1,m
imost=MIN(m,j+mah1)
jmost=MIN(m,j+mah2)
jp=j+1
ajj=a(j,0)
IF(ajj == zero)THEN
PRINT '(" Fails in LDUB_d:"/" Matrix requires pivoting or is singular")'
STOP
ENDIF
ajji=one/ajj
a(j,0)=ajji
DO i=jp,imost
aij=ajji*a(i,j-i)
a(i,j-i)=aij
a(i,jp-i:jmost-i)=a(i,jp-i:jmost-i)-aij*a(j,jp-j:jmost-j)
ENDDO
a(j,jp-j:jmost-j)=ajji*a(j,jp-j:jmost-j)
ENDDO
END SUBROUTINE DLDUB
!=============================================================================
SUBROUTINE L1UBB(a,b,m,mah1,mah2,mbh1,mbh2)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: l1ubb
!
! prgrmmr: R.J.Purser, 1996
!
! abstract: Form the [L]*[D]*[U] decomposition of asymmetric band-matrix
! [A] replace lower triangular elements of [A] by [D**-1]*[L]*[D],
! the upper by [U], replace matrix [B] by [D**-1]*[B].
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! A - input as band matrix, output as lower and upper triangulars with 1s
! implicitly assumed to lie on the main diagonal. The product of these
! triangular matrices is [D**-1]*[A], where [D] is a diagonal matrix.
! B - in as band matrix, out as same but premultiplied by diagonal [D**-1]
! M - Number of rows of A and B
! MAH1 - left half-width of fortran array A
! MAH2 - right half-width of fortran array A
! MBH1 - left half-width of fortran array B
! MBH2 - right half-width of fortran array B
!
! output argument list:
! A - input as band matrix, output as lower and upper triangulars with 1s
! implicitly assumed to lie on the main diagonal. The product of these
! triangular matrices is [D**-1]*[A], where [D] is a diagonal matrix.
! B - in as band matrix, out as same but premultiplied by diagonal [D**-1]
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m,mah1, mah2, mbh1, mbh2
REAL(r_kind) , INTENT(INOUT) :: a(m,-mah1:mah2), b(m,-mbh1:mbh2)
INTEGER(i_kind) :: j, imost, jmost, jleast, jp, i
REAL(r_kind) :: ajj, ajji, aij
DO j=1,m
imost=MIN(m,j+mah1)
jmost=MIN(m,j+mah2)
jleast=MAX(1,j-mah1)
jp=j+1
ajj=a(j,0)
IF(ajj == zero)STOP 'failure in L1UBB'
ajji=one/ajj
a(j,jleast-j:jmost-j) = ajji * a(j,jleast-j:jmost-j)
DO i=jp,imost
aij=a(i,j-i)
a(i,jp-i:jmost-i) = a(i,jp-i:jmost-i) - aij*a(j,jp-j:jmost-j)
ENDDO
a(j,0)=one
b(j,-mbh1:mbh2) = ajji * b(j,-mbh1:mbh2)
ENDDO
END SUBROUTINE L1UBB
!=============================================================================
SUBROUTINE DL1UBB(a,b,m,mah1,mah2,mbh1,mbh2)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: dl1ubb
!
! prgrmmr: R.J.Purser, 1996
!
! abstract: Form the [L]*[D]*[U] decomposition of asymmetric band-matrix
! [A] replace lower triangular elements of [A] by [D**-1]*[L]*[D],
! the upper by [U], replace matrix [B] by [D**-1]*[B].
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! A - input as band matrix, output as lower and upper triangulars with 1s
! implicitly assumed to lie on the main diagonal. The product of these
! triangular matrices is [D**-1]*[A], where [D] is a diagonal matrix.
! B - in as band matrix, out as same but premultiplied by diagonal [D**-1]
! M - Number of rows of A and B
! MAH1 - left half-width of fortran array A
! MAH2 - right half-width of fortran array A
! MBH1 - left half-width of fortran array B
! MBH2 - right half-width of fortran array B
!
! output argument list:
! A - input as band matrix, output as lower and upper triangulars with 1s
! implicitly assumed to lie on the main diagonal. The product of these
! triangular matrices is [D**-1]*[A], where [D] is a diagonal matrix.
! B - in as band matrix, out as same but premultiplied by diagonal [D**-1]
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: mah1, mah2, mbh1, mbh2
REAL(r_kind) , INTENT(INOUT) :: a(m,-mah1:mah2), b(m,-mbh1:mbh2)
INTEGER(i_kind) :: m,j, imost, jmost, jleast, jp, i
REAL(r_kind) :: ajj, ajji, aij
DO j=1,m
imost=MIN(m,j+mah1)
jmost=MIN(m,j+mah2)
jleast=MAX(1,j-mah1)
jp=j+1
ajj=a(j,0)
IF(ajj == zero)STOP 'failure in DL1UBB'
AJJI=one/AJJ
a(j,jleast-j:jmost-j) = ajji * a(j,jleast-j:jmost-j)
DO I=JP,IMOST
AIJ=A(I,J-I)
a(i,jp-i:jmost-i) = a(i,jp-i:jmost-i) - aij*a(j,jp-j:jmost-j)
ENDDO
A(J,0)=one
b(j,-mbh1:mbh2) = ajji * b(j,-mbh1:mbh2)
ENDDO
END SUBROUTINE DL1UBB
!=============================================================================
SUBROUTINE l1ueb(a,b,m,mah1,mah2,mbh1,mbh2)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: l1ueb
!
! prgrmmr: R.J.Purser, 1998
!
! abstract: Form the [L]*[D]*[U] decomposition of asymmetric band-matrix
! [A] replace all but row zero of the lower triangular
! elements of [A] by [D**-1]*[L]*[D], the upper by [U],
! replace matrix [B] by [D**-1]*[B].
! This is a special adaptation of L1UBB used to process quadarature weights
! for QEDBV etc in which the initial quadrature value is provided as input
! instead of being implicitly assumed zero (which is the case for QZDBV etc).
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! A - input as band matrix, output as lower and upper triangulars with 1s
! implicitly assumed to lie on the main diagonal. The product of these
! triangular matrices is [D**-1]*[A], where [D] is a diagonal matrix.
! B - in as band matrix, out as same but premultiplied by diagonal [D**-1]
! M - number of rows of B, one less than the rows of A (which has "row 0")
! MAH1 - left half-width of fortran array A
! MAH2 - right half-width of fortran array A
! MBH1 - left half-width of fortran array B
! MBH2 - right half-width of fortran array B
!
! output argument list:
! A - input as band matrix, output as lower and upper triangulars with 1s
! implicitly assumed to lie on the main diagonal. The product of these
! triangular matrices is [D**-1]*[A], where [D] is a diagonal matrix.
! B - in as band matrix, out as same but premultiplied by diagonal [D**-1]
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m,mah1, mah2, mbh1, mbh2
REAL(r_kind) , INTENT(INOUT) :: a(0:m,-mah1:mah2), b(m,-mbh1:mbh2)
INTEGER(i_kind) :: j, imost, jmost, jleast, jp, i
REAL(r_kind) :: ajj, ajji, aij
DO j=1,m
imost=MIN(m,j+mah1)
jmost=MIN(m,j+mah2)
jleast=MAX(0,j-mah1)
jp=j+1
ajj=a(j,0)
IF(ajj == zero)STOP 'failure in L1UEB'
ajji=one/ajj
a(j,jleast-j:jmost-j) = ajji * a(j,jleast-j:jmost-j)
DO i=jp,imost
aij=a(i,j-i)
a(i,jp-i:jmost-i) = a(i,jp-i:jmost-i) - aij*a(j,jp-j:jmost-j)
ENDDO
a(j,0)=one
b(j,-mbh1:mbh2) = ajji * b(j,-mbh1:mbh2)
ENDDO
END SUBROUTINE l1ueb
!=============================================================================
SUBROUTINE dl1ueb(a,b,m,mah1,mah2,mbh1,mbh2)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: dl1ueb
!
! prgrmmr: R.J.Purser, 1998
!
! abstract: Form the [L]*[D]*[U] decomposition of asymmetric band-matrix
! [A] replace all but row zero of the lower triangular
! elements of [A] by [D**-1]*[L]*[D], the upper by [U],
! replace matrix [B] by [D**-1]*[B].
! This is a special adaptation of L1UBB used to process quadarature weights
! for QEDBV etc in which the initial quadrature value is provided as input
! instead of being implicitly assumed zero (which is the case for QZDBV etc).
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! A - input as band matrix, output as lower and upper triangulars with 1s
! implicitly assumed to lie on the main diagonal. The product of these
! triangular matrices is [D**-1]*[A], where [D] is a diagonal matrix.
! B - in as band matrix, out as same but premultiplied by diagonal [D**-1]
! M - number of rows of B, one less than the rows of A (which has "row 0")
! MAH1 - left half-width of fortran array A
! MAH2 - right half-width of fortran array A
! MBH1 - left half-width of fortran array B
! MBH2 - right half-width of fortran array B
!
! output argument list:
! A - input as band matrix, output as lower and upper triangulars with 1s
! implicitly assumed to lie on the main diagonal. The product of these
! triangular matrices is [D**-1]*[A], where [D] is a diagonal matrix.
! B - in as band matrix, out as same but premultiplied by diagonal [D**-1]
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m,mah1, mah2, mbh1, mbh2
REAL(r_kind) , INTENT(INOUT) :: a(0:,-mah1:), b(:,-mbh1:)
INTEGER(i_kind) :: j, imost, jmost, jleast, jp, i
REAL(r_kind) :: ajj, ajji, aij
DO j=1,m
imost=MIN(m,j+mah1)
jmost=MIN(m,j+mah2)
jleast=MAX(0,j-mah1)
jp=j+1
ajj=a(j,0)
IF(ajj == zero)STOP 'failure in L1UEB_d'
ajji=one/ajj
a(j,jleast-j:jmost-j) = ajji * a(j,jleast-j:jmost-j)
DO i=jp,imost
aij=a(i,j-i)
a(i,jp-i:jmost-i) = a(i,jp-i:jmost-i) - aij*a(j,jp-j:jmost-j)
ENDDO
a(j,0)=one
b(j,-mbh1:mbh2) = ajji * b(j,-mbh1:mbh2)
ENDDO
END SUBROUTINE dl1ueb
!=============================================================================
SUBROUTINE L1LB(a,b,m,mah) ! Cholesky LU decomposition of Banded.
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: l1lb
!
! prgrmmr:
!
! abstract: Cholesky LU decomposition of Banded.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! A -
! M -
! MAH -
!
! output argument list:
! B -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m, mah
REAL(r_kind), INTENT(IN ) :: a(m,-mah:mah)
REAL(r_kind), INTENT( OUT) :: b(m,-mah:0)
INTEGER(i_kind) :: i, j,jmi
REAL(r_kind) :: s
CALL clib(b,m,m,mah,0)
DO j=1,m
s=a(j,0)-DOT_PRODUCT(b(j,-mah:-1),b(j,-mah:-1))
IF(s <= zero)THEN
PRINT '(" L1LB detects non-positivity at diagonal index",i5)',j
STOP
ENDIF
s=SQRT(s); b(j,0)=s; s=one/s
DO i=j+1,MIN(m,j+mah); jmi=j-i
b(i,jmi)=s*(a(i,jmi)-DOT_PRODUCT(b(i,-mah:jmi-1),b(j,-mah-jmi:-1)))
ENDDO
ENDDO
END SUBROUTINE L1LB
!=============================================================================
SUBROUTINE LDLB(a,b,d,m,mah) ! Modified Cholesky [L(D**-1)U, without sqrt]
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: ldlb
!
! prgrmmr:
!
! abstract: Modified Cholesky [L(D**-1)U, without sqrt]
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a -
! m -
! mah -
!
! output argument list:
! b -
! d -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m, mah
REAL(r_kind), INTENT(IN ) :: a(m,-mah:mah)
REAL(r_kind), INTENT( OUT) :: b(m,-mah:0)
REAL(r_kind), INTENT( OUT) :: d(m)
INTEGER(i_kind) :: i, j,k,jmi,lj,li
REAL(r_kind) :: s,t
CALL clib(b,m,m,mah,0); b(:,0)=one
DO j=1,m; lj=MAX(-mah,1-j)
s=a(j,0)
do k=lj,-1
s=s-b(j,k)**2*d(k+j)
enddo
IF(s <= zero)THEN
PRINT '(" LDLB detects non-positivity at diagonal index",i5)',j
STOP
ENDIF
d(j)=s; s=one/s
DO i=j+1,MIN(m,j+mah); jmi=j-i; li=MAX(-mah,1-i); lj=li-jmi
t=a(i,jmi)
do k=li,jmi-1
t=t-b(i,k)*b(j,k-jmi)*d(i+k)
enddo
b(i,jmi)=s*t
ENDDO
ENDDO
d=one/d
END SUBROUTINE LDLB
!=============================================================================
SUBROUTINE DLDLB(a,b,d,m,mah) ! Modified Cholesky [L(D**-1)U, without sqrt]
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: dl1lb
!
! prgrmmr:
!
! abstract: Modified Cholesky [L(D**-1)U, without sqrt]
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a -
! m -
! mah -
!
! output argument list:
! b -
! d -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m, mah
REAL(r_kind) , INTENT(IN ) :: a(m,-mah:mah)
REAL(r_kind) , INTENT( OUT) :: b(m,-mah:0)
REAL(r_kind) , INTENT( OUT) :: d(m)
INTEGER(i_kind) :: i, j,k,jmi,lj,li
REAL(r_kind) :: s,t
CALL clib_d(b,m,m,mah,0); b(:,0)=one
DO j=1,m; lj=MAX(-mah,1-j)
s=a(j,0)
do k=lj,-1
s=s-b(j,k)**2*d(k+j)
enddo
IF(s <= zero)THEN
PRINT '(" DLDLB detects non-positivity at diagonal index",i5)',j
STOP
ENDIF
d(j)=s; s=one/s
DO i=j+1,MIN(m,j+mah); jmi=j-i;
li=MAX(-mah,1-i);
lj=li-jmi;
t=a(i,jmi)
do k=li,jmi-1
t=t-b(i,k)*b(j,k-jmi)*d(i+k)
enddo
b(i,jmi)=s*t
ENDDO
ENDDO
d=one/d
END SUBROUTINE DLDLB
!=============================================================================
SUBROUTINE UDUB(a,b,d,m,mah) ! Modified reverse Cholesky [U(D**-1)U^t],
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: udub
!
! prgrmmr:
!
! abstract: Modified reverse Cholesky [U(D**-1)U^t],
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a -
! m -
! mah -
!
! output argument list:
! b -
! d -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m, mah
REAL(r_kind), INTENT(IN ) :: a(m,-mah:mah)
REAL(r_kind), INTENT( OUT) :: b(m,0:mah)
REAL(r_kind), INTENT( OUT) :: d(m)
REAL(r_kind), DIMENSION(m,-mah:mah ) :: at
REAL(r_kind), DIMENSION(m,-mah:0) :: bt
REAL(r_kind), DIMENSION(m) :: dt
at=a(m:1:-1,mah:-mah:-1); CALL ldlb(at,bt,dt,m,mah);
b=bt(m:1:-1,0:-mah:-1); d=dt(m:1:-1)
END SUBROUTINE UDUB
!=============================================================================
SUBROUTINE DUDUB(a,b,d,m,mah) ! Modified reverse Cholesky [U(D**-1)U^t],
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: dudub
!
! prgrmmr:
!
! abstract: Modified reverse Cholesky [U(D**-1)U^t],
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! a -
! m -
! mah -
!
! output argument list:
! b -
! d -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m, mah
REAL(r_kind), INTENT(IN ) :: a(m,-mah:mah)
REAL(r_kind), INTENT( OUT) :: b(m,0:mah)
REAL(r_kind), INTENT( OUT) :: d(m)
REAL(r_kind), DIMENSION(m,-mah:mah ) :: at
REAL(r_kind), DIMENSION(m,-mah:0) :: bt
REAL(r_kind), DIMENSION(m) :: dt
at=a(m:1:-1,mah:-mah:-1); CALL ldlb_d(at,bt,dt,m,mah);
b=bt(m:1:-1,0:-mah:-1); d=dt(m:1:-1)
END SUBROUTINE DUDUB
!=============================================================================
SUBROUTINE mulbv(a,v1,v2, m1,m2,mah1,mah2)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: mulbv
!
! prgrmmr: R.J.Purser, 1994
!
! abstract: MULtipication of a Banded matrix times a Vector.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! A - is the matrix
! V1 - the input vector
! M1 - the number of rows assumed for A and for V2
! M2 - the number of columns assumed for A and rows for V1
! MAH1 - the left half-bandwidth of fortran array A
! MAH2 - the right half-bandwidth of fortran array A
!
! output argument list:
! V2 - the output vector
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m1, m2, mah1, mah2
REAL(r_kind), INTENT(IN ) :: a(m1,-mah1:mah2), v1(m2)
REAL(r_kind), INTENT( OUT) :: v2(m1)
INTEGER(i_kind) :: j, i1,i2
v2 = zero
!=============================================================================
ENTRY madbv(a,v1,v2, m1,m2,mah1,mah2)
!=============================================================================
DO j=-mah1,mah2; i1=MAX(1,1-j); i2=MIN(m1,m2-j)
v2(i1:i2) = v2(i1:i2) + a(i1:i2,j)*v1(j+i1:j+i2)
ENDDO
RETURN
!=============================================================================
ENTRY msbbv(a,v1,v2, m1,m2,mah1,mah2)
!=============================================================================
DO j=-mah1,mah2; i1=MAX(1,1-j); i2=MIN(m1,m2-j)
v2(i1:i2) = v2(i1:i2) - a(i1:i2,j)*v1(j+i1:j+i2)
ENDDO
END SUBROUTINE mulbv
!=============================================================================
SUBROUTINE mulbx(a,v1,v2, m1,m2,mah1,mah2,my)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: mulbx
!
! prgrmmr: R.J.Purser, 1994
!
! abstract: MULtipication of a Banded matrix times parallel X-Vectors.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! A - is the matrix
! V1 - the array of input vectors
! M1 - the number of rows assumed for A and for V2
! M2 - the number of columns assumed for A and rows for V1
! MAH1 - the left half-bandwidth of fortran array A
! MAH2 - the right half-bandwidth of fortran array A
! MY - the number of parallel X-vectors
!
! output argument list:
! V2 - the array of output vectors
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m1, m2, mah1, mah2, my
REAL(r_kind), INTENT(IN ) :: a(m1,-mah1:mah2), v1(m2,my)
REAL(r_kind), INTENT( OUT) :: v2(m1,my)
INTEGER(i_kind) :: i,j
v2=zero
!=============================================================================
ENTRY madbx(a,v1,v2, m1,m2,mah1,mah2,my)
!=============================================================================
DO j=-mah1,mah2
DO i=MAX(1,1-j),MIN(m1,m2-j); v2(i,:)=v2(i,:)+a(i,j)*v1(i+j,:); ENDDO
ENDDO
RETURN
!=============================================================================
ENTRY msbbx(a,v1,v2, m1,m2,mah1,mah2,my)
!=============================================================================
DO j=-mah1,mah2
DO i=MAX(1,1-j),MIN(m1,m2-j); v2(i,:)=v2(i,:)-a(i,j)*v1(i+j,:); ENDDO
ENDDO
END SUBROUTINE mulbx
!=============================================================================
SUBROUTINE mulby(a,v1,v2, m1,m2,mah1,mah2,mx)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: mulby
!
! prgrmmr: R.J.Purser, 1994
!
! abstract: MULtipication of a Banded matrix times parallel Y-Vectors.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! A - is the matrix
! V1 - the array of input vectors
! M1 - the number of rows assumed for A and for V2
! M2 - the number of columns assumed for A and rows for V1
! MAH1 - the left half-bandwidth of fortran array A
! MAH2 - the right half-bandwidth of fortran array A
! MX - the length of each of the of parallel Y-vectors
!
! output argument list:
! V2 - the array of output vectors
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m1, m2, mah1, mah2, mx
REAL(r_kind), INTENT(IN ) :: a(m1,-mah1:mah2), v1(mx,m2)
REAL(r_kind), INTENT( OUT) :: v2(mx,m1)
INTEGER(i_kind) :: i,j
v2(1:mx,1:m1) = zero
ENTRY madby(a,v1,v2, m1,m2,mah1,mah2,mx)
DO j=-mah1,mah2
DO i=MAX(1,1-j),MIN(m1,m2-j); v2(:,i)=v2(:,i)+a(i,j)*v1(:,i+j); ENDDO
ENDDO
RETURN
ENTRY msbby(a,v1,v2, m1,m2,mah1,mah2,mx)
DO j=-mah1,mah2
DO i=MAX(1,1-j),MIN(m1,m2-j); v2(:,i)=v2(:,i)-a(i,j)*v1(:,i+j); ENDDO
ENDDO
END SUBROUTINE mulby
!=============================================================================
SUBROUTINE MULVB(v1,a,v2, m1,m2,mah1,mah2)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: mulvb
!
! prgrmmr: R.J.Purser, 1994
!
! abstract: MULtipication of a Vector times a Banded matrix.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! V1 - the input row-vector
! A - is the matrix
! M1 - the number of rows assumed for A and columns for V1
! M2 - the number of columns assumed for A and for V2
! MAH1 - the left half-bandwidth of fortran array A
! MAH2 - the right half-bandwidth of fortran array A
!
! output argument list:
! V2 - the output vector
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m1, m2, mah1, mah2
REAL(r_kind), INTENT(IN ) :: v1(m1), a(m1,-mah1:mah2)
REAL(r_kind), INTENT( OUT) :: v2(m2)
INTEGER(i_kind) :: j, i1,i2
v2=zero
!=============================================================================
ENTRY madvb(v1,a,v2, m1,m2,mah1,mah2)
!=============================================================================
DO j=-mah1,mah2; i1=MAX(1,1-j); i2=MIN(m1,m2-j)
v2(j+i1:j+i2)=v2(j+i1:j+i2)+v1(i1:i2)*a(i1:i2,j)
ENDDO
RETURN
!=============================================================================
ENTRY msbvb(v1,a,v2, m1,m2,mah1,mah2)
!=============================================================================
DO j=-mah1,mah2; i1=MAX(1,1-j); i2=MIN(m1,m2-j)
v2(j+i1:j+i2)=v2(j+i1:j+i2)-v1(i1:i2)*a(i1:i2,j)
ENDDO
END SUBROUTINE mulvb
!=============================================================================
SUBROUTINE mulxb(v1,a,v2, m1,m2,mah1,mah2,my)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: mulxb
!
! prgrmmr: R.J.Purser, 1994
!
! abstract: MULtipication of X-Vectors times Banded matrix.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! V1 - the array of input row-vectors
! A - is the matrix
! M1 - the number of rows assumed for A and columns for V1
! M2 - the number of columns assumed for A and V2
! MAH1 - the left half-bandwidth of fortran array A
! MAH2 - the right half-bandwidth of fortran array A
! MY - the number of parallel X-vectors
!
! output argument list:
! V2 - the array of output vectors
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m1, m2, mah1, mah2, my
REAL(r_kind), INTENT(IN ) :: v1(m1,my), a(m1,-mah1:mah2)
REAL(r_kind), INTENT( OUT) :: v2(m2,my)
INTEGER(i_kind) :: i,j
v2=zero
!=============================================================================
ENTRY madxb(v1,a,v2, m1,m2,mah1,mah2,my)
!=============================================================================
DO j=-mah1,mah2
DO i=MAX(1,1-j),MIN(m1,m2-j); v2(j+i,:)=v2(j+i,:)+v1(i,:)*a(i,j); ENDDO
ENDDO
RETURN
!=============================================================================
ENTRY msbxb(v1,a,v2, m1,m2,mah1,mah2,my)
!=============================================================================
DO j=-mah1,mah2
DO i=MAX(1,1-j),MIN(m1,m2-j); v2(j+i,:)=v2(j+i,:)-v1(i,:)*a(i,j); ENDDO
ENDDO
END SUBROUTINE mulxb
!=============================================================================
SUBROUTINE mulyb(v1,a,v2, m1,m2,mah1,mah2,mx)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: mulyb
!
! prgrmmr: R.J.Purser, 1994
!
! abstract: MULtipication of Y-Vectors times a Banded matrix.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! A - is the matrix
! V1 - the array of input row-vectors
! M1 - the number of rows assumed for A and colums for V1
! M2 - the number of columns assumed for A and V2
! MAH1 - the left half-bandwidth of fortran array A
! MAH2 - the right half-bandwidth of fortran array A
! MX - the length of each of the parallel Y-vectors
!
! output argument list:
! V2 - the array of output vectors
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m1, m2, mah1, mah2, mx
REAL(r_kind), INTENT(IN ) :: v1(mx,m1), a(m1,-mah1:mah2)
REAL(r_kind), INTENT( OUT) :: v2(mx,m2)
INTEGER(i_kind) :: i,j
v2=zero
ENTRY madyb(v1,a,v2, m1,m2,mah1,mah2,mx)
DO j=-mah1,mah2
DO i=MAX(1,1-j),MIN(m1,m2-j); v2(:,j+i)=v2(:,j+i)+v1(:,i)*a(i,j); ENDDO
ENDDO
RETURN
ENTRY msbyb(v1,a,v2, m1,m2,mah1,mah2,mx)
DO j=-mah1,mah2
DO i=MAX(1,1-j),MIN(m1,m2-j); v2(:,j+i)=v2(:,j+i)-v1(:,i)*a(i,j); ENDDO
ENDDO
END SUBROUTINE mulyb
!=============================================================================
SUBROUTINE mulbd(a,d,b,m1,m2,mah1,mah2)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: mulbd
!
! prgrmmr: R.J.Purser, 1994
!
! abstract: MULtipication of a Banded matrix times a Diagonal
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! A - is the input banded-matrix
! D - the diagonal matrix
! M1 - the number of rows assumed for A and for B
! M2 - number of columns assumed for A and B, number of elements of D
! MAH1 - the left half-bandwidth of arrays A and B
! MAH2 - the right half-bandwidth of arrays A and B
!
! output argument list:
! B - the output matrix
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m1, m2, mah1, mah2
REAL(r_kind), INTENT(IN ) :: d(m2)
REAL(r_kind), INTENT(INOUT) :: a(m1,-mah1:mah2),b(m1,-mah1:mah2)
INTEGER(i_kind) :: j, i1,i2
CALL clib(b,m1,m2,mah1,mah2)
DO j=-mah1,mah2; i1=MAX(1,1-j); i2=MIN(m1,m2-j)
b(i1:i2,j)=a(i1:i2,j)*d(j+i1:j+i2)
ENDDO
RETURN
!=============================================================================
ENTRY madbd(a,d,b,m1,m2,mah1,mah2)
!=============================================================================
DO j=-mah1,mah2; i1=MAX(1,1-j); i2=MIN(m1,m2-j)
b(i1:i2,j) = b(i1:i2,j)+a(i1:i2,j)*d(j+i1:j+i2)
ENDDO
RETURN
!=============================================================================
ENTRY msbbd(a,d,b,m1,m2,mah1,mah2)
!=============================================================================
DO j=-mah1,mah2; i1=MAX(1,1-j); i2=MIN(m1,m2-j)
b(i1:i2,j) = b(i1:i2,j)-a(i1:i2,j)*d(j+i1:j+i2)
ENDDO
END SUBROUTINE mulbd
!=============================================================================
SUBROUTINE muldb(d,a,b,m1,m2,mah1,mah2)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: muldb
!
! prgrmmr: R.J.Purser, 1994
!
! abstract: MULtipication of a Banded matrix times a Diagonal
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! D - the diagonal matrix
! A - is the input banded-matrix ! <-> if A and B are actually
! M1 - the number of rows assumed for A and for B
! M2 - number of columns assumed for A and B, number of elements of D
! MAH1 - the left half-bandwidth of arrays A and B
! MAH2 - the right half-bandwidth of arrays A and B
!
! output argument list:
! B - the output matrix ! <-> equivalent arrays.
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m1, m2, mah1, mah2
REAL(r_kind), INTENT(IN ) :: d(m1)
REAL(r_kind), INTENT(INOUT) :: a(m1,-mah1:mah2),b(m1,-mah1:mah2)
INTEGER(i_kind) :: j
CALL clib(b,m1,m2,mah1,mah2)
DO j=-mah1,mah2; b(:,j)=d(:)*a(:,j); ENDDO
END SUBROUTINE muldb
!=============================================================================
SUBROUTINE maddb(d,a,b,m1,mah1,mah2)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: maddb
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! d -
! M1 -
! MAH1 -
! MAH2 -
! a -
! b -
!
! output argument list:
! a -
! b -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m1, mah1, mah2
REAL(r_kind), INTENT(IN ) :: d(m1)
REAL(r_kind), INTENT(INOUT) :: a(m1,-mah1:mah2),b(m1,-mah1:mah2)
INTEGER(i_kind) :: j
DO j=-mah1,mah2; b(:,j)=b(:,j)+d(:)*a(:,j); ENDDO
END SUBROUTINE maddb
!=============================================================================
SUBROUTINE msbdb(d,a,b,m1,mah1,mah2)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: msbdb
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! d -
! M1 -
! MAH1 -
! MAH2 -
! a -
! b -
!
! output argument list:
! a -
! b -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m1, mah1, mah2
REAL(r_kind), INTENT(IN ) :: d(m1)
REAL(r_kind), INTENT(INOUT) :: a(m1,-mah1:mah2),b(m1,-mah1:mah2)
INTEGER(i_kind) :: j
DO j=-mah1,mah2; b(:,j)=b(:,j)-d(:)*a(:,j); ENDDO
END SUBROUTINE msbdb
!=============================================================================
SUBROUTINE udlbv(a,v, m,mah1,mah2)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: udlbv
!
! prgrmmr: R.J.Purser, 1994
!
! abstract: BACk-substitution step of linear inversion involving
! Banded matrix and Vector.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! A - encodes the (L)*(D**-1)*(U) factorization of the linear-system
! matrix, as supplied by subroutine LDUB
! V - input as right-hand-side vector, output as solution vector
! M - the number of rows assumed for A and for V
! MAH1 - the left half-bandwidth of fortran array A
! MAH2 - the right half-bandwidth of fortran array A
!
! output argument list:
! V - input as right-hand-side vector, output as solution vector
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ):: m, mah1, mah2
REAL(r_kind), INTENT(IN ):: a(m,-mah1:mah2)
REAL(r_kind), INTENT(INOUT):: v(m)
INTEGER(i_kind) :: i, j
REAL(r_kind) :: vj
DO j=1,m
vj=v(j)
DO i=j+1,MIN(m,j+mah1); v(i)=v(i)-a(i,j-i)*vj; ENDDO; v(j)=a(j,0)*vj
ENDDO
DO j=m,2,-1
vj=v(j)
DO i=MAX(1,j-mah2),j-1; v(i)=v(i)-a(i,j-i)*vj; ENDDO
ENDDO
END SUBROUTINE udlbv
!=============================================================================
SUBROUTINE udlbx(a,v, mx,mah1,mah2,my)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: udlbx
!
! prgrmmr: R.J.Purser, 1994
!
! abstract: BACk-substitution step of parallel linear inversion involving
! Banded matrix and X-Vectors.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! A - encodes the (L)*(D**-1)*(U) factorization of the linear-system
! matrix, as supplied by subroutine LDUB or, if N=NA, by LDUB
! V - input as right-hand-side vectors, output as solution vectors
! MX - the number of rows assumed for A and length of
! X-vectors stored in V
! MAH1 - the left half-bandwidth of fortran array A
! MAH2 - the right half-bandwidth of fortran array A
! MY - number of parallel X-vectors inverted
!
! output argument list:
! V - input as right-hand-side vectors, output as solution vectors
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: mx, mah1, mah2, my
REAL(r_kind), INTENT(IN ) :: a(mx,-mah1:mah2)
REAL(r_kind), INTENT(INOUT) :: v(mx,my)
INTEGER(i_kind) :: jx, ix
DO jx=1,mx
DO ix=jx+1,MIN(mx,jx+mah1); v(ix,:) = v(ix,:) - a(ix,jx-ix)*v(jx,:); ENDDO
v(jx,:) = a(jx,0) * v(jx,:)
ENDDO
DO jx=mx,2,-1
DO ix=MAX(1,jx-mah2),jx-1; v(ix,:) = v(ix,:) - a(ix,jx-ix)*v(jx,:); ENDDO
ENDDO
END SUBROUTINE udlbx
!=============================================================================
SUBROUTINE udlby(a,v, my,mah1,mah2,mx)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: udlby
!
! prgrmmr: R.J.Purser, 1994
!
! abstract: BACk-substitution step of parallel linear inversion involving
! Banded matrix and Y-Vectors.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! A - encodes the (L)*(D**-1)*(U) factorization of the linear-system
! matrix, as supplied by subroutine LDUB or, if N=NA, by LDUB
! V - input as right-hand-side vectors, output as solution vectors
! MY - the number of rows assumed for A and length of
! Y-vectors stored in V
! MAH1 - the left half-bandwidth of fortran array A
! MAH2 - the right half-bandwidth of fortran array A
! MX - number of parallel Y-vectors inverted
!
! output argument list:
! V - input as right-hand-side vectors, output as solution vectors
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: my, mah1, mah2, mx
REAL(r_kind), INTENT(IN ) :: a(my,-mah1:mah2)
REAL(r_kind), INTENT(INOUT) :: v(mx,my)
INTEGER(i_kind) :: iy, jy
DO jy=1,my
DO iy=jy+1,MIN(my,jy+mah1); v(:,iy) = v(:,iy)-a(iy,jy-iy)*v(:,jy); ENDDO
v(:,jy)=a(jy,0)*v(:,jy)
ENDDO
DO jy=my,2,-1
DO iy=MAX(1,jy-mah2),jy-1; v(:,iy)=v(:,iy)-a(iy,jy-iy)*v(:,jy); ENDDO
ENDDO
END SUBROUTINE udlby
!=============================================================================
SUBROUTINE udlvb(v,a, m,mah1,mah2)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: udlvb
!
! prgrmmr: R.J.Purser, 1994
!
! abstract: BACk-substitution step of linear inversion involving
! row-Vector and Banded matrix.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! A - encodes the (L)*(D**-1)*(U) factorization of the linear-system
! matrix, as supplied by subroutine LDUB
! M - the number of rows assumed for A and columns for V
! MAH1 - the left half-bandwidth of fortran array A
! MAH2 - the right half-bandwidth of fortran array A
!
! output argument list:
! V - input as right-hand-side row-vector, output as solution vector
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m, mah1, mah2
REAL(r_kind), INTENT(IN ) :: a(m,-mah1:mah2)
REAL(r_kind), INTENT(INOUT) :: v(m)
INTEGER(i_kind) :: i, j
REAL(r_kind) :: vi
DO i=1,m
vi=v(i)
DO j=i+1,MIN(m,i+mah2); v(j)=v(j)-vi*a(i,j-i); ENDDO
v(i)=vi*a(i,0)
ENDDO
DO i=m,2,-1
vi=v(i)
DO j=MAX(1,i-mah1),i-1; v(j)=v(j)-vi*a(i,j-i); ENDDO
ENDDO
END SUBROUTINE udlvb
!=============================================================================
SUBROUTINE udlxb(v,a, mx,mah1,mah2,my)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: udlxb
!
! prgrmmr: R.J.Purser, National Meteorological Center, 1994
!
! abstract: BACk-substitution step of parallel linear inversion involving
! Banded matrix and row-X-Vectors.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! V - input as right-hand-side vectors, output as solution vectors
! A - encodes the (L)*(D**-1)*(U) factorization of the linear-system
! matrix, as supplied by subroutine LDUB
! MX - the number of rows assumed for A and length of
! X-vectors stored in V
! MAH1 - the left half-bandwidth of fortran array A
! MAH2 - the right half-bandwidth of fortran array A
! MY - number of parallel X-vectors inverted
!
! output argument list:
! V - input as right-hand-side vectors, output as solution vectors
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: mx, mah1, mah2, my
REAL(r_kind), INTENT(IN ) :: a(mx,-mah1:mah2)
REAL(r_kind), INTENT(INOUT) :: v(mx,my)
INTEGER(i_kind) :: ix, jx
DO ix=1,mx
DO jx=ix+1,MIN(mx,ix+mah2); v(jx,:)=v(jx,:)-v(ix,:)*a(ix,jx-ix); ENDDO
v(ix,:)=v(ix,:)*a(ix,0)
ENDDO
DO ix=mx,2,-1
DO jx=MAX(1,ix-mah1),ix-1; v(jx,:)=v(jx,:)-v(ix,:)*a(ix,jx-ix); ENDDO
ENDDO
END SUBROUTINE udlxb
!=============================================================================
SUBROUTINE udlyb(v,a, my,mah1,mah2,mx)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: udlyb
!
! prgrmmr: R.J.Purser, National Meteorological Center, 1994
!
! abstract: BACk-substitution step of parallel linear inversion involving
! Banded matrix and row-Y-Vectors.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! V - input as right-hand-side vectors, output as solution vectors
! A - encodes the (L)*(D**-1)*(U) factorization of the linear-system
! matrix, as supplied by subroutine LDUB
! MY - the number of rows assumed for A and length of
! Y-vectors stored in V
! MAH1 - the left half-bandwidth of fortran array A
! MAH2 - the right half-bandwidth of fortran array A
! MX - number of parallel Y-vectors inverted
!
! output argument list:
! V - input as right-hand-side vectors, output as solution vectors
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: my, mah1, mah2, mx
REAL(r_kind), INTENT(IN ) :: a(my,-mah1:mah2)
REAL(r_kind), INTENT(INOUT) :: v(mx,my)
INTEGER(i_kind) :: iy, jy
DO iy=1,my
DO jy=iy+1,MIN(my,iy+mah2); v(:,jy)=v(:,jy)-v(:,iy)*a(iy,jy-iy); ENDDO
v(:,iy)=v(:,iy)*a(iy,0)
ENDDO
DO iy=my,2,-1
DO jy=MAX(1,iy-mah1),iy-1; v(:,jy)=v(:,jy)-v(:,iy)*a(iy,jy-iy); ENDDO
ENDDO
END SUBROUTINE udlyb
!=============================================================================
SUBROUTINE u1lbv(a,v, m,mah1,mah2)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: u1lbv
!
! prgrmmr: R.J.Purser, National Meteorological Center, 1996
!
! abstract: BACk-substitution step ((U**-1)*(L**-1)) of linear inversion
! involving special Banded matrix and right-Vector.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! A - encodes the [L]*[U] factorization of the linear-system
! matrix, as supplied by subroutine L1UBB
! V - input as right-hand-side vector, output as solution vector
! M - the number of rows assumed for A and for V
! MAH1 - the left half-bandwidth of fortran array A
! MAH2 - the right half-bandwidth of fortran array A
!
! output argument list:
! V - input as right-hand-side vector, output as solution vector
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m, mah1, mah2
REAL(r_kind), INTENT(IN ) :: a(m,-mah1:mah2)
REAL(r_kind), INTENT(INOUT) :: v(m)
INTEGER(i_kind) :: i, j
REAL(r_kind) :: vj
DO j=1,m
vj=v(j)
DO i=j+1,MIN(m,j+mah1); v(i)=v(i)-a(i,j-i)*vj; ENDDO
ENDDO
DO j=m,2,-1
vj=v(j)
DO i=MAX(1,j-mah2),j-1; v(i)=v(i)-a(i,j-i)*vj; ENDDO
ENDDO
END SUBROUTINE u1lbv
!=============================================================================
SUBROUTINE u1lbx(a,v, mx,mah1,mah2,my)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: u1lbx
!
! prgrmmr: R.J.Purser, National Meteorological Center, 1996
!
! abstract: Special BaCk-substitution step of parallel linear inversion
! involving Banded matrix and X-right-Vectors.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! A - encodes the [L]*[U] factorization of the linear-system
! matrix, as supplied by subroutine L1UBB
! V - input as right-hand-side vectors, output as solution vectors
! MX - the number of rows assumed for A and length of
! X-vectors stored in V
! MAH1 - the left half-bandwidth of fortran array A
! MAH2 - the right half-bandwidth of fortran array A
! MY - number of parallel X-vectors inverted
!
! output argument list:
! V - input as right-hand-side vectors, output as solution vectors
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: mx, mah1, mah2, my
REAL(r_kind), INTENT(IN ) :: a(mx,-mah1:mah2)
REAL(r_kind), INTENT(INOUT) :: v(mx,my)
INTEGER(i_kind) :: ix, jx
DO jx=1,mx
DO ix=jx+1,MIN(mx,jx+mah1); v(ix,:)=v(ix,:)-a(ix,jx-ix)*v(jx,:); ENDDO
ENDDO
DO jx=mx,2,-1
DO ix=MAX(1,jx-mah2),jx-1; v(ix,:)=v(ix,:)-a(ix,jx-ix)*v(jx,:); ENDDO
ENDDO
END SUBROUTINE u1lbx
!=============================================================================
SUBROUTINE u1lby(a,v, my,mah1,mah2,mx)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: u1lby
!
! prgrmmr: R.J.Purser, National Meteorological Center, 1996
!
! abstract: Special BaCk-substitution step of parallel linear inversion
! involving Banded matrix and Y-right-Vectors.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! A - encodes the [L]*[U] factorization of the linear-system
! matrix, as supplied by subroutine L1UBB
! V - input as right-hand-side vectors, output as solution vectors
! MY - the number of rows assumed for A and length of
! Y-vectors stored in V
! MAH1 - the left half-bandwidth of fortran array A
! MAH2 - the right half-bandwidth of fortran array A
! MX - number of parallel Y-vectors inverted
!
! output argument list:
! V - input as right-hand-side vectors, output as solution vectors
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: my, mah1, mah2, mx
REAL(r_kind), INTENT(IN ) :: a(my,-mah1:mah2)
REAL(r_kind), INTENT(INOUT) :: v(mx,my)
INTEGER(i_kind) :: iy, jy
DO jy=1,my
DO iy=jy+1,MIN(my,jy+mah1); v(:,iy)=v(:,iy)-a(iy,jy-iy)*v(:,jy); ENDDO
ENDDO
DO jy=my,2,-1
DO iy=MAX(1,jy-mah2),jy-1; v(:,iy)=v(:,iy)-a(iy,jy-iy)*v(:,jy); ENDDO
ENDDO
END SUBROUTINE u1lby
!=============================================================================
SUBROUTINE u1lvb(v,a, m,mah1,mah2)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: u1lvb
!
! prgrmmr: R.J.Purser, National Meteorological Center, 1996
!
! abstract: Special BaCk-substitution step of linear inversion involving
! left-Vector and Banded matrix.
!
! program history log:
! 2008-04-28 safford -- add subprogram doc block
!
! input argument list:
! V - input as right-hand-side row-vector, output as solution vector
! A - encodes the special [L]*[U] factorization of the linear-system
! matrix, as supplied by subroutine L1UBB
! M - the number of rows assumed for A and columns for V
! MAH1 - the left half-bandwidth of fortran array A
! MAH2 - the right half-bandwidth of fortran array A
!
! output argument list:
! V - input as right-hand-side row-vector, output as solution vector
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m, mah1, mah2
REAL(r_kind), INTENT(IN ) :: a(m,-mah1:mah2)
REAL(r_kind), INTENT(INOUT) :: v(m)
INTEGER(i_kind) :: i, j
REAL(r_kind) :: vi
DO i=1,m
vi=v(i)
DO j=i+1,MIN(m,i+mah2); v(j)=v(j)-vi*a(i,j-i); ENDDO
ENDDO
DO i=m,2,-1
vi=v(i)
DO j=MAX(1,i-mah1),i-1; v(j)=v(j)-vi*a(i,j-i); ENDDO
ENDDO
END SUBROUTINE u1lvb
!=============================================================================
SUBROUTINE u1lxb(v,a, mx,mah1,mah2,my)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: u1lxb
!
! prgrmmr: R.J.Purser, National Meteorological Center, 1996
!
! abstract: Special BaCk-substitution step of parallel linear inversion
! involving Banded matrix and X-left-Vectors.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! V - input as right-hand-side vectors, output as solution vectors
! A - encodes the special [L]*[U] factorization of the linear-system
! matrix, as supplied by subroutine L1UBB
! MX - the number of rows assumed for A and length of
! X-vectors stored in V
! MAH1 - the left half-bandwidth of fortran array A
! MAH2 - the right half-bandwidth of fortran array A
! MY - number of parallel X-vectors inverted
!
! output argument list:
! V - input as right-hand-side vectors, output as solution vectors
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: mx, mah1, mah2, my
REAL(r_kind), INTENT(IN ) :: a(mx,-mah1:mah2)
REAL(r_kind), INTENT(INOUT) :: v(mx,my)
INTEGER(i_kind) :: ix, jx
DO ix=1,mx
DO jx=ix+1,MIN(mx,ix+mah2); v(jx,:)=v(jx,:)-v(ix,:)*a(ix,jx-ix); ENDDO
ENDDO
DO ix=mx,2,-1
DO jx=MAX(1,ix-mah1),ix-1; v(jx,:)=v(jx,:)-v(ix,:)*a(ix,jx-ix); ENDDO
ENDDO
END SUBROUTINE u1lxb
!=============================================================================
SUBROUTINE u1lyb(v,a, my,mah1,mah2,mx)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: u1lyb
!
! prgrmmr: R.J.Purser, National Meteorological Center, 1996
!
! abstract: Special BaCk-substitution step of parallel linear inversion
! involving special Banded matrix and Y-left-Vectors.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! V - input as right-hand-side vectors, output as solution vectors
! A - encodes the [L]*[U] factorization of the linear-system
! matrix, as supplied by subroutine L1UBB
! MY - the number of rows assumed for A and length of
! Y-vectors stored in V
! MAH1 - the left half-bandwidth of fortran array A
! MAH2 - the right half-bandwidth of fortran array A
! MX - number of parallel Y-vectors inverted
!
! output argument list:
! V - input as right-hand-side vectors, output as solution vectors
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: my, mah1, mah2, mx
REAL(r_kind), INTENT(IN ) :: a(my,-mah1:mah2)
REAL(r_kind), INTENT(INOUT) :: v(mx,my)
INTEGER(i_kind) :: iy, jy
DO iy=1,my
DO jy=iy+1,MIN(my,iy+mah2); v(:,jy)=v(:,jy)-v(:,iy)*a(iy,jy-iy); ENDDO
ENDDO
DO iy=my,2,-1
DO jy=MAX(1,iy-mah1),iy-1; v(:,jy)=v(:,jy)-v(:,iy)*a(iy,jy-iy); ENDDO
ENDDO
END SUBROUTINE u1lyb
!=============================================================================
SUBROUTINE linbv(a,v,m,mah1,mah2)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: linbv
!
! prgrmmr: R.J.Purser, National Meteorological Center, 1994
!
! abstract: Solve LINear system with square Banded-matrix and vector V
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! A - system matrix on input, its [L]*[D**-1]*[U] factorization on exit
! V - vector of right-hand-sides on input, solution vector on exit
! M - order of matrix A
! MAH1 - left half-bandwidth of A
! MAH2 - right half-bandwidth of A
!
! output argument list:
! A - system matrix on input, its [L]*[D**-1]*[U] factorization on exit
! V - vector of right-hand-sides on input, solution vector on exit
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m, mah1, mah2
REAL(r_kind), INTENT(INOUT) :: a(m,-mah1:mah2), v(m)
CALL ldub(a,m,mah1,mah2)
CALL udlbv(a,v,m,mah1,mah2)
END SUBROUTINE linbv
!=============================================================================
SUBROUTINE wrtb(a,m1,m2,mah1,mah2)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: wrtb
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! m1 -
! m2 -
! mah1 -
! mah2 -
! a -
!
! output argument list:
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
INTEGER(i_kind), INTENT(IN ) :: m1, m2, mah1, mah2
REAL(r_kind), INTENT(IN ) :: a(m1,-mah1:mah2)
INTEGER(i_kind) :: i1, i2, i, j1, j2, j, nj1
DO i1=1,m1,20
i2=MIN(i1+19,m1)
PRINT '(7x,6(i2,10x))',(j,j=-mah1,mah2)
DO i=i1,i2
j1=MAX(-mah1,1-i)
j2=MIN(mah2,m2-i)
nj1=j1+mah1
IF(nj1==0)PRINT '(1x,i3,6(1x,e12.5))',i,(a(i,j),j=j1,j2)
IF(nj1==1)PRINT '(1x,i3,12x,5(1x,e12.5))',i,(a(i,j),j=j1,j2)
IF(nj1==2)PRINT '(1x,i3,24x,4(1x,e12.5))',i,(a(i,j),j=j1,j2)
IF(nj1==3)PRINT '(1x,i3,36x,3(1x,e12.5))',i,(a(i,j),j=j1,j2)
IF(nj1==4)PRINT '(1x,i3,48x,2(1x,e12.5))',i,(a(i,j),j=j1,j2)
IF(nj1==5)PRINT '(1x,i3,60x,1(1x,e12.5))',i,(a(i,j),j=j1,j2)
ENDDO
READ(*,*)
ENDDO
END SUBROUTINE wrtb
END MODULE MODULE_pmat2
module module_fitcons
!$$$ module documentation block
! . . . .
! module: module_fitcons
! prgmmr: purser
!
! abstract:
!
! program history log:
! 1994- - purser
! 2008-04-28 safford - add stander module documentation block
!
! subroutines included:
! setq
! lagw
! infit
!
! variable definitions:
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
!============================================================================
use kinds, only: i_kind
use kinds, only: r_kind => r_double
use constants, only: zero,one,two,three,five
implicit none
! set default to private
private
! set subroutines to public
public :: setq
public :: lagw
public :: infit
! set passed variables to public
public :: qco,dco,ldsig4,ldsig,nohp,sigc,nohm,no,sigb,noh,ico,wt,dwt
public :: nop,hunit2,nom,nnit,rcrit,q,hunit,dwt1,wt1,q1,hunit1
integer(i_kind),parameter :: noh=3, nohm=noh-1, nohp=noh+1,&
no=noh*2, nom=no-1, nop=no+1, nnit=7
real(r_kind),parameter :: sigc=three, sigb=two
real(r_kind),dimension(no) :: hunit,q,wt,dwt
real(r_kind),dimension(nom) :: hunit1,hunit2,q1,wt1,dwt1
real(r_kind),dimension(-noh:noh) :: qco
real(r_kind),dimension(-1-noh:noh):: ico,dco
real(r_kind) :: rcrit,ldsig,ldsig4
!============================================================================
contains
!============================================================================
SUBROUTINE setq(q,x,n)
!============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: setq
!
! prgrmmr: R.J.Purser, National Meteorological Center, 1994
!
! abstract: Precompute the N constant denominator factors of the
! N-point Lagrange polynomial interpolation formula.
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! X - The N abscissae.
! N - The number of points involved.
!
! output argument list:
! Q - The N denominator constants.
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
IMPLICIT NONE
INTEGER(i_kind), INTENT(in ) :: n
REAL(r_kind),DIMENSION(n),INTENT(in ) :: x
REAL(r_kind),DIMENSION(n),INTENT( out) :: q
!-----------------------------------------------------------------------------
INTEGER(i_kind) :: i,j
!=============================================================================
DO i=1,n
q(i)=one
DO j=1,n
IF(j /= i)q(i)=q(i)/(x(i)-x(j))
ENDDO
ENDDO
END SUBROUTINE setq
!============================================================================
SUBROUTINE lagw(x,xt,q,w,dw,n)
!============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: lagw
!
! prgrmmr: R.J.Purser, National Meteorological Center, 1994
!
! abstract: Construct the Lagrange weights and their derivatives when
! target abscissa is known and denominators Q have already
! been precomputed
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! X - Grid abscissae
! XT - Target abscissa
! Q - Q factors (denominators of the Lagrange weight formula)
! N - Number of grid points involved in the interpolation
!
! output argument list:
! W - Lagrange weights
! DW - Derivatives, dW/dX, of Lagrange weights W
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
IMPLICIT NONE
INTEGER(i_kind), INTENT(in ) :: n
REAL(r_kind), INTENT(in ) :: xt
REAL(r_kind),DIMENSION(n),INTENT(in ) :: x,q
REAL(r_kind),DIMENSION(n),INTENT( out) :: w,dw
!-----------------------------------------------------------------------------
REAL(r_kind),DIMENSION(n) :: sdit,d,di
INTEGER(i_kind) :: i,j
REAL(r_kind) :: p,s,sdil,sdir
!============================================================================
p=one ! ...will become product of all the d(i)=xt-x(i)
DO i=1,n
d(i)=xt-x(i)
p=p*d(i)
ENDDO
! test p to reveal whether any of the d(i) vanish:
IF(p == zero)THEN ! xt coincides with a grid point - use special code:
p=one ! p will become the product of the nonzero d(i),
s=zero ! s will become the corresponding sum of q(i)/d(i)
DO i=1,n
IF(d(i) == zero)THEN
j=i ! identify the grid index corresponding to present xt
w(j)=one ! interpolation weighted entirely to this one.
ELSE
w(i)=zero
p=p*d(i)
dw(i)=q(i)/d(i)
s=s+dw(i)
ENDIF
ENDDO
dw(j)=-s*p
DO i=1,n
IF(i /= j)dw(i)=dw(i)*p
ENDDO
ELSE ! xt is not a grid point - use generic code:
sdil=zero ! will become the sum of terms to the left.
sdir=zero ! will become the sum of terms to the right.
DO i=1,n
di(i)=one/d(i)
sdit(i)=sdil
sdil=sdil+di(i)
w(i)=q(i)*p*di(i)
ENDDO
DO i=n,1,-1
sdit(i)=sdit(i)+sdir
sdir=sdir+di(i)
dw(i)=w(i)*sdit(i)
ENDDO
ENDIF
END SUBROUTINE lagw
!============================================================================
subroutine infit
!============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: infit
!
! prgrmmr: R.J.Purser, National Meteorological Center, 1994
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block, rm unused vars
!
! input argument list:
!
! output argument list:
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
implicit none
integer(i_kind) :: i
real(r_kind) :: divq,divd
!============================================================================
! Initialize quantities that relate to interpolations:
do i=1,no; hunit(i)=i-noh; enddo
hunit1=hunit(:nom) ; hunit2=hunit(2:)
call setq(q,hunit,no) ; call setq(q1,hunit1,nom)
rcrit=SQRT(EPSILON(one))
!------------------------------------
! Initialize coefficients for quadrature, differencing and mdpt interpolation:
divq=967680_r_kind ; divd=1024_r_kind
qco(0)=862564_r_kind/divq ; dco(0)=1225_r_kind/divd ; ico(0)=1225_r_kind/(2*divd)
qco(1)= 57249_r_kind/divq ; dco(1)=-245_r_kind/(3*divd) ; ico(1)=-245_r_kind/(2*divd)
qco(2)= -5058_r_kind/divq ; dco(2)= 49_r_kind/(5*divd) ; ico(2)= 49_r_kind/(2*divd)
qco(3)= 367_r_kind/divq ; dco(3)= -five/(7*divd) ; ico(3)= -five/(2*divd)
qco(-1:-noh:-1) = qco(1:noh) ! complete the stencil of quadrature coeffs.
dco(-1:-nohp:-1) =-dco(0:noh) ! complete the stencil of difference coeffs
ico(-1:-nohp:-1) = ico(0:noh) ! complete the stencil of interpolation coeffs.
!------------------------------------
! Initial coefficients related to control of working grid resolution:
ldsig =log(sigc/sigb)
ldsig4=ldsig**4
end subroutine infit
end module module_fitcons
!=============================================================================
subroutine coefrf(sig,nu,n,m,bnf,lnf)
!=============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: coefrf
!
! prgrmmr: R.J.Purser, NCEP 2001
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! n, m -
! sig, nu -
!
! output argument list:
! bnf -
! lnf -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use module_pmat2
use kinds, only: i_kind
use kinds, only: r_kind => r_double
use constants, only: zero,half,one
implicit none
integer(i_kind), intent(IN ) :: n,m
real(r_kind), dimension(n), intent(IN ) :: sig,nu
real(r_kind), dimension(n), intent( OUT) :: bnf
real(r_kind), dimension(m,n), intent( OUT) :: lnf
!--------------------------------------------------------------------------
integer(i_kind), parameter :: irmax=6
real(r_kind), dimension(n,-m:m) :: s
real(r_kind), dimension(n,-m:0) :: sl
real(r_kind), dimension(n,-m:m,m) :: k,l
real(r_kind), dimension(n) :: eta
real(r_kind), dimension(irmax) :: bcofi,bcofh
integer(i_kind) :: i,i1,il,ir,ik
!--------------------------------------------------------------------------
! The coefficients bcofi are the reciprocals of the i=1 entries of TABLE 1
! of NCEP O.N. 431:
data bcofi/one, 12._r_kind, 90._r_kind, 560._r_kind, 3150._r_kind, 16632._r_kind/
!=============================================================================
bcofh=half/bcofi
do i=1,n
eta(i)=sig(i)*sqrt(nu(i))
enddo
k=zero
!-------------------------------------------------------------------------
! Set k(:, -1:1, 1) to be the K-matrix of (4.8)--(4.10) of NCEP O.N. 431:
!--------------------------------------------------------------------------
do i=1,n-1
k(i , 0,1)=k(i ,0,1)+eta(i+1)/eta(i )
k(i+1, 0,1)=k(i+1,0,1)+eta(i )/eta(i+1)
k(i , 1,1)=-one
k(i+1,-1,1)=-one
enddo
!-------------------------------------------------------------------------
! Set k(:, : , ir) to be the original K-matrix raised to the power of (ir):
!--------------------------------------------------------------------------
do ir=2,m
il=ir-1
call mulbb(k(:,-1:1,1),k(:,-il:il,il),k(:,-ir:ir,ir),n,n,1,1,il,il,ir,ir)
enddo
!-------------------------------------------------------------------------
! Pre- and post-multiply each of the m powers of K by the diagonal matrix,
! sigma, of NCEP O.N. 431, where the elements of sigma measure the smoothing
! scale of the quasi-Gaussian filter in grid-space units.
! Also, multiply each of the resulting banded matrices by .5*b_{1,ir} for
! the appropriate index, ir, corresponding to the power by which the original
! K was raised.
!--------------------------------------------------------------------------
do ir=1,m
call mulbd(k(:,-ir:ir,ir),sig,k(:,-ir:ir,ir),n,n,ir,ir)
call muldb(sig,k(:,-ir:ir,ir),k(:,-ir:ir,ir),n,n,ir,ir)
k(:,-ir:ir,ir)=k(:,-ir:ir,ir)*bcofh(ir)
enddo
s=zero
s(:,0)=one
do ir=1,m
l(:,-ir:ir,ir)=k(:,-ir:ir,ir)
s(:,-ir:ir)=s(:,-ir:ir)+l(:,-ir:ir,ir)
enddo
do i1=2,m
do ir=m,i1,-1
l(:,-ir:ir,ir)=zero
do ik=1,ir-i1+1
il=ir-ik
call madbb(k(:,-ik:ik,ik),l(:,-il:il,il),l(:,-ir:ir,ir), &
n,n,ik,ik,il,il,ir,ir)
enddo
l(:,-ir:ir,ir)=l(:,-ir:ir,ir)/i1
s(:,-ir:ir)=s(:,-ir:ir)+l(:,-ir:ir,ir)
enddo
enddo
call ldlb(s,sl,bnf,n,m)
do i1=1,m
do i=1,n
lnf(i1,i)=sl(i,-i1)
enddo
enddo
end subroutine coefrf
!============================================================================
subroutine ldlb1i(nol,lnf,bnf, &
ims,ime, &
its,ite )
!============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: ldlb1i
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! nol -
! ims, ime -
! its, ite -
! bnf -
! lnf -
!
! output argument list:
! bnf -
! lnf -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: i_kind
use kinds, only: r_kind => r_double
IMPLICIT NONE
INTEGER(i_kind), INTENT(IN ) :: nol
INTEGER(i_kind), INTENT(IN ) :: ims,ime
INTEGER(i_kind), INTENT(IN ) :: its,ite
REAL(r_kind), DIMENSION(ims:ime), &
INTENT(INOUT) :: bnf
REAL(r_kind), DIMENSION(nol, ims:ime), &
INTENT(INOUT) :: lnf
!----------------------------------------------------------------------------
INTEGER(i_kind) :: i,l,m,nola
real(r_kind) :: s
!============================================================================
do i=its,ite
nola=min(nol,i-its)
do l=nola,1,-1
s=lnf(l,i)
do m=l+1,nola
s=s-lnf(m,i)*bnf(i-m)*lnf(m-l,i-l)
enddo
lnf(l,i)=s/bnf(i-l)
enddo
s=bnf(i)
do l=1,nola
s=s-lnf(l,i)**2*bnf(i-l)
enddo
bnf(i)=s
enddo
end subroutine ldlb1i
!============================================================================
subroutine ldlb2i(nol,lnf,bnf, &
ims,ime, jms,jme, &
its,ite, jts,jte )
!============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: ldlb2i
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! nol -
! ims, ime, jms, jme -
! its, ite, jts, jte -
! bnf -
! lnf -
!
! output argument list:
! bnf -
! lnf -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: i_kind
use kinds, only: r_kind => r_double
IMPLICIT NONE
INTEGER(i_kind), INTENT(IN ) :: nol
INTEGER(i_kind), INTENT(IN ) :: ims,ime, jms,jme
INTEGER(i_kind), INTENT(IN ) :: its,ite, jts,jte
REAL(r_kind), DIMENSION(ims:ime, jms:jme), &
INTENT(INOUT) :: bnf
REAL(r_kind), DIMENSION(nol, ims:ime, jms:jme), &
INTENT(INOUT) :: lnf
!----------------------------------------------------------------------------
INTEGER(i_kind) :: i,j,l,m,nola
real(r_kind) :: s
!============================================================================
do j=jts,jte
do i=its,ite
nola=min(nol,i-its)
do l=nola,1,-1
s=lnf(l,i,j)
do m=l+1,nola
s=s-lnf(m,i,j)*bnf(i-m,j)*lnf(m-l,i-l,j)
enddo
lnf(l,i,j)=s/bnf(i-l,j)
enddo
s=bnf(i,j)
do l=1,nola
s=s-lnf(l,i,j)**2*bnf(i-l,j)
enddo
bnf(i,j)=s
enddo
enddo
end subroutine ldlb2i
!============================================================================
subroutine ldlb2j(nol,lnf,bnf, &
ims,ime, jms,jme, &
its,ite, jts,jte )
!============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: ldlb2
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! nol -
! ims,ime, jms,jme -
! its,ite, jts,jte -
! bnf -
! lnf -
!
! output argument list:
! bnf -
! lnf -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: i_kind
use kinds, only: r_kind => r_double
IMPLICIT NONE
INTEGER(i_kind), INTENT(IN ) :: nol
INTEGER(i_kind), INTENT(IN ) :: ims,ime, jms,jme
INTEGER(i_kind), INTENT(IN ) :: its,ite, jts,jte
REAL(r_kind), DIMENSION(ims:ime, jms:jme), &
INTENT(INOUT) :: bnf
REAL(r_kind), DIMENSION(nol, ims:ime, jms:jme), &
INTENT(INOUT) :: lnf
!----------------------------------------------------------------------------
INTEGER(i_kind) :: i,j,l,m,nola
real(r_kind) :: s
!============================================================================
do j=jts,jte
nola=min(nol,j-jts)
do i=its,ite
do l=nola,1,-1
s=lnf(l,i,j)
do m=l+1,nola
s=s-lnf(m,i,j)*bnf(i,j-m)*lnf(m-l,i,j-l)
enddo
lnf(l,i,j)=s/bnf(i,j-l)
enddo
s=bnf(i,j)
do l=1,nola
s=s-lnf(l,i,j)**2*bnf(i,j-l)
enddo
bnf(i,j)=s
enddo
enddo
end subroutine ldlb2j
!============================================================================
subroutine ldlb3i(nol,lnf,bnf, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
!============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: ldlb3i
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! nol -
! ims,ime, jms,jme, kms,kme -
! its,ite, jts,jte, kts,kte -
! bnf -
! lnf -
!
! output argument list:
! bnf -
! lnf -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: i_kind
use kinds, only: r_kind => r_double
IMPLICIT NONE
INTEGER(i_kind), INTENT(IN ) :: nol
INTEGER(i_kind), INTENT(IN ) :: ims,ime, jms,jme, kms,kme
INTEGER(i_kind), INTENT(IN ) :: its,ite, jts,jte, kts,kte
REAL(r_kind), DIMENSION(ims:ime, kms:kme, jms:jme), &
INTENT(INOUT) :: bnf
REAL(r_kind), DIMENSION(nol, ims:ime, kms:kme, jms:jme), &
INTENT(INOUT) :: lnf
!----------------------------------------------------------------------------
INTEGER(i_kind) :: i,j,k,l,m,nola
real(r_kind) :: s
!============================================================================
do j=jts,jte
do k=kts,kte
do i=its,ite
nola=min(nol,i-its)
do l=nola,1,-1
s=lnf(l,i,k,j)
do m=l+1,nola
s=s-lnf(m,i,k,j)*bnf(i-m,k,j)*lnf(m-l,i-l,k,j)
enddo
lnf(l,i,k,j)=s/bnf(i-l,k,j)
enddo
s=bnf(i,k,j)
do l=1,nola
s=s-lnf(l,i,k,j)**2*bnf(i-l,k,j)
enddo
bnf(i,k,j)=s
enddo
enddo
enddo
end subroutine ldlb3i
!============================================================================
subroutine ldlb3j(nol,lnf,bnf, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
!============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: ldlb3j
!
! prgrmmr: R.J.Purser, National Meteorological Center, 1994
!
! abstract:
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! nol -
! ims,ime, jms,jme, kms,kme -
! its,ite, jts,jte, kts,kte -
! bnf -
! lnf -
!
! output argument list:
! bnf -
! lnf -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: i_kind
use kinds, only: r_kind => r_double
IMPLICIT NONE
INTEGER(i_kind), INTENT(IN ) :: nol
INTEGER(i_kind), INTENT(IN ) :: ims,ime, jms,jme, kms,kme
INTEGER(i_kind), INTENT(IN ) :: its,ite, jts,jte, kts,kte
REAL(r_kind), DIMENSION(ims:ime, kms:kme, jms:jme), &
INTENT(INOUT) :: bnf
REAL(r_kind), DIMENSION(nol, ims:ime, kms:kme, jms:jme), &
INTENT(INOUT) :: lnf
!----------------------------------------------------------------------------
INTEGER(i_kind) :: i,j,k,l,m,nola
real(r_kind) :: s
!============================================================================
do j=jts,jte
nola=min(nol,j-jts)
do k=kts,kte
do i=its,ite
do l=nola,1,-1
s=lnf(l,i,k,j)
do m=l+1,nola
s=s-lnf(m,i,k,j)*bnf(i,k,j-m)*lnf(m-l,i,k,j-l)
enddo
lnf(l,i,k,j)=s/bnf(i,k,j-l)
enddo
s=bnf(i,k,j)
do l=1,nola
s=s-lnf(l,i,k,j)**2*bnf(i,k,j-l)
enddo
bnf(i,k,j)=s
enddo
enddo
enddo
end subroutine ldlb3j
SUBROUTINE hbnrf1i(a,nol,lnf,bnf, &
ims,ime, &
its,ite )
!============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: hbnrf1i
!
! prgrmmr:
!
! abstract: Horizontal basic inhomogeneous recursive filter,
! 1-dimensional, active index i
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! nol -
! ims,ime -
! its,ite -
! a -
! bnf -
! lnf -
!
! output argument list:
! a -
! bnf -
! lnf -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: i_kind
use kinds, only: r_kind => r_double
IMPLICIT NONE
INTEGER(i_kind), INTENT(IN ) :: nol
INTEGER(i_kind), INTENT(IN ) :: ims,ime
INTEGER(i_kind), INTENT(IN ) :: its,ite
REAL(r_kind), DIMENSION(ims:ime), &
INTENT(INOUT) :: a
REAL(r_kind), DIMENSION(ims:ime), &
INTENT(IN ) :: bnf
REAL(r_kind), DIMENSION(nol, ims:ime), &
INTENT(IN ) :: lnf
!----------------------------------------------------------------------------
INTEGER(i_kind) :: i,l,nola
!============================================================================
DO i=its+1,ite
nola=MIN(nol,i-its)
DO l=1,nola
a(i)=a(i)-lnf(l,i)*a(i-l)
ENDDO
ENDDO
DO i=its,ite
a(i)=bnf(i)*a(i)
ENDDO
DO i=ite-1,its,-1
nola=MIN(nol,ite-i)
DO l=1,nola
a(i)=a(i)-lnf(l,i+l)*a(i+l)
ENDDO
ENDDO
END SUBROUTINE hbnrf1i
SUBROUTINE hbnrf2i(a,nol,lnf,bnf, &
ims,ime, jms,jme, &
its,ite, jts,jte )
!============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: hbnrf2i
!
! prgrmmr:
!
! abstract: Horizontal basic inhomogeneous recursive filter,
! 2-dimensional, active index i
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! nol -
! ims,ime, jms,jme -
! its,ite, jts,jte -
! a -
! bnf -
! lnf -
!
! output argument list:
! a -
! bnf -
! lnf -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: i_kind
use kinds, only: r_kind => r_double
IMPLICIT NONE
INTEGER(i_kind), INTENT(IN ) :: nol
INTEGER(i_kind), INTENT(IN ) :: ims,ime, jms,jme
INTEGER(i_kind), INTENT(IN ) :: its,ite, jts,jte
REAL(r_kind), DIMENSION(ims:ime, jms:jme), &
INTENT(INOUT) :: a
REAL(r_kind), DIMENSION(ims:ime, jms:jme), &
INTENT(IN ) :: bnf
REAL(r_kind), DIMENSION(nol, ims:ime, jms:jme), &
INTENT(IN ) :: lnf
!----------------------------------------------------------------------------
INTEGER(i_kind) :: i,j,l,nola
!============================================================================
DO j=jts,jte
DO i=its+1,ite
nola=MIN(nol,i-its)
DO l=1,nola
a(i,j)=a(i,j)-lnf(l,i,j)*a(i-l,j)
ENDDO
ENDDO
DO i=its,ite
a(i,j)=bnf(i,j)*a(i,j)
ENDDO
DO i=ite-1,its,-1
nola=MIN(nol,ite-i)
DO l=1,nol
a(i,j)=a(i,j)-lnf(l,i+l,j)*a(i+l,j)
ENDDO
ENDDO
ENDDO
END SUBROUTINE hbnrf2i
SUBROUTINE hbnrf2j(a,nol,lnf,bnf, &
ims,ime, jms,jme, &
its,ite, jts,jte )
!============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: hbnrf1i
!
! prgrmmr:
!
! abstract: Horizontal basic inhomogeneous recursive filter,
! 2-dimensional, active index j
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! nol -
! ims,ime, jms,jme -
! its,ite, jts,jte -
! a -
! bnf -
! lnf -
!
! output argument list:
! a -
! bnf -
! lnf -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: i_kind
use kinds, only: r_kind => r_double
IMPLICIT NONE
INTEGER(i_kind), INTENT(IN ) :: nol
INTEGER(i_kind), INTENT(IN ) :: ims,ime, jms,jme
INTEGER(i_kind), INTENT(IN ) :: its,ite, jts,jte
REAL(r_kind), DIMENSION(ims:ime, jms:jme), &
INTENT(INOUT) :: a
REAL(r_kind), DIMENSION(ims:ime, jms:jme), &
INTENT(IN ) :: bnf
REAL(r_kind), DIMENSION(nol, ims:ime, jms:jme), &
INTENT(IN ) :: lnf
!----------------------------------------------------------------------------
INTEGER(i_kind) :: i,j,l,nola
!============================================================================
DO j=jts+1,jte
nola=MIN(nol,j-jts)
DO i=its,ite
DO l=1,nola
a(i,j)=a(i,j)-lnf(l,i,j)*a(i,j-l)
ENDDO
ENDDO
ENDDO
DO j=jts,jte
DO i=its,ite
a(i,j)=bnf(i,j)*a(i,j)
ENDDO
ENDDO
DO j=jte-1,jts,-1
nola=MIN(nol,jte-j)
DO i=its,ite
DO l=1,nola
a(i,j)=a(i,j)-lnf(l,i,j+l)*a(i,j+l)
ENDDO
ENDDO
ENDDO
END SUBROUTINE hbnrf2j
SUBROUTINE hbnrf3i(a,nol,lnf,bnf, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
!============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: hbnrf3i
!
! prgrmmr:
!
! abstract: Horizontal basic inhomogeneous recursive filter,
! 3-dimensional, active index i
!
! program history log:
! 2008-04-28 safford -- add subprogram doc block
!
! input argument list:
! nol -
! ims,ime, jms,jme, kms,kme -
! its,ite, jts,jte, kts,kte -
! a -
! bnf -
! lnf -
!
! output argument list:
! a -
! bnf -
! lnf -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: i_kind
use kinds, only: r_kind => r_double
IMPLICIT NONE
INTEGER(i_kind), INTENT(IN ) :: nol
INTEGER(i_kind), INTENT(IN ) :: ims,ime, jms,jme, kms,kme
INTEGER(i_kind), INTENT(IN ) :: its,ite, jts,jte, kts,kte
REAL(r_kind), DIMENSION(ims:ime, kms:kme, jms:jme), &
INTENT(INOUT) :: a
REAL(r_kind), DIMENSION(ims:ime, kms:kme, jms:jme), &
INTENT(IN ) :: bnf
REAL(r_kind), DIMENSION(nol, ims:ime, kms:kme, jms:jme), &
INTENT(IN ) :: lnf
!----------------------------------------------------------------------------
INTEGER(i_kind) :: i,j,k,l,nola
!============================================================================
DO j=jts,jte
DO k=kts,kte
DO i=its+1,ite
nola=MIN(nol,i-its)
DO l=1,nola
a(i,k,j)=a(i,k,j)-lnf(l,i,k,j)*a(i-l,k,j)
ENDDO
ENDDO
DO i=its,ite
a(i,k,j)=bnf(i,k,j)*a(i,k,j)
ENDDO
DO i=ite-1,its,-1
nola=MIN(nol,ite-i)
DO l=1,nola
a(i,k,j)=a(i,k,j)-lnf(l,i+l,k,j)*a(i+l,k,j)
ENDDO
ENDDO
ENDDO
ENDDO
END SUBROUTINE hbnrf3i
SUBROUTINE hbnrf3j(a,nol,lnf,bnf, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
!============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: hbnrf3j
!
! prgrmmr:
!
! abstract: Horizontal basic inhomogeneous recursive filter,
! 3-dimensional, active index j
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! nol -
! ims,ime, jms,jme, kms,kme -
! its,ite, jts,jte, kts,kte -
! a -
! bnf -
! lnf -
!
! output argument list:
! a -
! bnf -
! lnf -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: i_kind
use kinds, only: r_kind => r_double
IMPLICIT NONE
INTEGER(i_kind), INTENT(IN ) :: nol
INTEGER(i_kind), INTENT(IN ) :: ims,ime, jms,jme, kms,kme
INTEGER(i_kind), INTENT(IN ) :: its,ite, jts,jte, kts,kte
REAL(r_kind), DIMENSION(ims:ime, kms:kme, jms:jme), &
INTENT(INOUT) :: a
REAL(r_kind), DIMENSION(ims:ime, kms:kme, jms:jme), &
INTENT(IN ) :: bnf
REAL(r_kind), DIMENSION(nol, ims:ime, kms:kme, jms:jme), &
INTENT(IN ) :: lnf
!----------------------------------------------------------------------------
INTEGER(i_kind) :: i,j,k,l,nola
!============================================================================
DO j=jts+1,jte
nola=MIN(nol,j-jts)
DO k=kts,kte
DO i=its,ite
DO l=1,nola
a(i,k,j)=a(i,k,j)-lnf(l,i,k,j)*a(i,k,j-l)
ENDDO
ENDDO
ENDDO
ENDDO
DO j=jts,jte
DO k=kts,kte
DO i=its,ite
a(i,k,j)=bnf(i,k,j)*a(i,k,j)
ENDDO
ENDDO
ENDDO
DO j=jte-1,jts,-1
nola=MIN(nol,jte-j)
DO k=kts,kte
DO i=its,ite
DO l=1,nola
a(i,k,j)=a(i,k,j)-lnf(l,i,k,j+l)*a(i,k,j+l)
ENDDO
ENDDO
ENDDO
ENDDO
END SUBROUTINE hbnrf3j
SUBROUTINE vbnrf1k(a,nol,lnf,bnf, &
kms,kme, &
kts,kte )
!============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: vbnrf1k
!
! prgrmmr:
!
! abstract: Vertical bounded grid inhomogeneous recursive filter,
! 1-dimensional, active index k
!
! program history log:
! 2008-04-28 safford -- add subprogram doc block
!
! input argument list:
! nol -
! kms,kme -
! kts,kte -
! a -
! bnf -
! lnf -
!
! output argument list:
! a -
! bnf -
! lnf -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: i_kind
use kinds, only: r_kind => r_double
IMPLICIT NONE
INTEGER(i_kind), INTENT(IN ) :: nol
INTEGER(i_kind), INTENT(IN ) :: kms,kme
INTEGER(i_kind), INTENT(IN ) :: kts,kte
REAL(r_kind), DIMENSION(kms:kme), &
INTENT(INOUT) :: a
REAL(r_kind), DIMENSION(kms:kme), &
INTENT(IN ) :: bnf
REAL(r_kind), DIMENSION(nol, kms:kme), &
INTENT(IN ) :: lnf
!----------------------------------------------------------------------------
INTEGER(i_kind) :: k,l,nola
!============================================================================
DO k=kts+1,kte
nola=MIN(nol,k-kts)
DO l=1,nola
a(k)=a(k)-lnf(l,k)*a(k-l)
ENDDO
ENDDO
DO k=kts,kte
a(k)=bnf(k)*a(k)
ENDDO
DO k=kte-1,kts,-1
nola=MIN(nol,kte-k)
DO l=1,nola
a(k)=a(k)-lnf(l,k+l)*a(k+l)
ENDDO
ENDDO
END SUBROUTINE vbnrf1k
SUBROUTINE vbnrf2k(a,nol,lnf,bnf, &
ims,ime, kms,kme, &
its,ite, kts,kte )
!============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: vbnrf2k
!
! prgrmmr:
!
! abstract: Vertical bounded grid inhomogeneous recursive filter,
! 2-dimensional, active index k
!
! program history log:
! 2008-04-25 safford -- add subprogram doc block
!
! input argument list:
! nol -
! ims,ime, kms,kme -
! its,ite, kts,kte -
! a -
! bnf -
! lnf -
!
! output argument list:
! a -
! bnf -
! lnf -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: i_kind
use kinds, only: r_kind => r_double
IMPLICIT NONE
INTEGER(i_kind), INTENT(IN ) :: nol
INTEGER(i_kind), INTENT(IN ) :: ims,ime, kms,kme
INTEGER(i_kind), INTENT(IN ) :: its,ite, kts,kte
REAL(r_kind), DIMENSION(ims:ime, kms:kme), &
INTENT(INOUT) :: a
REAL(r_kind), DIMENSION(ims:ime, kms:kme), &
INTENT(IN ) :: bnf
REAL(r_kind), DIMENSION(nol, ims:ime, kms:kme), &
INTENT(IN ) :: lnf
!----------------------------------------------------------------------------
INTEGER(i_kind) :: i,k,l,nola
!============================================================================
DO k=kts+1,kte
nola=MIN(nol,k-kts)
DO i=its,ite
DO l=1,nola
a(i,k)=a(i,k)-lnf(l,i,k)*a(i,k-l)
ENDDO
ENDDO
ENDDO
DO k=kts,kte
DO i=its,ite
a(i,k)=bnf(i,k)*a(i,k)
ENDDO
ENDDO
DO k=kte-1,kts,-1
nola=MIN(nol,kte-k)
DO i=its,ite
DO l=1,nola
a(i,k)=a(i,k)-lnf(l,i,k+l)*a(i,k+l)
ENDDO
ENDDO
ENDDO
END SUBROUTINE vbnrf2k
SUBROUTINE vbnrf3k(a,nol,lnf,bnf, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
!============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: vbnrf3k
!
! prgrmmr:
!
! abstract: Vertical bounded grid inhomogeneous recursive filter,
! 3-dimensional, active index k
!
! program history log:
! 2008-04-28 safford -- add subprogram doc block
!
! input argument list:
! nol -
! ims,ime, jms,jme, kms,kme -
! its,ite, jts,jte, kts,kte -
! a -
! bnf -
! lnf -
!
! output argument list:
! a -
! bnf -
! lnf -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: i_kind
use kinds, only: r_kind => r_double
IMPLICIT NONE
INTEGER(i_kind), INTENT(IN ) :: nol
INTEGER(i_kind), INTENT(IN ) :: ims,ime, jms,jme, kms,kme
INTEGER(i_kind), INTENT(IN ) :: its,ite, jts,jte, kts,kte
REAL(r_kind), DIMENSION(ims:ime, kms:kme, jms:jme), &
INTENT(INOUT) :: a
REAL(r_kind), DIMENSION(ims:ime, kms:kme, jms:jme), &
INTENT(IN ) :: bnf
REAL(r_kind), DIMENSION(nol, ims:ime, kms:kme, jms:jme), &
INTENT(IN ) :: lnf
!----------------------------------------------------------------------------
INTEGER(i_kind) :: i,j,k,l,nola
!============================================================================
DO j=jts,jte
DO k=kts+1,kte
nola=MIN(nol,k-kts)
DO i=its,ite
DO l=1,nola
a(i,k,j)=a(i,k,j)-lnf(l,i,k,j)*a(i,k-l,j)
ENDDO
ENDDO
ENDDO
DO k=kts,kte
DO i=its,ite
a(i,k,j)=bnf(i,k,j)*a(i,k,j)
ENDDO
ENDDO
DO k=kte-1,kts,-1
nola=MIN(nol,kte-k)
DO i=its,ite
DO l=1,nola
a(i,k,j)=a(i,k,j)-lnf(l,i,k+l,j)*a(i,k+l,j)
ENDDO
ENDDO
ENDDO
ENDDO
END SUBROUTINE vbnrf3k
SUBROUTINE hbncij(a,hamp,nol,lnfi,bnfi,lnfj,bnfj, &
ims,ime, jms,jme, &
its,ite, jts,jte )
!============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: hbncij
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-28 safford -- add subprogram doc block
!
! input argument list:
! nol -
! ims,ime, jms,jme -
! its,ite, jts,jte -
! hamp,bnfi,bnfj -
! lnfi,lnfj -
! a -
!
! output argument list:
! a -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: i_kind
use kinds, only: r_kind => r_double
IMPLICIT NONE
INTEGER(i_kind), INTENT(IN ) :: nol
INTEGER(i_kind), INTENT(IN ) :: ims,ime, jms,jme
INTEGER(i_kind), INTENT(IN ) :: its,ite, jts,jte
REAL(r_kind), DIMENSION(ims:ime, jms:jme), &
INTENT(INOUT) :: a
REAL(r_kind), DIMENSION(ims:ime, jms:jme), &
INTENT(IN ) :: hamp,bnfi,bnfj
REAL(r_kind), DIMENSION(nol, ims:ime, jms:jme), &
INTENT(IN ) :: lnfi,lnfj
!----------------------------------------------------------------------------
INTEGER(i_kind) :: i,j
!============================================================================
DO j=jts,jte
DO i=its,ite
a(i,j)=hamp(i,j)*a(i,j)
ENDDO
ENDDO
!---------------
CALL hbnrf2i(a,nol,lnfi,bnfi, &
ims,ime, jms,jme, &
its,ite, jts,jte)
!----------
CALL hbnrf2j(a,nol,lnfj,bnfj, &
ims,ime, jms,jme, &
its,ite, jts,jte)
!----------
END SUBROUTINE hbncij
SUBROUTINE hbncji(a,hamp,nol,lnfi,bnfi,lnfj,bnfj, &
ims,ime, jms,jme, &
its,ite, jts,jte )
!============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: hbncji
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-28 safford -- add subprogram doc block
!
! input argument list:
! nol -
! ims,ime, jms,jme -
! its,ite, jts,jte -
! hamp,bnfi,bnfj -
! lnfi,lnfj -
! a -
!
! output argument list:
! a -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: i_kind
use kinds, only: r_kind => r_double
IMPLICIT NONE
INTEGER(i_kind), INTENT(IN ) :: nol
INTEGER(i_kind), INTENT(IN ) :: ims,ime, jms,jme
INTEGER(i_kind), INTENT(IN ) :: its,ite, jts,jte
REAL(r_kind), DIMENSION(ims:ime, jms:jme), &
INTENT(INOUT) :: a
REAL(r_kind), DIMENSION(ims:ime, jms:jme), &
INTENT(IN ) :: hamp,bnfi,bnfj
REAL(r_kind), DIMENSION(nol, ims:ime, jms:jme), &
INTENT(IN ) :: lnfi,lnfj
!----------------------------------------------------------------------------
INTEGER(i_kind) :: i,j
!============================================================================
CALL hbnrf2j(a,nol,lnfj,bnfj, &
ims,ime, jms,jme, &
its,ite, jts,jte)
!----------
CALL hbnrf2i(a,nol,lnfi,bnfi, &
ims,ime, jms,jme, &
its,ite, jts,jte)
!---------------
DO j=jts,jte
DO i=its,ite
a(i,j)=hamp(i,j)*a(i,j)
ENDDO
ENDDO
!---------------
END SUBROUTINE hbncji
SUBROUTINE hbncijk(a,hamp,nol,lnfi,bnfi,lnfj,bnfj,lnfk,bnfk, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
!============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: hbncijk
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-28 safford -- add subprogram doc block
!
! input argument list:
! nol -
! ims,ime, jms,jme, kms,kme -
! its,ite, jts,jte, kts,kte -
! hamp,bnfi,bnfj,bnfk -
! lnfi,lnfj,lnfk -
! a -
!
! output argument list:
! a -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: i_kind
use kinds, only: r_kind => r_double
IMPLICIT NONE
INTEGER(i_kind), INTENT(IN ) :: nol
INTEGER(i_kind), INTENT(IN ) :: ims,ime, jms,jme, kms,kme
INTEGER(i_kind), INTENT(IN ) :: its,ite, jts,jte, kts,kte
REAL(r_kind), DIMENSION(ims:ime, kms:kme, jms:jme), &
INTENT(INOUT) :: a
REAL(r_kind), DIMENSION(ims:ime, kms:kme, jms:jme), &
INTENT(IN ) :: hamp,bnfi,bnfj,bnfk
REAL(r_kind), DIMENSION(nol, ims:ime, kms:kme, jms:jme), &
INTENT(IN ) :: lnfi,lnfj,lnfk
!----------------------------------------------------------------------------
INTEGER(i_kind) :: i,j,k
!============================================================================
DO j=jts,jte
do k=kts,kte
DO i=its,ite
a(i,k,j)=hamp(i,k,j)*a(i,k,j)
ENDDO
enddo
ENDDO
!---------------
CALL hbnrf3i(a,nol,lnfi,bnfi, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte)
!----------
CALL hbnrf3j(a,nol,lnfj,bnfj, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte)
!----------
call vbnrf3k(a,nol,lnfk,bnfk, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte)
END SUBROUTINE hbncijk
SUBROUTINE hbnckji(a,hamp,nol,lnfi,bnfi,lnfj,bnfj,lnfk,bnfk, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
!============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: hbnckji
!
! prgrmmr:
!
! abstract:
!
! program history log:
! 2008-04-28 safford -- add subprogram doc block
!
! input argument list:
! nol -
! ims,ime, jms,jme, kms,kme -
! its,ite, jts,jte, kts,kte -
! hamp,bnfi,bnfj,bnfk -
! lnfi,lnfj,lnfk -
! a -
!
! output argument list:
! a -
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: i_kind
use kinds, only: r_kind => r_double
IMPLICIT NONE
INTEGER(i_kind), INTENT(IN ) :: nol
INTEGER(i_kind), INTENT(IN ) :: ims,ime, jms,jme, kms,kme
INTEGER(i_kind), INTENT(IN ) :: its,ite, jts,jte, kts,kte
REAL(r_kind), DIMENSION(ims:ime, kms:kme, jms:jme), &
INTENT(INOUT) :: a
REAL(r_kind), DIMENSION(ims:ime, kms:kme, jms:jme), &
INTENT(IN ) :: hamp,bnfi,bnfj,bnfk
REAL(r_kind), DIMENSION(nol, ims:ime, kms:kme, jms:jme), &
INTENT(IN ) :: lnfi,lnfj,lnfk
!----------------------------------------------------------------------------
INTEGER(i_kind) :: i,j,k
!============================================================================
call vbnrf3k(a,nol,lnfk,bnfk, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte)
!----------
CALL hbnrf3j(a,nol,lnfj,bnfj, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte)
!----------
CALL hbnrf3i(a,nol,lnfi,bnfi, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte)
!---------------
DO j=jts,jte
do k=kts,kte
DO i=its,ite
a(i,k,j)=hamp(i,k,j)*a(i,k,j)
ENDDO
enddo
ENDDO
!---------------
END SUBROUTINE hbnckji
!============================================================================
subroutine rfit(ng,sig,nu, ns,nw,ssig,snu,ins1,wts)
!============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: rfit
!
! prgrmmr: R. J. Purser, NCEP 2001
!
! abstract:
!
! program history log:
! 2008-04-28 safford -- add subprogram doc block, rm unused vars
!
! input argument list:
! ng
! sig,nu
! ins1
! wts
!
! output argument list:
! ins1
! wts
! ns,nw
! ssig,snu
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: i_kind
use kinds, only: r_kind => r_double
use constants, only: zero,one
use module_fitcons
implicit none
integer(i_kind), intent(IN ) :: ng
real(r_kind) , dimension(ng), intent(IN ) :: sig,nu
integer(i_kind), intent( OUT) :: ns,nw
real(r_kind) , dimension(ng), intent( OUT) :: ssig,snu
integer(i_kind), dimension(ng), intent(INOUT) :: ins1
real(r_kind) , dimension(no,ng),intent(INOUT) :: wts
!----------------------------------------------------------------------------
integer(i_kind) :: i,i1,im,k,l,is
real(r_kind) :: t
real(r_kind) , dimension(-nohm:ng+noh) :: dcdg
real(r_kind) , dimension(-noh:ng+noh) :: cofg,cofs
real(r_kind) , dimension(ng) :: dsdg,dhdg
!============================================================================
nw=0
do i=1,ng
dcdg(i)=one/sig(i)
if(sig(i) <= sigb)then
!----------------------------------------------------------------------------
! sig(i) below threshold; cleave to original grid spacing with ds/dg and
! dh/dg set accordingly:
!----------------------------------------------------------------------------
dsdg(i)=one ; dhdg(i)=zero
else
!----------------------------------------------------------------------------
! sig(i) exceeds basic threshold sigb, allowing working grid with coordinate
! s to differ from original grid with coordinate g. The formula for ds/dg
! is now <1 but tends smoothly to 1 again at the threshold value, sig=sigb.
! [The function for log(ds/dg) is based on the "hyper-hyperbola":
! y= (1+x**4)**(-4)-1, which rises very gradually from its base at x=y=0]
! Likewise, the perturbative component, dh/dg, is now < 0, but tends
! smoothly to 0 again at the threshold.
!----------------------------------------------------------------------------
t=ldsig-sqrt(sqrt(ldsig4+(ldsig-log(sigc*dcdg(i)))**4))
dsdg(i)=exp(t) ; dhdg(i)=dsdg(i)*t
endif
enddo
!----------------------------------------------------------------------------
! Apply mirror-symmetry to extrapolate beyond ends:
!----------------------------------------------------------------------------
do l=1,noh
dcdg(1-l)=dcdg(l); dcdg(ng+l)=dcdg(ng+1-l)
enddo
!----------------------------------------------------------------------------
! Integrate dc/dg wrt g to get c(g) at each of the points of the g-grid
! which is NOT staggered relative to the boundary
!----------------------------------------------------------------------------
cofg(0)=zero
do i=1,ng
cofg(i)=cofg(i-1)+dot_product(qco,dcdg(i-noh:i+noh))
enddo
do l=1,noh
cofg( -l)=-cofg( l)
cofg(ng+l)=-cofg(ng-l)+2*cofg(ng)
enddo
im=0
ns=0
!----------------------------------------------------------------------------
! loop over noncontiguous segments where it is numerically beneficial
! to employ a grid of relatively coarse resolution. The adoption of each
! alternative grid segment is subject to some conditions:
! 1) Each coarse-grid segment must span at least 5 points of the original grid
! 2) Each segment must shorten the tally of working grid points by at least 3.
! Subject to the above conditions, the coarse grid is blended smoothly
! with the original grid at internal thresholds and is designed to provide
! a resolution such that the smoothing scale, sigma, never exceeds the
! product, sigc*dg/ds, where sigc is a dimensionless parameter (e.g. sigc=3.)
! and dg/ds is the local working grid (s) spacing in units of the original
! grid (g) spacing.
!
! Each segment found is defined by its end points in the original grid,
! i1 and im. k is the counter for segments along this line.
! ns keeps count of the number of working grid (s-grid) points found so far.
!----------------------------------------------------------------------------
cofs(0)=zero
do k=1,ng
do i1=im+1,ng
if(i1< ng-3 .and. dhdg(i1) /= zero)exit
!----------------------------------------------------------------------------
! working s-grid continues to track the original g-grid; Set indices and
! weight for the trivial "interpolation" between these coincident grids:
!----------------------------------------------------------------------------
ns=ns+1
ins1(i1)=-ns
cofs(ns)=cofg(i1)
enddo
if(i1 > ng)exit
!----------------------------------------------------------------------------
! Having met the basic conditions for the start of a new segment in which
! the s-grid and g-grids may part company, seek the other end, im, of this
! possible segment:
!----------------------------------------------------------------------------
do im=i1+1,ng
if(dhdg(im) == zero)exit
enddo
im=im-1
if(im < i1+4)then
!----------------------------------------------------------------------------
! Segment too short to be viable; keep s-grid and g-grids synchronized:
!----------------------------------------------------------------------------
do i=i1,im
ns=ns+1
ins1(i)=-ns
cofs(ns)=cofg(i)
enddo
else
!----------------------------------------------------------------------------
! Segment long enough to be potentially viable. Call jfit to determine if
! the final condition is met, namely that the number of s-grid points
! in this segment is smaller than the g-grid tally by at least 3. If so,
! Fit an exact integer number of s-points into this segment and compute
! the indices and weights for the associated nontrivial interpolation
! from these (and neighboring) s-points back to the g-points of the segment:
!----------------------------------------------------------------------------
call jfit(ng,i1,im,ns,nw,cofg,dsdg,dhdg,cofs,ins1,wts)
if(ns<0) return
endif
enddo
if(ns<no .and. nw>0)then ! <- s-grid too short; use copy of g-grid instead
wts(:,1:nw)=zero
nw=0
do i=1,ng
ins1(i)=-i
cofs(i)=cofg(i)
enddo
ns=ng
endif
do l=1,noh
cofs( -l)=-cofs( l)
cofs(ns+l)=-cofs(ns-l)+2*cofs(ns)
enddo
do is=1,ns
ssig(is)=one/dot_product(dco,cofs(is-nohp:is+noh))
enddo
!----------------------------------------------------------------------------
! By applying adjoint-interpolation to the g-grid metric terms, obtain
! the corresponding metric terms for the new s-grid:
!----------------------------------------------------------------------------
call stogt(ns,ng,ins1,wts, snu,nu)
end subroutine rfit
!============================================================================
subroutine jfit(ng,ig1,igm,ns,iw,cofg,dsdg,dhdg,cofs,ins1,wts)
!============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: jfit
!
! prgrmmr: R. J. Purser, NCEP 2001
!
! abstract:
!
! program history log:
! 2008-04-28 safford -- add subprogram doc block, rm unused vars
!
! input argument list:
! ng,ig1,igm
! ns,iw
! dsdg,dhdg
! cofg
! cofs
! ins1
! wts
!
! output argument list:
! ns,iw
! cofs
! ins1
! wts
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: i_kind
use kinds, only: r_kind => r_double
use constants, only: zero,half,one,two
use module_fitcons
implicit none
integer(i_kind), intent(IN ) :: ng,ig1,igm
integer(i_kind), intent(INOUT) :: ns,iw
real(r_kind) , dimension(ng), intent(IN ) :: dsdg,dhdg
real(r_kind) , dimension(-noh:ng+noh), intent(IN ) :: cofg
real(r_kind) , dimension(-noh:ng+noh), intent(INOUT) :: cofs
integer(i_kind), dimension(ng), intent(INOUT) :: ins1
real(r_kind) , dimension(no,ng), intent(INOUT) :: wts
!----------------------------------------------------------------------------
real(r_kind) , dimension(-noh:ng+noh) :: sofg,dsdgt
real(r_kind) :: et,estar,destar,r,dr,sm
integer(i_kind) :: i,l,ie,iep,ie1,ien,ig0,is0,ism,init
!============================================================================
!----------------------------------------------------------------------------
! Form the definite integral sm, of ds/dg, within this segment:
!----------------------------------------------------------------------------
sm=sum(dsdg(ig1:igm))
!---------------------------------------------------------------------------
! Test whether it is worthwhile to allow s-grid to deviate from the original
! g-grid within this segment on the basis of the number of grid points that
! could potentially be eliminated (we require a saving > 3 per segment):
!---------------------------------------------------------------------------
if(sm > igm-ig1-2)then
!----------------------------------------------------------------------------
! This putative s-grid segment reduces the total number of grid points by an
! insufficient amount to justify the additional interpolations. Therefore,
! keep the original g-grid instead for this segment, and return:
!---------------------------------------------------------------------------
do i=ig1,igm
ns=ns+1
ins1(i)=-ns
cofs(ns)=cofg(i)
enddo
return
endif
!----------------------------------------------------------------------------
! s-grid segment achieves a worthwhile reduction of the number of points
! of the working grid. The tasks of the rest of this routine are to:
! (1) adjust the segment length in the s-metric to make it an integer;
! (2) find the s-coordinates of each g-grid points in this segment
! and hence the nontrivial interpolation indices and weights required
! to go from the s-grid to the g-grid (or adjoints going the other way);
! (3) use Newton iterations to find the accurate interpolation formulae
! that enable c(s) to be interpolated from the given c(g).
!----------------------------------------------------------------------------
ig0=ig1-1
is0=ns; ism=sm
!----------------------------------------------------------------------------
! Fractional remainder of sm, divided by the definite integral of dh/dg
! provides the adjustment factor that scales the perturbative component,
! dhdg, by exactly the amount that will make the segment integral of the
! perturbed grid-to-grid jacobian, dsdgt, the exact integer, ism:
!----------------------------------------------------------------------------
r=(sm-ism)/sum(dhdg(ig1:igm))
do i=ig1,igm
dsdgt(i)=dsdg(i)-r*dhdg(i)
enddo
!----------------------------------------------------------------------------
! Mirror-extrapolate adjusted ds/dg as an even-symmetry function at the
! ends of this segment. Note that the grid on which derivatives such as
! ds/dg reside is the one staggered wrt domain boundaries and segment
! end points. The indices of this grid go from ig1 to igm inside the
! segment. (The convention for the companion grid, NOT staggered wrt
! boundaries, is such that the two segment ends are denoted by indices,
! ig0=ig1-1 and igm.)
!----------------------------------------------------------------------------
do l=1,noh
dsdgt(ig1-l)=dsdgt(ig0 +l)
dsdgt(igm+l)=dsdgt(igm+1-l)
enddo
ism=is0+ism ! This integer also becomes (within round-off) the value, sofg(igm)
!----------------------------------------------------------------------------
! Set s(g) at both ends of the segment to be the appropriate integers:
!----------------------------------------------------------------------------
sofg(ig0)=is0; sofg(igm)=ism
!----------------------------------------------------------------------------
! Get s(g) inside the segment by performing a numerical quadrature of dsdgt:
!----------------------------------------------------------------------------
do i=ig1,igm
sofg(i)=sofg(i-1)+dot_product(qco,dsdgt(i-noh:i+noh))
enddo
!----------------------------------------------------------------------------
! Mirror-extrapolate s(g) as an odd-symmetry function at segment end points.
! Note that, being an inegral, s(g) resides on the grid NOT staggered wrt
! boundaries and segment end points.
!----------------------------------------------------------------------------
do l=1,noh
sofg(ig0-l)=two*is0-sofg(ig0+l)
sofg(igm+l)=two*ism-sofg(igm-l)
enddo
do i=ig1,igm
iw=iw+1 ; wts(:,iw)=zero
r=dot_product(ico,sofg(i-nohp:i+noh))+half
ie=r ! Take integer part...
ins1(i)=ie-nohm ! ...hence the index of the first point in the stencil...
r=r-ie ! ...then the fractional part to find interpolation weights:
call lagw(hunit1,r,q1,wt1,dwt1,nom) ! weights for left-biased stencil
wts(:nom,iw) = (one-r)*wt1 ! bias weight, 1-r
call lagw(hunit2,r,q1,wt1,dwt1,nom) ! weights for right-biased stencil
wts(2: ,iw) = wts(2: ,iw) +r*wt1 ! bias weight, r.
!----------------------------------------------------------------------------
! Exploit the mirror symmetries to confine the weight stencil to the
! domain interior, even though this may entail padding innermost end of
! the stencil with useless zeroes:
!----------------------------------------------------------------------------
L=1-INS1(I)
IF(L > 0)THEN ! FOLD LEFT OVERLAP OF L ELEMENTS BACK INSIDE:
WTS(1:L,IW) =WTS(L:1:-1,IW)+WTS(L+1:L*2,IW) ! FOLD INTO 1ST L
WTS(L+1:NO-L,IW) =WTS(L*2+1:NO,IW) ! SHIFT THE REST LEFT
WTS(NOP-L:NO,IW)=zero ! SET TRAILING L ELEMENTS TO ZERO
INS1(I)=1 ! RESET INDEX OF FIRST POINT OF STENCIL
ENDIF
l=ins1(i)+nom-ism
if(l > 0)then ! Fold right overlap of L elements back inside:
wts(nop-l:no,iw)=wts(no:nop-l:-1,iw)+wts(nop-l*2:no-l,iw) ! Fold last L
wts(l+1:no-l,iw)=wts(1:no-l*2,iw) ! Shift right
wts(1:l,iw)=zero ! Set first L elements to zero
ins1(i)=ism-nom ! reset index of first point of stencil
endif
enddo
ns=ism
!----------------------------------------------------------------------------
! Use Newton-Raphson iterations to locate the g-coordinates of all this
! segment's s-grid points. Then interpolate the function c to each of
! these s-grid points. (Note that, in the present context, the
! s- and g-grids are the ones NOT staggered wrt the domain boundaries.)
!----------------------------------------------------------------------------
ie=ig0
iloop: do i=is0+1,ism-1 ! Loop over s-grid target points interior to this segment
et=i
!----------------------------------------------------------------------------
! Find the g-grid interval containing this target:
!----------------------------------------------------------------------------
do iep=ie+1,igm-1; if(sofg(iep) > et)exit; enddo
do ie=iep-1,ig1,-1; if(sofg(ie) <= et)exit; enddo
ie1=ie-nohm; ien=ie+noh ! <-- Set the interpolation stencil range:
r=(et-sofg(ie))/(sofg(ie+1)-sofg(ie)) ! Linearly estimate interval fraction
!----------------------------------------------------------------------------
! Perform Newton-Raphson iterations to refine interval fraction, r:
!----------------------------------------------------------------------------
do init=1,nnit
call lagw(hunit,r,q,wt,dwt,no) ! Get Lagrange weights, wt and d(wt)/dg
estar =dot_product(wt, sofg(ie1:ien))-et ! <- Residual error, estar.
destar=dot_product(dwt,sofg(ie1:ien)) ! <- d(estar)/dg.
dr=-estar/destar ! <- Newton correction to r
r=r+dr ! <- Refined estimate, r
if(abs(dr) <= rcrit)then ! <- Converged enough yet?
wt=wt+dr*dwt ! <- Final refinement to wt
cofs(i)=dot_product(wt, cofg(ie1:ien)) ! <- Interpolate c(s)
cycle iloop
end if
enddo
! stop 'Too many Newton iterations' ! <- It never convergenced!
write(6,*)' Too many Newton iterations' ! <- It never convergenced!
ns=-1
return
enddo iloop
cofs(ism)=cofg(igm) ! <- End value directly
return
end subroutine jfit
!============================================================================
subroutine stog(ns,ng,ins1,wts, as,ag)
!============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: stog
!
! prgrmmr: R. J. Purser NCEP 2001
!
! abstract: Forward interpolation from s-grid to g-grid
!
! program history log:
! 2008-04-28 safford -- add subprogram doc block
!
! input argument list:
! ns,ng - sizes of s and g grids
! ins1 - array of 1st stencil indices (s-grid) for each target (g) point.
! wts - interpolation weights for each target (g-grid point).
! as - s-grid array of source data.
!
! output argument list:
! ag - g-grid array of interpolated target data.
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: i_kind
use kinds, only: r_kind => r_double
use constants, only: zero
implicit none
integer(i_kind), parameter :: noh=3,no=noh*2,nom=no-1
integer(i_kind), intent(IN ) :: ns,ng
integer(i_kind), dimension(ng), intent(IN ) :: ins1
real(r_kind) , dimension(no,ng),intent(IN ) :: wts
real(r_kind) , dimension(ns), intent(IN ) :: as
real(r_kind) , dimension(ng), intent( OUT) :: ag
!----------------------------------------------------------------------------
integer(i_kind) :: i,is,iw
!============================================================================
iw=0
ag=zero
do i=1,ng
is=ins1(i)
if(is>0)then
iw=iw+1
ag(i)=dot_product(wts(:,iw),as(is:is+nom))
else
ag(i)=as(-is)
endif
enddo
end subroutine stog
!============================================================================
subroutine stogt(ns,ng,ins1,wts, as,ag)
!============================================================================
!$$$ subprogram documentation block
! . . .
! subprogram: stogt
!
! prgrmmr: R. J. Purser NCEP 2001
!
! abstract: Perform the transpose of the operation defined by stog
!
! program history log:
! 2008-04-28 safford -- add subprogram doc block
!
! input argument list:
! ns,ng
! ins1
! wts
! ag
!
! output argument list:
! as
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$ end documentation block
use kinds, only: i_kind
use kinds, only: r_kind => r_double
use constants, only: zero
implicit none
integer(i_kind), parameter :: noh=3,no=noh*2,nom=no-1
integer(i_kind), intent(IN ) :: ns,ng
integer(i_kind), dimension(ng), intent(IN ) :: ins1
real(r_kind) , dimension(no,ng),intent(IN ) :: wts
real(r_kind) , dimension(ns), intent( OUT) :: as
real(r_kind) , dimension(ng), intent(IN ) :: ag
!----------------------------------------------------------------------------
integer(i_kind) :: i,is,iw
!============================================================================
iw=0
as=zero
do i=1,ng
is=ins1(i)
if(is>0)then
iw=iw+1
as(is:is+nom)=as(is:is+nom)+wts(:,iw)*ag(i)
else
as(-is)=as(-is)+ag(i)
endif
enddo
end subroutine stogt