slvblk Subroutine

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.

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.

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.

14 February 2007

Carl DeBoor

Carl DeBoor,
A Practical Guide to Splines,
Springer, 2001,
ISBN: 0387953663,
LC: QA1.A647.v27.

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.