Modular Analysis of Almost Block Diagonal Systems of Equations
1 Modular Analysis of Almost Block Diagonal Systems of Equations
Tarek M. A. El-Mistikawy
Department of Engineering Mathematics and Physics, Cairo University, Giza 12211, Egypt
Abstract
Almost block diagonal linear systems of equations can be exemplified by two modules. This makes it possible to construct all sequential forms of band and/or block elimination methods, six old and fourteen new. It allows easy assessment of the methods on the basis of their operation counts, storage needs, and admissibility of partial pivoting. It unveils a robust partial pivoting strategy- local pivoting . Extension of modular analysis to bordered systems is also included.
Keywords : almost block diagonal systems; sequential solution methods; modular analysis; partial pivoting; local pivoting; bordered systems
1. INTRODUCTION
Systems of equations with almost block diagonal (ABD) matrix of coefficients are frequently encountered in numerical solutions of sets of ordinary or partial differential equations. They occur, for example, when Keller’s box scheme is applied: to parabolic equations, and to elliptic equations with an alternating-direction-implicit approach. Several other situations in which ABD systems occur are described by Amodio et al. [1]. Sequential solution methods for ABD systems were developed, improved, and thoroughly studied to the extent that it was concluded in Ref. [1] that they needed little further study. Traditionally, sequential solution methods of ABD systems performed LU decompositions of the matrix of coefficients G through either band (scalar) elimination or block tridiagonal elimination. The famous COLROW algorithm [2], which is highly regarded for its performance, was incorporated in several applications [3,4,5,6,7]. It utilizes Lam’s alternating column/row pivoting [8] and Varah’s correspondingly alternating column/row scalar elimination [9]. The efficient block tridiagonal methods included Keller’s Block Tridiagonal Row elimination method [10, §5, case i], and El-Mistikawy’s Block Tridiagonal Column elimination method [11]. Both methods could apply a suitable form of Keller’s mixed pivoting strategy [10], which is more expensive than Lam’s. The present paper is intended to explore other variants of the LU decomposition of G . It does not follow the traditional approaches of treating the matrix of coefficients as a banded matrix, or casting it into a block tridiagonal form. It, rather, adopts a new approach, modular analysis , which offers a simple and unified way of expressing and assessing solution methods for ABD systems. The matrix of coefficients G (or, more specifically, its significant part containing the non-zero blocks) is disassembled into an ordered set of modules. (In fact, two different sets of modules are identified.) Each module is an entity that has the following components: a stem , two fins and , a head , and a tail . By arranging the modules in such a way that the head of a module is added to the tail of the next, the significant part of G can be re-assembled. The module exemplifies the matrix, but is much easier to analyze. 2 All possible methods of LU decomposition of G can be formulated as decompositions of . Fourteen new methods are thus discovered. The operation counts and storage needs are easily estimated, revealing the method with the best performance on each account. The validity and stability of the elimination methods are of primary concern to both numerical analysts and algorithm users. Validity means that division by a zero is never encountered, whereas stability guards against round-off-error growth. To insure validity and achieve stability, pivoting is called for [12]. Full pivoting is computationally expensive. It requires full two dimensional search for the pivots. Moreover, it destroys the banded form of the matrix of coefficients. Partial pivoting strategies, though potentially less stable, are considerably less expensive. Uni-directional (row or column) pivoting makes a minor change to the form of G by introducing few extraneous elements. Lam’s alternating pivoting [8], which involves alternating sequences of row pivoting and column pivoting, maintains the form of G . When G is nonsingular, Lam’s pivoting guarantees validity, and if followed by correspondingly alternating elimination it produces multipliers that are bounded by unity; thus enhancing stability. This approach was proposed by Varah [9] in his L G ~ U decomposition method. It was developed afterwards into a more efficient LU version that was adopted by the well-known COLROW solver [2]. The present approach of modular analysis shows that Lam’s pivoting (with Varah’s arrangement) applies to five other elimination methods. A more robust, though more expensive strategy- Local Pivoting , is hereinafter introduced. It performs full pivoting over the same segments of G (or ) to which Lam’s pivoting is applied.
2. PROBLEM DESCRIPTION
Consider the almost block diagonal system of equations gz G whose augmented matrix of coefficients g GG has the form shown in Fig. 1, where the blocks with leading character A , B , C , and D have m , n , m , and n columns, respectively. The blocks with leading character g have r columns indicating as many right-hand sides. The trailing character m or n (and subsequently p = m + n ) indicates the number of rows of a block; or the order of a square matrix such as an identity matrix I , a zero matrix , and a lower L or an upper U triangular matrix. Blanks indicate zero blocks. JJJ JJJJJ JJJJJ jjjjj jjjjj gnBnAn gmDmCmBmAm gnDnCnBnAn gmDmCmBmAm gnDnCnBnAn gmDmCmBmAm gnDnCnBnAn gmDmCm
11 11111 111 22211 12211 111 G Fig. 1: The augmented matrix of coefficients 3 The matrix of unknowns z is, likewise, written as ttJtJtjtjtttt znzmznzmznzmznzm z where the superscript t denotes the transpose. Such a system of equations results, for example, in the finite difference solution of p first order ODEs on a grid of J points; with m conditions at one side to be marked with j , and n conditions at the other side that is to be marked with Jj . Then, each column of the submatrix ttjtj znzm contains the p unknowns of the j th grid point corresponding to a right-hand side.
3. MODULAR ANALYSIS
The description of the decomposition methods for the augmented matrix of coefficients G can be made easy and concise through the introduction of modules of G . Two different modules can be identified as described below. This module, for j =1→ J -1, takes the form
111 11 jjjjj jjjjj jjjAjAjAj gmDmCmBmAm gnDnCnBnAn gmDmCm where, as a rule, the dotted line defines the partitioning to left and right entities. The dashed lines define the partitioning AjAj AjAjAj of Aj to the following components. The stem jj jjAj BnAn DmCm The head jjjAjAjAj gmDmCm The fins jjAj
BmAm jjj jAjAjAj gnDnCn gm
11
The module has the tail jjjAjAjAj gmDmCm which is defined through the head-tail relation jjjAjAj gmDmCm This makes it possible to construct the significant part of G by arranging the set of modules in such a way that the tail of Aj adds to the head of Aj . Minor adjustments need only to be invoked at both ends of G . Specifically, we define the truncated modules AA and AJAJAJ . (For convenience, we shall occasionally drop the subscript and/or superscript identifying a module and its components, as well as the subscript identifying its blocks. This should not cause any confusion.) The head of the module is yet to be defined. It is taken to be related to the other components of by (3.1) in order to allow for decompositions of having the form (3.2) The generic relations = , = , = , and = then hold; leading to = = ( )( )= ( ) = as defined in (3.1). Further, writing gives = = = and = We are now ready to present several old as well as new (marked with a square ■ ) elimination methods in terms of decompositions of the left module Aj . Several inflections of the blocks of G are involved and are defined in Appendix A. The sequence in which the blocks are manipulated for: decomposing the stem, processing the fins, and handling the head (evaluating the head and applying the head-tail relation to determine the tail of the succeeding module Aj ), is mentioned along with the equations (from Appendix A) involved. The correctness of the decompositions may be checked by carrying out the matrix multiplications, using the equalities of Appendix A, and comparing with the un-decomposed form of the module. 5 These methods perform scalar decomposition of both pivotal blocks Cm and Bn . The triangular matrices L and U appear explicitly, marked if unit diagonal with a circumflex ˆ . Four methods are involved. Scalar Row/Scalar Row (SRSR) Elimination (Case iii of [10, §5]) nDnCUnmDUmmBmA nLnAmL A ˆˆ Scalar Row/Scalar Column (SRSC) Elimination ■ nDnCnUmDUmmBmA LnnAmL A ˆˆ Scalar Column/Scalar Row (SCSR) Elimination (Method implemented by the COLROW algorithm) nDnCUnmDmUmBmA nLnALm A ˆˆ Scalar Column/Scalar Column (SCSC) Elimination (Case iv of [10, §5]) nDnCnUmDmUmBmA LnnALm A ˆˆ All four scalar elimination methods apply the following sequence of manipulations. Stem:
LmUm (A6), Dm ″(A19), An ′(A9), Bn (A2b), LnUn (A5) Fins: Am ′(A7), Bm (A1b), Bm ″(A17), Cn ′(A11), Dn ′(A13) Head: Cm (A3b), Dm (A4b) The difference is in how the decompositions of the pivotal blocks Cm and Bn are carried out, through row ( UL ˆ ) or column ( UL ˆ ) manipulation. 6 These methods perform scalar decomposition of Cm and block decomposition of Bn . It involves four methods. Scalar Row/Block Row (SRBR) Elimination ■ DnCnnBmDUmBmmA InnAmL A *ˆ Scalar Column/Block Row (SCBR) Elimination ■ DnCnnBDmmUBmmA InnALm A ˆ* Both SRBR and SCBR methods apply the following sequence of manipulations. Stem: mC (A6), An ′(A9), Dm ″(A19), Bn (A2b), nB (A5) Fins: Am ′(A7) , Bm (A1b), Bm ″(A17), Bm* (A18) Head: Cm (A3a), Dm (A4a) The difference is in how the decompositions of the pivotal block Cm is carried out, through row ( UL ˆ ) or column ( UL ˆ ) manipulation. Scalar Row/Block Column (SRBC) Elimination ■ DnCnInmDUmBmmA nBnAmL A ˆ Scalar Column/Block Column (SCBC) Elimination ■ DnCnInDmmUBmmA nBnALm A ˆ Both SRBC and SCBC methods apply the following sequence of manipulations. Stem: mC (A6), An ′(A9), Dm ″(A19), Bn (A2b), nB (A5) Fins: Am ′(A7), Bm (A1b), Cn ′(A11), Dn ′(A13), Cn° (A12),
Dn° (A14) Head: Cm (A3c), Dm (A4c) The difference is in how the decompositions of the pivotal block Cm is carried out, through row ( UL ˆ ) or column ( UL ˆ ) manipulation. 7 These methods perform block decomposition of Cm and scalar decomposition of Bn . It involves four methods. Block Row/Scalar Row (BRSR) Elimination ■ nDnCUnDmmCmBAm nLAnIm A ˆ Block Row/Scalar Column (BRSC) Elimination ■ nDnCnUDmmCmBAm LnAnIm A ˆ Both BRSR and BRSC methods apply the following sequence of manipulations. Stem: mC (A6), An ′(A9), An° (A10), Bn (A2c), nB (A5) Fins: Am ′(A7), Am° (A8), Bm (A1c), Bm ″(A17) Head: Cm (A3b), Dm (A4b) The difference is in how the decompositions of the pivotal block Bn is carried out, through row ( UL ˆ ) or column ( UL ˆ ) manipulation. Block Column/Scalar Row (BCSR) Elimination ■ nDnCUnDmImmBAm nLAnmC A *ˆ Block Column/Scalar Column (BCSC) Elimination ■ nDnCnUDmImmBAm LnAnmC A ˆ * Both BCSR and BCSC methods apply the following sequence of manipulations. Stem: mC (A6), Dm ″(A19), Dm* (A20), Bn (A2a), nB (A5) Fins: Bm (A1a), Bm ″(A17), Cn ′(A11), Dn ′(A13) Head: Cm (A3b), Dm (A4b) The difference is in how the decompositions of the pivotal block Bn is carried out, through row ( UL ˆ ) or column ( UL ˆ ) manipulation. 8 These methods perform block decomposition of both pivotal blocks Cm and Bn , in which the decomposed pivotal blocks mC (≡ LmUm ) and nB (≡ LnUn ) appear. It involves four methods.
Block Row/Block Row (BRBR) Elimination ■ DnCnnBDmmCBmAm InAnIm A * Stem: mC (A6), An ′(A9), An° (A10), Bn (A2c), nB (A5) Fins: Am ′(A7), Am° (A8), Bm (A1c), Bm ″(A17), Bm* (A18) Head: Cm (A3a), Dm (A4a) Block Row/Block Column (BRBC) Elimination ■ DnCnInDmmCBmAm nBAnIm A Stem: mC (A6), An ′(A9), An° (A10), Bn (A2c), nB (A5) Fins: Am ′(A7), Am° (A8), Bm (A1c), Cn ′(A11), Dn ′(A13), Cn° (A12),
Dn° (A14) Head: Cm (A3c), Dm (A4c) Block Column/Block Row (BCBR) Elimination ■ DnCnnBDmImBmAm InAnmC A ** Stem: mC (A6), Dm ″(A19), Dm* (A20), Bn (A2a), nB (A5) Fins: Bm (A1a), Bm ″(A17), Bm* (A18) Head: Cm (A3a), Dm (A4a) Block Column/Block Column (BCBC) Elimination ■ DnCnInDmImBmAm nBAnmC A * Stem: mC (A6), Dm ″(A19), Dm* (A20), Bn (A2a), nB (A5) Fins: Bm (A1a), Cn ′(A11), Dn ′(A13), Cn° (A12),
Dn° (A14) Head: Cm (A3c), Dm (A4c) 9 These methods perform the identity decomposition Ip or Ip , where the decomposed stem stands for any of the nonidentity (scalar or block) decompositions given above. This is an arbitrariness that can be put to advantage, as shown below. Aligned Block-Tridiagonal Row (ABTR) Elimination (
Case i of [10, §5] ) DnCnImBmAm Ip AA ** Traditionally ABTR elimination decomposes the stem A through scalar row manipulation. However, as indicated in §5 on pivoting, it would prove profitable to decompose the stem through an alternating manipulation method as in !C?R elimination, where ! and ? stand for S or B. Using SCSR elimination, the following sequence of manipulations applies. Stem: mULm ˆ (A6), Dm ″(A19), An ′(A9), Bn (A2b), L ˆ nUn (A5) Fins: Am ′(A7), Bm (A1b), Bm ″(A17), Bm* (A18), Am ″(A15), Am* (A16) Head: Cm (A3a), Dm (A4a) Aligned Block-Tridiagonal Column (ABTC) Elimination (
Case ii of [10, §5] ) DnCn DmCmIpBmAm eeAA The method introduces the two extraneous blocks e Cm and e Dm . Traditionally ABTC elimination (now, ABTCt) decomposes the stem A through scalar row manipulation as in SRSR elimination. Stem: L ˆ mUm (A6), Dm ″(A19), An ′(A9), Bn (A2b), L ˆ nUn (A5) Fins: Cn ′(A11), Dn ′(A13), Cn° (A12),
Dn° (A14), e Cm (A25b), e Dm (A26b) Head: Cm (A29), Dm (A30) A more efficient sequence (defining ABTCe) decomposes the stem through block column manipulation as in BCBC elimination. Stem: mC (A6), Dm ″(A19), Dm* (A20), Bn (A2a), nB (A5) Fins: Cn ′(A11), Dn ′(A13), Cn° (A12),
Dn° (A14), e Cm (A25a), e Dm (A26a) Head: Bm (A1a), Cm (A3c), Dm (A4c) Note that the use of (A3c) and (A4c) for evaluating the head is a result of substituting expressions (A25a) and (A26a) for e Cm and e Dm into (A29) and (A30), respectively. 11 The second module of G takes, for j =1→ J -2, the form
211 111 jjj jjj jjjj jjjjDjDjDj gmBmAm gnBnAn gmDmCmBm gnDnCnBn with the partitioning DjDj DjDjDj which defines the following components of Dj . The stem jj jjDj CmBm CnBn The head
21 11 jj jjDjDjDj gmBm gnBn
The fins jjDj AmAn0n and jj jjDjDjDj gmDm gnDn The tail jj jjDjDjDj gmBm gnBn The head-tail relation for Dj takes the form jj jjDjDj gmBm gnBn
11 so that, with minor end adjustments involving the use of the truncated modules
211 111 1110 gmBmAm gnBnAn gmDmCm D , JJJ JJJJ JJJJDJ gnBnAn gmDmCmBm gnDnCnBn , and
JJDJ gnBn the significant part of G can be constructed by arranging the modules in such a way that the tail of Dj adds to the head of Dj . The head of D is defined in accordance with relation (3.1) to allow for decompositions of the form (3.2). Such decompositions of D lead to all but the aligned block-tridiagonal elimination methods given above; e.g., Scalar Column/Scalar Column (SCSC) Elimination mDmU nDnCnUmA nALmmBLn D ˆˆ Stem: Ln nU ˆ (A5), Cn ′(A11), Bm ″(A17), Cm (A3b), Lm mU ˆ (A6) Fins: Dn ′(A13), Dm (A4b), Dm ″(A19), An ′(A9), Am ′(A7) Head: Bn (A2b), Bm (A1b) Block Row/Block Row (BRBR) Elimination ■ * DmmC DnCnnBAmAnImBmIn D Stem: nB (A5), Bm ″(A17), Bm* (A18), Cm (A3a), mC (A6) Fins: Dm (A4a), An ′(A9), Am ′(A7), An° (A10),
Am° (A8) Head: Bn (A2c), Bm (A1c) Decomposition of D leads to two other block tridiagonal elimination methods, that cannot be generated from A . 12 Displaced Block-Tridiagonal Row (DBTR) Elimination ■ DmDnAmBm AnBn Ip
DeeD The method introduces the two extraneous blocks e Bn and e Bm . An efficient sequence of its calculations decomposes the stem D through block row manipulation as in BRBR elimination. Stem: nB (A5), Bm ″(A17), Bm *(A18), Cm (A3a), mC (A6) Fins: An ′(A9), Am ′(A7), An° (A10),
Am° (A8) Head: Dm (A4a), Bn (A2c), Bm (A1c) Note that e Bm and e Bn are not evaluated, since they need not to be stored. Their contribution to the head is estimated through substituting their expressions of (A23) and (A24) into (A27) and (A28) leading to (A2c) and (A1c), respectively. Displaced Block-Tridiagonal Column (DBTC) Elimination ** DmDnIpAmAn0n DD The method, as described in Ref. [11], decomposes the stem D through scalar column manipulation as in SCSC elimination However, as with ABTR elimination, it would prove profitable to decompose the stem through an alternating manipulation method as in !C?R elimination, where ! and ? stand for S or B. Using SCSR elimination, the following sequence of manipulations applies. Stem: Ln nU ˆ (A5), Cn ′(A11), Bm ″(A17), Cm (A3b), L ˆ mUm (A6) Fins: Dn ′(A13), Dm (A4b), Dm ″(A19), Dm *(A20), Dn ″(A21), Dn *(A22) Head: Bn (A2a), Bm (A1a) We observe that a decomposition LUG that is based on any of the elimination methods described in §3.1 and §3.2 can be constructed from the corresponding decomposition jjj of the module. The significant part of L (or U ) is obtained by arranging the j ’s (or j ’s), sequentially, in order. The procedure for solving the matrix equation gz G , which can be described in terms of manipulation of the augmented matrix g GG , can, similarly, be described in terms of 13 manipulation of the augmented module jjj . The manipulation of G applies a forward sweep which corresponds to the decomposition LUG Z UL , that is followed by a backward sweep which corresponds to the decomposition z I UU . Similarly, the manipulation of j applies a forward sweep involving two steps. The first step performs the decomposition jjj . The second step evaluates the head jjj , then applies the head-tail relation to determine the tail of j . In a backward sweep, two steps are applied to jjjj leading to the solution module j z ( ttjtjAj znzmz , ttjtjDj zmznz ). With j z known, the first step uses ʓ j ( ʓ Aj = Aj z , ʓ Dj = j zn ) in the back substitution relation jj ʓ j j to contract j to jjj . The second step solves jjj z for j z which is equivalent to the decomposition jjj zIp .
4. OPERATION COUNTS AND STORAGE NEEDS
The operation counts are measured by the number of multiplications (mul) with the understanding that the number of additions is comparable. The storage needs are measured by the number of locations (loc) required to store arrays calculated in the forward sweep for use in the backward sweep; provided that the elements of G are not stored but are generated when needed, as is usually done in the numerical solution of a set of ODEs, for example. The modules introduced in §3 allow easy evaluation of the elimination methods. Per module (i.e., per grid point), each method requires, for the manipulation of G , as many operations as it requires to manipulate , and as many storage locations as it requires to store . The manipulation of jjj involves manipulating j which requires the decomposition of the stem j , the processing of the fins j and j , and the handling of the head j ; as well as manipulating j which requires the evaluation of j then j in the forward sweep, and j then j z in the backward sweep. All methods require pp (mul) for decomposing the stem. To illustrate, a scalar decomposition of the stem Aj involves the use of Eqs. (A6,19,9,2,5) which requires mm , nmm , nmm , mn , nn (mul), respectively, totaling pp (mul); whereas a scalar decomposition of the stem Dj involves the use of Eqs. (A5,11,17,3,6) which requires nn , mnn , mnn , nm , mm (mul), respectively, totaling, again, pp (mul). Scalar/Scalar elimination methods apply such scalar decomposition of j . 14 BRBR\BCBC\BCBR\BRBC elimination applies block decomposition of Aj or Dj which involves the use of Eqs. (A6,9\19\19\9,10\20\20\10, 2,5) or (A5,17\11\17\11, 18\12\18\12,3,6), requiring, respectively, mm , nmm , nmm , mn , nn (mul) or nn , mnn , mnn , nm , mm (mul), totaling in either case pp (mul). ABT\DBT row or column elimination applies an identity decomposition involving either scalar or block decomposition of Aj \ Dj . All methods, which introduce no extraneous elements, require pmn (mul) for evaluating the head jjj in preparation for evaluating the tail j with the help of the head-tail relation, which involves no multiplications. Moreover, all such methods require rp (mul) to handle the right module j : The determination of the head jjj in the forward sweep requires pmr (mul), while the back substitution relation jj ʓ j j in the backward sweep requires pnr (mul); adding up to rp (mul). The solution of jjj for j in the forward sweep, and the solution of jjj z for j z in the backward sweep, together require rp (mul), irrespective of the decomposition jjj of the stem. These methods differ only in the operation counts for processing the fins j and j , with BCR elimination requiring the least count pmn (mul). The only two methods, ABTC and DBTR elimination, which introduce extraneous nonzero elements through elimination, perform additional calculations for evaluating the extraneous blocks and/or the blocks they contribute to. As for the storage needs, all methods require pr (loc) to store j . They differ in the number of locations needed for storing the j ’s presented in §3, with DBTC elimination requiring the least number pn (loc). Note that, in methods involving scalar elimination, square blocks need to be reserved for storing the triangular blocks Um and/or Un . Table (1B) of Appendix B contains the above information, allowing for clear comparison among the methods. It is noted that some methods show preference to m > n and others to m < n . Taking m > n , a fair comparison is possible by interchanging m and n for the methods which prefer m < n . (When solving a set of ODEs, this amounts to assigning j =1 to the side with n boundary conditions.) The four methods SCSR, SRSC, BCBR, and BRBC whose operation counts are not biased toward m or n can be applied on two processors of a parallel machine; one operating from j =1 and the other from j = J , with equal efficiency, considerably reducing the calculation time. The savings achieved by BCBR elimination in operations and by DBTC elimination in storage are of leading order significance when p >>1, in the two distinguished limits m ~ n ~ p /2 and ( m , n )~( p ,1). In comparison with the scalar elimination methods, savings of at least ~ p (mul) and ~ p (loc) are achieved in the former limit, while savings of at least ~ p (mul) and ~ p (loc) are achieved in the latter limit. 15
5. PIVOTING
The role of pivoting is to guarantee validity and stability. Partial pivoting (aiming at evading division by a zero) is all one needs to insure validity. It is also generally accepted that partial pivoting (aiming further at binding the multipliers) leads to stability [12]. Partial pivoting strategies involve either one-dimensional search for the pivot (e.g., row pivoting, column pivoting, and Lam’s alternating pivoting) or the more expensive two-dimensional search (e.g., Keller’s mixed pivoting [10] and the to-be-introduced local pivoting). Lam’s alternating pivoting [8] is given, here, special attention. It involves one-dimensional search for the pivot, but unlike uni-directional (row or column) pivoting, it introduces no extraneous elements and, furthermore, it guarantees validity. Expressed in terms of the present notation, Lam’s strategy applies column pivoting to
DmCm A in order to form and decompose a nonsingular pivotal block Cm , and applies row pivoting to tttD BmBn in order to form and decompose a nonsingular pivotal block Bn . These are valid processes since A is of rank m and D is of rank n , as can be shown following the reasoning of Keller [10, §5, Theorem]. It is emphasized here that applying column interchanges within A is mandatory for guaranteed validity. Row interchanges among the m linearly independent rows of A do not guarantee validity as the first m columns may be linearly dependent. Likewise, row interchanges within D are mandatory. As advocated by Varah [9], Lam’s alternating pivoting should be followed by correspondingly alternating elimination. Step by step, column pivoting is to be followed by column elimination and row pivoting by row elimination. This enhances stability since the multipliers are then bounded by unity. Obviously, SCSR, SCBR, BCSR and BCBR elimination methods can benefit from this. So can do ABTR and DBTC elimination methods- contrary to the common belief- if the decomposition of the stem is carried out as in SCSR, SCBR, BCSR or BCBR elimination. To enhance stability further we introduce the Local Pivoting Strategy which applies full pivoting (maximum pivot strategy) to the segments A and D . Only those methods which inflect Dm to at least mD and Bm to at least mB can apply local pivoting. (Note that Keller’s mixed pivoting [10], if interpreted as applying full pivoting to A and row pivoting to D , is midway between Lam’s and local pivoting.) Table (1B) indicates the admissibility of Lam’s and/or local pivoting by an asterisk *.
6. BORDERED AND DOUBLE-BORDERED SYSTEMS
Extension of modular analysis to bordered (BABD) and double-bordered (DBABD) almost block diagonal systems of equations is possible. For example, for a DBABD system with right-bordered array of q columns tttJtJtjtjtt HkHnHmHnHmHnHm and lower-bordered array of k (≤ n + q ) rows gkHkBkAkBkAkBkAk JJjj
16 we define the A-module jjjjjj jjjjjj jjjjjj jjjjAjAj AjAjAj gkHkBkAkBkAk gmHmDmCmBmAm gnHnDnCnBnAn gmHmDmCm
Its tail jjjj jjjjAj gkHkBkAk gmHmDmCm having the hidden segment jj gkHk which does not appear in Aj , is determined through the head-tail relation jjjj jjjjAjAj gkHkBkAk gmHmDmCm All twenty elimination methods given in §3 are applicable. Handling the new bordered arrays requires pq (loc) but different numbers of (mul). The methods that do not introduce extraneous elements require pqrpqqp (mul), when k = q .
6. CONCLUSION
Using the novel approach of modular analysis, we have analyzed the sequential solution methods for almost block diagonal systems of equations. Two modules have been identified and have made it possible to express and assess all possible band and block elimination methods on the basis of their operation counts, storage needs, and admissibility of partial pivoting. Modular analysis has also been extended to bordered almost block diagonal systems. Implementation of the four distinguished methods SCSR, BCSR, BCBR and DBTC within the COLROW platform [13], which was designed to give SCSR elimination its best performance, showed that the other three methods could outperform SCSR, in some cases (when m >> n ). FORTRAN codes and related materials are given in Appendix C. REFERENCES [1] Amodio P, Cash JR, Roussos G, Wright RW, Fairweather G, Gladwell I, Kraut GL, Paprzycki M. Almost block diagonal linear systems: sequential and parallel solution techniques, and applications.
Numer Linear Algebr , 7, 275-317, 2000. [2] Diaz JC, Fairweather G, Keast P. FORTRAN packages for solving certain almost block diagonal linear systems by modified alternate row and column elimination.
ACM T Math Software , 9, 358-375, 1983. [3] Cash JR, Wright RW. A deferred correction method for nonlinear two-point boundary value problems.
SIAM J Sci Stat Comp , 12, 971-989, 1991 17 [4] Cash JR, Moore G, Wright RW. An automatic continuation strategy for the solution of singularly perturbed linear boundary value problems.
J Comput Phys , 122, 266-279, 1995. [5] Enright WH, Muir PH. Runge-Kutta software with defect control for boundary value ODES.
SIAM J Sci Comput ,17, 479-497, 1996. [6] Keast P, Muir PH. Algorithm 688: EPDCOL: A more efficient PDECOL code.
ACM T Math Software , 17, 153-166, 1991. [7] Wright RW, Cash JR, Moore G. Mesh selection for stiff two-point boundary value problems.
Numer Algorithms , 7, 205-224, 1994. [8] Lam DC.
Implementation of the box scheme and model analysis of diffusion-convection equations . PhD Thesis, University of Waterloo: Waterloo, Ontario, Canada; 1974. [9] Varah JM. Alternate row and column elimination for solving certain linear systems.
SIAM J Numer Anal , 13:71-75, 1976. [10] Keller HB. Accurate difference methods for nonlinear two-point boundary value problems.
SIAM J Numer Anal , 11, 305-320, 1974. [11] El-Mistikawy TMA. Solution of Keller’s box equations for direct and inverse boundary-layer problems.
AIAA J , 32, 1538-1541, 1994. [12] Duff IS, Erisman AM, Reid JK.
Direct Methods for Sparse Matrices . Clarendon Press: Oxford, UK; 1986. [13] El-Mistikawy TMA. Modular Analysis of Sequential Solution Methods for Almost Block Diagonal Systems of Equations.
Advances in Numerical Analysis , 2013, 563872, 2013.
Appendix A: Inflections of Blocks of G In §3, two modules, A and S , of the matrix of coefficients G are introduced and decomposed to generate the elimination methods. The process involves inflections of the blocks of G , which proceed for a block E , say, according to the following scheme. * EEEE EEE The following equalities are to be used to determine the blocks with underscored leading character. The number of multiplications involved is given between wiggly brackets. When a ± or ∓ sign appears, the upper (or lower) sign corresponds to the triangular matrix U (or L ) involved being unit diagonal. Blocks- E * DmAmmDmAAmDmmBBm cba { nm } (A1) * DmAnmDnAAnDmnBBn cba { mn } (A2) CnBmnCmBCnBmmCCm cba * { nm } (A3) DnBmnDmBDnBmmDDm cba * { mn } (A4) 18 Blocks- E (Decomposed blocks) )( nUnLnBBn { nn } )( mUmLmCCm { mm } (A5,6) Blocks-and EE AmUmmA { mmm } mALmmA { mmm } (A7,8) AnUmnA { nmm } nALmnA { nmm } (A9,10) CnnCLn { mnn } nCnCUn { mnn } (A11,12) DnnDLn { nnn } nDnDUn { nnn } (A13,14) -Blocks*and EE nABmmAmA * { nm } mALmmA * { mmm } (A15,16) BmUnmB { mnn } mBLnmB * { mnn } (A17,18) DmmDLm { nmm } mDmDUm * { nmm } (A19,20) * DmnCnDnD { mn } nDnDUn * { nnn } (A21,22) Blocks- e E and related Blocks- E * BmAmmB e { nm } * BmAnnB e { mn } (A23,24) CnDmmC e * { nm } CnmDmCUm e { nm + }2 mmm (A25a,b) DnDmmD e * { mn } DnmDmDUm e { mn + nmm } (A26a,b) DmAmDnBmmBBm e { mn + nm } (A27) DmAnDnBnnBBn e { n + mn } (A28) BmCnAmCmmCCm e { m + nm } (A29) BmDnAmDmmDDm e { nm + mn } (A30) 19 Appendix B: Assessment Table
The operation counts, storage needs, and admissibility of Lam's and local pivoting of all twenty elimination methods are contained in Table (B1) given in this appendix. Note that the methods involving scalar elimination reserve square blocks for storing the triangular blocks Um and/or Un . The expressions given in the footnotes of Table (B1) should, therefore, replace the corresponding table entries. Table (B1) : Assessment of Elimination Methods Method Operation Counts (mul) Storage Needs (loc) Pivoting rppmnpp pnpr Lam’s Local SRSR nmnm nmp a SCSC nmnm nmp a SCSR nmnm nmp a SRSC ■ nmnm nmp a SRBR ■ mm mnp a SCBR ■ mm mnp a SRBC ■ mnm mnp b SCBC ■ mnm mnp b BRSR ■ nnm nmp a BRSC ■ nnm nmp a BCSR ■ nn nmp c BCSC ■ nn nmp c BRBR ■ m mnp BCBC ■ n mn BCBR ■ pn BRBC ■ nm pm ABTR pm p DBTC pn DBTR ■ mrnpm )( p ABTCe pmrnpmn pm ABTCt pmrmnmp pm a mnp b pm c pn Appendix C: Supplementary Materials
The following supplementary materials include listings of four FORTRAN codes: SCSR.FOR; a Fortran program for the SCSR elimination method, BCSR.FOR; a Fortran program for the BCSR elimination method, BCBR.FOR; a Fortran program for the BCBR elimination method, DBTC.FOR; a Fortran program for the DBTC elimination method. SCSR.FOR is essentially the COLROW algorithm, without FLAG statements and with renaming of some variables. BCSR.FOR, BCBR.FOR, and DBTC.FOR are modifications of SCSR.FOR which apply BCSR, BCBR, and DBTC eliminations, instead of SCSR elimination. In the four codes, the integer variables: LM is the number of rows in the top block, LN is the number of rows in the bottom block, KK is the number of repetitions of the repeated block, and MM is the number of times the solution is to be carried out. They appear in PARAMETER statements in the main program and in the two subroutines, where they currently assume the values LM=10, LN=1, KK=10, MM=1000000. A change in the values of LM and LN should be applied to all PARAMETER statements, keeping LM+LN=11. The supplementary materials also include two input data listings: TESTA; containing the top and bottom blocks, and TESTB; containing the repeated block. The data in TESTA and TESTB correspond to a coefficient matrix whose entries are chosen randomly. The right hand side column is calculated such that the solution column is (1.0 1.1 1.2 … 1.8 1.9 2.0) t repeated KK+1 times; irrespective of the values of LM and LN. A change in the values of LM and LN can be effected by positioning the line BOTTOM BLOCK in TESTA; according to the new values. Currently, it is positioned such that LM=10, LN=1. SCSR.FOR Fortran program
CCCCC C SCSR CCCCC IMPLICIT REAL*8(A-H,O-Z) PARAMETER(LM=10,LN=1,KK=10,MM=1000000, , LM1=LM+1,LN1=LN+1,KK1=KK+1, , LC=LM+LN,LC1=LC+1,LA=2*LC,LA1=LA+1,N=LC*KK1) DIMENSION TOPBLK(LM,LC),ARRAY(LC,LA,KK),BOTBLK(LN,LC) DIMENSION IPIVOT(N),B(N) DIMENSION A1(LM,LC1),A2(LN,LC1),A3(LC,LA1) OPEN(1,FILE='TESTA') READ(1,1000) ((A1(I,J),J=1,LC1),I=1,LM) READ(1,1000) ((A2(I,J),J=1,LC1),I=1,LN) CLOSE(1) OPEN(2,FILE='TESTB') READ(2,2000) ((A3(I,J),J=1,LA1),I=1,LC) CLOSE(2) DO 4 M=1,MM DO 1 I=1,LM B(I)=A1(I,LC1) DO 1 J=1,LC J1 = J + 1 JN = J - LM JN1 = JN + 1 INCRJ=INCR+J IPVT = JN ROWMAX = DABS(ARRAY(JN,J,K)) DO 120 I=JN1,LC TEMPIV = DABS(ARRAY(I,J,K)) IF (TEMPIV.LE.ROWMAX) GO TO 120 IPVT = I ROWMAX = TEMPIV 120 CONTINUE INCRP=INCRM+IPVT IPIVOT(INCRJ)=INCRP IF (IPVT.EQ.JN) GO TO 140 DO 130 L=J,LA SWAP = ARRAY(IPVT,L,K) ARRAY(IPVT,L,K) = ARRAY(JN,L,K) ARRAY(JN,L,K) = SWAP 130 CONTINUE SWAP = B(INCRP) B(INCRP) = B(INCRJ) B(INCRJ) = SWAP 140 CONTINUE ROWPIV = ARRAY(JN,J,K) DO 150 I=JN1,LC ARRAY(I,J,K) = ARRAY(I,J,K)/ROWPIV 150 CONTINUE DO 160 L=J1,LA ROWMLT = ARRAY(JN,L,K) DO 160 I=JN1,LC 160 ARRAY(I,L,K) = ARRAY(I,L,K) - ROWMLT*ARRAY(I,J,K) ROWMLT = B(INCRJ) DO 170 I=JN1,LC INCRI=INCRM+I 170 B(INCRI) = B(INCRI) - ROWMLT*ARRAY(I,J,K) 180 CONTINUE DO 310 I=LN1,LC IM = I + LM IM1 = IM + 1 IMC = IM - LC INCRN = INCR + IM IN=I-LN I1 = I + 1 IPVT = IM COLMAX = DABS(ARRAY(I,IPVT,K)) DO 190 J=IM1,LA TEMPIV = DABS(ARRAY(I,J,K)) IF (TEMPIV.LE.COLMAX) GO TO 190 IPVT = J COLMAX = TEMPIV 190 CONTINUE IPIVOT(INCRN) = INCR + IPVT IF (IPVT.EQ.IM) GO TO 240 DO 200 L=I,LC SWAP = ARRAY(L,IPVT,K) ARRAY(L,IPVT,K) = ARRAY(L,IM,K) ARRAY(L,IM,K) = SWAP 200 CONTINUE IPC = IPVT - LC IF (K.EQ.KK) GO TO 220 DO 210 L=1,LC SWAP = ARRAY(L,IPC,K1) ARRAY(L,IPC,K1) = ARRAY(L,IMC,K1) ARRAY(L,IMC,K1) = SWAP 210 CONTINUE GO TO 240 220 CONTINUE DO 230 L=1,LN SWAP = BOTBLK(L,IPC) BOTBLK(L,IPC) = BOTBLK(L,IMC) BOTBLK(L,IMC) = SWAP 230 CONTINUE 240 CONTINUE CII = ARRAY(I,IM,K) DO 270 J=IM1,LA CIJ = ARRAY(I,J,K)/CII ARRAY(I,J,K) = CIJ DO 250 L=I1,LC 250 ARRAY(L,J,K) = ARRAY(L,J,K) - ARRAY(L,IM,K)*CIJ JC = J - LC IF (K.LT.KK) THEN DO 260 L=1,LC 260 ARRAY(L,JC,K1) = ARRAY(L,JC,K1) -ARRAY(L,IMC,K1)*CIJ ELSE DO 265 L=1,LN 265 BOTBLK(L,JC) = BOTBLK(L,JC) -BOTBLK(L,IMC)*CIJ ENDIF 270 CONTINUE BI=B(INCRN)/CII B(INCRN)=BI DO 280 L=I1,LC LINCRM=L+INCRM 280 B(LINCRM)=B(LINCRM)-ARRAY(L,IM,K)*BI IF (K.LT.KK) THEN DO 290 L=1,LC LINCRB=L+INCRB 290 B(LINCRB)=B(LINCRB)-ARRAY(L,IN,K1)*BI ELSE DO 295 L=1,LN LINCRB=L+INCRB 295 B(LINCRB)=B(LINCRB)-BOTBLK(L,IN)*BI ENDIF 310 CONTINUE INCR = INCRC 320 CONTINUE IF (LN.EQ.1) GO TO 400 INCRM=INCR+LM LC2=LC-1 DO 390 J=LM1,LC2 J1 = J + 1 JN2 = J - LM JN21 = JN2 + 1 INCRJ=INCR+J IPVT = JN2 ROWMAX = DABS(BOTBLK(JN2,J)) DO 330 I=JN21,LN TEMPIV = DABS(BOTBLK(I,J)) IF (TEMPIV.LE.ROWMAX) GO TO 330 IPVT = I ROWMAX = TEMPIV 330 CONTINUE INCRP=INCRM+IPVT IPIVOT(INCRJ)=INCRP IF (IPVT.EQ.JN2) GO TO 350 DO 340 L=J,LC SWAP = BOTBLK(IPVT,L) BOTBLK(IPVT,L) = BOTBLK(JN2,L) BOTBLK(JN2,L) = SWAP 340 CONTINUE SWAP = B(INCRP) B(INCRP) = B(INCRJ) B(INCRJ) = SWAP 350 CONTINUE ROWPIV = BOTBLK(JN2,J) DO 360 I=JN21,LN 360 BOTBLK(I,J) = BOTBLK(I,J)/ROWPIV DO 370 L=J1,LC ROWMLT = BOTBLK(JN2,L) DO 370 I=JN21,LN 370 BOTBLK(I,L) = BOTBLK(I,L) - ROWMLT*BOTBLK(I,J) ROWMLT = B(INCRJ) DO 380 I=JN21,LN INCRI=INCRM+I 380 B(INCRI) = B(INCRI) - ROWMLT*BOTBLK(I,J) 390 CONTINUE 400 CONTINUE RETURN END SUBROUTINE CRSLVE(TOPBLK,ARRAY,BOTBLK,IPIVOT,B) IMPLICIT REAL*8(A-H,O-Z) PARAMETER(LM=10,LN=1,KK=10, , LM1=LM+1,LN1=LN+1,KK1=KK+1, , LC=LM+LN,LC1=LC+1,LA=2*LC,N=LC*KK1) DIMENSION TOPBLK(LM,LC),ARRAY(LC,LA,KK),BOTBLK(LN,LC) DIMENSION IPIVOT(N),B(N) INCR=KK*LC INCRM=INCR+LM DO 210 L=1,LN J = LC1 - L INCRJ = INCR + J LN1L = LN1 - L B(INCRJ) = B(INCRJ)/BOTBLK(LN1L,J) IF (L.EQ.LN) GO TO 200 BINCRJ = B(INCRJ) LNL = LN - L DO 190 I=1,LNL INCRI = INCRM + I B(INCRI) = B(INCRI) - BOTBLK(I,J)*BINCRJ 190 CONTINUE 200 CONTINUE 210 CONTINUE DO 300 LK=1,KK K = KK1 - LK INCR = INCR - LC DO 240 L1=LN1,LC I = LC + LN1 - L1 IM = I + LM IM1 = IM + 1 INCRN = INCR + IM DOTPRD = B(INCRN) DO 220 J=IM1,LA INCRJ = INCR + J DOTPRD = DOTPRD - ARRAY(I,J,K)*B(INCRJ) 220 CONTINUE B(INCRN) = DOTPRD IPVTN = IPIVOT(INCRN) IF (INCRN.EQ.IPVTN) GO TO 230 SWAP = B(INCRN) B(INCRN) = B(IPVTN) B(IPVTN) = SWAP 230 CONTINUE 240 CONTINUE INCRM = INCR + LM DO 260 J=LC1,LA INCRJ = INCR + J BINCRJ = B(INCRJ) DO 250 I=1,LN INCRI = INCRM + I B(INCRI) = B(INCRI) - ARRAY(I,J,K)*BINCRJ 250 CONTINUE 260 CONTINUE DO 290 L=1,LN J = LC1 - L INCRJ = INCR + J LN1L = LN1 - L B(INCRJ) = B(INCRJ)/ARRAY(LN1L,J,K) IF (L.EQ.LN) GO TO 280 BINCRJ = B(INCRJ) LNL = LN - L DO 270 I=1,LNL INCRI = INCRM + I B(INCRI) = B(INCRI) - ARRAY(I,J,K)*BINCRJ 270 CONTINUE 280 CONTINUE 290 CONTINUE 300 CONTINUE DO 330 L=1,LM I = LM1 - L I1 = I + 1 DOTPRD = B(I) DO 310 J=I1,LC DOTPRD = DOTPRD - TOPBLK(I,J)*B(J) 310 CONTINUE B(I) = DOTPRD IPVTI = IPIVOT(I) IF (I.EQ.IPVTI) GO TO 320 SWAP = B(I) B(I) = B(IPVTI) B(IPVTI) = SWAP 320 CONTINUE 330 CONTINUE RETURN END
BCSR.FOR Fortran program
CCCCC C BCSR CCCCC IMPLICIT REAL*8(A-H,O-Z) PARAMETER(LM=10,LN=1,KK=10,MM=1000000, , LM1=LM+1,LN1=LN+1,KK1=KK+1, , LC=LM+LN,LC1=LC+1,LA=2*LC,LA1=LA+1,N=LC*KK1) DIMENSION TOPBLK(LM,LC),ARRAY(LC,LA,KK),BOTBLK(LN,LC) DIMENSION IPIVOT(N),B(N) DIMENSION A1(LM,LC1),A2(LN,LC1),A3(LC,LA1) OPEN(1,FILE='TESTA') READ(1,1000) ((A1(I,J),J=1,LC1),I=1,LM) READ(1,1000) ((A2(I,J),J=1,LC1),I=1,LN) CLOSE(1) OPEN(2,FILE='TESTB') READ(2,2000) ((A3(I,J),J=1,LA1),I=1,LC) CLOSE(2) DO 4 M=1,MM DO 1 I=1,LM B(I)=A1(I,LC1) DO 1 J=1,LC 1 TOPBLK(I,J)=A1(I,J) DO 2 I=1,LN B(I+N-LN)=A2(I,LC1) DO 2 J=1,LC 2 BOTBLK(I,J)=A2(I,J) DO 3 K=1,KK DO 3 I=1,LC B(I+K*LC-LN)=A3(I,LA1) DO 3 J=1,LA 3 ARRAY(I,J,K)=A3(I,J) CALL CRDCMP(TOPBLK,ARRAY,BOTBLK,IPIVOT,B) CALL CRSLVE(TOPBLK,ARRAY,BOTBLK,IPIVOT,B) 4 CONTINUE OPEN(3,FILE='TESTC') WRITE (3,3000) (B(I),I=1,N) CLOSE(3) STOP 1000 FORMAT(5X/(11F6.0,F8.0)) 2000 FORMAT(5X/(22F6.0,F8.0)) 3000 FORMAT(3E25.15) END SUBROUTINE CRDCMP(TOPBLK,ARRAY,BOTBLK,IPIVOT,B) IMPLICIT REAL*8(A-H,O-Z) PARAMETER(LM=10,LN=1,KK=10, , LM1=LM+1,LN1=LN+1,KK1=KK+1, , LC=LM+LN,LC1=LC+1,LA=2*LC,N=LC*KK1) DIMENSION TOPBLK(LM,LC),ARRAY(LC,LA,KK),BOTBLK(LN,LC) DIMENSION IPIVOT(N),B(N) LB=LC+LM LB1=LB+1 DO 80 I=1,LM I1 = I + 1 IPVT = I COLMAX = DABS(TOPBLK(I,I)) DO 30 J=I1,LC TEMPIV = DABS(TOPBLK(I,J)) IF (TEMPIV.LE.COLMAX) GO TO 30 IPVT = J COLMAX = TEMPIV 30 CONTINUE IPIVOT(I) = IPVT IF (IPVT.EQ.I) GO TO 60 DO 40 L=1,LM SWAP = TOPBLK(L,IPVT) TOPBLK(L,IPVT) = TOPBLK(L,I) TOPBLK(L,I) = SWAP 40 CONTINUE DO 50 L=1,LC SWAP = ARRAY(L,IPVT,1) ARRAY(L,IPVT,1) = ARRAY(L,I,1) ARRAY(L,I,1) = SWAP 50 CONTINUE 60 CONTINUE CII = TOPBLK(I,I) DO 70 J=I1,LC CIJ = TOPBLK(I,J)/CII TOPBLK(I,J) = CIJ DO 70 L=I1,LM 70 TOPBLK(L,J) = TOPBLK(L,J) - TOPBLK(L,I)*CIJ BI=B(I)/CII B(I)=BI DO 75 L=I1,LM 75 B(L)=B(L)-TOPBLK(L,I)*BI 80 CONTINUE DO 110 L=LM,1,-1 L1=L-1 BL=B(L) DO 90 I=1,L1 CIL=TOPBLK(I,L) DO 85 J=LM1,LC 85 TOPBLK(I,J)=TOPBLK(I,J)-CIL*TOPBLK(L,J) 90 B(I)=B(I)-CIL*BL DO 100 I=1,LC IM=I+LM CIL=ARRAY(I,L,1) DO 95 J=LM1,LC 95 ARRAY(I,J,1) = ARRAY(I,J,1) - CIL*TOPBLK(L,J) 100 B(IM)=B(IM)-CIL*BL 110 CONTINUE INCR = 0 DO 320 K=1,KK K1 = K + 1 INCRM=INCR+LM INCRC=INCR+LC INCRB=INCR+LB DO 180 J=LM1,LC J1 = J + 1 JN = J - LM JN1 = JN + 1 INCRJ=INCR+J IPVT = JN ROWMAX = DABS(ARRAY(JN,J,K)) DO 120 I=JN1,LC TEMPIV = DABS(ARRAY(I,J,K)) IF (TEMPIV.LE.ROWMAX) GO TO 120 IPVT = I ROWMAX = TEMPIV 120 CONTINUE INCRP=INCRM+IPVT IPIVOT(INCRJ)=INCRP IF (IPVT.EQ.JN) GO TO 140 DO 130 L=J,LA SWAP = ARRAY(IPVT,L,K) ARRAY(IPVT,L,K) = ARRAY(JN,L,K) ARRAY(JN,L,K) = SWAP 130 CONTINUE SWAP = B(INCRP) B(INCRP) = B(INCRJ) B(INCRJ) = SWAP 140 CONTINUE ROWPIV = ARRAY(JN,J,K) DO 150 I=JN1,LC ARRAY(I,J,K) = ARRAY(I,J,K)/ROWPIV 150 CONTINUE DO 160 L=J1,LA ROWMLT = ARRAY(JN,L,K) DO 160 I=JN1,LC 160 ARRAY(I,L,K) = ARRAY(I,L,K) - ROWMLT*ARRAY(I,J,K) ROWMLT = B(INCRJ) DO 170 I=JN1,LC INCRI=INCRM+I 170 B(INCRI) = B(INCRI) - ROWMLT*ARRAY(I,J,K) 180 CONTINUE DO 260 I=LN1,LC IM = I + LM IM1 = IM + 1 IMC = IM - LC INCRN = INCR + IM IN=I-LN I1 = I + 1 IPVT = IM COLMAX = DABS(ARRAY(I,IPVT,K)) DO 190 J=IM1,LA TEMPIV = DABS(ARRAY(I,J,K)) IF (TEMPIV.LE.COLMAX) GO TO 190 IPVT = J COLMAX = TEMPIV 190 CONTINUE IPIVOT(INCRN) = INCR + IPVT IF (IPVT.EQ.IM) GO TO 240 DO 200 L=LN1,LC SWAP = ARRAY(L,IPVT,K) ARRAY(L,IPVT,K) = ARRAY(L,IM,K) ARRAY(L,IM,K) = SWAP 200 CONTINUE IPC = IPVT - LC IF (K.EQ.KK) GO TO 220 DO 210 L=1,LC SWAP = ARRAY(L,IPC,K1) ARRAY(L,IPC,K1) = ARRAY(L,IMC,K1) ARRAY(L,IMC,K1) = SWAP 210 CONTINUE GO TO 240 220 CONTINUE DO 230 L=1,LN SWAP = BOTBLK(L,IPC) BOTBLK(L,IPC) = BOTBLK(L,IMC) BOTBLK(L,IMC) = SWAP 230 CONTINUE 240 CONTINUE CII = ARRAY(I,IM,K) DO 250 J=IM1,LA CIJ = ARRAY(I,J,K)/CII ARRAY(I,J,K) = CIJ DO 250 L=I1,LC 250 ARRAY(L,J,K) = ARRAY(L,J,K) - ARRAY(L,IM,K)*CIJ BI=B(INCRN)/CII B(INCRN)=BI DO 255 L=I1,LC LINCRM=L+INCRM 255 B(LINCRM)=B(LINCRM)-ARRAY(L,IM,K)*BI 260 CONTINUE DO 300 L=LM,1,-1 L1=L-1 L1N=L1+LN BL=B(L+INCRC) DO 270 I=LN1,L1N IINCRM=I+INCRM CIL=ARRAY(I,L+LC,K) DO 265 J=LB1,LA 265 ARRAY(I,J,K)=ARRAY(I,J,K)-CIL*ARRAY(L+LN,J,K) 270 B(IINCRM)=B(IINCRM)-CIL*BL IF(K.LT.KK) THEN DO 280 I=1,LC IINCRB=I+INCRB CIL=ARRAY(I,L,K1) DO 275 J=LM1,LC 275 ARRAY(I,J,K1) = ARRAY(I,J,K1) - CIL*ARRAY(L+LN,J+LC,K) 280 B(IINCRB)=B(IINCRB)-CIL*BL ELSE DO 290 I=1,LN IINCRB=I+INCRB CIL=BOTBLK(I,L) DO 285 J=LM1,LC 285 BOTBLK(I,J) = BOTBLK(I,J) - CIL*ARRAY(L+LN,J+LC,K) 290 B(IINCRB)=B(IINCRB)-CIL*BL ENDIF 300 CONTINUE INCR = INCRC 320 CONTINUE IF (LN.EQ.1) GO TO 400 INCRM=INCR+LM LC2=LC-1 DO 390 J=LM1,LC2 J1 = J + 1 JN2 = J - LM JN21 = JN2 + 1 INCRJ=INCR+J IPVT = JN2 ROWMAX = DABS(BOTBLK(JN2,J)) DO 330 I=JN21,LN TEMPIV = DABS(BOTBLK(I,J)) IF (TEMPIV.LE.ROWMAX) GO TO 330 IPVT = I ROWMAX = TEMPIV 330 CONTINUE INCRP=INCRM+IPVT IPIVOT(INCRJ)=INCRP IF (IPVT.EQ.JN2) GO TO 350 DO 340 L=J,LC SWAP = BOTBLK(IPVT,L) BOTBLK(IPVT,L) = BOTBLK(JN2,L) BOTBLK(JN2,L) = SWAP 340 CONTINUE SWAP = B(INCRP) B(INCRP) = B(INCRJ) B(INCRJ) = SWAP 350 CONTINUE ROWPIV = BOTBLK(JN2,J) DO 360 I=JN21,LN 360 BOTBLK(I,J) = BOTBLK(I,J)/ROWPIV DO 370 L=J1,LC ROWMLT = BOTBLK(JN2,L) DO 370 I=JN21,LN 370 BOTBLK(I,L) = BOTBLK(I,L) - ROWMLT*BOTBLK(I,J) ROWMLT = B(INCRJ) DO 380 I=JN21,LN INCRI=INCRM+I 380 B(INCRI) = B(INCRI) - ROWMLT*BOTBLK(I,J) 390 CONTINUE 400 CONTINUE RETURN END SUBROUTINE CRSLVE(TOPBLK,ARRAY,BOTBLK,IPIVOT,B) IMPLICIT REAL*8(A-H,O-Z) PARAMETER(LM=10,LN=1,KK=10, , LM1=LM+1,LN1=LN+1,KK1=KK+1, , LC=LM+LN,LC1=LC+1,LA=2*LC,N=LC*KK1) DIMENSION TOPBLK(LM,LC),ARRAY(LC,LA,KK),BOTBLK(LN,LC) DIMENSION IPIVOT(N),B(N) LB=LC+LM LB1=LB+1 INCR=KK*LC INCRM=INCR+LM DO 210 L=1,LN J = LC1 - L INCRJ = INCR + J LN1L = LN1 - L B(INCRJ) = B(INCRJ)/BOTBLK(LN1L,J) IF (L.EQ.LN) GO TO 200 BINCRJ = B(INCRJ) LNL = LN - L DO 190 I=1,LNL INCRI = INCRM + I B(INCRI) = B(INCRI) - BOTBLK(I,J)*BINCRJ 190 CONTINUE
200 CONTINUE 210 CONTINUE DO 300 LK=1,KK K = KK1 - LK INCR = INCR - LC DO 225 L1=LN1,LC I = LC + LN1 - L1 IM = I + LM INCRN = INCR + IM DOTPRD = B(INCRN) DO 220 J=LB1,LA INCRJ = INCR + J DOTPRD = DOTPRD - ARRAY(I,J,K)*B(INCRJ) 220 CONTINUE B(INCRN) = DOTPRD 225 CONTINUE DO 240 L1=LN1,LC I = LC + LN1 - L1 IM = I + LM INCRN = INCR + IM IPVTN = IPIVOT(INCRN) IF (INCRN.EQ.IPVTN) GO TO 230 SWAP = B(INCRN) B(INCRN) = B(IPVTN) B(IPVTN) = SWAP 230 CONTINUE 240 CONTINUE INCRM = INCR + LM DO 260 J=LC1,LA INCRJ = INCR + J BINCRJ = B(INCRJ) DO 250 I=1,LN INCRI = INCRM + I B(INCRI) = B(INCRI) - ARRAY(I,J,K)*BINCRJ 250 CONTINUE 260 CONTINUE DO 290 L=1,LN J = LC1 - L INCRJ = INCR + J LN1L = LN1 - L B(INCRJ) = B(INCRJ)/ARRAY(LN1L,J,K) IF (L.EQ.LN) GO TO 280 BINCRJ = B(INCRJ) LNL = LN - L DO 270 I=1,LNL INCRI = INCRM + I B(INCRI) = B(INCRI) - ARRAY(I,J,K)*BINCRJ 270 CONTINUE 280 CONTINUE 290 CONTINUE 300 CONTINUE DO 315 L=1,LM I = LM1 - L DOTPRD = B(I) DO 310 J=LM1,LC DOTPRD = DOTPRD - TOPBLK(I,J)*B(J) 310 CONTINUE B(I) = DOTPRD 315 CONTINUE DO 330 L=1,LM I = LM1 - L IPVTI = IPIVOT(I) IF (I.EQ.IPVTI) GO TO 320 SWAP = B(I) B(I) = B(IPVTI) B(IPVTI) = SWAP 320 CONTINUE 330 CONTINUE RETURN END
BCBR.FOR Fortran program
CCCCC C BCBR CCCCC IMPLICIT REAL*8(A-H,O-Z) PARAMETER(LM=10,LN=1,KK=10,MM=1000000, , LM1=LM+1,LN1=LN+1,KK1=KK+1, , LC=LM+LN,LC1=LC+1,LA=2*LC,LA1=LA+1,N=LC*KK1) DIMENSION TOPBLK(LM,LC),ARRAY(LC,LA,KK),BOTBLK(LN,LC) DIMENSION IPIVOT(N),B(N) DIMENSION A1(LM,LC1),A2(LN,LC1),A3(LC,LA1) OPEN(1,FILE='TESTA') READ(1,1000) ((A1(I,J),J=1,LC1),I=1,LM) READ(1,1000) ((A2(I,J),J=1,LC1),I=1,LN) CLOSE(1) OPEN(2,FILE='TESTB') READ(2,2000) ((A3(I,J),J=1,LA1),I=1,LC) CLOSE(2) DO 4 M=1,MM DO 1 I=1,LM B(I)=A1(I,LC1) DO 1 J=1,LC 1 TOPBLK(I,J)=A1(I,J) DO 2 K=1,KK DO 2 I=1,LC B(I+K*LC-LN)=A3(I,LA1) DO 2 J=1,LA 2 ARRAY(I,J,K)=A3(I,J) DO 3 I=1,LN B(I+N-LN)=A2(I,LC1) DO 3 J=1,LC 3 BOTBLK(I,J)=A2(I,J) CALL CRDCMP(TOPBLK,ARRAY,BOTBLK,IPIVOT,B) CALL CRSLVE(TOPBLK,ARRAY,BOTBLK,IPIVOT,B) 4 CONTINUE OPEN(3,FILE='TESTC') WRITE (3,3000) (B(I),I=1,N) CLOSE(3) STOP 1000 FORMAT(5X/(11F6.0,F8.0)) 2000 FORMAT(5X/(22F6.0,F8.0)) 3000 FORMAT(3E25.15) END SUBROUTINE CRDCMP(TOPBLK,ARRAY,BOTBLK,IPIVOT,B) IMPLICIT REAL*8(A-H,O-Z) PARAMETER(LM=10,LN=1,KK=10, , LM1=LM+1,LN1=LN+1,KK1=KK+1, , LC=LM+LN,LC1=LC+1,LA=2*LC,N=LC*KK1) DIMENSION TOPBLK(LM,LC),ARRAY(LC,LA,KK),BOTBLK(LN,LC) DIMENSION IPIVOT(N),B(N) LB=LC+LM LB1=LB+1 DO 80 I=1,LM I1 = I + 1 IPVT = I COLMAX = DABS(TOPBLK(I,I)) DO 30 J=I1,LC TEMPIV = DABS(TOPBLK(I,J)) IF (TEMPIV.LE.COLMAX) GO TO 30 IPVT = J COLMAX = TEMPIV 30 CONTINUE IPIVOT(I) = IPVT IF (IPVT.EQ.I) GO TO 60 DO 40 L=1,LM SWAP = TOPBLK(L,IPVT) TOPBLK(L,IPVT) = TOPBLK(L,I) TOPBLK(L,I) = SWAP 40 CONTINUE DO 50 L=1,LC SWAP = ARRAY(L,IPVT,1) ARRAY(L,IPVT,1) = ARRAY(L,I,1) ARRAY(L,I,1) = SWAP 50 CONTINUE 60 CONTINUE CII = TOPBLK(I,I) DO 70 J=I1,LC CIJ = TOPBLK(I,J)/CII TOPBLK(I,J) = CIJ DO 70 L=I1,LM 70 TOPBLK(L,J) = TOPBLK(L,J) - TOPBLK(L,I)*CIJ BI=B(I)/CII B(I)=BI DO 75 L=I1,LM 75 B(L)=B(L)-TOPBLK(L,I)*BI 80 CONTINUE DO 110 L=LM,1,-1 L1=L-1 BL=B(L) DO 90 I=1,L1 CIL=TOPBLK(I,L) DO 85 J=LM1,LC 85 TOPBLK(I,J)=TOPBLK(I,J)-CIL*TOPBLK(L,J) 90 B(I)=B(I)-CIL*BL DO 100 I=1,LC IM=I+LM CIL=ARRAY(I,L,1) DO 95 J=LM1,LC 95 ARRAY(I,J,1) = ARRAY(I,J,1) - CIL*TOPBLK(L,J) 100 B(IM)=B(IM)-CIL*BL 110 CONTINUE INCR = 0 DO 320 K=1,KK K1 = K + 1 INCRM=INCR+LM INCRC=INCR+LC INCRB=INCR+LB DO 170 J=LM1,LC J1 = J + 1 JN = J - LM JN1 = JN + 1 INCRJ=INCR+J IPVT = JN ROWMAX = DABS(ARRAY(JN,J,K)) DO 120 I=JN1,LC TEMPIV = DABS(ARRAY(I,J,K)) IF (TEMPIV.LE.ROWMAX) GO TO 120 IPVT = I ROWMAX = TEMPIV 120 CONTINUE INCRP=INCRM+IPVT IPIVOT(INCRJ)=INCRP IF (IPVT.EQ.JN) GO TO 140 DO 130 L=LM1,LA SWAP = ARRAY(IPVT,L,K) ARRAY(IPVT,L,K) = ARRAY(JN,L,K) ARRAY(JN,L,K) = SWAP 130 CONTINUE SWAP = B(INCRP) B(INCRP) = B(INCRJ) B(INCRJ) = SWAP 140 CONTINUE ROWPIV = ARRAY(JN,J,K) DO 150 I=JN1,LC ARRAY(I,J,K) = ARRAY(I,J,K)/ROWPIV 150 CONTINUE DO 160 L=J1,LC ROWMLT = ARRAY(JN,L,K) DO 160 I=JN1,LC 160 ARRAY(I,L,K) = ARRAY(I,L,K) - ROWMLT*ARRAY(I,J,K) 170 CONTINUE L2=LC DO 180 L=LN,2,-1 L1=L2-1 DO 175 J=L1,1,-1 ANJ=ARRAY(L,J,K) DO 175 I=LN1,LC 175 ARRAY(I,J,K)=ARRAY(I,J,K)-ARRAY(I,L2,K)*ANJ 180 L2=L1 DO 182 L=1,LN LLM=L+LM BL=B(L+INCRM) DO 182 I=LN1,LC AIN=ARRAY(I,LLM,K) DO 181 J=LC1,LA 181 ARRAY(I,J,K)=ARRAY(I,J,K)-AIN*ARRAY(L,J,K) IINCRM=I+INCRM 182 B(IINCRM)=B(IINCRM)-AIN*BL DO 260 I=LN1,LC IM = I + LM IM1 = IM + 1 IMC = IM - LC INCRN = INCR + IM IN=I-LN I1 = I + 1 IPVT = IM COLMAX = DABS(ARRAY(I,IPVT,K)) DO 190 J=IM1,LA TEMPIV = DABS(ARRAY(I,J,K)) IF (TEMPIV.LE.COLMAX) GO TO 190 IPVT = J COLMAX = TEMPIV 190 CONTINUE IPIVOT(INCRN) = INCR + IPVT IF (IPVT.EQ.IM) GO TO 240 DO 200 L=LN1,LC SWAP = ARRAY(L,IPVT,K) ARRAY(L,IPVT,K) = ARRAY(L,IM,K) ARRAY(L,IM,K) = SWAP 200 CONTINUE IPC = IPVT - LC IF (K.EQ.KK) GO TO 220 DO 210 L=1,LC SWAP = ARRAY(L,IPC,K1) ARRAY(L,IPC,K1) = ARRAY(L,IMC,K1) ARRAY(L,IMC,K1) = SWAP 210 CONTINUE GO TO 240 220 CONTINUE DO 230 L=1,LN SWAP = BOTBLK(L,IPC) BOTBLK(L,IPC) = BOTBLK(L,IMC) BOTBLK(L,IMC) = SWAP 230 CONTINUE 240 CONTINUE CII = ARRAY(I,IM,K) DO 250 J=IM1,LA CIJ = ARRAY(I,J,K)/CII ARRAY(I,J,K) = CIJ DO 250 L=I1,LC 250 ARRAY(L,J,K) = ARRAY(L,J,K) - ARRAY(L,IM,K)*CIJ BI=B(INCRN)/CII B(INCRN)=BI DO 255 L=I1,LC LINCRM=L+INCRM 255 B(LINCRM)=B(LINCRM)-ARRAY(L,IM,K)*BI 260 CONTINUE DO 300 L=LM,1,-1 L1=L-1 L1N=L1+LN BL=B(L+INCRC) DO 270 I=LN1,L1N IINCRM=I+INCRM CIL=ARRAY(I,L+LC,K) DO 265 J=LB1,LA 265 ARRAY(I,J,K)=ARRAY(I,J,K)-CIL*ARRAY(L+LN,J,K) 270 B(IINCRM)=B(IINCRM)-CIL*BL IF(K.LT.KK) THEN DO 280 I=1,LC IINCRB=I+INCRB CIL=ARRAY(I,L,K1) DO 275 J=LM1,LC 275 ARRAY(I,J,K1) = ARRAY(I,J,K1) - CIL*ARRAY(L+LN,J+LC,K) 280 B(IINCRB)=B(IINCRB)-CIL*BL ELSE DO 290 I=1,LN IINCRB=I+INCRB CIL=BOTBLK(I,L) DO 285 J=LM1,LC 285 BOTBLK(I,J) = BOTBLK(I,J) - CIL*ARRAY(L+LN,J+LC,K) 290 B(IINCRB)=B(IINCRB)-CIL*BL ENDIF 300 CONTINUE INCR = INCRC 320 CONTINUE IF (LN.EQ.1) GO TO 400 INCRM=INCR+LM LC2=LC-1 DO 390 J=LM1,LC2 J1 = J + 1 JN = J - LM JN1 = JN + 1 INCRJ=INCR+J IPVT = JN ROWMAX = DABS(BOTBLK(JN,J)) DO 330 I=JN1,LN TEMPIV = DABS(BOTBLK(I,J)) IF (TEMPIV.LE.ROWMAX) GO TO 330 IPVT = I ROWMAX = TEMPIV 330 CONTINUE INCRP=INCRM+IPVT IPIVOT(INCRJ)=INCRP IF (IPVT.EQ.JN) GO TO 350 DO 340 L=LM1,LC SWAP = BOTBLK(IPVT,L) BOTBLK(IPVT,L) = BOTBLK(JN,L) BOTBLK(JN,L) = SWAP 340 CONTINUE SWAP = B(INCRP) B(INCRP) = B(INCRJ) B(INCRJ) = SWAP 350 CONTINUE ROWPIV = BOTBLK(JN,J) DO 360 I=JN1,LN 360 BOTBLK(I,J) = BOTBLK(I,J)/ROWPIV DO 370 L=J1,LC ROWMLT = BOTBLK(JN,L) DO 370 I=JN1,LN 370 BOTBLK(I,L) = BOTBLK(I,L) - ROWMLT*BOTBLK(I,J) ROWMLT = B(INCRJ) DO 380 I=JN1,LN INCRI=INCRM+I 380 B(INCRI) = B(INCRI) - ROWMLT*BOTBLK(I,J) 390 CONTINUE 400 CONTINUE RETURN END SUBROUTINE CRSLVE(TOPBLK,ARRAY,BOTBLK,IPIVOT,B) IMPLICIT REAL*8(A-H,O-Z) PARAMETER(LM=10,LN=1,KK=10, , LM1=LM+1,LN1=LN+1,KK1=KK+1, , LC=LM+LN,LC1=LC+1,LA=2*LC,N=LC*KK1) DIMENSION TOPBLK(LM,LC),ARRAY(LC,LA,KK),BOTBLK(LN,LC) DIMENSION IPIVOT(N),B(N) LB=LC+LM LB1=LB+1 INCR=KK*LC INCRM=INCR+LM DO 210 L=1,LN J = LC1 - L INCRJ = INCR + J LN1L = LN1 - L B(INCRJ) = B(INCRJ)/BOTBLK(LN1L,J) IF (L.EQ.LN) GO TO 200 BINCRJ = B(INCRJ) LNL = LN - L DO 190 I=1,LNL INCRI = INCRM + I B(INCRI) = B(INCRI) - BOTBLK(I,J)*BINCRJ 190 CONTINUE 200 CONTINUE
210 CONTINUE DO 300 LK=1,KK K = KK1 - LK INCR = INCR - LC DO 225 L1=LN1,LC I = LC + LN1 - L1 IM = I + LM INCRN = INCR + IM DOTPRD = B(INCRN) DO 220 J=LB1,LA INCRJ = INCR + J DOTPRD = DOTPRD - ARRAY(I,J,K)*B(INCRJ) 220 CONTINUE B(INCRN) = DOTPRD 225 CONTINUE DO 240 L1=LN1,LC I = LC + LN1 - L1 IM = I + LM INCRN = INCR + IM IPVTN = IPIVOT(INCRN) IF (INCRN.EQ.IPVTN) GO TO 230 SWAP = B(INCRN) B(INCRN) = B(IPVTN) B(IPVTN) = SWAP 230 CONTINUE 240 CONTINUE INCRM = INCR + LM DO 260 J=LC1,LA INCRJ = INCR + J BINCRJ = B(INCRJ) DO 250 I=1,LN INCRI = INCRM + I B(INCRI) = B(INCRI) - ARRAY(I,J,K)*BINCRJ 250 CONTINUE 260 CONTINUE DO 295 L=2,LN L1=L-1 L1LM=L1+LM BINCRL=B(INCRM+L1) DO 295 I=L,LN INCRI=INCRM+I 295 B(INCRI)=B(INCRI)-BINCRL*ARRAY(I,L1LM,K) DO 290 L=1,LN J = LC1 - L INCRJ = INCR + J LN1L = LN1 - L B(INCRJ) = B(INCRJ)/ARRAY(LN1L,J,K) IF (L.EQ.LN) GO TO 280 BINCRJ = B(INCRJ) LNL = LN - L DO 270 I=1,LNL INCRI = INCRM + I B(INCRI) = B(INCRI) - ARRAY(I,J,K)*BINCRJ 270 CONTINUE 280 CONTINUE 290 CONTINUE 300 CONTINUE DO 315 L=1,LM I = LM1 - L DOTPRD = B(I) DO 310 J=LM1,LC DOTPRD = DOTPRD - TOPBLK(I,J)*B(J) 310 CONTINUE B(I) = DOTPRD 315 CONTINUE DO 330 L=1,LM I = LM1 - L IPVTI = IPIVOT(I) IF (I.EQ.IPVTI) GO TO 320 SWAP = B(I) B(I) = B(IPVTI) B(IPVTI) = SWAP 320 CONTINUE 330 CONTINUE RETURN END
DBTC.FOR Fortran program
CCCCC C DBTC CCCCC IMPLICIT REAL*8(A-H,O-Z) PARAMETER(LM=10,LN=1,KK=10,MM=1000000, , LM1=LM+1,LN1=LN+1,KK1=KK+1, , LC=LM+LN,LC1=LC+1,LA=2*LC,LA1=LA+1,N=LC*KK1) DIMENSION TOPBLK(LM,LC),ARRAY(LC,LA,KK),BOTBLK(LN,LC) DIMENSION IPIVOT(N),B(N) DIMENSION A1(LM,LC1),A2(LN,LC1),A3(LC,LA1) OPEN(1,FILE='TESTA') READ(1,1000) ((A1(I,J),J=1,LC1),I=1,LM) READ(1,1000) ((A2(I,J),J=1,LC1),I=1,LN) CLOSE(1) OPEN(2,FILE='TESTB') READ(2,2000) ((A3(I,J),J=1,LA1),I=1,LC) CLOSE(2) DO 4 M=1,MM DO 1 I=1,LM B(I)=A1(I,LC1) DO 1 J=1,LC 1 TOPBLK(I,J)=A1(I,J) DO 2 K=1,KK DO 2 I=1,LC B(I+K*LC-LN)=A3(I,LA1) DO 2 J=1,LA 2 ARRAY(I,J,K)=A3(I,J) DO 3 I=1,LN B(I+N-LN)=A2(I,LC1) DO 3 J=1,LC 3 BOTBLK(I,J)=A2(I,J) CALL CRDCMP(TOPBLK,ARRAY,BOTBLK,IPIVOT,B) CALL CRSLVE(TOPBLK,ARRAY,BOTBLK,IPIVOT,B) 4 CONTINUE OPEN(3,FILE='TESTC') WRITE (3,3000) (B(I),I=1,N) CLOSE(3) STOP 1000 FORMAT(5X/(11F6.0,F8.0)) 2000 FORMAT(5X/(22F6.0,F8.0)) 3000 FORMAT(3E25.15) END SUBROUTINE CRDCMP(TOPBLK,ARRAY,BOTBLK,IPIVOT,B) IMPLICIT REAL*8(A-H,O-Z) PARAMETER(LM=10,LN=1,KK=10, , LM1=LM+1,LN1=LN+1,KK1=KK+1, , LC=LM+LN,LC1=LC+1,LA=2*LC,N=LC*KK1) DIMENSION TOPBLK(LM,LC),ARRAY(LC,LA,KK),BOTBLK(LN,LC) DIMENSION IPIVOT(N),B(N) LB=LC+LM LB1=LB+1 DO 80 I=1,LM I1 = I + 1 IPVT = I COLMAX = DABS(TOPBLK(I,I)) DO 30 J=I1,LC TEMPIV = DABS(TOPBLK(I,J)) IF (TEMPIV.LE.COLMAX) GO TO 30 IPVT = J COLMAX = TEMPIV 30 CONTINUE IPIVOT(I) = IPVT IF (IPVT.EQ.I) GO TO 60 DO 40 L=1,LM SWAP = TOPBLK(L,IPVT) TOPBLK(L,IPVT) = TOPBLK(L,I) TOPBLK(L,I) = SWAP 40 CONTINUE DO 50 L=1,LC SWAP = ARRAY(L,IPVT,1) ARRAY(L,IPVT,1) = ARRAY(L,I,1) ARRAY(L,I,1) = SWAP 50 CONTINUE 60 CONTINUE CII = TOPBLK(I,I) DO 70 J=I1,LC CIJ = TOPBLK(I,J)/CII TOPBLK(I,J) = CIJ DO 70 L=I1,LM 70 TOPBLK(L,J) = TOPBLK(L,J) - TOPBLK(L,I)*CIJ BI=B(I)/CII B(I)=BI DO 75 L=I1,LM 75 B(L)=B(L)-TOPBLK(L,I)*BI 80 CONTINUE DO 110 L=LM,1,-1 L1=L-1 BL=B(L) DO 90 I=1,L1 CIL=TOPBLK(I,L) DO 85 J=LM1,LC 85 TOPBLK(I,J)=TOPBLK(I,J)-CIL*TOPBLK(L,J) 90 B(I)=B(I)-CIL*BL DO 100 I=1,LC IM=I+LM CIL=ARRAY(I,L,1) DO 95 J=LM1,LC 95 ARRAY(I,J,1) = ARRAY(I,J,1) - CIL*TOPBLK(L,J) 100 B(IM)=B(IM)-CIL*BL 110 CONTINUE INCR = 0 DO 320 K=1,KK K1 = K + 1 INCRM=INCR+LM INCRC=INCR+LC INCRB=INCR+LB DO 180 J=LM1,LC J1 = J + 1 JN = J - LM JN1 = JN + 1 INCRJ=INCR+J IPVT = JN ROWMAX = DABS(ARRAY(JN,J,K)) DO 120 I=JN1,LC TEMPIV = DABS(ARRAY(I,J,K)) IF (TEMPIV.LE.ROWMAX) GO TO 120 IPVT = I ROWMAX = TEMPIV 120 CONTINUE INCRP=INCRM+IPVT IPIVOT(INCRJ)=INCRP IF (IPVT.EQ.JN) GO TO 140 DO 130 L=J,LA SWAP = ARRAY(IPVT,L,K) ARRAY(IPVT,L,K) = ARRAY(JN,L,K) ARRAY(JN,L,K) = SWAP 130 CONTINUE SWAP = B(INCRP) B(INCRP) = B(INCRJ) B(INCRJ) = SWAP 140 CONTINUE ROWPIV = ARRAY(JN,J,K) DO 150 I=JN1,LC ARRAY(I,J,K) = ARRAY(I,J,K)/ROWPIV 150 CONTINUE DO 160 L=J1,LA ROWMLT = ARRAY(JN,L,K) DO 160 I=JN1,LC 160 ARRAY(I,L,K) = ARRAY(I,L,K) - ROWMLT*ARRAY(I,J,K) ROWMLT = B(INCRJ) DO 170 I=JN1,LC INCRI=INCRM+I 170 B(INCRI) = B(INCRI) - ROWMLT*ARRAY(I,J,K) 180 CONTINUE DO 260 I=LN1,LC IM = I + LM IM1 = IM + 1 IMC = IM - LC INCRN = INCR + IM IN=I-LN I1 = I + 1 IPVT = IM COLMAX = DABS(ARRAY(I,IPVT,K)) DO 190 J=IM1,LA TEMPIV = DABS(ARRAY(I,J,K)) IF (TEMPIV.LE.COLMAX) GO TO 190 IPVT = J COLMAX = TEMPIV 190 CONTINUE IPIVOT(INCRN) = INCR + IPVT IF (IPVT.EQ.IM) GO TO 240 DO 200 L=1,LC CC SWAP = ARRAY(L,IPVT,K) ARRAY(L,IPVT,K) = ARRAY(L,IM,K) ARRAY(L,IM,K) = SWAP 200 CONTINUE IPC = IPVT - LC IF (K.EQ.KK) GO TO 220 DO 210 L=1,LC SWAP = ARRAY(L,IPC,K1) ARRAY(L,IPC,K1) = ARRAY(L,IMC,K1) ARRAY(L,IMC,K1) = SWAP 210 CONTINUE GO TO 240 220 CONTINUE DO 230 L=1,LN SWAP = BOTBLK(L,IPC) BOTBLK(L,IPC) = BOTBLK(L,IMC) BOTBLK(L,IMC) = SWAP 230 CONTINUE 240 CONTINUE CII = ARRAY(I,IM,K) DO 254 J=IM1,LA CIJ = ARRAY(I,J,K)/CII ARRAY(I,J,K) = CIJ DO 250 L=I1,LC 250 ARRAY(L,J,K) = ARRAY(L,J,K) - ARRAY(L,IM,K)*CIJ 254 CONTINUE BI=B(INCRN)/CII B(INCRN)=BI DO 255 L=I1,LC LINCRM=L+INCRM 255 B(LINCRM)=B(LINCRM)-ARRAY(L,IM,K)*BI 260 CONTINUE DO 270 L=LM,1,-1 L1=L-1 L1N=L1+LN BL=B(L+INCRC) DO 270 I=LN1,L1N IINCRM=I+INCRM CIL=ARRAY(I,L+LC,K) DO 265 J=LB1,LA 265 ARRAY(I,J,K)=ARRAY(I,J,K)-CIL*ARRAY(L+LN,J,K) 270 B(IINCRM)=B(IINCRM)-CIL*BL DO 259 I=LN1,LC IM=I+LM IN=I-LN DO 2540 J=LB1,LA CIJ = ARRAY(I,J,K) DO 251 L=1,LN 251 ARRAY(L,J,K) = ARRAY(L,J,K) - ARRAY(L,IM,K)*CIJ JC = J - LC IF (K.LT.KK) THEN DO 252 L=1,LC 252 ARRAY(L,JC,K1) = ARRAY(L,JC,K1) -ARRAY(L,IN,K1)*CIJ ELSE DO 253 L=1,LN 253 BOTBLK(L,JC) = BOTBLK(L,JC) -BOTBLK(L,IN)*CIJ ENDIF 2540 CONTINUE BI=B(I+INCRM) DO 256 L=1,LN LINCRM=L+INCRM 256 B(LINCRM)=B(LINCRM)-ARRAY(L,IM,K)*BI IF (K.LT.KK) THEN DO 257 L=1,LC LINCRB=L+INCRB 257 B(LINCRB)=B(LINCRB)-ARRAY(L,IN,K1)*BI ELSE DO 258 L=1,LN LINCRB=L+INCRB 258 B(LINCRB)=B(LINCRB)-BOTBLK(L,IN)*BI ENDIF 259 CONTINUE DO 310 L=1,LN LL = LC1 - L INCRL = INCR + LL LN1L = LN1 - L LNL = LN - L ALL=ARRAY(LN1L,LL,K) BINCRL = B(INCRL)/ALL B(INCRL)=BINCRL DO 301 J=LB1,LA ALJ=ARRAY(LN1L,J,K)/ALL ARRAY(LN1L,J,K)=ALJ DO 301 I=1,LNL 301 ARRAY(I,J,K)=ARRAY(I,J,K)-ARRAY(I,LL,K)*ALJ DO 302 I=1,LNL INCRI = INCRM + I 302 B(INCRI) = B(INCRI) - ARRAY(I,LL,K)*BINCRL 310 CONTINUE INCR = INCRC 320 CONTINUE IF (LN.EQ.1) GO TO 400 INCRM=INCR+LM LC2=LC-1 DO 390 J=LM1,LC2 J1 = J + 1 JN2 = J - LM JN21 = JN2 + 1 INCRJ=INCR+J IPVT = JN2 ROWMAX = DABS(BOTBLK(JN2,J)) DO 330 I=JN21,LN TEMPIV = DABS(BOTBLK(I,J)) IF (TEMPIV.LE.ROWMAX) GO TO 330 IPVT = I ROWMAX = TEMPIV 330 CONTINUE INCRP=INCRM+IPVT IPIVOT(INCRJ)=INCRP IF (IPVT.EQ.JN2) GO TO 350 DO 340 L=J,LC SWAP = BOTBLK(IPVT,L) BOTBLK(IPVT,L) = BOTBLK(JN2,L) BOTBLK(JN2,L) = SWAP 340 CONTINUE SWAP = B(INCRP) B(INCRP) = B(INCRJ) B(INCRJ) = SWAP 350 CONTINUE ROWPIV = BOTBLK(JN2,J) DO 360 I=JN21,LN 360 BOTBLK(I,J) = BOTBLK(I,J)/ROWPIV DO 370 L=J1,LC ROWMLT = BOTBLK(JN2,L) DO 370 I=JN21,LN 370 BOTBLK(I,L) = BOTBLK(I,L) - ROWMLT*BOTBLK(I,J) ROWMLT = B(INCRJ) DO 380 I=JN21,LN INCRI=INCRM+I 380 B(INCRI) = B(INCRI) - ROWMLT*BOTBLK(I,J) 390 CONTINUE 400 CONTINUE RETURN END SUBROUTINE CRSLVE(TOPBLK,ARRAY,BOTBLK,IPIVOT,B) IMPLICIT REAL*8(A-H,O-Z) PARAMETER(LM=10,LN=1,KK=10, , LM1=LM+1,LN1=LN+1,KK1=KK+1, , LC=LM+LN,LC1=LC+1,LA=2*LC,N=LC*KK1) DIMENSION TOPBLK(LM,LC),ARRAY(LC,LA,KK),BOTBLK(LN,LC) DIMENSION IPIVOT(N),B(N) LB=LC+LM LB1=LB+1 INCR=KK*LC INCRM=INCR+LM DO 210 L=1,LN J = LC1 - L INCRJ = INCR + J LN1L = LN1 - L B(INCRJ) = B(INCRJ)/BOTBLK(LN1L,J) IF (L.EQ.LN) GO TO 200 BINCRJ = B(INCRJ) LNL = LN - L DO 190 I=1,LNL INCRI = INCRM + I B(INCRI) = B(INCRI) - BOTBLK(I,J)*BINCRJ 190 CONTINUE 200 CONTINUE 210 CONTINUE DO 300 LK=1,KK K = KK1 - LK INCR = INCR - LC DO 225 L=1,LC IC = LC1-L IB = LB1 - L INCRI = INCR + IB DOTPRD = B(INCRI) DO 220 J=LB1,LA INCRJ = INCR + J DOTPRD = DOTPRD - ARRAY(IC,J,K)*B(INCRJ) 220 CONTINUE B(INCRI) = DOTPRD 225 CONTINUE DO 240 L=1,LM IB = LB1-L INCRI = INCR + IB IPVTI = IPIVOT(INCRI) IF (INCRI.EQ.IPVTI) GO TO 230 SWAP = B(INCRI) B(INCRI) = B(IPVTI) B(IPVTI) = SWAP 230 CONTINUE 240 CONTINUE 300 CONTINUE DO 315 L=1,LM I = LM1 - L DOTPRD = B(I) DO 310 J=LM1,LC DOTPRD = DOTPRD - TOPBLK(I,J)*B(J) 310 CONTINUE B(I) = DOTPRD 315 CONTINUE DO 330 L=1,LM I = LM1 - L IPVTI = IPIVOT(I) IF (I.EQ.IPVTI) GO TO 320 SWAP = B(I) B(I) = B(IPVTI) B(IPVTI) = SWAP 320 CONTINUE 330 CONTINUE RETURN END TESTA input data
TOP BLOCK .10 -.22 -.24 -.42 -.37 -.77 -.99 -.96 -.89 .85 -.28 -6.412 -.63 .09 -.10 -.07 .51 -.02 .01 -.52 .07 .48 -.54 -0.968 .32 -.29 .02 -.81 .29 .00 -.05 -.91 .00 .00 .69 -0.869 -.25 -.09 -.91 -.17 -.46 -.92 -.14 .98 -.34 .70 -.53 -2.586 .76 -.90 -.64 -.08 .95 .15 .15 -.46 -.48 .93 -.39 0.034 -.06 -.72 -.91 -.14 .36 -.69 -.40 -.93 -.61 -.97 -.12 -8.059 -.21 .54 -.53 .97 .91 .58 -.32 .27 .33 .72 -.20 4.662 -.57 .04 .26 -.04 .69 -.65 -.57 .83 -.42 -.56 -.18 -1.956 .89 -.62 -.07 -.63 .28 -.54 -.29 .52 .67 .00 -.68 -0.847 .10 -.01 -.25 -.22 .06 .81 .11 .56 .05 .63 -.43 2.357 BOTTOM BLOCK .88 .48 .52 -.87 .71 .51 .52 -.33 -.46 -.33 .85 3.176