Mercurial > octave-nkf
comparison libcruft/lapack/sgttrf.f @ 7789:82be108cc558
First attempt at single precision tyeps
* * *
corrections to qrupdate single precision routines
* * *
prefer demotion to single over promotion to double
* * *
Add single precision support to log2 function
* * *
Trivial PROJECT file update
* * *
Cache optimized hermitian/transpose methods
* * *
Add tests for tranpose/hermitian and ChangeLog entry for new transpose code
author | David Bateman <dbateman@free.fr> |
---|---|
date | Sun, 27 Apr 2008 22:34:17 +0200 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
7788:45f5faba05a2 | 7789:82be108cc558 |
---|---|
1 SUBROUTINE SGTTRF( N, DL, D, DU, DU2, IPIV, INFO ) | |
2 * | |
3 * -- LAPACK routine (version 3.1) -- | |
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. | |
5 * November 2006 | |
6 * | |
7 * .. Scalar Arguments .. | |
8 INTEGER INFO, N | |
9 * .. | |
10 * .. Array Arguments .. | |
11 INTEGER IPIV( * ) | |
12 REAL D( * ), DL( * ), DU( * ), DU2( * ) | |
13 * .. | |
14 * | |
15 * Purpose | |
16 * ======= | |
17 * | |
18 * SGTTRF computes an LU factorization of a real tridiagonal matrix A | |
19 * using elimination with partial pivoting and row interchanges. | |
20 * | |
21 * The factorization has the form | |
22 * A = L * U | |
23 * where L is a product of permutation and unit lower bidiagonal | |
24 * matrices and U is upper triangular with nonzeros in only the main | |
25 * diagonal and first two superdiagonals. | |
26 * | |
27 * Arguments | |
28 * ========= | |
29 * | |
30 * N (input) INTEGER | |
31 * The order of the matrix A. | |
32 * | |
33 * DL (input/output) REAL array, dimension (N-1) | |
34 * On entry, DL must contain the (n-1) sub-diagonal elements of | |
35 * A. | |
36 * | |
37 * On exit, DL is overwritten by the (n-1) multipliers that | |
38 * define the matrix L from the LU factorization of A. | |
39 * | |
40 * D (input/output) REAL array, dimension (N) | |
41 * On entry, D must contain the diagonal elements of A. | |
42 * | |
43 * On exit, D is overwritten by the n diagonal elements of the | |
44 * upper triangular matrix U from the LU factorization of A. | |
45 * | |
46 * DU (input/output) REAL array, dimension (N-1) | |
47 * On entry, DU must contain the (n-1) super-diagonal elements | |
48 * of A. | |
49 * | |
50 * On exit, DU is overwritten by the (n-1) elements of the first | |
51 * super-diagonal of U. | |
52 * | |
53 * DU2 (output) REAL array, dimension (N-2) | |
54 * On exit, DU2 is overwritten by the (n-2) elements of the | |
55 * second super-diagonal of U. | |
56 * | |
57 * IPIV (output) INTEGER array, dimension (N) | |
58 * The pivot indices; for 1 <= i <= n, row i of the matrix was | |
59 * interchanged with row IPIV(i). IPIV(i) will always be either | |
60 * i or i+1; IPIV(i) = i indicates a row interchange was not | |
61 * required. | |
62 * | |
63 * INFO (output) INTEGER | |
64 * = 0: successful exit | |
65 * < 0: if INFO = -k, the k-th argument had an illegal value | |
66 * > 0: if INFO = k, U(k,k) is exactly zero. The factorization | |
67 * has been completed, but the factor U is exactly | |
68 * singular, and division by zero will occur if it is used | |
69 * to solve a system of equations. | |
70 * | |
71 * ===================================================================== | |
72 * | |
73 * .. Parameters .. | |
74 REAL ZERO | |
75 PARAMETER ( ZERO = 0.0E+0 ) | |
76 * .. | |
77 * .. Local Scalars .. | |
78 INTEGER I | |
79 REAL FACT, TEMP | |
80 * .. | |
81 * .. Intrinsic Functions .. | |
82 INTRINSIC ABS | |
83 * .. | |
84 * .. External Subroutines .. | |
85 EXTERNAL XERBLA | |
86 * .. | |
87 * .. Executable Statements .. | |
88 * | |
89 INFO = 0 | |
90 IF( N.LT.0 ) THEN | |
91 INFO = -1 | |
92 CALL XERBLA( 'SGTTRF', -INFO ) | |
93 RETURN | |
94 END IF | |
95 * | |
96 * Quick return if possible | |
97 * | |
98 IF( N.EQ.0 ) | |
99 $ RETURN | |
100 * | |
101 * Initialize IPIV(i) = i and DU2(I) = 0 | |
102 * | |
103 DO 10 I = 1, N | |
104 IPIV( I ) = I | |
105 10 CONTINUE | |
106 DO 20 I = 1, N - 2 | |
107 DU2( I ) = ZERO | |
108 20 CONTINUE | |
109 * | |
110 DO 30 I = 1, N - 2 | |
111 IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN | |
112 * | |
113 * No row interchange required, eliminate DL(I) | |
114 * | |
115 IF( D( I ).NE.ZERO ) THEN | |
116 FACT = DL( I ) / D( I ) | |
117 DL( I ) = FACT | |
118 D( I+1 ) = D( I+1 ) - FACT*DU( I ) | |
119 END IF | |
120 ELSE | |
121 * | |
122 * Interchange rows I and I+1, eliminate DL(I) | |
123 * | |
124 FACT = D( I ) / DL( I ) | |
125 D( I ) = DL( I ) | |
126 DL( I ) = FACT | |
127 TEMP = DU( I ) | |
128 DU( I ) = D( I+1 ) | |
129 D( I+1 ) = TEMP - FACT*D( I+1 ) | |
130 DU2( I ) = DU( I+1 ) | |
131 DU( I+1 ) = -FACT*DU( I+1 ) | |
132 IPIV( I ) = I + 1 | |
133 END IF | |
134 30 CONTINUE | |
135 IF( N.GT.1 ) THEN | |
136 I = N - 1 | |
137 IF( ABS( D( I ) ).GE.ABS( DL( I ) ) ) THEN | |
138 IF( D( I ).NE.ZERO ) THEN | |
139 FACT = DL( I ) / D( I ) | |
140 DL( I ) = FACT | |
141 D( I+1 ) = D( I+1 ) - FACT*DU( I ) | |
142 END IF | |
143 ELSE | |
144 FACT = D( I ) / DL( I ) | |
145 D( I ) = DL( I ) | |
146 DL( I ) = FACT | |
147 TEMP = DU( I ) | |
148 DU( I ) = D( I+1 ) | |
149 D( I+1 ) = TEMP - FACT*D( I+1 ) | |
150 IPIV( I ) = I + 1 | |
151 END IF | |
152 END IF | |
153 * | |
154 * Check for a zero on the diagonal of U. | |
155 * | |
156 DO 40 I = 1, N | |
157 IF( D( I ).EQ.ZERO ) THEN | |
158 INFO = I | |
159 GO TO 50 | |
160 END IF | |
161 40 CONTINUE | |
162 50 CONTINUE | |
163 * | |
164 RETURN | |
165 * | |
166 * End of SGTTRF | |
167 * | |
168 END |