slvblk Subroutine

subroutine slvblk ( bloks, integs, nbloks, b, ipivot, x, iflag )

  !*****************************************************************************80
  !
  !! SLVBLK solves the almost block diagonal linear system A * x = b.  
  !
  !  Discussion:
  !
  !    Such almost block diagonal matrices arise naturally in piecewise 
  !    polynomial interpolation or approximation and in finite element 
  !    methods for two-point boundary value problems.  The PLU factorization 
  !    method is implemented here to take advantage of the special structure 
  !    of such systems for savings in computing time and storage requirements.
  !
  !    SLVBLK relies on several auxiliary programs:
  !
  !    FCBLOK (BLOKS,INTEGS,NBLOKS,IPIVOT,SCRTCH,IFLAG)  
  !    factors the matrix A.
  !
  !    SBBLOK (BLOKS,INTEGS,NBLOKS,IPIVOT,B,X)
  !    solves the system A*X=B once A is factored.
  !
  !    DTBLOK (BLOKS,INTEGS,NBLOKS,IPIVOT,IFLAG,DETSGN,DETLOG) 
  !    computes the determinant of A once it has been factored.
  !
  !  Block structure of A:
  !
  !    The NBLOKS blocks are stored consecutively in the array BLOKS.
  !
  !    The first block has its (1,1)-entry at BLOKS(1), and, if the I-th
  !    block has its (1,1)-entry at BLOKS(INDEX(I)), then
  !
  !      INDEX(I+1) = INDEX(I) + NROW(I) * NCOL(I).
  !
  !    The blocks are pieced together to give the interesting part of A
  !    as follows.  For I=1,2,..., NBLOKS-1, the (1,1)-entry of the next
  !    block (the (I+1)st block) corresponds to the (LAST+1,LAST+1)-entry
  !    of the current I-th block.  Recall LAST = INTEGS(3,I) and note that
  !    this means that
  !
  !    A: every block starts on the diagonal of A.
  !
  !    B: the blocks overlap (usually). the rows of the (I+1)st block
  !       which are overlapped by the I-th block may be arbitrarily 
  !       defined initially.  They are overwritten during elimination.
  !
  !    The right hand side for the equations in the I-th block are stored
  !    correspondingly as the last entries of a piece of B of length NROW
  !    (= INTEGS(1,I)) and following immediately in B the corresponding
  !    piece for the right hand side of the preceding block, with the right 
  !    hand side for the first block starting at B(1).  In this, the right 
  !    hand side for an equation need only be specified once on input, 
  !    in the first block in which the equation appears.
  !
  !  Example:
  !
  !    The test driver for this package contains an example, a linear
  !    system of order 11, whose nonzero entries are indicated in the
  !    following diagram by their row and column index modulo 10.  Next to it
  !    are the contents of the INTEGS arrray when the matrix is taken to
  !    be almost block diagonal with NBLOKS = 5, and below it are the five
  !    blocks.
  !
  !                        NROW1 = 3, NCOL1 = 4
  !             11 12 13 14
  !             21 22 23 24   NROW2 = 3, NCOL2 = 3
  !             31 32 33 34
  !    LAST1 = 2      43 44 45
  !                   53 54 55            NROW3 = 3, NCOL3 = 4
  !          LAST2 = 3         66 67 68 69   NROW4 = 3, NCOL4 = 4
  !                            76 77 78 79      NROW5 = 4, NCOL5 = 4
  !                            86 87 88 89
  !                   LAST3 = 1   97 98 99 90
  !                      LAST4 = 1   08 09 00 01
  !                                  18 19 10 11
  !                         LAST5 = 4
  !
  !    Actual input to BLOKS shown by rows of blocks of A.
  !    The ** items are arbitrary.
  !
  !    11 12 13 14  / ** ** **  / 66 67 68 69  / ** ** ** **  / ** ** ** **
  !    21 22 23 24 /  43 44 45 /  76 77 78 79 /  ** ** ** ** /  ** ** ** **
  !    31 32 33 34/   53 54 55/   86 87 88 89/   97 98 99 90/   08 09 00 01
  !                                                             18 19 10 11
  !
  !    INDEX = 1      INDEX = 13  INDEX = 22     INDEX = 34     INDEX = 46
  !
  !    Actual right hand side values with ** for arbitrary values:
  !
  !      B1 B2 B3 ** B4 B5 B6 B7 B8 ** ** B9 ** ** B10 B11
  !
  !    It would have been more efficient to combine block 3 with block 4.
  !
  !  Modified:
  !
  !    14 February 2007
  !
  !  Author:
  !
  !    Carl DeBoor
  !
  !  Reference:
  !
  !    Carl DeBoor,
  !    A Practical Guide to Splines,
  !    Springer, 2001,
  !    ISBN: 0387953663,
  !    LC: QA1.A647.v27.
  !
  !  Parameters:
  !
  !    Input/output, real ( kind = 8 ) BLOKS(*), a one-dimenional array, 
  !    of length sum ( INTEGS(1,1:NBLOKS) * INTEGS(2,1:NBLOKS) ).
  !    On input, contains the blocks of the almost block diagonal matrix A.  
  !    The array INTEGS describes the block structure.
  !    On output, contains correspondingly the PLU factorization
  !    of A, if IFLAG /= 0.  Certain entries in BLOKS are arbitrary, 
  !    where the blocks overlap.
  !
  !    Input, integer ( kind = 4 ) INTEGS(3,NBLOKS), description of the block 
  !    structure of A.
  !    integs(1,I) = number of rows of block I = nrow;
  !    integs(2,I) = number of colums of block I = ncol;
  !    integs(3,I) = number of elimination steps in block I = last.
  !    The linear system is of order n = sum ( integs(3,i), i=1,...,nbloks ),
  !    but the total number of rows in the blocks is
  !    nbrows=sum( integs(1,i) ; i = 1,...,nbloks)
  !
  !    Input, integer ( kind = 4 ) NBLOKS, the number of blocks.
  !
  !    Input, real ( kind = 8 ) B(NBROWS), the right hand side.  Certain entries 
  !    are arbitrary, corresponding to rows of the blocks which overlap.  See 
  !    the block structure in the example.
  !
  !    Output, integer ( kind = 4 ) IPIVOT(NBROWS), the pivoting sequence used.
  !
  !    Output, real ( kind = 8 ) X(N), the computed solution, if iflag /= 0.
  !
  !    Output, integer ( kind = 4 ) IFLAG.
  !    = (-1)**(number of interchanges during factorization) if A is invertible;
  !    = 0 if A is singular.
  !
  implicit none

  integer ( kind = 4 ) nbloks

  real ( kind = 8 ) b(*)
  real ( kind = 8 ) bloks(*)
  integer ( kind = 4 ) iflag
  integer ( kind = 4 ) integs(3,nbloks)
  integer ( kind = 4 ) ipivot(*)
  real ( kind = 8 ) x(*)
  !
  !  In the call to FCBLOK, X is used for temporary storage.
  !
  call fcblok ( bloks, integs, nbloks, ipivot, x, iflag )

  if ( iflag == 0 ) then
     return
  end if

  call sbblok ( bloks, integs, nbloks, ipivot, b, x )

  return
end subroutine slvblk