      subroutine slvblk ( bloks, integs, nbloks, b, ipivot, x, iflag )
c    this program solves  the  linear system  a*x = b  where a is an
c  almost block diagonal matrix.  such almost block diagonal matrices
c  arise naturally in piecewise polynomial interpolation or approx-
c  imation and in finite element methods for two-point boundary value
c  problems.  the plu factorization method is implemented here to take
c  advantage of the special structure of such systems for savings in
c  computing time and storage requirements.
c
c                  parameters
c  bloks   a one-dimenional array, of length
c                   sum( integs(1,i)*integs(2,i) ; i = 1,nbloks )
c          on input, contains the blocks of the almost block diagonal
c          matrix  a  .  the array integs (see below and the example)
c          describes the block structure.
c          on output, contains correspondingly the plu factorization
c          of  a  (if iflag .ne. 0).  certain of the entries into bloks
c          are arbitrary (where the blocks overlap).
c  integs  integer array description of the block structure of  a .
c            integs(1,i) = no. of rows of block i        =  nrow
c            integs(2,i) = no. of colums of block i      =  ncol
c            integs(3,i) = no. of elim. steps in block i =  last
c                          i  = 1,2,...,nbloks
c          the linear system is of order
c                n  =  sum ( integs(3,i) , i=1,...,nbloks ),
c          but the total number of rows in the blocks is
c              nbrows = sum( integs(1,i) ; i = 1,...,nbloks)
c  nbloks  number of blocks
c  b       right side of the linear system, array of length nbrows.
c          certain of the entries are arbitrary, corresponding to
c          rows of the blocks which overlap (see block structure and
c          the example below).
c  ipivot  on output, integer array containing the pivoting sequence
c          used. length is nbrows
c  x       on output, contains the computed solution (if iflag .ne. 0)
c          length is n.
c  iflag   on output, integer
c            = (-1)**(no. of interchanges during factorization)
c                   if  a  is invertible
c            = 0    if  a  is singular
c
c                   auxiliary programs
c  fcblok (bloks,integs,nbloks,ipivot,scrtch,iflag)  factors the matrix
c           a , and is used for this purpose in slvblk. its arguments
c          are as in slvblk, except for
c              scrtch = a work array of length max(integs(1,i)).
c
c  sbblok (bloks,integs,nbloks,ipivot,b,x)  solves the system a*x = b
c          once  a  is factored. this is done automatically by slvblk
c          for one right side b, but subsequent solutions may be
c          obtained for additional b-vectors. the arguments are all
c          as in slvblk.
c
c  dtblok (bloks,integs,nbloks,ipivot,iflag,detsgn,detlog) computes the
c          determinant of  a  once slvblk or fcblok has done the fact-
c          orization.the first five arguments are as in slvblk.
c              detsgn  = sign of the determinant
c              detlog  = natural log of the determinant
c
c             ------ block structure of  a  ------
c  the nbloks blocks are stored consecutively in the array  bloks .
c  the first block has its (1,1)-entry at bloks(1), and, if the i-th
c  block has its (1,1)-entry at bloks(index(i)), then
c         index(i+1) = index(i)  +  nrow(i)*ncol(i) .
c    the blocks are pieced together to give the interesting part of  a
c  as follows.  for i = 1,2,...,nbloks-1, the (1,1)-entry of the next
c  block (the (i+1)st block ) corresponds to the (last+1,last+1)-entry
c  of the current i-th block.  recall last = integs(3,i) and note that
c  this means that
c      a. every block starts on the diagonal of  a .
c      b. the blocks overlap (usually). the rows of the (i+1)st block
c         which are overlapped by the i-th block may be arbitrarily de-
c         fined initially. they are overwritten during elimination.
c    the right side for the equations in the i-th block are stored cor-
c  respondingly as the last entries of a piece of  b  of length  nrow
c  (= integs(1,i)) and following immediately in  b  the corresponding
c  piece for the right side of the preceding block, with the right side
c  for the first block starting at  b(1) . in this, the right side for
c  an equation need only be specified once on input, in the first block
c  in which the equation appears.
c
c             ------ example and test driver ------
c    the test driver for this package contains an example, a linear
c  system of order 11, whose nonzero entries are indicated in the fol-
c  lowing schema by their row and column index modulo 10. next to it
c  are the contents of the  integs  arrray when the matrix is taken to
c  be almost block diagonal with  nbloks = 5, and below it are the five
c  blocks.
c
c                      nrow1 = 3, ncol1 = 4
c           11 12 13 14
c           21 22 23 24   nrow2 = 3, ncol2 = 3
c           31 32 33 34
c  last1 = 2      43 44 45
c                 53 54 55            nrow3 = 3, ncol3 = 4
c        last2 = 3         66 67 68 69   nrow4 = 3, ncol4 = 4
c                          76 77 78 79      nrow5 = 4, ncol5 = 4
c                          86 87 88 89
c                 last3 = 1   97 98 99 90
c                    last4 = 1   08 09 00 01
c                                18 19 10 11
c                       last5 = 4
c
c         actual input to bloks shown by rows of blocks of  a .
c      (the ** items are arbitrary, this storage is used by slvblk)
c
c  11 12 13 14  / ** ** **  / 66 67 68 69  / ** ** ** **  / ** ** ** **
c  21 22 23 24 /  43 44 45 /  76 77 78 79 /  ** ** ** ** /  ** ** ** **
c  31 32 33 34/   53 54 55/   86 87 88 89/   97 98 99 90/   08 09 00 01
c                                                           18 19 10 11
c
c  index = 1      index = 13  index = 22     index = 34     index = 46
c
c         actual right side values with ** for arbitrary values
c  b1 b2 b3 ** b4 b5 b6 b7 b8 ** ** b9 ** ** b10 b11
c
c  (it would have been more efficient to combine block 3 with block 4)
c
      integer integs(3,nbloks),ipivot(1),iflag
      real bloks(1),b(1),x(1)
c     in the call to fcblok,  x  is used for temporary storage.
      call fcblok(bloks,integs,nbloks,ipivot,x,iflag)
      if (iflag .eq. 0)                 return
      call sbblok(bloks,integs,nbloks,ipivot,b,x)
                                        return
      end
      subroutine fcblok ( bloks, integs, nbloks, ipivot, scrtch, iflag )
calls subroutines  f a c t r b  and  s h i f t b .
c
c   f c b l o k  supervises the plu factorization with pivoting of
c  scaled rows of the almost block diagonal matrix stored in the arrays
c   b l o k s  and  i n t e g s .
c
c   factrb = subprogram which carries out steps 1,...,last of gauss
c            elimination (with pivoting) for an individual block.
c   shiftb = subprogram which shifts the remaining rows to the top of
c            the next block
c
c parameters
c    bloks   an array that initially contains the almost block diagonal
c            matrix  a  to be factored, and on return contains the com-
c            puted factorization of  a .
c    integs  an integer array describing the block structure of  a .
c    nbloks  the number of blocks in  a .
c    ipivot  an integer array of dimension  sum (integs(1,n) ; n=1,
c            ...,nbloks) which, on return, contains the pivoting stra-
c            tegy used.
c    scrtch  work area required, of length  max (integs(1,n) ; n=1,
c            ...,nbloks).
c    iflag   output parameter;
c            = 0  in case matrix was found to be singular.
c            otherwise,
c            = (-1)**(number of row interchanges during factorization)
c
      integer integs(3,nbloks),ipivot(1),iflag, i,index,indexb,indexn,
     *        last,ncol,nrow
      real bloks(1),scrtch(1)
      iflag = 1
      indexb = 1
      indexn = 1
      i = 1
c                        loop over the blocks.  i  is loop index
   10    index = indexn
         nrow = integs(1,i)
         ncol = integs(2,i)
         last = integs(3,i)
c        carry out elimination on the i-th block until next block
c        enters, i.e., for columns 1,...,last  of i-th block.
         call factrb(bloks(index),ipivot(indexb),scrtch,nrow,ncol,last,
     *            iflag)
c         check for having reached a singular block or the last block
         if (iflag .eq. 0 .or. i .eq. nbloks)
     *                                  return
         i = i+1
         indexn = nrow*ncol + index
c              put the rest of the i-th block onto the next block
         call shiftb(bloks(index),ipivot(indexb),nrow,ncol,last,
     *            bloks(indexn),integs(1,i),integs(2,i))
         indexb = indexb + nrow
                                        go to 10
      end
      subroutine factrb ( w, ipivot, d, nrow, ncol, last, iflag )
c  adapted from p.132 of 'element.numer.analysis' by conte-de boor
c
c  constructs a partial plu factorization, corresponding to steps 1,...,
c   l a s t   in gauss elimination, for the matrix  w  of order
c   ( n r o w ,  n c o l ), using pivoting of scaled rows.
c
c  parameters
c    w       contains the (nrow,ncol) matrix to be partially factored
c            on input, and the partial factorization on output.
c    ipivot  an integer array of length nrow containing a record of the
c            pivoting strategy used; row ipivot(i) is used during the
c            i-th elimination step, i=1,...,last.
c    d       a work array of length nrow used to store row sizes
c            temporarily.
c    nrow    number of rows of w.
c    ncol    number of columns of w.
c    last    number of elimination steps to be carried out.
c    iflag   on output, equals iflag on input times (-1)**(number of
c            row interchanges during the factorization process), in
c            case no zero pivot was encountered.
c            otherwise, iflag = 0 on output.
c
      integer ipivot(nrow),ncol,last,iflag, i,ipivi,ipivk,j,k,kp1
      real w(nrow,ncol),d(nrow), awikdi,colmax,ratio,rowmax
c  initialize ipivot, d
      do 10 i=1,nrow
         ipivot(i) = i
         rowmax = 0.
         do 9 j=1,ncol
    9       rowmax = amax1(rowmax, abs(w(i,j)))
         if (rowmax .eq. 0.)            go to 999
   10    d(i) = rowmax
c gauss elimination with pivoting of scaled rows, loop over k=1,.,last
      k = 1
c        as pivot row for k-th step, pick among the rows not yet used,
c        i.e., from rows ipivot(k),...,ipivot(nrow), the one whose k-th
c        entry (compared to the row size) is largest. then, if this row
c        does not turn out to be row ipivot(k), redefine ipivot(k) ap-
c        propriately and record this interchange by changing the sign
c        of  i f l a g .
   11    ipivk = ipivot(k)
         if (k .eq. nrow)               go to 21
         j = k
         kp1 = k+1
         colmax = abs(w(ipivk,k))/d(ipivk)
c              find the (relatively) largest pivot
         do 15 i=kp1,nrow
            ipivi = ipivot(i)
            awikdi = abs(w(ipivi,k))/d(ipivi)
            if (awikdi .le. colmax)     go to 15
               colmax = awikdi
               j = i
   15       continue
         if (j .eq. k)                  go to 16
         ipivk = ipivot(j)
         ipivot(j) = ipivot(k)
         ipivot(k) = ipivk
         iflag = -iflag
   16    continue
c        if pivot element is too small in absolute value, declare
c        matrix to be noninvertible and quit.
         if (abs(w(ipivk,k))+d(ipivk) .le. d(ipivk))
     *                                  go to 999
c        otherwise, subtract the appropriate multiple of the pivot
c        row from remaining rows, i.e., the rows ipivot(k+1),...,
c        ipivot(nrow), to make k-th entry zero. save the multiplier in
c        its place.
         do 20 i=kp1,nrow
            ipivi = ipivot(i)
            w(ipivi,k) = w(ipivi,k)/w(ipivk,k)
            ratio = -w(ipivi,k)
            do 20 j=kp1,ncol
   20          w(ipivi,j) = ratio*w(ipivk,j) + w(ipivi,j)
         k = kp1
c        check for having reached the next block.
         if (k .le. last)               go to 11
                                        return
c     if  last  .eq. nrow , check now that pivot element in last row
c     is nonzero.
   21 if( abs(w(ipivk,nrow))+d(ipivk) .gt. d(ipivk) )
     *                                  return
c                   singularity flag set
  999 iflag = 0
                                        return
      end
      subroutine shiftb ( ai, ipivot, nrowi, ncoli, last,
     *                    ai1, nrowi1, ncoli1 )
c  shifts the rows in current block, ai, not used as pivot rows, if
c  any, i.e., rows ipivot(last+1),...,ipivot(nrowi), onto the first
c  mmax = nrow-last rows of the next block, ai1, with column last+j of
c  ai  going to column j , j=1,...,jmax=ncoli-last. the remaining col-
c  umns of these rows of ai1 are zeroed out.
c
c                             picture
c
c       original situation after         results in a new block i+1
c       last = 2 columns have been       created and ready to be
c       done in factrb (assuming no      factored by next factrb call.
c       interchanges of rows)
c                   1
c              x  x 1x  x  x           x  x  x  x  x
c                   1
c              0  x 1x  x  x           0  x  x  x  x
c  block i          1                       ---------------
c  nrowi = 4   0  0 1x  x  x           0  0 1x  x  x  0  01
c  ncoli = 5        1                       1             1
c  last = 2    0  0 1x  x  x           0  0 1x  x  x  0  01
c  -------------------------------          1             1   new
c                   1x  x  x  x  x          1x  x  x  x  x1  block
c                   1                       1             1   i+1
c  block i+1        1x  x  x  x  x          1x  x  x  x  x1
c  nrowi1= 5        1                       1             1
c  ncoli1= 5        1x  x  x  x  x          1x  x  x  x  x1
c  -------------------------------          1-------------1
c                   1
c
      integer ipivot(nrowi),last, ip,j,jmax,jmaxp1,m,mmax
      real ai(nrowi,ncoli),ai1(nrowi1,ncoli1)
      mmax = nrowi - last
      jmax = ncoli - last
      if (mmax .lt. 1 .or. jmax .lt. 1) return
c              put the remainder of block i into ai1
      do 10 m=1,mmax
         ip = ipivot(last+m)
         do 10 j=1,jmax
   10       ai1(m,j) = ai(ip,last+j)
      if (jmax .eq. ncoli1)             return
c              zero out the upper right corner of ai1
      jmaxp1 = jmax + 1
      do 20 j=jmaxp1,ncoli1
         do 20 m=1,mmax
   20       ai1(m,j) = 0.
                                        return
      end
      subroutine sbblok ( bloks, integs, nbloks, ipivot, b, x )
calls subroutines  s u b f o r  and  s u b b a k .
c
c  supervises the solution (by forward and backward substitution) of
c  the linear system  a*x = b  for x, with the plu factorization of  a
c  already generated in  f c b l o k .  individual blocks of equations
c  are solved via  s u b f o r  and  s u b b a k .
c
c parameters
c    bloks, integs, nbloks, ipivot    are as on return from fcblok.
c    b       the right side, stored corresponding to the storage of
c            the equations. see comments in  s l v b l k  for details.
c    x       solution vector
c
      integer integs(3,nbloks),ipivot(1), i,index,indexb,indexx,j,last,
     *        nbp1,ncol,nrow
      real bloks(1),b(1),x(1)
c
c      forward substitution pass
c
      index = 1
      indexb = 1
      indexx = 1
      do 20 i=1,nbloks
         nrow = integs(1,i)
         last = integs(3,i)
         call subfor(bloks(index),ipivot(indexb),nrow,last,b(indexb),
     *               x(indexx))
         index = nrow*integs(2,i) + index
         indexb = indexb + nrow
   20    indexx = indexx + last
c
c     back substitution pass
c
      nbp1 = nbloks + 1
      do 30 j=1,nbloks
         i = nbp1 - j
         nrow = integs(1,i)
         ncol = integs(2,i)
         last = integs(3,i)
         index = index - nrow*ncol
         indexb = indexb - nrow
         indexx = indexx - last
   30    call subbak(bloks(index),ipivot(indexb),nrow,ncol,last,
     *               x(indexx))
                                        return
      end
      subroutine subfor ( w, ipivot, nrow, last, b, x )
c  carries out the forward pass of substitution for the current block,
c  i.e., the action on the right side corresponding to the elimination
c  carried out in  f a c t r b  for this block.
c     at the end, x(j) contains the right side of the transformed
c  ipivot(j)-th equation in this block, j=1,...,nrow. then, since
c  for i=1,...,nrow-last, b(nrow+i) is going to be used as the right
c  side of equation  i  in the next block (shifted over there from
c  this block during factorization), it is set equal to x(last+i) here.
c
c parameters
c    w, ipivot, nrow, last  are as on return from factrb.
c    b(j)   is expected to contain, on input, the right side of j-th
c           equation for this block, j=1,...,nrow.
c    b(nrow+j)   contains, on output, the appropriately modified right
c           side for equation j in next block, j=1,...,nrow-last.
c    x(j)   contains, on output, the appropriately modified right
c           side of equation ipivot(j) in this block, j=1,...,last (and
c           even for j=last+1,...,nrow).
c
      integer ipivot(nrow), ip,jmax,k
c     dimension b(nrow + nrow-last)
      real w(nrow,last),b(1),x(nrow)
      ip = ipivot(1)
      x(1) = b(ip)
      if (nrow .eq. 1)                  go to 99
      do 15 k=2,nrow
         ip = ipivot(k)
         jmax = amin0(k-1,last)
         sum = 0.
         do 14 j=1,jmax
   14       sum = w(ip,j)*x(j) + sum
   15    x(k) = b(ip) - sum
c
c     transfer modified right sides of equations ipivot(last+1),...,
c     ipivot(nrow) to next block.
      nrowml = nrow - last
      if (nrowml .eq. 0)                go to 99
      lastp1 = last+1
      do 25 k=lastp1,nrow
   25    b(nrowml+k) = x(k)
   99                                   return
      end
      subroutine subbak ( w, ipivot, nrow, ncol, last, x )
c  carries out backsubstitution for current block.
c
c parameters
c    w, ipivot, nrow, ncol, last  are as on return from factrb.
c    x(1),...,x(ncol)  contains, on input, the right side for the
c            equations in this block after backsubstitution has been
c            carried up to but not including equation ipivot(last).
c            means that x(j) contains the right side of equation ipi-
c            vot(j) as modified during elimination, j=1,...,last, while
c            for j .gt. last, x(j) is already a component of the solut-
c            ion vector.
c    x(1),...,x(ncol) contains, on output, the components of the solut-
c            ion corresponding to the present block.
c
      integer ipivot(nrow),last,  ip,j,k,kp1
      real w(nrow,ncol),x(ncol), sum
      k = last
      ip = ipivot(k)
      sum = 0.
      if (k .eq. ncol)                  go to 4
      kp1 = k+1
    2    do 3 j=kp1,ncol
    3       sum = w(ip,j)*x(j) + sum
    4    x(k) = (x(k) - sum)/w(ip,k)
         if (k .eq. 1)                  return
         kp1 = k
         k = k-1
         ip = ipivot(k)
         sum = 0.
                                        go to 2
      end
      subroutine dtblok ( bloks, integs, nbloks, ipivot, iflag,
     *                    detsgn, detlog )
c  computes the determinant of an almost block diagonal matrix whose
c  plu factorization has been obtained previously in fcblok.
c  *** the logarithm of the determinant is computed instead of the
c  determinant itself to avoid the danger of overflow or underflow
c  inherent in this calculation.
c
c parameters
c    bloks, integs, nbloks, ipivot, iflag  are as on return from fcblok.
c            in particular, iflag = (-1)**(number of interchanges dur-
c            ing factorization) if successful, otherwise iflag = 0.
c    detsgn  on output, contains the sign of the determinant.
c    detlog  on output, contains the natural logarithm of the determi-
c            nant if determinant is not zero. otherwise contains 0.
c
      integer integs(3,nbloks),ipivot(1),iflag, i,indexp,ip,k,last
      real bloks(1),detsgn,detlog
c
      detsgn = iflag
      detlog = 0.
      if (iflag .eq. 0)                 return
      index = 0
      indexp = 0
      do 2 i=1,nbloks
         nrow = integs(1,i)
         last = integs(3,i)
         do 1 k=1,last
            ip = index + nrow*(k-1) + ipivot(indexp+k)
            detlog = detlog + alog(abs(bloks(ip)))
    1       detsgn = detsgn*sign(1.,bloks(ip))
         index = nrow*integs(2,i) + index
    2    indexp = indexp + nrow
                                        return
      end