pmat.f90

Go to the documentation of this file.

!! Utility routines for various linear inversions and Cholesky.
!!
!! Originally, these routines were copies of the purely "inversion" members
!! of pmat1.f90 (a most extensive collection of matrix routines -- not just
!! inversions). As well as having both single and double precision versions
!! of each routine, these versions also make provision for a more graceful
!! termination in cases where the system matrix is detected to be 
!! essentially singular (and therefore noninvertible). This provision takes
!! the form of an optional "failure flag", FF, which is normally returned
!! as .FALSE., but is returned as .TRUE. when inversion fails.
!! In Sep 2012, these routines were collected together into pmat.f90 so
!! that all the main matrix routines could be in the same library, pmat.a.
!! 
!! @author R. J. Purser
module pmat
use pkind, only: spi,sp,dp,spc,dpc
use pietc, only: t,f
implicit none
private
public:: ldum,udlmm,inv,l1lm,ldlm,invl,invu
interface swpvv;  module procedure sswpvv,dswpvv,cswpvv;        end interface
interface ldum
   module procedure sldum,dldum,cldum,sldumf,dldumf,cldumf;     end interface
interface udlmm
   module procedure sudlmm,dudlmm,cudlmm,sudlmv,dudlmv,cudlmv;  end interface
interface inv
   module procedure                                                           &
sinvmt, dinvmt, cinvmt, slinmmt, dlinmmt, clinmmt, slinmvt, dlinmvt, clinmvt, &
sinvmtf,dinvmtf,cinvmtf,slinmmtf,dlinmmtf,clinmmtf,slinmvtf,dlinmvtf,clinmvtf,&
iinvf
                                                               end interface
interface l1lm;   module procedure sl1lm,dl1lm,sl1lmf,dl1lmf;  end interface
interface ldlm;   module procedure sldlm,dldlm,sldlmf,dldlmf;  end interface
interface invl;   module procedure sinvl,dinvl,slinlv,dlinlv;  end interface
interface invu;   module procedure sinvu,dinvu,slinuv,dlinuv;  end interface

contains

subroutine sswpvv(d,e)!                                                [swpvv]
real(sp),    intent(inout) :: d(:), e(:)
real(sp)                   :: tv(size(d))
tv = d; d = e; e = tv
end subroutine sswpvv

subroutine dswpvv(d,e)!                                                [swpvv]
real(dp), intent(inout) :: d(:), e(:)
real(dp)                :: tv(size(d))
tv = d; d = e; e = tv
end subroutine dswpvv

subroutine cswpvv(d,e)!                                                [swpvv]
complex(dpc),intent(inout) :: d(:), e(:)
complex(dpc)               :: tv(size(d))
tv = d; d = e; e = tv
end subroutine cswpvv

subroutine sinvmt(a)!                                                    [inv]
real(sp),dimension(:,:),intent(INOUT):: a
logical                              :: ff
call sinvmtf(a,ff)
if(ff)stop 'In sinvmt; Unable to invert matrix'
end subroutine sinvmt

subroutine dinvmt(a)!                                                    [inv]
real(dp),dimension(:,:),intent(inout):: a
logical                              :: ff
call dinvmtf(a,ff)
if(ff)stop 'In dinvmt; Unable to invert matrix'
end subroutine dinvmt

subroutine cinvmt(a)!                                                    [inv]
complex(dpc),dimension(:,:),intent(inout):: a
logical                                  :: ff
call cinvmtf(a,ff)
if(ff)stop 'In cinvmt; Unable to invert matrix'
end subroutine cinvmt

subroutine sinvmtf(a,ff)!                                                [inv]
use pietc_s, only: u1
real(sp),dimension(:,:),intent(inout):: a
logical,                intent(  out):: ff 
integer(spi)                     :: m,i,j,jp,l
real(sp)                         :: d
integer(spi),dimension(size(a,1)):: ipiv
m=size(a,1)
if(m /= size(a,2))stop 'In sinvmtf; matrix passed to sinvmtf is not square'
! Perform a pivoted L-D-U decomposition on matrix a:
call sldumf(a,ipiv,d,ff)
if(ff)then
   print '(" In sinvmtf; failed call to sldumf")'
   return
endif
! Invert upper triangular portion U in place:
do i=1,m; a(i,i)=u1/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 sswpvv(a(:,j),a(:,l)); enddo
end subroutine sinvmtf

subroutine dinvmtf(a,ff)!                                                [inv]
real(dp),dimension(:,:),intent(inout):: a
logical,                intent(  out):: ff
integer(spi)                         :: m,i,j,jp,l
real(dp)                             :: d
integer(spi), dimension(size(a,1))   :: ipiv
m=size(a,1)
if(m /= size(a,2))stop 'In inv; matrix passed to dinvmtf is not square'
! Perform a pivoted L-D-U decomposition on matrix a:
call dldumf(a,ipiv,d,ff)
if(ff)then
   print '(" In dinvmtf; failed call to dldumf")'
   return
endif
! Invert upper triangular portion U in place:
do i=1,m; a(i,i)=1_dp/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 dswpvv(a(:,j),a(:,l)); enddo
end subroutine dinvmtf

subroutine cinvmtf(a,ff)!                                                [inv]
use pietc, only: c1
complex(dpc),dimension(:,:),intent(INOUT):: a
logical,                    intent(  OUT):: ff
integer(spi)                     :: m,i,j,jp,l
complex(dpc)                     :: d
integer(spi),dimension(size(a,1)):: ipiv
m=size(a,1)
if(m /= size(a,2))stop 'In inv; matrix passed to cinvmtf is not square'
! Perform a pivoted L-D-U decomposition on matrix a:
call cldumf(a,ipiv,d,ff)
if(ff)then
   print '(" In cinvmtf; failed call to cldumf")'
   return
endif
! Invert upper triangular portion U in place:
do i=1,m; a(i,i)=c1/a(i,i); enddo
do i=1,m-1
   do j=i+1,m; a(i,j)=-a(j,j)*sum(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)-sum(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)+sum(a(jp:m,j)*a(i,jp:m)); enddo
   do i=jp,m; a(i,j)=sum(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 cswpvv(a(:,j),a(:,l)); enddo
end subroutine cinvmtf

subroutine slinmmt(a,b)!                                                 [inv]
real(sp),dimension(:,:),intent(inout):: a,b
logical                              :: ff
call slinmmtf(a,b,ff)
if(ff)stop 'In slinmmt; unable to invert linear system'
end subroutine slinmmt

subroutine dlinmmt(a,b)!                                                 [inv]
real(dp),dimension(:,:),intent(inout):: a,b
logical                              :: ff
call dlinmmtf(a,b,ff)
if(ff)stop 'In dlinmmt; unable to invert linear system'
end subroutine dlinmmt

subroutine clinmmt(a,b)!                                                 [inv]
complex(dpc),dimension(:,:),intent(inout):: a,b
logical                                  :: ff
call clinmmtf(a,b,ff)
if(ff)stop 'In clinmmt; unable to invert linear system'
end subroutine clinmmt

subroutine slinmmtf(a,b,ff)!                                             [inv]
real(sp),   dimension(:,:),intent(inout):: a,b
logical,                   intent(  out):: ff
integer(spi),dimension(size(a,1))       :: ipiv
integer(spi)                            :: m
real(sp)                                :: d
m=size(a,1)
if(m /= size(a,2))stop 'In inv; matrix passed to slinmmtf is not square'
if(m /= size(b,1))&
     stop 'In inv; matrix and vectors in slinmmtf have unmatched sizes'
call sldumf(a,ipiv,d,ff)
if(ff)then
   print '("In slinmmtf; failed call to sldumf")'
   return
endif
call sudlmm(a,b,ipiv)
end subroutine slinmmtf

subroutine dlinmmtf(a,b,ff)!                                             [inv]
real(dp),dimension(:,:),   intent(inout):: a,b
logical,                   intent(  out):: ff
integer(spi),dimension(size(a,1)):: ipiv
integer(spi):: m 
real(dp)    :: d
m=size(a,1)
if(m /= size(a,2))stop 'In inv; matrix passed to dlinmmtf is not square'
if(m /= size(b,1))&
     stop 'In inv; matrix and vectors in dlinmmtf have unmatched sizes'
call dldumf(a,ipiv,d,ff)
if(ff)then
   print '("In dlinmmtf; failed call to dldumf")'
   return
endif
call dudlmm(a,b,ipiv)
end subroutine dlinmmtf

subroutine clinmmtf(a,b,ff)!                                             [inv]
complex(dpc),dimension(:,:),intent(INOUT):: a,b
logical,                    intent(  OUT):: ff
integer(spi),dimension(size(a,1)):: ipiv
integer(spi)                     :: m 
complex(dpc)                     :: d
m=size(a,1)
if(m /= size(a,2))stop 'In inv; matrix passed to dlinmmtf is not square'
if(m /= size(b,1))&
     stop 'In inv; matrix and vectors in dlinmmtf have unmatched sizes'
call cldumf(a,ipiv,d,ff)
if(ff)then
   print '("In clinmmtf; failed call to cldumf")'
   return
endif
call cudlmm(a,b,ipiv)
end subroutine clinmmtf

subroutine slinmvt(a,b)!                                                 [inv]
real(sp),dimension(:,:),intent(inout):: a
real(sp),dimension(:),  intent(inout):: b
logical:: ff
call slinmvtf(a,b,ff)
if(ff)stop 'In slinmvt; matrix singular, unable to continue'
end subroutine slinmvt

subroutine dlinmvt(a,b)!                                                 [inv]
real(dp),dimension(:,:),intent(inout):: a
real(dp),dimension(:),  intent(inout):: b
logical                              :: ff
call dlinmvtf(a,b,ff)
if(ff)stop 'In dlinmvt; matrix singular, unable to continue'
end subroutine dlinmvt

subroutine clinmvt(a,b)!                                                 [inv]
complex(dpc),   dimension(:,:),intent(inout):: a
complex(dpc),   dimension(:),  intent(inout):: b
logical                                     :: ff
call clinmvtf(a,b,ff)
if(ff)stop 'In clinmvt; matrix singular, unable to continue'
end subroutine clinmvt

subroutine slinmvtf(a,b,ff)!                                             [inv]
real(sp),dimension(:,:),intent(inout):: a
real(sp),dimension(:),  intent(inout):: b
logical,                intent(  out):: ff
integer(spi),dimension(size(a,1))    :: ipiv
real(sp)                             :: d
if(size(a,1) /= size(a,2).or. size(a,1) /= size(b))&
     stop 'In inv; In slinmvtf; incompatible array dimensions'
call sldumf(a,ipiv,d,ff)
if(ff)then
   print '("In slinmvtf; failed call to sldumf")'
   return
endif
call sudlmv(a,b,ipiv) 
end subroutine slinmvtf

subroutine dlinmvtf(a,b,ff)!                                             [inv]
real(dp),dimension(:,:),intent(inout):: a
real(dp),dimension(:),  intent(inout):: b
logical,                intent(  out):: ff
integer(spi), dimension(size(a,1))   :: ipiv
real(dp)                             :: d
if(size(a,1) /= size(a,2).or. size(a,1) /= size(b))&
     stop 'In inv; incompatible array dimensions passed to dlinmvtf'
call dldumf(a,ipiv,d,ff)
if(ff)then
   print '("In dlinmvtf; failed call to dldumf")'
   return
endif
call dudlmv(a,b,ipiv)
end subroutine dlinmvtf

subroutine clinmvtf(a,b,ff)!                                             [inv]
complex(dpc),dimension(:,:),intent(inout):: a
complex(dpc),dimension(:),  intent(inout):: b
logical,                    intent(  out):: ff
integer, dimension(size(a,1))            :: ipiv
complex(dpc)                             :: d
if(size(a,1) /= size(a,2).or. size(a,1) /= size(b))&
     stop 'In inv; incompatible array dimensions passed to clinmvtf'
call cldumf(a,ipiv,d,ff)
if(ff)then
   print '("In clinmvtf; failed call to cldumf")'
   return
endif
call cudlmv(a,b,ipiv)
end subroutine clinmvtf

subroutine iinvf(imat,ff)!                                               [inv]
integer(spi),dimension(:,:),intent(INOUT):: imat
logical,                    intent(  OUT):: ff
real(dp),parameter                           :: eps=1.e-6_dp
real(dp),dimension(size(imat,1),size(imat,1)):: dmat
integer(spi)                                 :: m,i,j
m=size(imat,1)
if(m /= size(imat,2))stop 'In inv; matrix passed to iinvf is not square'
dmat=imat; call inv(dmat,ff)
if(.not.ff)then
   do j=1,m
      do i=1,m
         imat(i,j)=nint(dmat(i,j)); if(abs(dmat(i,j)-imat(i,j))>eps)ff=t
      enddo
   enddo
endif
end subroutine iinvf

subroutine sldum(a,ipiv,d)!                                             [ldum]
real(sp),    intent(inout) :: a(:,:) 
real(sp),    intent(  out) :: d
integer(spi),intent(  out) :: ipiv(:)
logical:: ff
call sldumf(a,ipiv,d,ff)
if(ff)stop 'In sldum; matrix singular, unable to continue'
end subroutine sldum

subroutine dldum(a,ipiv,d)!                                             [ldum]
real(dp),    intent(inout) :: a(:,:) 
real(dp),    intent(  out) :: d
integer(spi),intent(  out) :: ipiv(:)
logical:: ff
call dldumf(a,ipiv,d,ff)
if(ff)stop 'In dldum; matrix singular, unable to continue'
end subroutine dldum

subroutine cldum(a,ipiv,d)!                                             [ldum]
complex(dpc),intent(inout) :: a(:,:) 
complex(dpc),intent(out  ) :: d
integer(spi),intent(out  ) :: ipiv(:)
logical:: ff
call cldumf(a,ipiv,d,ff)
if(ff)stop 'In cldum; matrix singular, unable to continue'
end subroutine cldum

subroutine sldumf(a,ipiv,d,ff)!                                         [ldum]
use pietc_s,only: u0,u1
real(sp),    intent(inout) :: a(:,:) 
real(sp),    intent(  out) :: d
integer(spi),intent(  out) :: ipiv(:)
logical,     intent(  out) :: ff
integer(spi):: m,i, j, jp, ibig, jm
real(sp)    :: s(size(a,1)),  aam, aa, abig,  ajj, ajji, aij
ff=f
m=size(a,1)
do i=1,m
  aam=u0
  do j=1,m
    aa=abs(a(i,j))
    if(aa > aam)aam=aa
  enddo
  if(aam == u0)then
    print '("In sldumf; row ",i6," of matrix vanishes")',i
    ff=t
    return
  endif
  s(i)=u1/aam
enddo
d=1_sp
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 sswpvv(a(j,:),a(ibig,:))
      s(ibig)=s(j)
   endif
   ajj=a(j,j)
   if(ajj == u0)then
      jm=j-1
      print '(" failure in sldumf:"/" matrix singular, rank=",i3)',jm
      ff=t
      return
   endif
   ajji=u1/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 sldumf

subroutine dldumf(a,ipiv,d,ff)!                                         [ldum]
use pietc, only: u0,u1
real(dp),    intent(inout) :: a(:,:) 
real(dp),    intent(  out) :: d
integer,     intent(  out) :: ipiv(:)
logical(spi),intent(  out) :: ff
integer(spi)               :: m,i, j, jp, ibig, jm
real(dp)                   :: s(size(a,1)),  aam, aa, abig,  ajj, ajji, aij
ff=f
m=size(a,1)
do i=1,m
   aam=u0
   do j=1,m
      aa=abs(a(i,j))
      if(aa > aam)aam=aa
   enddo
   if(aam == u0)then
      print '("In dldumf;  row ",i6," of matrix vanishes")',i
      ff=t
      return
   endif
   s(i)=u1/aam
enddo
d=u1
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 dswpvv(a(j,:),a(ibig,:))
      s(ibig)=s(j)
   endif
   ajj=a(j,j)
   if(ajj == u0)then
      jm=j-1
      print '(" Failure in dldumf:"/" matrix singular, rank=",i3)',jm
      ff=t
      return
   endif
   ajji=u1/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 dldumf

subroutine cldumf(a,ipiv,d,ff)!                                         [ldum]
use pietc, only: u0,u1,c0,c1
complex(dpc), intent(inout)  :: a(:,:) 
complex(dpc), intent(  out)  :: d
integer(spi), intent(  out)  :: ipiv(:)
logical,      intent(  out)  :: ff
integer(spi)                 :: m,i, j, jp, ibig, jm
complex(dpc)                 :: ajj, ajji, aij
real(dp)                     :: aam,aa,abig
real(dp),dimension(size(a,1)):: s
ff=f
m=size(a,1)
do i=1,m
   aam=u0
   do j=1,m
      aa=abs(a(i,j))
      if(aa > aam)aam=aa
   enddo
   if(aam == u0)then
      print '("In cldumf;  row ",i6," of matrix vanishes")',i
      ff=t
      return
   endif
   s(i)=u1/aam
enddo
d=c1
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 cswpvv(a(j,:),a(ibig,:))
      s(ibig)=s(j)
   endif
   ajj=a(j,j)
   if(ajj == c0)then
      jm=j-1
      print '(" Failure in cldumf:"/" matrix singular, rank=",i3)',jm
      ff=t
      return
   endif
   ajji=c1/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 cldumf

subroutine sudlmm(a,b,ipiv)!                                           [udlmm]
use pietc_s, only: u1
integer(spi),dimension(:),  intent(in)    :: ipiv 
real(sp),    dimension(:,:),intent(in)    :: a 
real(sp),    dimension(:,:),intent(inout) :: b 
integer(spi):: m,i, k, l
real(sp)    :: s,aiii
m=size(a,1)
do k=1,size(b,2) !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=u1/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 sudlmm

subroutine dudlmm(a,b,ipiv)!                                           [udlmm]
use pietc, only: u1
integer(spi),dimension(:),  intent(in   ) :: ipiv 
real(dp),    dimension(:,:),intent(in   ) :: a 
real(dp),    dimension(:,:),intent(inout) :: b 
integer(spi):: m,i, k, l
real(dp)    :: s,aiii
m=size(a,1)
do k=1, size(b,2)!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=u1/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 cudlmm(a,b,ipiv)!                                           [udlmm]
use pietc, only: c1
integer(spi),dimension(:),  intent(in   ) :: ipiv 
complex(dpc),dimension(:,:),intent(in   ) :: a 
complex(dpc),dimension(:,:),intent(inout) :: b 
integer(spi):: m,i, k, l
complex(dpc):: s,aiii
m=size(a,1)
do k=1, size(b,2)!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=c1/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 cudlmm

subroutine sudlmv(a,b,ipiv)!                                           [udlmv]
use pietc_s, only: u1
integer(spi),dimension(:),  intent(in   ):: ipiv 
real(sp),    dimension(:,:),intent(in   ):: a 
real(sp),    dimension(:),  intent(inout):: b 
integer(spi):: m,i, l
real(sp)    :: 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=u1/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 sudlmv

subroutine dudlmv(a,b,ipiv)!                                           [udlmv]
use pietc, only: u1
integer(spi),dimension(:),  intent(in   ) :: ipiv(:) 
real(dp),    dimension(:,:),intent(in   ) :: a(:,:) 
real(dp),    dimension(:),  intent(inout) :: b(:) 
integer(spi):: m,i, l
real(dp)    :: 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=u1/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 cudlmv(a,b,ipiv)!                                           [udlmv]
use pietc, only: c1
integer(spi),dimension(:),  intent(in   ) :: ipiv(:) 
complex(dpc),dimension(:,:),intent(in   ) :: a(:,:) 
complex(dpc),dimension(:),  intent(inout) :: b(:) 
integer(spi):: m,i, l
complex(dpc):: 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=c1/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 cudlmv

subroutine sl1lm(a,b) !                                                 [l1lm]
real(sp),intent(in   ):: a(:,:)
real(sp),intent(inout):: b(:,:)
logical:: ff
call sl1lmf(a,b,ff)
if(ff)stop 'In sl1lm; matrix singular, unable to continue'
end subroutine sl1lm

subroutine dl1lm(a,b) !                                                 [l1lm]
real(dp),intent(in   ):: a(:,:)
real(dp),intent(inout):: b(:,:)
logical:: ff
call dl1lmf(a,b,ff)
if(ff)stop 'In dl1lm; matrix singular, unable to continue'
end subroutine dl1lm

subroutine sl1lmf(a,b,ff)!                                              [L1Lm] 
use pietc_s, only: u0
real(sp),intent(in   ):: a(:,:)
real(sp),intent(inout):: b(:,:)
logical, intent(  out):: ff
integer(spi):: m,j, jm, jp, i
real(sp)    :: s, bjji
m=size(a,1)
ff=f
do j=1,m
   jm=j-1
   jp=j+1
   s = a(j,j) - sum(b(j,1:jm)*b(j,1:jm))
   ff=(s <= u0)
   if(ff)then
      print '("sL1Lmf detects nonpositive a, rank=",i6)',jm
      return
   endif
   b(j,j)=sqrt(s)
   bjji=1_sp/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) = u0
enddo
end subroutine sl1lmf

subroutine dl1lmf(a,b,ff) !                                             [L1Lm]
use pietc, only: u0,u1
real(dp),intent(in   ) :: a(:,:) 
real(dp),intent(inout) :: b(:,:) 
logical, intent(  out) :: ff
integer(spi):: m,j, jm, jp, i
real(dp)    :: s, bjji
m=size(a,1)
ff=f
do j=1,m
  jm=j-1
  jp=j+1
  s = a(j,j) - sum(b(j,1:jm)*b(j,1:jm))
  ff=(s <= u0)
  if(ff)then
     print '("dL1LMF detects nonpositive A, rank=",i6)',jm
     return
  endif
  b(j,j)=sqrt(s)
  bjji=u1/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) = u0
enddo
end subroutine dl1lmf

subroutine sldlm(a,b,d)!                                                [LdLm]
real(sp),intent(in   ):: a(:,:)
real(sp),intent(inout):: b(:,:)
real(sp),intent(  out):: d(:)
logical:: ff
call sldlmf(a,b,d,ff)
if(ff)stop 'In sldlm; matrix singular, unable to continue'
end subroutine sldlm

subroutine dldlm(a,b,d)!                                                [LdLm]
real(dp),intent(in   ):: a(:,:)
real(dp),intent(inout):: b(:,:)
real(dp),intent(  out):: d(:)
logical:: ff
call dldlmf(a,b,d,ff)
if(ff)stop 'In dldlm; matrix singular, unable to continue'
end subroutine dldlm

subroutine sldlmf(a,b,d,ff) !                                           [LDLM]
use pietc_s, only: u0,u1
real(sp), intent(in   ):: a(:,:)
real(sp), intent(inout):: b(:,:)
real(sp), intent(  out):: d(:)
logical,  intent(  out):: ff
integer(spi):: m,j, jm, jp, i
real(sp)    :: bjji
m=size(a,1)
ff=f
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) = u1
  ff=(d(j) == u0)
  if(ff)then
     print '("In sldlmf; singularity of matrix detected")'
     print '("Rank of matrix: ",i6)',jm
     return
  endif
  bjji=u1/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
  b(1:jm,j)=u0
enddo
end subroutine sldlmf

subroutine dldlmf(a,b,d,ff) !                                           [LDLM]
use pietc, only: u0,u1
real(dp), intent(IN   ) :: a(:,:)
real(dp), intent(INOUT) :: b(:,:)
real(dp), intent(  OUT) :: d(:)
logical,  intent(  OUT) :: ff
integer(spi):: m,j, jm, jp, i
real(dp)    :: bjji
m=size(a,1)
ff=f
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) = 1
  ff=(d(j) == u0)
  if(ff)then
     print '("In dldlmf; singularity of matrix detected")'
     print '("Rank of matrix: ",i6)',jm
     return
  endif
  bjji=u1/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
  b(1:jm,j)=u0
enddo
end subroutine dldlmf

subroutine sinvu(a)!                                                     [invu]
real(sp),dimension(:,:),intent(inout):: a
a=transpose(a); call sinvl(a); a=transpose(a)
end subroutine sinvu

subroutine dinvu(a)!                                                     [invu]
real(dp),dimension(:,:),intent(inout):: a
a=transpose(a); call dinvl(a); a=transpose(a)
end subroutine dinvu

subroutine sinvl(a)!                                                     [invl]
use pietc_s, only: u0,u1
real(sp), intent(inout) :: a(:,:) 
integer(spi):: m,j, i
m=size(a,1)
do j=m,1,-1
   a(1:j-1,j) = u0
   a(j,j)=u1/a(j,j)
   do i=j+1,m
      a(i,j)=-a(i,i)*sum(a(j:i-1,j)*a(i,j:i-1))
   enddo
enddo
end subroutine sinvl

subroutine dinvl(a)!                                                     [invl]
use pietc, only: u0,u1
real(dp), intent(inout) :: a(:,:) 
integer(spi):: m,j, i
m=size(a,1)
do j=m,1,-1
   a(1:j-1,j) = u0
   a(j,j)=u1/a(j,j)
   do i=j+1,m
      a(i,j)=-a(i,i)*sum(a(j:i-1,j)*a(i,j:i-1))
   enddo
enddo
end subroutine dinvl

subroutine slinlv(a,u)!                                                  [invl]
real(sp),intent(in   ) :: a(:,:)
real(sp),intent(inout) :: u(:)
integer(spi):: i
if(size(a,1) /= size(a,2) .or. size(a,1) /= size(u))&
     stop 'In slinlv; incompatible array dimensions'
do i=1,size(u); u(i)=(u(i) - sum(u(:i-1)*a(i,:i-1)))/a(i,i); enddo
end subroutine slinlv

subroutine dlinlv(a,u)!                                                  [invl]
real(dp),intent(in   ) :: a(:,:)
real(dp),intent(inout) :: u(:)
integer(spi):: i
if(size(a,1) /= size(a,2) .or. size(a,1) /= size(u))&
     stop 'In dlinlv; incompatible array dimensions'
do i=1,size(u); u(i)=(u(i) - sum(u(:i-1)*a(i,:i-1)))/a(i,i); enddo
end subroutine dlinlv

subroutine slinuv(a,u)!                                                  [invu]
real(sp),intent(in   ) :: a(:,:)
real(sp),intent(inout) :: u(:)
integer(spi):: i
if(size(a,1) /= size(a,2) .or. size(a,1) /= size(u))&
     stop 'In linuv; incompatible array dimensions'
do i=size(u),1,-1; u(i)=(u(i) - sum(a(i+1:,i)*u(i+1:)))/a(i,i); enddo
end subroutine slinuv

subroutine dlinuv(a,u)!                                                  [invu]
real(dp), intent(in   ) :: a(:,:)
real(dp), intent(inout) :: u(:)
integer(spi)            :: i
if(size(a,1) /= size(a,2) .or. size(a,1) /= size(u))&
     stop 'In dlinuv; incompatible array dimensions'
do i=size(u),1,-1; u(i)=(u(i) - sum(a(i+1:,i)*u(i+1:)))/a(i,i); enddo
end subroutine dlinuv

end module pmat