Mercurial > octave

comparison libcruft/lapack/zgetrf.f @ 2329:30c606bec7a8

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1996-07-19 01:29:05 by jwe] Initial revision
author jwe
date Fri, 19 Jul 1996 01:29:55 +0000
parents
children 15cddaacbc2d
comparison
equal deleted inserted replaced
2328:b44c3b2a5fce 2329:30c606bec7a8
1 SUBROUTINE ZGETRF( M, N, A, LDA, IPIV, INFO )
2 *
3 * -- LAPACK routine (version 2.0) --
4 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
5 * Courant Institute, Argonne National Lab, and Rice University
6 * September 30, 1994
7 *
8 * .. Scalar Arguments ..
9 INTEGER INFO, LDA, M, N
10 * ..
11 * .. Array Arguments ..
12 INTEGER IPIV( * )
13 COMPLEX*16 A( LDA, * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * ZGETRF computes an LU factorization of a general M-by-N matrix A
20 * using partial pivoting with row interchanges.
21 *
22 * The factorization has the form
23 * A = P * L * U
24 * where P is a permutation matrix, L is lower triangular with unit
25 * diagonal elements (lower trapezoidal if m > n), and U is upper
26 * triangular (upper trapezoidal if m < n).
27 *
28 * This is the right-looking Level 3 BLAS version of the algorithm.
29 *
30 * Arguments
31 * =========
32 *
33 * M (input) INTEGER
34 * The number of rows of the matrix A. M >= 0.
35 *
36 * N (input) INTEGER
37 * The number of columns of the matrix A. N >= 0.
38 *
39 * A (input/output) COMPLEX*16 array, dimension (LDA,N)
40 * On entry, the M-by-N matrix to be factored.
41 * On exit, the factors L and U from the factorization
42 * A = P*L*U; the unit diagonal elements of L are not stored.
43 *
44 * LDA (input) INTEGER
45 * The leading dimension of the array A. LDA >= max(1,M).
46 *
47 * IPIV (output) INTEGER array, dimension (min(M,N))
48 * The pivot indices; for 1 <= i <= min(M,N), row i of the
49 * matrix was interchanged with row IPIV(i).
50 *
51 * INFO (output) INTEGER
52 * = 0: successful exit
53 * < 0: if INFO = -i, the i-th argument had an illegal value
54 * > 0: if INFO = i, U(i,i) is exactly zero. The factorization
55 * has been completed, but the factor U is exactly
56 * singular, and division by zero will occur if it is used
57 * to solve a system of equations.
58 *
59 * =====================================================================
60 *
61 * .. Parameters ..
62 COMPLEX*16 ONE
63 PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
64 * ..
65 * .. Local Scalars ..
66 INTEGER I, IINFO, J, JB, NB
67 * ..
68 * .. External Subroutines ..
69 EXTERNAL XERBLA, ZGEMM, ZGETF2, ZLASWP, ZTRSM
70 * ..
71 * .. External Functions ..
72 INTEGER ILAENV
73 EXTERNAL ILAENV
74 * ..
75 * .. Intrinsic Functions ..
76 INTRINSIC MAX, MIN
77 * ..
78 * .. Executable Statements ..
79 *
80 * Test the input parameters.
81 *
82 INFO = 0
83 IF( M.LT.0 ) THEN
84 INFO = -1
85 ELSE IF( N.LT.0 ) THEN
86 INFO = -2
87 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
88 INFO = -4
89 END IF
90 IF( INFO.NE.0 ) THEN
91 CALL XERBLA( 'ZGETRF', -INFO )
92 RETURN
93 END IF
94 *
95 * Quick return if possible
96 *
97 IF( M.EQ.0 .OR. N.EQ.0 )
98 $ RETURN
99 *
100 * Determine the block size for this environment.
101 *
102 NB = ILAENV( 1, 'ZGETRF', ' ', M, N, -1, -1 )
103 IF( NB.LE.1 .OR. NB.GE.MIN( M, N ) ) THEN
104 *
105 * Use unblocked code.
106 *
107 CALL ZGETF2( M, N, A, LDA, IPIV, INFO )
108 ELSE
109 *
110 * Use blocked code.
111 *
112 DO 20 J = 1, MIN( M, N ), NB
113 JB = MIN( MIN( M, N )-J+1, NB )
114 *
115 * Factor diagonal and subdiagonal blocks and test for exact
116 * singularity.
117 *
118 CALL ZGETF2( M-J+1, JB, A( J, J ), LDA, IPIV( J ), IINFO )
119 *
120 * Adjust INFO and the pivot indices.
121 *
122 IF( INFO.EQ.0 .AND. IINFO.GT.0 )
123 $ INFO = IINFO + J - 1
124 DO 10 I = J, MIN( M, J+JB-1 )
125 IPIV( I ) = J - 1 + IPIV( I )
126 10 CONTINUE
127 *
128 * Apply interchanges to columns 1:J-1.
129 *
130 CALL ZLASWP( J-1, A, LDA, J, J+JB-1, IPIV, 1 )
131 *
132 IF( J+JB.LE.N ) THEN
133 *
134 * Apply interchanges to columns J+JB:N.
135 *
136 CALL ZLASWP( N-J-JB+1, A( 1, J+JB ), LDA, J, J+JB-1,
137 $ IPIV, 1 )
138 *
139 * Compute block row of U.
140 *
141 CALL ZTRSM( 'Left', 'Lower', 'No transpose', 'Unit', JB,
142 $ N-J-JB+1, ONE, A( J, J ), LDA, A( J, J+JB ),
143 $ LDA )
144 IF( J+JB.LE.M ) THEN
145 *
146 * Update trailing submatrix.
147 *
148 CALL ZGEMM( 'No transpose', 'No transpose', M-J-JB+1,
149 $ N-J-JB+1, JB, -ONE, A( J+JB, J ), LDA,
150 $ A( J, J+JB ), LDA, ONE, A( J+JB, J+JB ),
151 $ LDA )
152 END IF
153 END IF
154 20 CONTINUE
155 END IF
156 RETURN
157 *
158 * End of ZGETRF
159 *
160 END