c  Part IV, Problem 2.2
c  Diffusive decays of toroidal fields
c==========================================
c  This file contains the FORTRAN-77 subroutines:
c
c  matxvec  <-- matrix-vector multiply
c  tribu    <-- triangularized banded matrix
c  rhsbu    <-- reduces RHS given a triangularized matrix
c  fixbc    <-- implements nul-Dirichlet boundary conditions
c
c  See file README.txt for additional information
c
c******************************************
      subroutine matxvec(n,mhb,mb,gk,a,b)
c******************************************
c     matrix-vector multiply for shifted-stored matrix:
c
c     {b}=[gk]{a}
c
c     nm is maximum matrix dimension
c     nb is maximum full bandwidth
c
c     INPUT:
c     n   is actual matrix dimension
c     mhb is actual half bandwidth
c     mb  is actual full bandwidth
c     gk  is rectangular array of dimension (mb,n)
c     a   is vector of dimension (n) by which matrix is multiplied
c
c     OUTPUT:
c     b   is vector of dimension (n) that results from [gk]{b}
c==========================================
c     if   (i,j)       is position in original squared matrix,
c     then (i+mhb-j,j) is position in shifted matrix
c******************************************
      parameter(nb=3,nm=200)
      dimension gk(nb,nm),a(nm),b(nm)

      do 1 j=1,n
	 sum=0.
	 jpmb=j+mhb
	 do 2 i=max(1,jpmb-n),min(mb,jpmb-1)
	    jshift=jpmb-i
	    sum=sum+gk(i,jshift)*a(jshift)
    2    continue
	 b(j)=sum
    1 continue

      return
      end

c******************************************
      subroutine fixbc(n,mhb,mb,gk)
c******************************************
c     hardwiring of Dirichlet boundary conditions
c     at 1-D mesh endpoints
c     SPECIALIZED FOR NUL BOUNDARY CONDITION
c
c     The matrix system is
c
c     [gk]{a}={f}
c
c     nm is maximum matrix dimension
c     nb is maximum full bandwidth
c
c     INPUT:
c     n   is actual matrix dimension
c     mhb is actual half bandwidth
c     mb  is actual full bandwidth
c     gk  is rectangular array of dimension (mb,n)
c
c     OUTPUT:
c     gk  is modified form of system matrix
c==========================================
c     if   (i,j)       is position in original squared matrix,
c     then (i+mhb-j,j) is position in shifted matrix
c******************************************
      parameter(nb=3,nm=200)
      dimension gk(nb,nm)
      dimension uebc(2),iuebc(2)

c---------- safety check on array sizes
      if(mb.gt.nb)then
	 write(*,10)
         stop ' *FAIL*'
      endif
   10 format(1x,'FIXBC: bandwith mb too large: set nb=mb')
      if(n.gt.nm)then
	 write(*,11)
         stop ' *FAIL*'
      endif
   11 format(1x,'FIXBC: mesh size n too large: set nm=n')

      iuebc(1)=1
      iuebc(2)=n
      uebc(1) =0.
      uebc(2) =0.
      numbc   =2
      
c---------- 1. modify system matrix
      do 1 ibc=1,numbc
         jj=iuebc(ibc)
c --> zero matrix "line" (shifted version)
         do 2 i=2,mhb
            ishp=mhb-i+1
            ishm=mhb+i-1
            jrp=jj+i-1
            jrm=jj-i+1
            if(jrm.ge.1) gk(ishm,jrm)=0.
            if(jrp.le.n) gk(ishp,jrp)=0.
    2    continue
c --> zero matrix column
         do 3 i=1,mb
            gk(i,jj)=0.
    3    continue
c --> set corresponding diagonal element of system matrix to Real(one)
	 gk(mhb,jj)=1.
    1 continue

      return
      end

c******************************************
      subroutine tribu(n,mhb,mb,gk)
c******************************************
c     triangularized shifted matrix
c     (modified alg., after Golub & Van Loan, Matrix Computations, p. 152)
c     NO PIVOTING IS CARRIED OUT, so no antisymmetric matrices please
c
c     nm is maximum matrix dimension
c     nb is maximum full bandwidth
c
c     INPUT:
c     n   is actual matrix dimension
c     mhb is actual half bandwidth
c     mb  is actual full bandwidth
c     gk  is rectangular array of dimension (mb,n)
c
c     OUTPUT:
c     gk  is triangularized form of input matrix
c==========================================
c     if   (i,j)       is position in original squared matrix,
c     then (i+mhb-j,j) is position in shifted matrix
c******************************************
      parameter(nb=3,nm=200)
      dimension gk(nb,nm)

      do 1 k=1,n-1
	 if(gk(mhb,k).eq.0.)goto 99
	 c=1./gk(mhb,k)
	 do 2 i=mhb+1,mb
            gk(i,k)=c*gk(i,k)
    2    continue
    5    continue
         do 3 j=1,min(mhb-1,n-k)
	    kpj=k+j
	    c=gk(mhb-j,kpj)
	    if(c.ne.0.)then
CDIR$ IVDEP
	    do 4 i=mhb+1,min(mb,n-k+mhb)
	       imj=i-j
	       gk(imj,kpj)=gk(imj,kpj)
     +                          -c*gk(i,k)
    4       continue
	    end if
    3    continue
    1 continue
      return
   99 write(*,100) k,gk(mhb,k)
      stop 'tribu'
  100 format(//,' TRIBU: set of equations may be singular..',/,
     +'   diagonal term of equation',i5,'  at tribu is equal to',
     +1pe15.8)
      end

c******************************************
      subroutine rhsbu(n,mhb,mb,gk,gf,u)
c******************************************
c     forward reduction & backsubstitution
c     (modified alg., after Golub & Van Loan, p. 152)
c     nb is maximum matrix dimension
c     nm is maximum full bandwidth
c
c     INPUT:
c     n   is actual matrix dimension
c     mhb is actual half bandwidth
c     mb  is actual full bandwidth
c     gk  is triangularized rectangular matrix of dimension (mb,n)
c     gf  is RHS vector
c
c     OUTPUT:
c     u   is solution vector     
c     gf  is modified on output
c******************************************
      parameter(nb=3,nm=200)
      dimension gk(nb,nm),gf(nm),u(nm)

c---------- forward reduction
      do 1 j=1,n
	 if(gk(mhb,j).eq.0.)goto 99
	 mhbmj=mhb-j
	 do 2 i=j+1,min(n,j+mhb-1)
	    ishift=i+mhbmj
	    gf(i)=gf(i)-gk(ishift,j)*gf(j)
    2    continue
    1 continue

c---------- backsubstitution 

      do 3 j=n,1,-1
	 u(j)=gf(j)/gk(mhb,j)
	 mhbmj=mhb-j
	 do 4 i=max(1,j-mhb+1),j-1
	    ishift=i+mhbmj
            gf(i)=gf(i)-gk(ishift,j)*u(j)
    4    continue
    3 continue
      return
   99 write(*,100) j,gk(mhb,j)
      stop 'rhsbu'
  100 format(//,' RHSBU: set of equations are singular.',/,
     +'   diagonal term of equation',i5,'  at rhsbu is equal to',
     +1pe15.8)
      end