KTH framework for Nek5000 toolboxes; testing version  0.0.1
dgetrf.f
Go to the documentation of this file.
1  SUBROUTINE dgetrf( 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 * March 31, 1993
7 *
8 * .. Scalar Arguments ..
9  INTEGER INFO, LDA, M, N
10 * ..
11 * .. Array Arguments ..
12  INTEGER IPIV( * )
13  DOUBLE PRECISION A( LDA, * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DGETRF 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) DOUBLE PRECISION 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  DOUBLE PRECISION ONE
63  parameter( one = 1.0d+0 )
64 * ..
65 * .. Local Scalars ..
66  INTEGER I, IINFO, J, JB, NB
67 * ..
68 * .. External Subroutines ..
69  EXTERNAL dgemm, dgetf2, dlaswp, dtrsm, 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( 'DGETRF', -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, 'DGETRF', ' ', m, n, -1, -1 )
103  IF( nb.LE.1 .OR. nb.GE.min( m, n ) ) THEN
104 *
105 * Use unblocked code.
106 *
107  CALL dgetf2( 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 dgetf2( 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 dlaswp( 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 dlaswp( 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 dtrsm( '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 dgemm( '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 DGETRF
159 *
160  END
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
Definition: dgemm.f:3
subroutine dgetf2(M, N, A, LDA, IPIV, INFO)
Definition: dgetf2.f:2
subroutine dgetrf(M, N, A, LDA, IPIV, INFO)
Definition: dgetrf.f:2
subroutine dlaswp(N, A, LDA, K1, K2, IPIV, INCX)
Definition: dlaswp.f:2
subroutine dtrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
Definition: dtrsm.f:3
subroutine xerbla(SRNAME, INFO)
Definition: xerbla.f:2