cwidth Subroutine

This routine is a variation of the theme in the algorithm
by Martin and Wilkinson.  It solves the linear system

  A * X = B

of NEQU equations in case A is almost block diagonal with all
blocks having NCOLS columns using no more storage than it takes to
store the interesting part of A.  Such systems occur in the determination
of the B-spline coefficients of a spline approximation.

The interesting part of A is taken to consist of NBLOKS
consecutive blocks, with the I-th block made up of NROWI = INTEGS(1,I)
consecutive rows and NCOLS consecutive columns of A, and with
the first LASTI = INTEGS(2,I) columns to the left of the next block.
These blocks are stored consecutively in the work array W.

For example, here is an 11th order matrix and its arrangement in
the work array W.  (The interesting entries of A are indicated by
their row and column index modulo 10.)

               ---   A   ---                          ---   W   ---

                  NROW1=3
       11 12 13 14                                     11 12 13 14
       21 22 23 24                                     21 22 23 24
       31 32 33 34      NROW2=2                        31 32 33 34
LAST1=2      43 44 45 46                               43 44 45 46
             53 54 55 56         NROW3=3               53 54 55 56
      LAST2=3         66 67 68 69                      66 67 68 69
                      76 77 78 79                      76 77 78 79
                      86 87 88 89   NROW4=1            86 87 88 89
               LAST3=1   97 98 99 90   NROW5=2         97 98 99 90
                  LAST4=1   08 09 00 01                08 09 00 01
                            18 19 10 11                18 19 10 11
                     LAST5=4

For this interpretation of A as an almost block diagonal matrix,
we have NBLOKS = 5, and the INTEGS array is

                      I = 1   2   3   4   5
                K =
INTEGS(K,I) =      1      3   2   3   1   2
                   2      2   3   1   1   4

Gauss elimination with scaled partial pivoting is used, but
multipliers are not saved in order to save storage.  Rather, the
right hand side is operated on during elimination.  The two parameters
IPVTEQ and LASTEQ are used to keep track of the action.  IPVTEQ 
is the index of the variable to be eliminated next, from equations 
IPVTEQ+1,...,LASTEQ, using equation IPVTEQ, possibly after an 
interchange, as the pivot equation.

The entries in the pivot column are always in column
1 of W.  This is accomplished by putting the entries in rows
IPVTEQ+1,...,LASTEQ revised by the elimination of the IPVTEQ-th
variable one to the left in W.  In this way, the columns of the
equations in a given block, as stored in W, will be aligned with
those of the next block at the moment when these next equations 
become involved in the elimination process.

Thus, for the above example, the first elimination steps proceed
as follows.

*11 12 13 14    11 12 13 14    11 12 13 14    11 12 13 14
*21 22 23 24   *22 23 24       22 23 24       22 23 24
*31 32 33 34   *32 33 34      *33 34          33 34
 43 44 45 46    43 44 45 46   *43 44 45 46   *44 45 46        
 53 54 55 56    53 54 55 56   *53 54 55 56   *54 55 56
 66 67 68 69    66 67 68 69    66 67 68 69    66 67 68 69

In all other respects, the procedure is standard, including the
scaled partial pivoting.

14 February 2007

Carl DeBoor

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

Roger Martin, James Wilkinson,
Solution of Symmetric and Unsymmetric Band Equations and
the Calculation of Eigenvectors of Band Matrices,
Numerische Mathematik,
Volume 9, Number 4, December 1976, pages 279-301.

Input/output, real ( kind = 8 ) W(NEQU,NCOLS), on input, contains
the interesting part of the almost block diagonal coefficient matrix 
A.  The array INTEGS describes the storage scheme.  On output, W 
contains the upper triangular factor U of the LU factorization of a 
possibly permuted version of A.  In particular, the determinant of 
A could now be found as
  IFLAG * W(1,1) * W(2,1) * ... * W(NEQU,1).

Input/output, real ( kind = 8 ) B(NEQU); on input, the right hand
side of the linear system.  On output, B has been overwritten by
other information.

Input, integer ( kind = 4 ) NEQU, the number of equations.

Input, integer ( kind = 4 ) NCOLS, the block width, that is, the number of 
columns in each block.

Input, integer ( kind = 4 ) INTEGS(2,NEQU), describes the block structure 
of A.
INTEGS(1,I) = number of rows in block I = NROW.
INTEGS(2,I) = number of elimination steps in block I = overhang over 
next block = LAST.

Input, integer ( kind = 4 ) NBOKS, the number of blocks.

Workspace, real D(NEQU), used to contain row sizes.  If storage is 
scarce, the array X could be used in the calling sequence for D.

Output, real ( kind = 8 ) X(NEQU), the computed solution, if 
IFLAG is nonzero.

Output, integer ( kind = 4 ) IFLAG, error flag.
= (-1)**(number of interchanges during elimination) if A is invertible;
= 0 if A is singular.