Mercurial > octave-nkf

comparison libcruft/lapack/zgttrf.f @ 7034:68db500cb558

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-10-16 18:54:19 by jwe]
author jwe
date Tue, 16 Oct 2007 18:54:23 +0000
parents 57077d0ddc8e
children
comparison
equal deleted inserted replaced
7033:f0142f2afdc6 7034:68db500cb558
1 SUBROUTINE ZGTTRF( N, DL, D, DU, DU2, IPIV, INFO ) 1 SUBROUTINE ZGTTRF( N, DL, D, DU, DU2, IPIV, INFO )
2 * 2 *
3 * -- LAPACK routine (version 2.0) -- 3 * -- LAPACK routine (version 3.1) --
4 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 * Courant Institute, Argonne National Lab, and Rice University 5 * November 2006
6 * September 30, 1994
7 * 6 *
8 * .. Scalar Arguments .. 7 * .. Scalar Arguments ..
9 INTEGER INFO, N 8 INTEGER INFO, N
10 * .. 9 * ..
11 * .. Array Arguments .. 10 * .. Array Arguments ..
27 * 26 *
28 * Arguments 27 * Arguments
29 * ========= 28 * =========
30 * 29 *
31 * N (input) INTEGER 30 * N (input) INTEGER
32 * The order of the matrix A. N >= 0. 31 * The order of the matrix A.
33 * 32 *
34 * DL (input/output) COMPLEX*16 array, dimension (N-1) 33 * DL (input/output) COMPLEX*16 array, dimension (N-1)
35 * On entry, DL must contain the (n-1) subdiagonal elements of 34 * On entry, DL must contain the (n-1) sub-diagonal elements of
36 * A. 35 * A.
36 *
37 * On exit, DL is overwritten by the (n-1) multipliers that 37 * On exit, DL is overwritten by the (n-1) multipliers that
38 * define the matrix L from the LU factorization of A. 38 * define the matrix L from the LU factorization of A.
39 * 39 *
40 * D (input/output) COMPLEX*16 array, dimension (N) 40 * D (input/output) COMPLEX*16 array, dimension (N)
41 * On entry, D must contain the diagonal elements of A. 41 * On entry, D must contain the diagonal elements of A.
42 *
42 * On exit, D is overwritten by the n diagonal elements of the 43 * On exit, D is overwritten by the n diagonal elements of the
43 * upper triangular matrix U from the LU factorization of A. 44 * upper triangular matrix U from the LU factorization of A.
44 * 45 *
45 * DU (input/output) COMPLEX*16 array, dimension (N-1) 46 * DU (input/output) COMPLEX*16 array, dimension (N-1)
46 * On entry, DU must contain the (n-1) superdiagonal elements 47 * On entry, DU must contain the (n-1) super-diagonal elements
47 * of A. 48 * of A.
49 *
48 * On exit, DU is overwritten by the (n-1) elements of the first 50 * On exit, DU is overwritten by the (n-1) elements of the first
49 * superdiagonal of U. 51 * super-diagonal of U.
50 * 52 *
51 * DU2 (output) COMPLEX*16 array, dimension (N-2) 53 * DU2 (output) COMPLEX*16 array, dimension (N-2)
52 * On exit, DU2 is overwritten by the (n-2) elements of the 54 * On exit, DU2 is overwritten by the (n-2) elements of the
53 * second superdiagonal of U. 55 * second super-diagonal of U.
54 * 56 *
55 * IPIV (output) INTEGER array, dimension (N) 57 * IPIV (output) INTEGER array, dimension (N)
56 * The pivot indices; for 1 <= i <= n, row i of the matrix was 58 * The pivot indices; for 1 <= i <= n, row i of the matrix was
57 * interchanged with row IPIV(i). IPIV(i) will always be either 59 * interchanged with row IPIV(i). IPIV(i) will always be either
58 * i or i+1; IPIV(i) = i indicates a row interchange was not 60 * i or i+1; IPIV(i) = i indicates a row interchange was not
59 * required. 61 * required.
60 * 62 *
61 * INFO (output) INTEGER 63 * INFO (output) INTEGER
62 * = 0: successful exit 64 * = 0: successful exit
63 * < 0: if INFO = -i, the i-th argument had an illegal value 65 * < 0: if INFO = -k, the k-th argument had an illegal value
64 * > 0: if INFO = i, U(i,i) is exactly zero. The factorization 66 * > 0: if INFO = k, U(k,k) is exactly zero. The factorization
65 * has been completed, but the factor U is exactly 67 * has been completed, but the factor U is exactly
66 * singular, and division by zero will occur if it is used 68 * singular, and division by zero will occur if it is used
67 * to solve a system of equations. 69 * to solve a system of equations.
68 * 70 *
69 * ===================================================================== 71 * =====================================================================
70 * 72 *
73 * .. Parameters ..
74 DOUBLE PRECISION ZERO
75 PARAMETER ( ZERO = 0.0D+0 )
76 * ..
71 * .. Local Scalars .. 77 * .. Local Scalars ..
72 INTEGER I 78 INTEGER I
73 COMPLEX*16 FACT, TEMP, ZDUM 79 COMPLEX*16 FACT, TEMP, ZDUM
74 * .. 80 * ..
75 * .. Intrinsic Functions ..
76 INTRINSIC ABS, DBLE, DIMAG
77 * ..
78 * .. External Subroutines .. 81 * .. External Subroutines ..
79 EXTERNAL XERBLA 82 EXTERNAL XERBLA
80 * .. 83 * ..
81 * .. Parameters .. 84 * .. Intrinsic Functions ..
82 COMPLEX*16 CZERO 85 INTRINSIC ABS, DBLE, DIMAG
83 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) )
84 * .. 86 * ..
85 * .. Statement Functions .. 87 * .. Statement Functions ..
86 DOUBLE PRECISION CABS1 88 DOUBLE PRECISION CABS1
87 * .. 89 * ..
88 * .. Statement Function definitions .. 90 * .. Statement Function definitions ..
100 * Quick return if possible 102 * Quick return if possible
101 * 103 *
102 IF( N.EQ.0 ) 104 IF( N.EQ.0 )
103 $ RETURN 105 $ RETURN
104 * 106 *
105 * Initialize IPIV(i) = i 107 * Initialize IPIV(i) = i and DU2(i) = 0
106 * 108 *
107 DO 10 I = 1, N 109 DO 10 I = 1, N
108 IPIV( I ) = I 110 IPIV( I ) = I
109 10 CONTINUE 111 10 CONTINUE
112 DO 20 I = 1, N - 2
113 DU2( I ) = ZERO
114 20 CONTINUE
110 * 115 *
111 DO 20 I = 1, N - 1 116 DO 30 I = 1, N - 2
112 IF( DL( I ).EQ.CZERO ) THEN 117 IF( CABS1( D( I ) ).GE.CABS1( DL( I ) ) ) THEN
113 *
114 * Subdiagonal is zero, no elimination is required.
115 *
116 IF( D( I ).EQ.CZERO .AND. INFO.EQ.0 )
117 $ INFO = I
118 IF( I.LT.N-1 )
119 $ DU2( I ) = CZERO
120 ELSE IF( CABS1( D( I ) ).GE.CABS1( DL( I ) ) ) THEN
121 * 118 *
122 * No row interchange required, eliminate DL(I) 119 * No row interchange required, eliminate DL(I)
123 * 120 *
124 FACT = DL( I ) / D( I ) 121 IF( CABS1( D( I ) ).NE.ZERO ) THEN
125 DL( I ) = FACT 122 FACT = DL( I ) / D( I )
126 D( I+1 ) = D( I+1 ) - FACT*DU( I ) 123 DL( I ) = FACT
127 IF( I.LT.N-1 ) 124 D( I+1 ) = D( I+1 ) - FACT*DU( I )
128 $ DU2( I ) = CZERO 125 END IF
129 ELSE 126 ELSE
130 * 127 *
131 * Interchange rows I and I+1, eliminate DL(I) 128 * Interchange rows I and I+1, eliminate DL(I)
132 * 129 *
133 FACT = D( I ) / DL( I ) 130 FACT = D( I ) / DL( I )
134 D( I ) = DL( I ) 131 D( I ) = DL( I )
135 DL( I ) = FACT 132 DL( I ) = FACT
136 TEMP = DU( I ) 133 TEMP = DU( I )
137 DU( I ) = D( I+1 ) 134 DU( I ) = D( I+1 )
138 D( I+1 ) = TEMP - FACT*D( I+1 ) 135 D( I+1 ) = TEMP - FACT*D( I+1 )
139 IF( I.LT.N-1 ) THEN 136 DU2( I ) = DU( I+1 )
140 DU2( I ) = DU( I+1 ) 137 DU( I+1 ) = -FACT*DU( I+1 )
141 DU( I+1 ) = -FACT*DU( I+1 ) 138 IPIV( I ) = I + 1
139 END IF
140 30 CONTINUE
141 IF( N.GT.1 ) THEN
142 I = N - 1
143 IF( CABS1( D( I ) ).GE.CABS1( DL( I ) ) ) THEN
144 IF( CABS1( D( I ) ).NE.ZERO ) THEN
145 FACT = DL( I ) / D( I )
146 DL( I ) = FACT
147 D( I+1 ) = D( I+1 ) - FACT*DU( I )
142 END IF 148 END IF
143 IPIV( I ) = IPIV( I ) + 1 149 ELSE
150 FACT = D( I ) / DL( I )
151 D( I ) = DL( I )
152 DL( I ) = FACT
153 TEMP = DU( I )
154 DU( I ) = D( I+1 )
155 D( I+1 ) = TEMP - FACT*D( I+1 )
156 IPIV( I ) = I + 1
144 END IF 157 END IF
145 20 CONTINUE
146 IF( D( N ).EQ.CZERO .AND. INFO.EQ.0 ) THEN
147 INFO = N
148 RETURN
149 END IF 158 END IF
159 *
160 * Check for a zero on the diagonal of U.
161 *
162 DO 40 I = 1, N
163 IF( CABS1( D( I ) ).EQ.ZERO ) THEN
164 INFO = I
165 GO TO 50
166 END IF
167 40 CONTINUE
168 50 CONTINUE
150 * 169 *
151 RETURN 170 RETURN
152 * 171 *
153 * End of ZGTTRF 172 * End of ZGTTRF
154 * 173 *