pmat2.f90

Go to the documentation of this file.

 
 
module pmat2
!============================================================================
use    pkind, only: spi,sp,dp,dpc
implicit none
private
public:: avco 
public:: dfco 
public:: dfco2 
public:: clipb 
public:: cad1b 
public:: csb1b 
public:: cad2b 
public:: csb2b 
public:: ldub 
public:: ldltb 
public:: udlb 
public:: l1ubb 
public:: l1ueb 
public:: ltdlbv 
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 
real(dp),parameter:: zero=0
 
interface avco;   module procedure avco,  davco,  tavco;   end interface
interface dfco;   module procedure dfco,  ddfco,  tdfco;   end interface
interface dfco2;  module procedure dfco2, ddfco2, tdfco2;  end interface
interface clipb;  module procedure clib,  clib_d, clib_c;  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 ldub;   module procedure ldub,  dldub;           end interface
interface ldltb;  module procedure ldltb, dldltb;          end interface
interface l1ubb;  module procedure l1ubb, dl1ubb;          end interface
interface l1ueb;  module procedure l1ueb, dl1ueb;          end interface
interface ltdlbv; module procedure ltdlbv,dltdlbv;         end interface
interface udlb;   module procedure udlb,  dudlb;           end interface
interface udlbv;  module procedure udlbv, dudlbv;          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 avco(na,nb,za,zb,z0,a,b) !                                   [AVCO]
use pietc, only: u0,u1
use pmat,  only: inv
implicit none
integer(spi),intent(in ):: na,nb
real(sp),    intent(in ):: za(na),zb(nb),z0
real(sp),    intent(out):: a(na),b(nb)
!-----------------------------------------------------------------------------
integer(spi)                    :: na1,nab,i
real(sp), dimension(na+nb,na+nb):: w
real(sp), dimension(na)         :: za0,pa
real(sp), dimension(nb)         :: zb0,pb
real(sp), dimension(na+nb)      :: ab
!=============================================================================
na1=na+1;     nab=na+nb
za0=za-z0;    zb0=zb-z0
pa=u1;        pb=-u1
w=u0;         ab=u0
w(1,1:na)=u1; ab(1)=u1
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 davco(na,nb,za,zb,z0,a,b) !                                  [AVCO]
use pietc, only: u0,u1
use pmat, only: inv
implicit none
integer(spi),intent(IN ):: na,nb
real(dp),    intent(IN ):: za(na),zb(nb),z0
real(dp),    intent(OUT):: a(na),b(nb)
!-----------------------------------------------------------------------------
integer(spi)                   :: na1,nab,i
real(dp),dimension(na+nb,na+nb):: w
real(dp),dimension(na)         :: za0,pa
real(dp),dimension(nb)         :: zb0,pb
real(dp),dimension(na+nb)      :: ab
!=============================================================================
na1=na+1;     nab=na+nb
za0=za-z0;    zb0=zb-z0
pa=u1;        pb=-u1
w=u0;         ab=u0
w(1,1:na)=u1; ab(1)=u1
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 davco
 
subroutine tavco(xa,xb,a,b)!                                            [AVCO]
implicit none
real(dp),dimension(:),intent(IN ):: xa,xb
real(dp),dimension(:),intent(OUT):: a,b
!-----------------------------------------------------------------------------
integer(spi):: na,nb
!=============================================================================
na=size(xa); if(na /= size(a))stop 'In tavco; sizes of a and xa different'
nb=size(xb); if(nb /= size(b))stop 'In tavco; sizes of b and xb different'
call davco(na,nb,xa,xb,zero,a,b)
end subroutine tavco
 
subroutine dfco(na,nb,za,zb,z0,a,b)!                                    [DFCO]
use pietc_s, only: u0,u1
use pmat,    only: inv
implicit none
integer(spi),intent(IN )        :: na,nb
real(sp),    intent(IN )        :: za(na),zb(nb),z0
real(sp),    intent(OUT)        :: a(na),b(nb)
!-----------------------------------------------------------------------------
integer(spi):: na1,nab,i
real(sp), dimension(na+nb,na+nb):: w
real(sp), dimension(na)         :: za0,pa
real(sp), dimension(nb)         :: zb0,pb
real(sp), dimension(na+nb)      :: ab
!=============================================================================
na1=na+1; nab=na+nb
za0=za-z0; zb0=zb-z0
pa=u1;     pb=-u1
w=u0;         ab=u0
w(1,1:na)=u1; ab(1)=u1
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 ddfco(na,nb,za,zb,z0,a,b) ! Real(dp) version of              [DFCO]
use pietc, only: u0,u1
use pmat,  only: inv
implicit none
integer(spi),intent(in) :: na,nb
real(dp),    intent(in) :: za(na),zb(nb),z0
real(dp),    intent(out):: a(na),b(nb)
!-----------------------------------------------------------------------------
integer(spi)                    :: na1,nab,i
real(dp), dimension(na+nb,na+nb):: w
real(dp), dimension(na)         :: za0,pa
real(dp), dimension(nb)         :: zb0,pb
real(dp), dimension(na+nb)      :: ab
!=============================================================================
na1=na+1; nab=na+nb
za0=za-z0; zb0=zb-z0
pa=u1;     pb=-u1
w=u0;         ab=u0
w(1,1:na)=u1; ab(1)=u1
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 ddfco 
 
subroutine tdfco(xa,xb,a,b)!                                            [DFCO]
implicit none
real(dp),dimension(:),intent(IN ):: xa,xb
real(dp),dimension(:),intent(OUT):: a,b
!-----------------------------------------------------------------------------
integer(spi):: na,nb
!=============================================================================
na=size(xa); if(na /= size(a))stop 'In tdfco; sizes of a and xa different'
nb=size(xb); if(nb /= size(b))stop 'In tdfco; sizes of b and xb different'
call ddfco(na,nb,xa,xb,zero,a,b)
end subroutine tdfco
 
subroutine dfco2(na,nb,za,zb,z0,a,b)!                                  [DFCO2] 
use pietc_s, only: u0,u1
use pmat, only: inv
implicit none
integer(spi), intent(IN ):: na,nb
real(sp),     intent(IN ):: za(na),zb(nb),z0
real(sp),     intent(OUT):: a(na),b(nb)
!-----------------------------------------------------------------------------
integer(spi)                    :: na1,nab,i
real(sp), dimension(na+nb,na+nb):: w
real(sp), dimension(na)         :: za0,pa
real(sp), dimension(nb)         :: zb0,pb
real(sp), dimension(na+nb)      :: ab
!=============================================================================
na1=na+1; nab=na+nb
za0=za-z0; zb0=zb-z0
pa=u1;     pb=-u1
w=u0;         ab=u0
w(1,1:na)=u1; ab(1)=u1
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 ddfco2(na,nb,za,zb,z0,a,b) ! Real(dp) version of            [DFCO2]
use pietc, only: u0,u1
use pmat, only: inv
implicit none
integer(spi),intent(IN )           :: na,nb
real(dp),    intent(IN )           :: za(na),zb(nb),z0
real(dp),    intent(OUT)           :: a(na),b(nb)
!-----------------------------------------------------------------------------
integer(spi)                    :: na1,nab,i
real(dp), dimension(na+nb,na+nb):: w
real(dp), dimension(na)         :: za0,pa
real(dp), dimension(nb)         :: zb0,pb
real(dp), dimension(na+nb)      :: ab
!=============================================================================
na1=na+1; nab=na+nb
za0=za-z0; zb0=zb-z0
pa=u1;     pb=-u1
w=u0;         ab=u0
w(1,1:na)=u1; ab(1)=u1
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 ddfco2 
 
subroutine tdfco2(xa,xb,a,b)!                                          [DFCO2]
!=============================================================================
real(dp),dimension(:),intent(IN ):: xa,xb
real(dp),dimension(:),intent(OUT):: a,b
!-----------------------------------------------------------------------------
integer(spi):: na,nb
!=============================================================================
na=size(xa); if(na /= size(a))stop 'In tdfco2; sizes of a and xa different'
nb=size(xb); if(nb /= size(b))stop 'In tdfco2; sizes of b and xb different'
call ddfco2(na,nb,xa,xb,zero,a,b)
end subroutine tdfco2
 
pure subroutine clib(m1,m2,mah1,mah2,a)!                               [CLIPB]
use pietc_s, only: u0
implicit none
integer(spi), intent(IN   ) :: m1, m2, mah1, mah2
real(sp),     intent(INOUT) :: a(m1,-mah1:mah2)
integer(spi):: j
do j=1,mah1; a(1:min(m1,j),-j)=u0; enddo
do j=m2-m1+1,mah2; a(max(1,m2-j+1):m1,j)=u0; enddo
end subroutine clib
 
pure subroutine clib_d(m1,m2,mah1,mah2,a)!                             [CLIPB]
use pietc, only: u0
implicit none
integer(spi),intent(IN   ) :: m1, m2, mah1, mah2
real(dp),    intent(INOUT) :: a(m1,-mah1:mah2)
integer(spi):: j
do j=1,mah1; a(1:min(m1,j),-j)=u0; enddo
do j=m2-m1+1,mah2; a(max(1,m2-j+1):m1,j)=u0; enddo
end subroutine clib_d
 
pure subroutine clib_c(m1,m2,mah1,mah2,a)!                              [CLIPB]
use pietc, only: c0
implicit none
integer(spi), intent(IN   ) :: m1, m2, mah1, mah2
complex(dpc), intent(INOUT) :: a(m1,-mah1:mah2)
integer(spi):: j
do j=1,mah1; a(1:min(m1,j),-j)=c0; enddo
do j=m2-m1+1,mah2; a(max(1,m2-j+1):m1,j)=c0; enddo
end subroutine clib_c
 
subroutine cad1b(m1,mah1,mah2,mirror2,a)!                              [CAD1B]
use pietc_s, only: u0
implicit none
integer(spi),intent(IN   ):: m1,mah1,mah2,mirror2
real(sp),    intent(INOUT):: a(0:m1-1,-mah1:mah2)
!-----------------------------------------------------------------------------
integer(spi):: 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)=u0              ! zero the exterior part
   enddo
enddo
end subroutine cad1b
 
subroutine csb1b(m1,mah1,mah2,mirror2,a)!                              [CSB1B]
use pietc_s, only: u0
implicit none
integer(spi),intent(IN   ):: m1,mah1,mah2,mirror2
real(sp),    intent(INOUT):: a(0:m1-1,-mah1:mah2)
!-----------------------------------------------------------------------------
integer(spi):: i,i2,jm,jp,jpmax
!=============================================================================
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)=u0              ! zero the exterior part
   enddo
enddo
end subroutine csb1b
 
subroutine cad2b(m1,m2,mah1,mah2,mirror2,a)!                           [CAD2B]
use pietc_s, only: u0
implicit none
integer(spi),intent(IN   ):: m1,m2,mah1,mah2,mirror2
real(sp),    intent(INOUT):: a(1-m1:0,m1-m2-mah1:m1-m2+mah2)
!-----------------------------------------------------------------------------
integer(spi):: i,i2,jm,jp,jmmin,nah1,nah2
!=============================================================================
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)=u0              ! zero the exterior part
   enddo
enddo
end subroutine cad2b
 
subroutine csb2b(m1,m2,mah1,mah2,mirror2,a)!                           [CSB2B]
use pietc_s, only: u0
implicit none
integer(spi),intent(IN   ):: m1,m2,mah1,mah2,mirror2
real(sp),    intent(INOUT):: a(1-m1:0,m1-m2-mah1:m1-m2+mah2)
!-----------------------------------------------------------------------------
integer(spi):: i,i2,jm,jp,jmmin,nah1,nah2
!=============================================================================
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)=u0              ! zero the exterior part
   enddo
enddo
end subroutine csb2b
 
subroutine ldub(m,mah1,mah2,a)!                                         [LDUB]
use pietc_s, only: u0,u1
implicit none
integer(spi),intent(IN   ):: m,mah1, mah2 
real(sp),    intent(INOUT):: a(m,-mah1:mah2) 
!-----------------------------------------------------------------------------
integer(spi):: j, imost, jmost, jp, i
real(sp)    :: 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 == u0)then
      print '(" Failure in LDUB:"/" Matrix requires pivoting or is singular")'
      stop
   endif
   ajji=u1/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,1:jmost-j)
   enddo
   a(j,1:jmost-j)=ajji*a(j,1:jmost-j)
enddo
end subroutine ldub
 
subroutine dldub(m,mah1,mah2,a) ! Real(dp) version of                   [LDUB]
use pietc, only: u0,u1
implicit none
integer(spi),intent(IN   ):: m,mah1, mah2 
real(dp),    intent(INOUT):: a(m,-mah1:mah2) 
!-----------------------------------------------------------------------------
integer(spi):: j, imost, jmost, jp, i
real(dp)    :: 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 == u0)then
    print '(" Fails in LDUB_d:"/" Matrix requires pivoting or is singular")'
    stop
  endif
  ajji=u1/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,1:jmost-j)
  enddo
  a(j,1:jmost-j)=ajji*a(j,1:jmost-j)
enddo
end subroutine dldub
 
subroutine ldltb(m,mah1,a) ! Real(sp) version of                       [LDLTB]
use pietc_s, only: u0,u1
integer(spi),intent(IN   ):: m,mah1
real(sp),    intent(INOUT):: a(m,-mah1:0) 
!-----------------------------------------------------------------------------
integer(spi):: j, imost, jp, i,k
real(sp)    :: ajj, ajji, aij
!=============================================================================
do j=1,m
  imost=min(m,j+mah1)
  jp=j+1
  ajj=a(j,0)
  if(ajj == u0)then
    print '(" Fails in LDLTB:"/" Matrix requires pivoting or is singular")'
    stop
  endif
  ajji=u1/ajj
  a(j,0)=ajji
  do i=jp,imost
    aij=a(i,j-i)
    a(i,j-i)=ajji*aij
    do k=jp,i
       a(i,k-i)=a(i,k-i)-aij*a(k,j-k)
    enddo
  enddo
enddo
end subroutine ldltb
 
subroutine dldltb(m,mah1,a) ! Real(dp) version of                   [LDLTB]
use pietc, only: u0,u1
integer(spi),intent(IN   ) :: m,mah1
real(dp),    intent(INOUT) :: a(m,-mah1:0) 
!-----------------------------------------------------------------------------
integer(spi):: j, imost, jp, i,k
real(dp)    :: ajj, ajji, aij
!=============================================================================
do j=1,m
   imost=min(m,j+mah1)
   jp=j+1
   ajj=a(j,0)
   if(ajj == u0)then
      print '(" Fails in LDLTB_d:"/" Matrix requires pivoting or is singular")'
      stop
   endif
   ajji=u1/ajj
   a(j,0)=ajji
   do i=jp,imost
      aij=a(i,j-i)
      a(i,j-i)=ajji*aij
      do k=jp,i
         a(i,k-i)=a(i,k-i)-aij*a(k,j-k)
      enddo
   enddo
enddo
end subroutine dldltb
 
subroutine udlb(m,mah1,mah2,a) ! Reversed-index version of ldub        [UDLB]
implicit none
integer(spi),                    intent(IN   ) :: m,mah1,mah2
real(sp),dimension(m,-mah1:mah2),intent(INOUT) :: a(m,-mah1:mah2)
!-----------------------------------------------------------------------------
real(sp),dimension(m,-mah2:mah1):: at
!=============================================================================
at=a(m:1:-1,mah2:-mah1:-1); call ldub(m,mah2,mah1,at)
a=at(m:1:-1,mah1:-mah2:-1)
end subroutine udlb 
 
subroutine dudlb(m,mah1,mah2,a) ! real(dp) version of udlb              [UDLB]
implicit none
integer(spi),                    intent(IN   ) :: m,mah1,mah2
real(dp),dimension(m,-mah1:mah2),intent(INOUT) :: a(m,-mah1:mah2)
!-----------------------------------------------------------------------------
real(dp),dimension(m,-mah2:mah1):: at
!=============================================================================
at=a(m:1:-1,mah2:-mah1:-1); call dldub(m,mah2,mah1,at)
a=at(m:1:-1,mah1:-mah2:-1)
end subroutine dudlb 
 
subroutine l1ubb(m,mah1,mah2,mbh1,mbh2,a,b)!                           [L1UBB] 
use pietc_s, only: u0,u1
implicit none
integer(spi), intent(IN   ) ::  m,mah1, mah2, mbh1, mbh2 
real(sp),     intent(INOUT) :: a(m,-mah1:mah2), b(m,-mbh1:mbh2)
!-----------------------------------------------------------------------------
integer(spi):: j, imost, jmost, jleast, jp, i
real(sp)    :: 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 == u0)stop 'In L1UBB; zero element found in diagonal factor'
   ajji=u1/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)=u1
   b(j,-mbh1:mbh2) = ajji * b(j,-mbh1:mbh2)
enddo
end subroutine l1ubb
 
subroutine dl1ubb(m,mah1,mah2,mbh1,mbh2,a,b) ! Real(dp) version of     [L1UBB]
use pietc, only: u0,u1
implicit none
integer(spi),intent(IN   ) :: m,mah1, mah2, mbh1, mbh2 
real(dp),    intent(INOUT) :: a(m,-mah1:mah2), b(m,-mbh1:mbh2)
!-----------------------------------------------------------------------------
integer(spi):: j, imost, jmost, jleast, jp, i
real(dp)    :: 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 == u0)stop 'In L1UBB_d; zero element found in diagonal factor'
   ajji=u1/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)=u1
   b(j,-mbh1:mbh2) = ajji * b(j,-mbh1:mbh2)
enddo
end subroutine dl1ubb
 
subroutine l1ueb(m,mah1,mah2,mbh1,mbh2,a,b)!                           [L1UEB]
use pietc_s, only: u0,u1
implicit none
integer(spi),intent(IN   ) :: m,mah1, mah2, mbh1, mbh2 
real(sp),    intent(INOUT) :: a(0:m,-mah1:mah2), b(m,-mbh1:mbh2)
!-----------------------------------------------------------------------------
integer(spi):: j, imost, jmost, jleast, jp, i
real(sp)    :: 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 == u0)stop 'In L1UEB; zero element found in diagonal factor'
   ajji=u1/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)=u1
   b(j,-mbh1:mbh2) = ajji * b(j,-mbh1:mbh2)
enddo
end subroutine l1ueb
 
subroutine dl1ueb(m,mah1,mah2,mbh1,mbh2,a,b) ! Real(dp) version of     [L1UEB]
use pietc, only: u0,u1
implicit none
integer(spi),intent(IN   ):: m,mah1, mah2, mbh1, mbh2 
real(dp),    intent(INOUT):: a(0:,-mah1:), b(:,-mbh1:)
!-----------------------------------------------------------------------------
integer(spi):: j, imost, jmost, jleast, jp, i
real(dp)    :: 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 == u0)stop 'In L1UEB_D; zero element found in diagonal factor'
   ajji=u1/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)=u1
   b(j,-mbh1:mbh2) = ajji * b(j,-mbh1:mbh2)
enddo
end subroutine dl1ueb
 
subroutine udlbv(m,mah1,mah2,a,v)!                                    [UDLBV]
implicit none
integer(spi),intent(IN   ):: m, mah1, mah2
real(sp),    intent(IN   ):: a(m,-mah1:mah2)
real(sp),    intent(INOUT):: v(m)
!-----------------------------------------------------------------------------
integer(spi):: i, j
real(sp)    :: 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 dudlbv(m,mah1,mah2,a,v)!                                    [udlbv]
implicit none
integer(spi),intent(IN   ) :: m, mah1, mah2
real(dp),    intent(IN   ) :: a(m,-mah1:mah2)
real(dp),    intent(INOUT) :: v(m)
!-----------------------------------------------------------------------------
integer(spi):: i, j
real(dp)    :: 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 dudlbv
 
subroutine ltdlbv(m,mah1,a,v)!                                        [ltdlbv]
implicit none
integer(spi),intent(IN   ) :: m, mah1
real(sp),    intent(IN   ) :: a(m,-mah1:0)
real(sp),    intent(INOUT) :: v(m)
!-----------------------------------------------------------------------------
integer(spi):: i, j
real(sp)    :: 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-mah1),j-1; v(i)=v(i)-a(j,i-j)*vj; enddo
enddo
end subroutine ltdlbv
 
 
subroutine dltdlbv(m,mah1,a,v)!                                       [ltdlbv]
implicit none
integer(spi),intent(IN   ) :: m, mah1
real(dp),    intent(IN   ) :: a(m,-mah1:0)
real(dp),    intent(INOUT) :: v(m)
!-----------------------------------------------------------------------------
integer(spi):: i, j
real(dp)    :: 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-mah1),j-1; v(i)=v(i)-a(j,i-j)*vj; enddo
enddo
end subroutine dltdlbv
 
subroutine udlbx(mx,mah1,mah2,my,a,v)!                                [UDLBX]
implicit none
integer(spi),intent(IN   ) :: mx, mah1, mah2, my
real(sp),    intent(IN   ) :: a(mx,-mah1:mah2)
real(sp),    intent(INOUT) :: v(mx,my)
!-----------------------------------------------------------------------------
integer(spi):: 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(my,mah1,mah2,mx,a,v)!                                [UDLBY]
implicit none
integer(spi),intent(IN   ) :: my, mah1, mah2, mx
real(sp),    intent(IN   ) :: a(my,-mah1:mah2)
real(sp),    intent(INOUT) :: v(mx,my)
!-----------------------------------------------------------------------------
integer(spi):: 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(m,mah1,mah2,v,a)!                                    [UDLVB]
implicit none
integer(spi), intent(IN   ) :: m, mah1, mah2
real(sp),     intent(IN   ) :: a(m,-mah1:mah2)
real(sp),     intent(INOUT) :: v(m)
!-----------------------------------------------------------------------------
integer(spi):: i, j
real(sp)    :: 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(mx,mah1,mah2,my,v,a)!                                [UDLXB]
implicit none
integer(spi),intent(IN   ) :: mx, mah1, mah2, my
real(sp),    intent(IN   ) :: a(mx,-mah1:mah2)
real(sp),    intent(INOUT) :: v(mx,my)
!-----------------------------------------------------------------------------
integer(spi):: 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(my,mah1,mah2,mx,v,a)!                                [UDLYB]
implicit none
integer(spi),intent(IN   ) :: my, mah1, mah2, mx
real(sp),    intent(IN   ) :: a(my,-mah1:mah2)
real(sp),    intent(INOUT) :: v(mx,my)
!-----------------------------------------------------------------------------
integer(spi):: 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(m,mah1,mah2,a,v)!                                    [U1LBV]
implicit none
integer(spi),intent(IN   ) :: m, mah1, mah2
real(sp),    intent(IN   ) :: a(m,-mah1:mah2)
real(sp),    intent(INOUT) :: v(m)
!-----------------------------------------------------------------------------
integer(spi):: i, j
real(sp)    :: 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(mx,mah1,mah2,my,a,v)!                                [U1LBX]
implicit none
integer(spi),intent(IN   ) :: mx, mah1, mah2, my
real(sp),    intent(IN   ) :: a(mx,-mah1:mah2)
real(sp),    intent(INOUT) :: v(mx,my)
!-----------------------------------------------------------------------------
integer(spi):: 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(my,mah1,mah2,mx,a,v)!                                [U1LBY]
implicit none
integer(spi),intent(IN   ) :: my, mah1, mah2, mx
real(sp),    intent(IN   ) :: a(my,-mah1:mah2)
real(sp),    intent(INOUT) :: v(mx,my)
!-----------------------------------------------------------------------------
integer(spi):: 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(m,mah1,mah2,v,a)!                                    [U1LVB]
implicit none
integer(spi),intent(IN   ) :: m, mah1, mah2
real(sp),    intent(IN   ) :: a(m,-mah1:mah2)
real(sp),    intent(INOUT) :: v(m)
!-----------------------------------------------------------------------------
integer(spi):: i, j
real(sp)    :: 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(mx,mah1,mah2,my,v,a)!                                [U1LXB]
implicit none
integer(spi),intent(IN   ) :: mx, mah1, mah2, my
real(sp),    intent(IN   ) :: a(mx,-mah1:mah2)
real(sp),    intent(INOUT) :: v(mx,my)
!-----------------------------------------------------------------------------
integer(spi):: 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(my,mah1,mah2,mx,v,a)!                                [U1LYB]
implicit none
integer(spi),intent(IN   ) :: my, mah1, mah2, mx
real(sp),    intent(IN   ) :: a(my,-mah1:mah2)
real(sp),    intent(INOUT) :: v(mx,my)
!-----------------------------------------------------------------------------
integer(spi):: 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(m,mah1,mah2,a,v)!                                     [LINBV]
implicit none
integer(spi),intent(IN   ) :: m, mah1, mah2
real(sp),    intent(INOUT) :: a(m,-mah1:mah2), v(m)
!=============================================================================
call ldub(m,mah1,mah2,a)
call udlbv(m,mah1,mah2,a,v)
end subroutine linbv
 
subroutine wrtb(m1,m2,mah1,mah2,a)!                                     [WRTB]
implicit none
integer(spi),intent(IN) :: m1, m2, mah1, mah2
real(sp),    intent(IN) :: a(m1,-mah1:mah2)
!-----------------------------------------------------------------------------
integer(spi):: 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 pmat2