Routine: DDBTRF()  File: SRC\ddbtrf.f

 
 
# lines: 336
  # code: 336
  # comment: 0
  # blank:0
# Variables:18
# Callers:0
# Callings:0
# Words:164
# Keywords:59
 

 

 
..
     .. Local Scalars ..
     ..
     .. Local Arrays ..
     ..
     .. External Functions ..
     ..
     .. External Subroutines ..
     ..
     .. Intrinsic Functions ..
     ..
     .. Executable Statements ..
     KV is the number of superdiagonals in the factor U
     Test the input parameters.
     Quick return if possible
     Determine the block size for this environment
     The block size must not exceed the limit set by the size of the
     local arrays WORK13 and WORK31.
        Use unblocked code
        Use blocked code
        Zero the superdiagonal elements of the work array WORK13
        Zero the subdiagonal elements of the work array WORK31
        JU is the index of the last column affected by the current
        stage of the factorization
           The active part of the matrix is partitioned
              A11   A12   A13
              A21   A22   A23
              A31   A32   A33
           Here A11, A21 and A31 denote the current block of JB columns
           which is about to be factorized. The number of rows in the
           partitioning are JB, I2, I3 respectively, and the numbers
           of columns are JB, J2, J3. The superdiagonal elements of A13
           and the subdiagonal elements of A31 lie outside the band.
           J2 and J3 are computed after JU has been updated.
           Factorize the current block of JB columns
              Find pivot and test for singularity. KM is the number of
              subdiagonal elements in the current column.
                 Compute multipliers
 
Display dynamic version Find AutoScroll Reload FontSize: - + Hide Comments Hide Blanks Frame FullScreen MailPrint

  
001        SUBROUTINE DDBTRF( M , N , KL , KU , AB , LDAB , INFO )
002  
003  *     Written by Andrew J. Cleary , University of Tennessee.
004  *     August , 1996.
005  *     Modified from DGBTRF :
006  *     -- LAPACK routine(preliminary version) --
007  *     Univ. of Tennessee , Univ. of California Berkeley , NAG Ltd. ,
008  *     Courant Institute , Argonne National Lab , and Rice University
009  *     August 6 , 1991
010  
011  *     .. Scalar Arguments ..
012        INTEGER INFO , KL , KU , LDAB , M , N
013  *     ..
014  *     .. Array Arguments ..
015        DOUBLE PRECISION AB( LDAB , * )
016  *     ..
017  
018  *     Purpose
019  *     === ====
020  
021  *     Ddbtrf computes an LU factorization of a real m - by - n band matrix A
022  *     without using partial pivoting or row interchanges.
023  
024  *     This is the blocked version of the algorithm , calling Level 3 BLAS.
025  
026  *     Arguments
027  *     === ======
028  
029  *     M(input) INTEGER
030  *     The number of rows of the matrix A. M >= 0.
031  
032  *     N(input) INTEGER
033  *     The number of columns of the matrix A. N >= 0.
034  
035  *     KL(input) INTEGER
036  *     The number of subdiagonals within the band of A. KL >= 0.
037  
038  *     KU(input) INTEGER
039  *     The number of superdiagonals within the band of A. KU >= 0.
040  
041  *     AB(input / output) REAL array , dimension(LDAB , N)
042  *     On entry , the matrix A in band storage , in rows KL + 1 to
043  *     2*KL + KU + 1 ; rows 1 to KL of the array need not be set.
044  *     The j - th column of A is stored in the j - th column of the
045  *     array AB as follows :
046  *     AB(kl + ku + 1 + i - j , j) = A(i , j) for max(1 , j - ku) <= i <= min(m , j + kl)
047  
048  *     On exit , details of the factorization : U is stored as an
049  *     upper triangular band matrix with KL + KU superdiagonals in
050  *     rows 1 to KL + KU + 1 , and the multipliers used during the
051  *     factorization are stored in rows KL + KU + 2 to 2*KL + KU + 1.
052  *     See below for further details.
053  
054  *     LDAB(input) INTEGER
055  *     The leading dimension of the array AB. LDAB >= 2*KL + KU + 1.
056  
057  *     INFO(output) INTEGER
058  *     = 0 : successful exit
059  *     < 0 : if INFO = - i , the i - th argument had an illegal value
060  *     > 0 : if INFO = + i , U(i , i) is exactly zero. The factorization
061  *     has been completed , but the factor U is exactly
062  *     singular , and division by zero will occur if it is used
063  *     to solve a system of equations.
064  
065  *     Further Details
066  *     === ============
067  
068  *     The band storage scheme is illustrated by the following example , when
069  *     M = N = 6 , KL = 2 , KU = 1 :
070  
071  *     On entry : On exit :
072  
073  *     * a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56
074  *     a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66
075  *     a21 a32 a43 a54 a65 * m21 m32 m43 m54 m65 *
076  *     a31 a42 a53 a64 * * m31 m42 m53 m64 * *
077  
078  *     Array elements marked * are not used by the routine.
079  
080  *     === ==================================================================
081  
082  *     .. Parameters ..
083        DOUBLE PRECISION ONE , ZERO
084        PARAMETER( ONE = 1.0D + 0 )
085        PARAMETER( ZERO = 0.0D + 0 )
086        INTEGER NBMAX , LDWORK
087        PARAMETER( NBMAX = 64 , LDWORK = NBMAX + 1 )
088        CALL DSCAL( KM , ONE / AB( KV + 1 , JJ ) , AB( KV + 2 , JJ ) ,
089       $1 )
090  
091  *     Update trailing submatrix within the band and within
092  *     the current block. JM is the index of the last column
093  *     which needs to be updated.
094  
095        JM = MIN( JU , J + JB - 1 )
096        IF( JM.GT.JJ ) THEN
097            CALL DGER( KM , JM - JJ , - ONE , AB( KV + 2 , JJ ) , 1 ,
098       $    AB( KV , JJ + 1 ) , LDAB - 1 ,
099       $    AB( KV + 1 , JJ + 1 ) , LDAB - 1 )
100        END IF
101        END IF
102  
103  *     Copy current column of A31 into the work array WORK31
104  
105        NW = MIN( JJ - J + 1 , I3 )
106        IF( NW.GT.0 )
107       $    CALL DCOPY( NW , AB( KV + KL + 1 - JJ + J , JJ ) , 1 ,
108       $    WORK31( 1 , JJ - J + 1 ) , 1 )
109     80 CONTINUE
110        IF( J + JB.LE.N ) THEN
111  
112  *         Apply the row interchanges to the other blocks.
113  
114            J2 = MIN( JU - J + 1 , KV ) - JB
115            J3 = MAX( 0 , JU - J - KV + 1 )
116  
117  *         Update the relevant part of the trailing submatrix
118  
119            IF( J2.GT.0 ) THEN
120  
121  *             Update A12
122  
123                CALL DTRSM( 'Left' , 'Lower' , 'No transpose' , 'Unit' ,
124       $        JB , J2 , ONE , AB( KV + 1 , J ) , LDAB - 1 ,
125       $        AB( KV + 1 - JB , J + JB ) , LDAB - 1 )
126  
127                IF( I2.GT.0 ) THEN
128  
129  *                 Update A22
130  
131                    CALL DGEMM( 'No transpose' , 'No transpose' , I2 , J2 ,
132       $            JB , - ONE , AB( KV + 1 + JB , J ) , LDAB - 1 ,
133       $            AB( KV + 1 - JB , J + JB ) , LDAB - 1 , ONE ,
134       $            AB( KV + 1 , J + JB ) , LDAB - 1 )
135                END IF
136  
137                IF( I3.GT.0 ) THEN
138  
139  *                 Update A32
140  
141                    CALL DGEMM( 'No transpose' , 'No transpose' , I3 , J2 ,
142       $            JB , - ONE , WORK31 , LDWORK ,
143       $            AB( KV + 1 - JB , J + JB ) , LDAB - 1 , ONE ,
144       $            AB( KV + KL + 1 - JB , J + JB ) , LDAB - 1 )
145                END IF
146            END IF
147  
148            IF( J3.GT.0 ) THEN
149  
150  *             Copy the lower triangle of A13 into the work array
151  *             WORK13
152  
153                DO 130 JJ = 1 , J3
154                    DO 120 II = JJ , JB
155                        WORK13( II , JJ ) = AB( II - JJ + 1 , JJ + J + KV - 1 )
156    120             CONTINUE
157    130         CONTINUE
158  
159  *             Update A13 in the work array
160  
161                CALL DTRSM( 'Left' , 'Lower' , 'No transpose' , 'Unit' ,
162       $        JB , J3 , ONE , AB( KV + 1 , J ) , LDAB - 1 ,
163       $        WORK13 , LDWORK )
164  
165                IF( I2.GT.0 ) THEN
166  
167  *                 Update A23
168  
169                    CALL DGEMM( 'No transpose' , 'No transpose' , I2 , J3 ,
170       $            JB , - ONE , AB( KV + 1 + JB , J ) , LDAB - 1 ,
171       $            WORK13 , LDWORK , ONE , AB( 1 + JB , J + KV ) ,
172       $            LDAB - 1 )
173                END IF
174  
175                IF( I3.GT.0 ) THEN
176  
177  *                 Update A33
178  
179                    CALL DGEMM( 'No transpose' , 'No transpose' , I3 , J3 ,
180       $            JB , - ONE , WORK31 , LDWORK , WORK13 ,
181       $            LDWORK , ONE , AB( 1 + KL , J + KV ) , LDAB - 1 )
182                END IF
183  
184  *             Copy the lower triangle of A13 back into place
185  
186                DO 150 JJ = 1 , J3
187                    DO 140 II = JJ , JB
188                        AB( II - JJ + 1 , JJ + J + KV - 1 ) = WORK13( II , JJ )
189    140             CONTINUE
190    150         CONTINUE
191            END IF
192        ELSE
193        END IF
194  
195  *     copy the upper triangle of A31 back into place
196  
197        DO 170 JJ = J + JB - 1 , J , - 1
198  
199  *         Copy the current column of A31 back into place
200  
201            NW = MIN( I3 , JJ - J + 1 )
202            IF( NW.GT.0 )
203       $        CALL DCOPY( NW , WORK31( 1 , JJ - J + 1 ) , 1 ,
204       $        AB( KV + KL + 1 - JJ + J , JJ ) , 1 )
205    170 CONTINUE
206    180 CONTINUE
207        END IF
208  
209        RETURN
210  
211  *     End of DDBTRF
212  
213        END