      subroutine cwidth ( w,b,nequ,ncols,integs,nbloks, d, x,iflag )
c  this program is a variation of the theme in the algorithm bandet1
c  by martin and wilkinson (numer.math. 9(1976)279-307). it solves
c  the linear system
c                           a*x  =  b
c  of  nequ  equations in case  a  is almost block diagonal with all
c  blocks having  ncols  columns using no more storage than it takes to
c  store the interesting part of  a . such systems occur in the determ-
c  ination of the b-spline coefficients of a spline approximation.
c
c                           parameters
c  w     on input, a two-dimensional array of size (nequ,ncols) contain-
c        ing the interesting part of the almost block diagonal coeffici-
c        ent matrix  a (see description and example below). the array
c        integs  describes the storage scheme.
c        on output, w  contains the upper triangular factor  u  of the
c        lu factorization of a possibly permuted version of  a . in par-
c        ticular, the determinant of  a  could now be found as
c            iflag*w(1,1)*w(2,1)* ... * w(nequ,1)  .
c  b     on input, the right side of the linear system, of length  nequ.
c        the contents of  b  are changed during execution.
c  nequ  number of equations in system
c  ncols block width, i.e., number of columns in each block.
c  integs  integer array, of size (2,nequ), describing the block struct-
c        ure of  a .
c           integs(1,i) = no. of rows in block i               =  nrow
c           integs(2,i) = no. of elimination steps in block i
c                       = overhang over next block             =  last
c  nbloks  number of blocks
c  d     work array, to contain row sizes . if storage is scarce, the
c        array  x  could be used in the calling sequence for  d .
c  x     on output, contains computed solution (if iflag .ne. 0), of
c        length  nequ .
c  iflag  on output, integer
c         = (-1)**(no.of interchanges during elimination)
c                if  a  is invertible
c         =  0   if  a  is singular
c
c        ------  block structure of  a  ------
c     the interesting part of  a  is taken to consist of  nbloks  con-
c  secutive blocks, with the i-th block made up of  nrowi = integs(1,i)
c  consecutive rows and  ncols  consecutive columns of  a , and with
c  the first  lasti = integs(2,i) columns to the left of the next block.
c  these blocks are stored consecutively in the workarray  w .
c     for example, here is an 11th order matrix and its arrangement in
c  the workarray  w . (the interesting entries of  a  are indicated by
c  their row and column index modulo 10.)
c
c                  ---   a   ---                          ---   w   ---
c
c                     nrow1=3
c          11 12 13 14                                     11 12 13 14
c          21 22 23 24                                     21 22 23 24
c          31 32 33 34      nrow2=2                        31 32 33 34
c   last1=2      43 44 45 46                               43 44 45 46
c                53 54 55 56         nrow3=3               53 54 55 56
c         last2=3         66 67 68 69                      66 67 68 69
c                         76 77 78 79                      76 77 78 79
c                         86 87 88 89   nrow4=1            86 87 88 89
c                  last3=1   97 98 99 90   nrow5=2         97 98 99 90
c                     last4=1   08 09 00 01                08 09 00 01
c                               18 19 10 11                18 19 10 11
c                        last5=4
c
c  for this interpretation of  a  as an almost block diagonal matrix,
c  we have  nbloks = 5 , and the integs array is
c
c                        i=  1   2   3   4   5
c                  k=
c  integs(k,i) =      1      3   2   3   1   2
c                     2      2   3   1   1   4
c
c       --------  method  --------
c     gauss elimination with scaled partial pivoting is used, but mult-
c  ipliers are  n o t  s a v e d  in order to save storage. rather, the
c  right side is operated on during elimination.
c     the two parameters
c                  i p v t e q   and  l a s t e q
c  are used to keep track of the action.  ipvteq is the index of the
c  variable to be eliminated next, from equations  ipvteq+1,...,lasteq,
c  using equation  ipvteq (possibly after an interchange) as the pivot
c  equation. the entries in the pivot column are  a l w a y s  in column
c  1 of  w . this is accomplished by putting the entries in rows
c  ipvteq+1,...,lasteq  revised by the elimination of the  ipvteq-th
c  variable one to the left in  w . in this way, the columns of the
c  equations in a given block (as stored in  w ) will be aligned with
c  those of the next block at the moment when these next equations be-
c  come involved in the elimination process.
c     thus, for the above example, the first elimination steps proceed
c  as follows.
c
c  *11 12 13 14    11 12 13 14    11 12 13 14    11 12 13 14
c  *21 22 23 24   *22 23 24       22 23 24       22 23 24
c  *31 32 33 34   *32 33 34      *33 34          33 34
c   43 44 45 46    43 44 45 46   *43 44 45 46   *44 45 46        etc.
c   53 54 55 56    53 54 55 56   *53 54 55 56   *54 55 56
c   66 67 68 69    66 67 68 69    66 67 68 69    66 67 68 69
c        .              .              .              .
c
c     in all other respects, the procedure is standard, including the
c  scaled partial pivoting.
c
      integer iflag,integs(2,nbloks),ncols,nequ,   i,ii,icount,ipvteq
     *                     ,istar,j,lastcl,lasteq,lasti,nexteq,nrowad
      real b(nequ),d(nequ),w(nequ,ncols),x(nequ),   awi1od,colmax
     *                                             ,rowmax,temp
      iflag = 1
      ipvteq = 0
      lasteq = 0
c                                 the i-loop runs over the blocks
      do 50 i=1,nbloks
c
c        the equations for the current block are added to those current-
c        ly involved in the elimination process, by increasing  lasteq
c        by  integs(1,i) after the rowsize of these equations has been
c        recorded in the array  d .
c
         nrowad = integs(1,i)
         do 10 icount=1,nrowad
            nexteq = lasteq + icount
            rowmax = 0.
            do 5 j=1,ncols
    5          rowmax = amax1(rowmax,abs(w(nexteq,j)))
            if (rowmax .eq. 0.)         go to 999
   10       d(nexteq) = rowmax
         lasteq = lasteq + nrowad
c
c        there will be  lasti = integs(2,i)  elimination steps before
c        the equations in the next block become involved. further,
c        l a s t c l  records the number of columns involved in the cur-
c        rent elimination step. it starts equal to  ncols  when a block
c        first becomes involved and then drops by one after each elim-
c        ination step.
c
         lastcl = ncols
         lasti = integs(2,i)
         do 30 icount=1,lasti
            ipvteq = ipvteq + 1
            if (ipvteq .lt. lasteq)     go to 11
            if ( abs(w(ipvteq,1))+d(ipvteq) .gt. d(ipvteq) )
     *                                  go to 50
                                        go to 999
c
c        determine the smallest  i s t a r  in  (ipvteq,lasteq)  for
c        which  abs(w(istar,1))/d(istar)  is as large as possible, and
c        interchange equations  ipvteq  and  istar  in case  ipvteq
c        .lt. istar .
c
   11       colmax = abs(w(ipvteq,1))/d(ipvteq)
            istar = ipvteq
            ipvtp1 = ipvteq + 1
            do 13 ii=ipvtp1,lasteq
               awi1od = abs(w(ii,1))/d(ii)
               if (awi1od .le. colmax)  go to 13
               colmax = awi1od
               istar = ii
   13          continue
            if ( abs(w(istar,1))+d(istar) .eq. d(istar) )
     *                                  go to 999
            if (istar .eq. ipvteq)      go to 16
            iflag = -iflag
            temp = d(istar)
            d(istar) = d(ipvteq)
            d(ipvteq) = temp
            temp = b(istar)
            b(istar) = b(ipvteq)
            b(ipvteq) = temp
            do 14 j=1,lastcl
               temp = w(istar,j)
               w(istar,j) = w(ipvteq,j)
   14          w(ipvteq,j) = temp
c
c        subtract the appropriate multiple of equation  ipvteq  from
c        equations  ipvteq+1,...,lasteq to make the coefficient of the
c        ipvteq-th unknown (presently in column 1 of  w ) zero, but
c        store the new coefficients in  w  one to the left from the old.
c
   16       do 20 ii=ipvtp1,lasteq
               ratio = w(ii,1)/w(ipvteq,1)
               do 18 j=2,lastcl
   18             w(ii,j-1) = w(ii,j) - ratio*w(ipvteq,j)
               w(ii,lastcl) = 0.
   20          b(ii) = b(ii) - ratio*b(ipvteq)
   30       lastcl = lastcl - 1
   50    continue
c
c  at this point,  w  and  b  contain an upper triangular linear system
c  equivalent to the original one, with  w(i,j) containing entry
c  (i, i-1+j ) of the coefficient matrix. solve this system by backsub-
c  stitution, taking into account its block structure.
c
c                          i-loop over the blocks, in reverse order
      i = nbloks
   59    lasti = integs(2,i)
         jmax = ncols - lasti
         do 70 icount=1,lasti
            sum = 0.
            if (jmax .eq. 0)            go to 61
            do 60 j=1,jmax
   60          sum = sum + x(ipvteq+j)*w(ipvteq,j+1)
   61       x(ipvteq) = (b(ipvteq)-sum)/w(ipvteq,1)
            jmax = jmax + 1
   70       ipvteq = ipvteq - 1
         i = i - 1
         if (i .gt. 0)                  go to 59
                                        return
  999 iflag = 0
                                        return
      end