subroutine simpin1(wgts,wgtsx1,wgtsx11,iwgts, &
iflag,x1in,nin,iord,lbig,x1grid,n1grid,nf, &
ndx1,ndx11,ixi,tl,alocal,blocal)
!$$$ subprogram documentation block
! . . . .
! subprogram: simpin1
! prgmmr: parrish org: np22 date: 2000-01-01
!
! abstract: This routine computes interpolation weights for 1-d
! interpolation
!
! program history log:
! 2000-01-01 parrish
! 2004-06-22 treadon - update documentation
! 2008-04-12 safford - rm unused vars
!
! input argument list:
! x1in - coordinates of interpolatees
! nin - number of interpolatees
! iord - order of interpolation (1=linear, 2=quadratic, etc)
! lbig - number of interpolating points (=iord+1)
! x1grid - coordinates of interpolator grid
! (can be monotonic increasing or decreasing)
! n1grid - dimension of interpolator grid
! nf - =1, then return interpolation weights for function f
! ndx1 - =1, then return interpolation weights for df/dx1
! ndx11 - =1, then return interpolation weights for d2f/(dx1*dx1)
!
! output argument list:
! wgts - interpolation weights for function ( wgts(nin,lbig) )
! wgtsx1 - interpolation weights for df/dx1
! wgtsx11 - interpolation weights for d2f/(dx1*dx1)
! note: if any of these 3 are not asked for, they
! can be dummy arguments
! iwgts - absolute grid addresses ( iwgts(nin,lbig) )
! iflag - flag for each interpolatee (iflag(nin))
! =0, then point too close to edge of domain, but weights
! are computed anyway, in case it is considered
! OK to extrapolate.
! =1, then weights computed
!
! attributes:
! language: f90
! machine: ibm rs/6000 sp
!
!$$$ end documentation block
use kinds, only: r_kind,i_kind
use constants, only: zero, half, one
implicit none
integer(i_kind),intent(in ) :: nin,iord,lbig,n1grid,nf,ndx1,ndx11
real(r_kind) ,intent(in ) :: x1in(nin),x1grid(n1grid)
real(r_kind) ,intent( out) :: wgts(nin,lbig)
integer(i_kind),intent( out) :: iwgts(nin,lbig)
real(r_kind) ,intent( out) :: wgtsx1(nin,lbig),wgtsx11(nin,lbig)
integer(i_kind),intent(in ) :: ixi(0:iord)
real(r_kind) ,intent(in ) :: tl(lbig,lbig,n1grid,2),alocal(n1grid,2),blocal(n1grid,2)
real(r_kind) dx1gridi(-3:n1grid+3)
integer(i_kind) i1ref(nin)
integer(i_kind) ix1sign(nin)
real(r_kind) x1p(nin)
integer(i_kind) ix1temp(nin)
integer(i_kind) iflag(nin)
real(r_kind) z0(max(64,nin/lbig),lbig)
integer(i_kind), allocatable, dimension(:) :: ix1grid
real(r_kind) dxa,dxb,dxc,dxd,dxe,dxf,dxg,x1ref,xboundleft,xboundright
real(r_kind) dxfine,dxminmin,dxthis,dxfinei,dxmax,dxmin
integer(i_kind) i1ref0,i1fine,l,iflip,nn,j,k,ip,nend,nstart,lp,&
nthis,nmaxright,n,nminleft,iximx,ntl
integer(i_kind) ii,n1fine,iximax,iximn,i
ntl=max(64,nin/lbig)
! Find and mark all points that are outside interval of interpolation
iximx=maxval(ixi) ; iximn=minval(ixi)
iximax=max(abs(iximx),abs(iximn))
! write(6,*)' iximx,mn=',iximx,iximn
! write(6,*)' iximax=',iximax
! write(6,*)' ixi=',ixi
nminleft=abs(iximn)+1 ; nmaxright=n1grid-iximax
if(iord==1) then
xboundleft=x1grid(1)
xboundright=x1grid(n1grid)
else
xboundleft=x1grid(1+iximax)-.49999_r_kind*(x1grid(1+iximax)-x1grid(iximax))
xboundright=x1grid(n1grid-iximax)+ &
.49999_r_kind*(x1grid(n1grid-iximax+1)-x1grid(n1grid-iximax))
end if
! write(6,*) ' xboundleft=',xboundleft
! write(6,*) ' xboundright=',xboundright
! do i=1,iximax+1
! write(6,*) ' x1grid(',i,')=',x1grid(i)
! end do
! do i=n1grid,n1grid-iximax,-1
! write(6,*)' x1grid(',i,')=',x1grid(i)
! end do
if(x1grid(n1grid)>x1grid(1)) then
iflag=1
do n=1,nin
if(x1in(n)<=xboundleft.or.x1in(n)>=xboundright) iflag(n)=0
end do
end if
if(x1grid(n1grid)<x1grid(1)) then
iflag=1
do n=1,nin
if(x1in(n)>=xboundleft.or.x1in(n)<=xboundright) iflag(n)=0
end do
end if
! Set up uniform fine grid to use in finding interpolation coordinates
dxmax=-huge(dxmax) ; dxmin=huge(dxmin)
do i=1,n1grid-1
dxthis=x1grid(i+1)-x1grid(i)
dx1gridi(i)=one/dxthis
dxmax=max(dxthis,dxmax) ; dxmin=min(dxthis,dxmin)
end do
dx1gridi(-3:0)=dx1gridi(1)
dx1gridi(n1grid:n1grid+3)=dx1gridi(n1grid-1)
if(dxmax*dxmin<=zero) then
write(6,*)' INTERPOLATION GRID NOT MONOTONIC IN SIMPIN1'
call stop2(999)
end if
dxminmin=min(abs(dxmax),abs(dxmin))
dxfine=sign(dxminmin,dxmax)
dxfinei=one/dxfine
n1fine=1+ceiling((x1grid(n1grid)-x1grid(1))/dxfine)
! write(6,*)' in simpin1, n1grid,n1fine=',n1grid,n1fine
allocate (ix1grid(n1fine))
ii=1
do i=1,n1fine
x1ref=x1grid(1)+(i-1)*dxfine
if(dxfinei*(x1ref-x1grid(ii+1))>=-.001_r_kind) ii=min(ii+1,n1grid-1)
ix1grid(i)=ii
! write(6,*)' x1fine,x1grid=',x1ref,x1grid(ix1grid(i))
end do
! Now get i1ref, index of interval containing point
i1ref=-1
do n=1,nin
i1fine=1+nint((x1in(n)-x1grid(1))*dxfinei)
i1fine=max(1,min(i1fine,n1fine))
i1ref0=max(4,min(ix1grid(i1fine),n1grid-4))
dxa=(x1in(n)-x1grid(i1ref0-3))*dx1gridi(i1ref0-3)
dxb=(x1in(n)-x1grid(i1ref0-2))*dx1gridi(i1ref0-2)
dxc=(x1in(n)-x1grid(i1ref0-1 ))*dx1gridi(i1ref0-1)
dxd=(x1in(n)-x1grid(i1ref0 ))*dx1gridi(i1ref0 )
dxe=(x1in(n)-x1grid(i1ref0+1 ))*dx1gridi(i1ref0+1)
dxf=(x1in(n)-x1grid(i1ref0+2))*dx1gridi(i1ref0+2)
dxg=(x1in(n)-x1grid(i1ref0+3))*dx1gridi(i1ref0+3)
if(dxa<=one) i1ref(n)=i1ref0-3
if(dxb>=zero.and.dxb<=one) i1ref(n)=i1ref0-2
if(dxc>=zero.and.dxc<=one) i1ref(n)=i1ref0-1
if(dxd>=zero.and.dxd<=one) i1ref(n)=i1ref0
if(dxe>=zero.and.dxe<=one) i1ref(n)=i1ref0+1
if(dxf>=zero.and.dxf<=one) i1ref(n)=i1ref0+2
if(dxg>=zero) i1ref(n)=i1ref0+3
i1ref(n)=max(nminleft,min(i1ref(n),nmaxright))
end do
deallocate (ix1grid)
! write(6,*)' max,min(i1ref)=',maxval(i1ref),minval(i1ref)
! i1ref now has index of interval containing point
! adjust to be closest one to interpolating point
ix1sign=1
do n=1,nin
dxa=(x1in(n)-x1grid(i1ref(n)))*dx1gridi(i1ref(n))
if(dxa>half.and.i1ref(n)<nmaxright) then
i1ref(n)=i1ref(n)+1
ix1sign(n)=-1
end if
end do
! get interpolation indices
l=0
do i=0,iord
ix1temp=0
l=l+1
do n=1,nin
ix1temp(n)=i1ref(n)+ix1sign(n)*ixi(i)
iwgts(n,l)=ix1temp(n)
end do
! write(6,*)' max,min(ix1temp)=',maxval(ix1temp),minval(ix1temp)
end do
! write(6,*)' max,min(iwgts)=',maxval(iwgts),minval(iwgts)
! write(6,*)' both of above should be >= 1 and <= ',n1grid
! check that we got all points
! dx1max=-huge(dx1max)
! dx1min=huge(dx1min)
! numinside=0
! numgtp5=0
! numoutside=0
! halfpeps=half+epsilon(x)
! do n=1,nin
! if(iflag(n)/=0) then
! numinside=numinside+1
! dx1top=x1in(n)-x1grid(i1ref(n))
! if(dx1top*dx1gridi(i1ref(n))<zero) then
! dx1this=(x1in(n)-x1grid(i1ref(n)-1)) &
! *dx1gridi(i1ref(n)-1)-one
! else
! dx1this=dx1top*dx1gridi(i1ref(n))
! end if
! if(abs(dx1this)>halfpeps) then
! numgtp5=numgtp5+1
! end if
! dx1max=max(dx1this,dx1max)
! dx1min=min(dx1this,dx1min)
! else
! numoutside=numoutside+1
! end if
! end do
! write(6,*)' numinside,outside=',numinside,numoutside
! write(6,*)' numgtp5=',numgtp5
! write(6,*)' dx1min,max=',dx1min,dx1max
! Get taylor matrices and invert
nstart=1
nend=min(nin,ntl)
do while (nstart<=nend)
nthis=nend-nstart+1
if(nf==1) then ! get weights for interpolating function
! get taylor vector(s)
lp=0
x1p(nstart:nend)=one
do ip=0,iord
lp=lp+1
z0(1:nthis,lp)=x1p(nstart:nend)
do n=1,nthis
nn=n+nstart-1
iflip=1
if(ix1sign(nn)<0) iflip=2
x1p(nn)=x1p(nn)*(alocal(i1ref(nn),iflip)* &
(x1in(nn)-x1grid(max(1,i1ref(nn)))) &
+blocal(i1ref(nn),iflip))
end do
end do
! get interpolation weights
do k=1,lbig
wgts(nstart:nend,k)=zero
do j=1,lbig
do n=1,nthis
nn=n+nstart-1
iflip=1
if(ix1sign(nn)<0) iflip=2
wgts(nn,k)=wgts(nn,k)+z0(n,j)*tl(j,k,i1ref(nn),iflip)
end do
end do
end do
end if
if(ndx1==1) then ! get weights for df/dx1
! get taylor vector(s)
lp=0
x1p(nstart:nend)=one
do ip=0,iord
lp=lp+1
do n=1,nthis
nn=n+nstart-1
iflip=1
if(ix1sign(nn)<0) iflip=2
z0(n,lp)=ip*x1p(nn)*alocal(i1ref(nn),iflip)
end do
if(ip>0) then
do n=1,nthis
nn=n+nstart-1
iflip=1
if(ix1sign(nn)<0) iflip=2
x1p(nn)=x1p(nn)*(alocal(i1ref(nn),iflip)* &
(x1in(nn)-x1grid(max(1,i1ref(nn)))) &
+blocal(i1ref(nn),iflip))
end do
end if
end do
! get interpolation weights
do k=1,lbig
wgtsx1(nstart:nend,k)=zero
do j=1,lbig
do n=1,nthis
nn=n+nstart-1
iflip=1
if(ix1sign(nn)<0) iflip=2
wgtsx1(nn,k)=wgtsx1(nn,k)+z0(n,j)*tl(j,k,i1ref(nn),iflip)
end do
end do
end do
end if
if(ndx11==1) then ! get weights for d2f/(dx1*dx1)
! get taylor vector(s)
lp=0
x1p(nstart:nend)=one
do ip=0,iord
lp=lp+1
do n=1,nthis
nn=n+nstart-1
iflip=1
if(ix1sign(nn)<0) iflip=2
z0(n,lp)=ip*(ip-1)*x1p(nn)*alocal(i1ref(nn),iflip)**2
end do
if(ip>1) then
do n=1,nthis
nn=n+nstart-1
iflip=1
if(ix1sign(nn)<0) iflip=2
x1p(nn)=x1p(nn)*(alocal(i1ref(nn),iflip)* &
(x1in(nn)-x1grid(max(1,i1ref(nn)))) &
+blocal(i1ref(nn),iflip))
end do
end if
end do
! get interpolation weights
do k=1,lbig
wgtsx11(nstart:nend,k)=zero
do j=1,lbig
do n=1,nthis
nn=n+nstart-1
iflip=1
if(ix1sign(nn)<0) iflip=2
wgtsx11(nn,k)=wgtsx11(nn,k)+z0(n,j)*tl(j,k,i1ref(nn),iflip)
end do
end do
end do
end if
nstart=nstart+ntl
nend=min(nend+ntl,nin)
end do
return
end subroutine simpin1
subroutine vinvmm(b,a,m,nb,na,ninv,ninv0)
!$$$ subprogram documentation block
! . . . .
! subprogram: vinvmm
! prgmmr: purser org: np20 date: 1998-01-01
!
! abstract: This routine is a vector version of purserlib routine
! invmm. The purpose of this routine is to invert a
! collection of matrices, possibly in place (a=b), using
! the l-u decomposition method
!
! program history log:
! 1998-01-01 purser
! 2000-01-01 parrish
! 2004-06-22 treadon - update documentation
!
! input argument list:
! b - input matrices
! m - order of matrices
! nb - leading dimension of b
! na - leading dimension of a
! ninv - number of matrices to invert
! ninv0 - maximum number of matrices
!
! output argument list:
! a - output inverses of b (can have a=b in calling program)
!
! attributes:
! language: f90
! machine: ibm rs/6000 sp
!
!$$$ end documentation block
use kinds, only: r_kind,i_kind
use constants, only: zero, one
implicit none
integer(i_kind),intent(in ) :: m,nb,na,ninv,ninv0
real(r_kind) ,intent(inout) :: a(ninv0,na,*),b(ninv0,nb,*)
integer(i_kind) ipiv(ninv,m)
real(r_kind) s(ninv),d(ninv)
integer(i_kind) i,j,k,n
do j=1,m
do i=1,m
a(1:ninv,i,j)=b(1:ninv,i,j)
end do
end do
call vlufm(a,ipiv,d,m,na,ninv,ninv0,s)
! invert u in place:
do i=1,m
a(1:ninv,i,i)=one/a(1:ninv,i,i)
end do
do i=1,m-1
do j=i+1,m
s=zero
do k=i,j-1
s(1:ninv)=s(1:ninv)-a(1:ninv,i,k)*a(1:ninv,k,j)
end do
a(1:ninv,i,j)=a(1:ninv,j,j)*s(1:ninv)
end do
end do
! invert l in place assuming implicitly diagonal elements of unity
do j=1,m-1
do i=j+1,m
s(1:ninv)=-a(1:ninv,i,j)
do k=j+1,i-1
s(1:ninv)=s(1:ninv)-a(1:ninv,i,k)*a(1:ninv,k,j)
end do
a(1:ninv,i,j)=s(1:ninv)
end do
end do
! form the product of u**-1 and l**-1 in place
do j=1,m-1
do i=1,j
s(1:ninv)=a(1:ninv,i,j)
do k=j+1,m
s(1:ninv)=s(1:ninv)+a(1:ninv,i,k)*a(1:ninv,k,j)
end do
a(1:ninv,i,j)=s(1:ninv)
end do
do i=j+1,m
s=zero
do k=i,m
s(1:ninv)=s(1:ninv)+a(1:ninv,i,k)*a(1:ninv,k,j)
end do
a(1:ninv,i,j)=s(1:ninv)
end do
end do
! permute columns according to ipiv
do j=m-1,1,-1
do i=1,m
do n=1,ninv
s(n)=a(n,i,j)
a(n,i,j)=a(n,i,ipiv(n,j))
a(n,i,ipiv(n,j))=s(n)
end do
end do
end do
return
end subroutine vinvmm
subroutine vlufm(a,ipiv,d,m,na,ninv,ninv0,s)
!$$$ subprogram documentation block
! . . . .
! subprogram: vlufm
! prgmmr: purser org: np20 date: 1998-01-01
!
! abstract: This routine is a vector version of purserlib routine
! lufm. This routine is a pivot routine used by matrix
! inversion routine vinvmm.
!
! program history log:
! 1998-01-01 purser
! 2000-01-01 parrish
! 2004-06-22 treadon - update documentation
!
! input argument list:
! a - original matrices
! m - order of matrices
! na - leading dimension of matrices
! ninv - number of matrices to invert
! ninv0 - maximum number of matrices
! s - work space of dimension ninv
!
! output argument list:
! a - permuted matrices
! ipiv - record of permutations
! d - signs of determinants of matrices
!
! attributes:
! language: f90
! machine: ibm rs/6000 sp
!
!$$$ end documentation block
use kinds, only: r_kind,i_kind
use constants, only: zero, one
implicit none
integer(i_kind),intent(in ) :: m,na,ninv,ninv0
real(r_kind) ,intent(inout) :: a(ninv0,na,*)
integer(i_kind),intent( out) :: ipiv(ninv,m)
real(r_kind) ,intent( out) :: d(ninv)
real(r_kind) ,intent(inout) :: s(ninv)
real(r_kind) ajj(ninv)
integer(i_kind) i,j,jm,jp,k,n
real(r_kind) aa
d=one
ipiv(1:ninv,m)=m
do j=1,m-1
jp=j+1
ajj(1:ninv)=abs(a(1:ninv,j,j))
ipiv(1:ninv,j)=j
do i=jp,m
do n=1,ninv
aa=abs(a(n,i,j))
if(aa>ajj(n))then
ipiv(n,j)=i
ajj(n)=aa
end if
end do
end do
! swap rows, recording changed sign of determinant
do n=1,ninv
if(ipiv(n,j)/=j) d(n)=-d(n)
end do
do k=1,m
do n=1,ninv
s(n)=a(n,j,k)
a(n,j,k)=a(n,ipiv(n,j),k)
a(n,ipiv(n,j),k)=s(n)
end do
end do
do n=1,ninv
ajj(n)=a(n,j,j)
if(ajj(n)==zero)then
jm=j-1
write(6,100)jm
100 format(' failure in lufact:',/,' matrix singular, rank=',i3)
endif
end do
ajj=one/ajj
do i=jp,m
a(1:ninv,i,j)=ajj(1:ninv)*a(1:ninv,i,j)
do k=jp,m
a(1:ninv,i,k)=a(1:ninv,i,k)-a(1:ninv,i,j)*a(1:ninv,j,k)
end do
end do
end do
return
end subroutine vlufm