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