subroutine slvblk ( bloks, integs, nbloks, b, ipivot, x, iflag ) c this program solves the linear system a*x = b where a is an c almost block diagonal matrix. such almost block diagonal matrices c arise naturally in piecewise polynomial interpolation or approx- c imation and in finite element methods for two-point boundary value c problems. the plu factorization method is implemented here to take c advantage of the special structure of such systems for savings in c computing time and storage requirements. c c parameters c bloks a one-dimenional array, of length c sum( integs(1,i)*integs(2,i) ; i = 1,nbloks ) c on input, contains the blocks of the almost block diagonal c matrix a . the array integs (see below and the example) c describes the block structure. c on output, contains correspondingly the plu factorization c of a (if iflag .ne. 0). certain of the entries into bloks c are arbitrary (where the blocks overlap). c integs integer array description of the block structure of a . c integs(1,i) = no. of rows of block i = nrow c integs(2,i) = no. of colums of block i = ncol c integs(3,i) = no. of elim. steps in block i = last c i = 1,2,...,nbloks c the linear system is of order c n = sum ( integs(3,i) , i=1,...,nbloks ), c but the total number of rows in the blocks is c nbrows = sum( integs(1,i) ; i = 1,...,nbloks) c nbloks number of blocks c b right side of the linear system, array of length nbrows. c certain of the entries are arbitrary, corresponding to c rows of the blocks which overlap (see block structure and c the example below). c ipivot on output, integer array containing the pivoting sequence c used. length is nbrows c x on output, contains the computed solution (if iflag .ne. 0) c length is n. c iflag on output, integer c = (-1)**(no. of interchanges during factorization) c if a is invertible c = 0 if a is singular c c auxiliary programs c fcblok (bloks,integs,nbloks,ipivot,scrtch,iflag) factors the matrix c a , and is used for this purpose in slvblk. its arguments c are as in slvblk, except for c scrtch = a work array of length max(integs(1,i)). c c sbblok (bloks,integs,nbloks,ipivot,b,x) solves the system a*x = b c once a is factored. this is done automatically by slvblk c for one right side b, but subsequent solutions may be c obtained for additional b-vectors. the arguments are all c as in slvblk. c c dtblok (bloks,integs,nbloks,ipivot,iflag,detsgn,detlog) computes the c determinant of a once slvblk or fcblok has done the fact- c orization.the first five arguments are as in slvblk. c detsgn = sign of the determinant c detlog = natural log of the determinant c c ------ block structure of a ------ c the nbloks blocks are stored consecutively in the array bloks . c the first block has its (1,1)-entry at bloks(1), and, if the i-th c block has its (1,1)-entry at bloks(index(i)), then c index(i+1) = index(i) + nrow(i)*ncol(i) . c the blocks are pieced together to give the interesting part of a c as follows. for i = 1,2,...,nbloks-1, the (1,1)-entry of the next c block (the (i+1)st block ) corresponds to the (last+1,last+1)-entry c of the current i-th block. recall last = integs(3,i) and note that c this means that c a. every block starts on the diagonal of a . c b. the blocks overlap (usually). the rows of the (i+1)st block c which are overlapped by the i-th block may be arbitrarily de- c fined initially. they are overwritten during elimination. c the right side for the equations in the i-th block are stored cor- c respondingly as the last entries of a piece of b of length nrow c (= integs(1,i)) and following immediately in b the corresponding c piece for the right side of the preceding block, with the right side c for the first block starting at b(1) . in this, the right side for c an equation need only be specified once on input, in the first block c in which the equation appears. c c ------ example and test driver ------ c the test driver for this package contains an example, a linear c system of order 11, whose nonzero entries are indicated in the fol- c lowing schema by their row and column index modulo 10. next to it c are the contents of the integs arrray when the matrix is taken to c be almost block diagonal with nbloks = 5, and below it are the five c blocks. c c nrow1 = 3, ncol1 = 4 c 11 12 13 14 c 21 22 23 24 nrow2 = 3, ncol2 = 3 c 31 32 33 34 c last1 = 2 43 44 45 c 53 54 55 nrow3 = 3, ncol3 = 4 c last2 = 3 66 67 68 69 nrow4 = 3, ncol4 = 4 c 76 77 78 79 nrow5 = 4, ncol5 = 4 c 86 87 88 89 c last3 = 1 97 98 99 90 c last4 = 1 08 09 00 01 c 18 19 10 11 c last5 = 4 c c actual input to bloks shown by rows of blocks of a . c (the ** items are arbitrary, this storage is used by slvblk) c c 11 12 13 14 / ** ** ** / 66 67 68 69 / ** ** ** ** / ** ** ** ** c 21 22 23 24 / 43 44 45 / 76 77 78 79 / ** ** ** ** / ** ** ** ** c 31 32 33 34/ 53 54 55/ 86 87 88 89/ 97 98 99 90/ 08 09 00 01 c 18 19 10 11 c c index = 1 index = 13 index = 22 index = 34 index = 46 c c actual right side values with ** for arbitrary values c b1 b2 b3 ** b4 b5 b6 b7 b8 ** ** b9 ** ** b10 b11 c c (it would have been more efficient to combine block 3 with block 4) c integer integs(3,nbloks),ipivot(1),iflag real bloks(1),b(1),x(1) c in the call to fcblok, x is used for temporary storage. call fcblok(bloks,integs,nbloks,ipivot,x,iflag) if (iflag .eq. 0) return call sbblok(bloks,integs,nbloks,ipivot,b,x) return end subroutine fcblok ( bloks, integs, nbloks, ipivot, scrtch, iflag ) calls subroutines f a c t r b and s h i f t b . c c f c b l o k supervises the plu factorization with pivoting of c scaled rows of the almost block diagonal matrix stored in the arrays c b l o k s and i n t e g s . c c factrb = subprogram which carries out steps 1,...,last of gauss c elimination (with pivoting) for an individual block. c shiftb = subprogram which shifts the remaining rows to the top of c the next block c c parameters c bloks an array that initially contains the almost block diagonal c matrix a to be factored, and on return contains the com- c puted factorization of a . c integs an integer array describing the block structure of a . c nbloks the number of blocks in a . c ipivot an integer array of dimension sum (integs(1,n) ; n=1, c ...,nbloks) which, on return, contains the pivoting stra- c tegy used. c scrtch work area required, of length max (integs(1,n) ; n=1, c ...,nbloks). c iflag output parameter; c = 0 in case matrix was found to be singular. c otherwise, c = (-1)**(number of row interchanges during factorization) c integer integs(3,nbloks),ipivot(1),iflag, i,index,indexb,indexn, * last,ncol,nrow real bloks(1),scrtch(1) iflag = 1 indexb = 1 indexn = 1 i = 1 c loop over the blocks. i is loop index 10 index = indexn nrow = integs(1,i) ncol = integs(2,i) last = integs(3,i) c carry out elimination on the i-th block until next block c enters, i.e., for columns 1,...,last of i-th block. call factrb(bloks(index),ipivot(indexb),scrtch,nrow,ncol,last, * iflag) c check for having reached a singular block or the last block if (iflag .eq. 0 .or. i .eq. nbloks) * return i = i+1 indexn = nrow*ncol + index c put the rest of the i-th block onto the next block call shiftb(bloks(index),ipivot(indexb),nrow,ncol,last, * bloks(indexn),integs(1,i),integs(2,i)) indexb = indexb + nrow go to 10 end subroutine factrb ( w, ipivot, d, nrow, ncol, last, iflag ) c adapted from p.132 of 'element.numer.analysis' by conte-de boor c c constructs a partial plu factorization, corresponding to steps 1,..., c l a s t in gauss elimination, for the matrix w of order c ( n r o w , n c o l ), using pivoting of scaled rows. c c parameters c w contains the (nrow,ncol) matrix to be partially factored c on input, and the partial factorization on output. c ipivot an integer array of length nrow containing a record of the c pivoting strategy used; row ipivot(i) is used during the c i-th elimination step, i=1,...,last. c d a work array of length nrow used to store row sizes c temporarily. c nrow number of rows of w. c ncol number of columns of w. c last number of elimination steps to be carried out. c iflag on output, equals iflag on input times (-1)**(number of c row interchanges during the factorization process), in c case no zero pivot was encountered. c otherwise, iflag = 0 on output. c integer ipivot(nrow),ncol,last,iflag, i,ipivi,ipivk,j,k,kp1 real w(nrow,ncol),d(nrow), awikdi,colmax,ratio,rowmax c initialize ipivot, d do 10 i=1,nrow ipivot(i) = i rowmax = 0. do 9 j=1,ncol 9 rowmax = amax1(rowmax, abs(w(i,j))) if (rowmax .eq. 0.) go to 999 10 d(i) = rowmax c gauss elimination with pivoting of scaled rows, loop over k=1,.,last k = 1 c as pivot row for k-th step, pick among the rows not yet used, c i.e., from rows ipivot(k),...,ipivot(nrow), the one whose k-th c entry (compared to the row size) is largest. then, if this row c does not turn out to be row ipivot(k), redefine ipivot(k) ap- c propriately and record this interchange by changing the sign c of i f l a g . 11 ipivk = ipivot(k) if (k .eq. nrow) go to 21 j = k kp1 = k+1 colmax = abs(w(ipivk,k))/d(ipivk) c find the (relatively) largest pivot do 15 i=kp1,nrow ipivi = ipivot(i) awikdi = abs(w(ipivi,k))/d(ipivi) if (awikdi .le. colmax) go to 15 colmax = awikdi j = i 15 continue if (j .eq. k) go to 16 ipivk = ipivot(j) ipivot(j) = ipivot(k) ipivot(k) = ipivk iflag = -iflag 16 continue c if pivot element is too small in absolute value, declare c matrix to be noninvertible and quit. if (abs(w(ipivk,k))+d(ipivk) .le. d(ipivk)) * go to 999 c otherwise, subtract the appropriate multiple of the pivot c row from remaining rows, i.e., the rows ipivot(k+1),..., c ipivot(nrow), to make k-th entry zero. save the multiplier in c its place. do 20 i=kp1,nrow ipivi = ipivot(i) w(ipivi,k) = w(ipivi,k)/w(ipivk,k) ratio = -w(ipivi,k) do 20 j=kp1,ncol 20 w(ipivi,j) = ratio*w(ipivk,j) + w(ipivi,j) k = kp1 c check for having reached the next block. if (k .le. last) go to 11 return c if last .eq. nrow , check now that pivot element in last row c is nonzero. 21 if( abs(w(ipivk,nrow))+d(ipivk) .gt. d(ipivk) ) * return c singularity flag set 999 iflag = 0 return end subroutine shiftb ( ai, ipivot, nrowi, ncoli, last, * ai1, nrowi1, ncoli1 ) c shifts the rows in current block, ai, not used as pivot rows, if c any, i.e., rows ipivot(last+1),...,ipivot(nrowi), onto the first c mmax = nrow-last rows of the next block, ai1, with column last+j of c ai going to column j , j=1,...,jmax=ncoli-last. the remaining col- c umns of these rows of ai1 are zeroed out. c c picture c c original situation after results in a new block i+1 c last = 2 columns have been created and ready to be c done in factrb (assuming no factored by next factrb call. c interchanges of rows) c 1 c x x 1x x x x x x x x c 1 c 0 x 1x x x 0 x x x x c block i 1 --------------- c nrowi = 4 0 0 1x x x 0 0 1x x x 0 01 c ncoli = 5 1 1 1 c last = 2 0 0 1x x x 0 0 1x x x 0 01 c ------------------------------- 1 1 new c 1x x x x x 1x x x x x1 block c 1 1 1 i+1 c block i+1 1x x x x x 1x x x x x1 c nrowi1= 5 1 1 1 c ncoli1= 5 1x x x x x 1x x x x x1 c ------------------------------- 1-------------1 c 1 c integer ipivot(nrowi),last, ip,j,jmax,jmaxp1,m,mmax real ai(nrowi,ncoli),ai1(nrowi1,ncoli1) mmax = nrowi - last jmax = ncoli - last if (mmax .lt. 1 .or. jmax .lt. 1) return c put the remainder of block i into ai1 do 10 m=1,mmax ip = ipivot(last+m) do 10 j=1,jmax 10 ai1(m,j) = ai(ip,last+j) if (jmax .eq. ncoli1) return c zero out the upper right corner of ai1 jmaxp1 = jmax + 1 do 20 j=jmaxp1,ncoli1 do 20 m=1,mmax 20 ai1(m,j) = 0. return end subroutine sbblok ( bloks, integs, nbloks, ipivot, b, x ) calls subroutines s u b f o r and s u b b a k . c c supervises the solution (by forward and backward substitution) of c the linear system a*x = b for x, with the plu factorization of a c already generated in f c b l o k . individual blocks of equations c are solved via s u b f o r and s u b b a k . c c parameters c bloks, integs, nbloks, ipivot are as on return from fcblok. c b the right side, stored corresponding to the storage of c the equations. see comments in s l v b l k for details. c x solution vector c integer integs(3,nbloks),ipivot(1), i,index,indexb,indexx,j,last, * nbp1,ncol,nrow real bloks(1),b(1),x(1) c c forward substitution pass c index = 1 indexb = 1 indexx = 1 do 20 i=1,nbloks nrow = integs(1,i) last = integs(3,i) call subfor(bloks(index),ipivot(indexb),nrow,last,b(indexb), * x(indexx)) index = nrow*integs(2,i) + index indexb = indexb + nrow 20 indexx = indexx + last c c back substitution pass c nbp1 = nbloks + 1 do 30 j=1,nbloks i = nbp1 - j nrow = integs(1,i) ncol = integs(2,i) last = integs(3,i) index = index - nrow*ncol indexb = indexb - nrow indexx = indexx - last 30 call subbak(bloks(index),ipivot(indexb),nrow,ncol,last, * x(indexx)) return end subroutine subfor ( w, ipivot, nrow, last, b, x ) c carries out the forward pass of substitution for the current block, c i.e., the action on the right side corresponding to the elimination c carried out in f a c t r b for this block. c at the end, x(j) contains the right side of the transformed c ipivot(j)-th equation in this block, j=1,...,nrow. then, since c for i=1,...,nrow-last, b(nrow+i) is going to be used as the right c side of equation i in the next block (shifted over there from c this block during factorization), it is set equal to x(last+i) here. c c parameters c w, ipivot, nrow, last are as on return from factrb. c b(j) is expected to contain, on input, the right side of j-th c equation for this block, j=1,...,nrow. c b(nrow+j) contains, on output, the appropriately modified right c side for equation j in next block, j=1,...,nrow-last. c x(j) contains, on output, the appropriately modified right c side of equation ipivot(j) in this block, j=1,...,last (and c even for j=last+1,...,nrow). c integer ipivot(nrow), ip,jmax,k c dimension b(nrow + nrow-last) real w(nrow,last),b(1),x(nrow) ip = ipivot(1) x(1) = b(ip) if (nrow .eq. 1) go to 99 do 15 k=2,nrow ip = ipivot(k) jmax = amin0(k-1,last) sum = 0. do 14 j=1,jmax 14 sum = w(ip,j)*x(j) + sum 15 x(k) = b(ip) - sum c c transfer modified right sides of equations ipivot(last+1),..., c ipivot(nrow) to next block. nrowml = nrow - last if (nrowml .eq. 0) go to 99 lastp1 = last+1 do 25 k=lastp1,nrow 25 b(nrowml+k) = x(k) 99 return end subroutine subbak ( w, ipivot, nrow, ncol, last, x ) c carries out backsubstitution for current block. c c parameters c w, ipivot, nrow, ncol, last are as on return from factrb. c x(1),...,x(ncol) contains, on input, the right side for the c equations in this block after backsubstitution has been c carried up to but not including equation ipivot(last). c means that x(j) contains the right side of equation ipi- c vot(j) as modified during elimination, j=1,...,last, while c for j .gt. last, x(j) is already a component of the solut- c ion vector. c x(1),...,x(ncol) contains, on output, the components of the solut- c ion corresponding to the present block. c integer ipivot(nrow),last, ip,j,k,kp1 real w(nrow,ncol),x(ncol), sum k = last ip = ipivot(k) sum = 0. if (k .eq. ncol) go to 4 kp1 = k+1 2 do 3 j=kp1,ncol 3 sum = w(ip,j)*x(j) + sum 4 x(k) = (x(k) - sum)/w(ip,k) if (k .eq. 1) return kp1 = k k = k-1 ip = ipivot(k) sum = 0. go to 2 end subroutine dtblok ( bloks, integs, nbloks, ipivot, iflag, * detsgn, detlog ) c computes the determinant of an almost block diagonal matrix whose c plu factorization has been obtained previously in fcblok. c *** the logarithm of the determinant is computed instead of the c determinant itself to avoid the danger of overflow or underflow c inherent in this calculation. c c parameters c bloks, integs, nbloks, ipivot, iflag are as on return from fcblok. c in particular, iflag = (-1)**(number of interchanges dur- c ing factorization) if successful, otherwise iflag = 0. c detsgn on output, contains the sign of the determinant. c detlog on output, contains the natural logarithm of the determi- c nant if determinant is not zero. otherwise contains 0. c integer integs(3,nbloks),ipivot(1),iflag, i,indexp,ip,k,last real bloks(1),detsgn,detlog c detsgn = iflag detlog = 0. if (iflag .eq. 0) return index = 0 indexp = 0 do 2 i=1,nbloks nrow = integs(1,i) last = integs(3,i) do 1 k=1,last ip = index + nrow*(k-1) + ipivot(indexp+k) detlog = detlog + alog(abs(bloks(ip))) 1 detsgn = detsgn*sign(1.,bloks(ip)) index = nrow*integs(2,i) + index 2 indexp = indexp + nrow return end