cwidth.f

Go to the documentation of this file.

      subroutine cwidth ( w,b,nequ,ncols,integs,nbloks, d, x,iflag )    

c*********************************************************************72
c
cc CWIDTH solves an almost block diagonal linear system.
c
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     
      implicit none

      integer nbloks

      integer iflag,integs(2,nbloks),ncols,nequ,   i,ii,icount,ipvteq   
     &  ,ipvtp1,istar,j,jmax,lastcl,lasteq,lasti,nexteq,nrowad   
      double precision b(nequ),d(nequ),w(nequ,ncols),x(nequ),
     &   awi1od,colmax,ratio,rowmax,sum,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.0d+00 
            do 5 j=1,ncols    
    5          rowmax = dmax1(rowmax,dabs(w(nexteq,j)))
            if (rowmax .eq. 0.0d+00)         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 ( dabs(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 = dabs(w(ipvteq,1))/d(ipvteq) 
            istar = ipvteq    
            ipvtp1 = ipvteq + 1     
            do 13 ii=ipvtp1,lasteq  
               awi1od = dabs(w(ii,1))/d(ii)
               if (awi1od .le. colmax)  go to 13
               colmax = awi1od
               istar = ii     
   13          continue 
            if ( dabs(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.0d+00   
   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.0d+00   
            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