Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
ElmerCSC
GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/mathlibs/src/lapack/cgetrf.f
5177 views
1
SUBROUTINE CGETRF( M, N, A, LDA, IPIV, INFO )
2
*
3
* -- LAPACK routine (version 3.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 A( LDA, * )
14
* ..
15
*
16
* Purpose
17
* =======
18
*
19
* CGETRF 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 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 ONE
63
PARAMETER ( ONE = ( 1.0E+0, 0.0E+0 ) )
64
* ..
65
* .. Local Scalars ..
66
INTEGER I, IINFO, J, JB, NB
67
* ..
68
* .. External Subroutines ..
69
EXTERNAL CGEMM, CGETF2, CLASWP, CTRSM, XERBLA
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( 'CGETRF', -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, 'CGETRF', ' ', M, N, -1, -1 )
103
IF( NB.LE.1 .OR. NB.GE.MIN( M, N ) ) THEN
104
*
105
* Use unblocked code.
106
*
107
CALL CGETF2( 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 CGETF2( 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 CLASWP( 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 CLASWP( 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 CTRSM( '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 CGEMM( '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 CGETRF
159
*
160
END
161
162