Routine: ZDBTRF()  File: SRC\zdbtrf.f

 
 
# lines: 339
  # code: 339
  # comment: 0
  # blank:0
# Variables:20
# Callers:0
# Callings:0
# Words:166
# Keywords:60
 

 

 
..
     .. 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 ZDBTRF( M , N , KL , KU , AB , LDAB , INFO )
002  
003  *     Written by Andrew J. Cleary , University of Tennessee.
004  *     August , 1996.
005  *     Modified from ZGBTRF :
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        COMPLEX*16 AB( LDAB , * )
016  *     ..
017  
018  *     Purpose
019  *     === ====
020  
021  *     Zdbtrf 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        COMPLEX*16 CONE , CZERO
087        PARAMETER( CONE =( 1.0D + 0 , 0.0D + 0 ) )
088        PARAMETER( CZERO =( 0.0D + 0 , 0.0D + 0 ) )
089        INTEGER NBMAX , LDWORK
090        PARAMETER( NBMAX = 64 , LDWORK = NBMAX + 1 )
091        CALL ZSCAL( KM , ONE / AB( KV + 1 , JJ ) , AB( KV + 2 , JJ ) ,
092       $1 )
093  
094  *     Update trailing submatrix within the band and within
095  *     the current block. JM is the index of the last column
096  *     which needs to be updated.
097  
098        JM = MIN( JU , J + JB - 1 )
099        IF( JM.GT.JJ ) THEN
100            CALL ZGERU( KM , JM - JJ , - CONE , AB( KV + 2 , JJ ) , 1 ,
101       $    AB( KV , JJ + 1 ) , LDAB - 1 ,
102       $    AB( KV + 1 , JJ + 1 ) , LDAB - 1 )
103        END IF
104        END IF
105  
106  *     Copy current column of A31 into the work array WORK31
107  
108        NW = MIN( JJ - J + 1 , I3 )
109        IF( NW.GT.0 )
110       $    CALL ZCOPY( NW , AB( KV + KL + 1 - JJ + J , JJ ) , 1 ,
111       $    WORK31( 1 , JJ - J + 1 ) , 1 )
112     80 CONTINUE
113        IF( J + JB.LE.N ) THEN
114  
115  *         Apply the row interchanges to the other blocks.
116  
117            J2 = MIN( JU - J + 1 , KV ) - JB
118            J3 = MAX( 0 , JU - J - KV + 1 )
119  
120  *         Update the relevant part of the trailing submatrix
121  
122            IF( J2.GT.0 ) THEN
123  
124  *             Update A12
125  
126                CALL ZTRSM( 'Left' , 'Lower' , 'No transpose' , 'Unit' ,
127       $        JB , J2 , CONE , AB( KV + 1 , J ) , LDAB - 1 ,
128       $        AB( KV + 1 - JB , J + JB ) , LDAB - 1 )
129  
130                IF( I2.GT.0 ) THEN
131  
132  *                 Update A22
133  
134                    CALL ZGEMM( 'No transpose' , 'No transpose' , I2 , J2 ,
135       $            JB , - CONE , AB( KV + 1 + JB , J ) , LDAB - 1 ,
136       $            AB( KV + 1 - JB , J + JB ) , LDAB - 1 , CONE ,
137       $            AB( KV + 1 , J + JB ) , LDAB - 1 )
138                END IF
139  
140                IF( I3.GT.0 ) THEN
141  
142  *                 Update A32
143  
144                    CALL ZGEMM( 'No transpose' , 'No transpose' , I3 , J2 ,
145       $            JB , - CONE , WORK31 , LDWORK ,
146       $            AB( KV + 1 - JB , J + JB ) , LDAB - 1 , CONE ,
147       $            AB( KV + KL + 1 - JB , J + JB ) , LDAB - 1 )
148                END IF
149            END IF
150  
151            IF( J3.GT.0 ) THEN
152  
153  *             Copy the lower triangle of A13 into the work array
154  *             WORK13
155  
156                DO 130 JJ = 1 , J3
157                    DO 120 II = JJ , JB
158                        WORK13( II , JJ ) = AB( II - JJ + 1 , JJ + J + KV - 1 )
159    120             CONTINUE
160    130         CONTINUE
161  
162  *             Update A13 in the work array
163  
164                CALL ZTRSM( 'Left' , 'Lower' , 'No transpose' , 'Unit' ,
165       $        JB , J3 , CONE , AB( KV + 1 , J ) , LDAB - 1 ,
166       $        WORK13 , LDWORK )
167  
168                IF( I2.GT.0 ) THEN
169  
170  *                 Update A23
171  
172                    CALL ZGEMM( 'No transpose' , 'No transpose' , I2 , J3 ,
173       $            JB , - CONE , AB( KV + 1 + JB , J ) , LDAB - 1 ,
174       $            WORK13 , LDWORK , CONE , AB( 1 + JB , J + KV ) ,
175       $            LDAB - 1 )
176                END IF
177  
178                IF( I3.GT.0 ) THEN
179  
180  *                 Update A33
181  
182                    CALL ZGEMM( 'No transpose' , 'No transpose' , I3 , J3 ,
183       $            JB , - CONE , WORK31 , LDWORK , WORK13 ,
184       $            LDWORK , CONE , AB( 1 + KL , J + KV ) , LDAB - 1 )
185                END IF
186  
187  *             Copy the lower triangle of A13 back into place
188  
189                DO 150 JJ = 1 , J3
190                    DO 140 II = JJ , JB
191                        AB( II - JJ + 1 , JJ + J + KV - 1 ) = WORK13( II , JJ )
192    140             CONTINUE
193    150         CONTINUE
194            END IF
195        ELSE
196        END IF
197  
198  *     copy the upper triangle of A31 back into place
199  
200        DO 170 JJ = J + JB - 1 , J , - 1
201  
202  *         Copy the current column of A31 back into place
203  
204            NW = MIN( I3 , JJ - J + 1 )
205            IF( NW.GT.0 )
206       $        CALL ZCOPY( NW , WORK31( 1 , JJ - J + 1 ) , 1 ,
207       $        AB( KV + KL + 1 - JJ + J , JJ ) , 1 )
208    170 CONTINUE
209    180 CONTINUE
210        END IF
211  
212        RETURN
213  
214  *     End of ZDBTRF
215  
216        END