1: *> \brief \b DLAORHR_COL_GETRFNP
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DLAORHR_COL_GETRFNP + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaorhr_col_getrfnp.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaorhr_col_getrfnp.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaorhr_col_getrfnp.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DLAORHR_COL_GETRFNP( M, N, A, LDA, D, INFO )
22: *
23: * .. Scalar Arguments ..
24: * INTEGER INFO, LDA, M, N
25: * ..
26: * .. Array Arguments ..
27: * DOUBLE PRECISION A( LDA, * ), D( * )
28: * ..
29: *
30: *
31: *> \par Purpose:
32: * =============
33: *>
34: *> \verbatim
35: *>
36: *> DLAORHR_COL_GETRFNP computes the modified LU factorization without
37: *> pivoting of a real general M-by-N matrix A. The factorization has
38: *> the form:
39: *>
40: *> A - S = L * U,
41: *>
42: *> where:
43: *> S is a m-by-n diagonal sign matrix with the diagonal D, so that
44: *> D(i) = S(i,i), 1 <= i <= min(M,N). The diagonal D is constructed
45: *> as D(i)=-SIGN(A(i,i)), where A(i,i) is the value after performing
46: *> i-1 steps of Gaussian elimination. This means that the diagonal
47: *> element at each step of "modified" Gaussian elimination is
48: *> at least one in absolute value (so that division-by-zero not
49: *> not possible during the division by the diagonal element);
50: *>
51: *> L is a M-by-N lower triangular matrix with unit diagonal elements
52: *> (lower trapezoidal if M > N);
53: *>
54: *> and U is a M-by-N upper triangular matrix
55: *> (upper trapezoidal if M < N).
56: *>
57: *> This routine is an auxiliary routine used in the Householder
58: *> reconstruction routine DORHR_COL. In DORHR_COL, this routine is
59: *> applied to an M-by-N matrix A with orthonormal columns, where each
60: *> element is bounded by one in absolute value. With the choice of
61: *> the matrix S above, one can show that the diagonal element at each
62: *> step of Gaussian elimination is the largest (in absolute value) in
63: *> the column on or below the diagonal, so that no pivoting is required
64: *> for numerical stability [1].
65: *>
66: *> For more details on the Householder reconstruction algorithm,
67: *> including the modified LU factorization, see [1].
68: *>
69: *> This is the blocked right-looking version of the algorithm,
70: *> calling Level 3 BLAS to update the submatrix. To factorize a block,
71: *> this routine calls the recursive routine DLAORHR_COL_GETRFNP2.
72: *>
73: *> [1] "Reconstructing Householder vectors from tall-skinny QR",
74: *> G. Ballard, J. Demmel, L. Grigori, M. Jacquelin, H.D. Nguyen,
75: *> E. Solomonik, J. Parallel Distrib. Comput.,
76: *> vol. 85, pp. 3-31, 2015.
77: *> \endverbatim
78: *
79: * Arguments:
80: * ==========
81: *
82: *> \param[in] M
83: *> \verbatim
84: *> M is INTEGER
85: *> The number of rows of the matrix A. M >= 0.
86: *> \endverbatim
87: *>
88: *> \param[in] N
89: *> \verbatim
90: *> N is INTEGER
91: *> The number of columns of the matrix A. N >= 0.
92: *> \endverbatim
93: *>
94: *> \param[in,out] A
95: *> \verbatim
96: *> A is DOUBLE PRECISION array, dimension (LDA,N)
97: *> On entry, the M-by-N matrix to be factored.
98: *> On exit, the factors L and U from the factorization
99: *> A-S=L*U; the unit diagonal elements of L are not stored.
100: *> \endverbatim
101: *>
102: *> \param[in] LDA
103: *> \verbatim
104: *> LDA is INTEGER
105: *> The leading dimension of the array A. LDA >= max(1,M).
106: *> \endverbatim
107: *>
108: *> \param[out] D
109: *> \verbatim
110: *> D is DOUBLE PRECISION array, dimension min(M,N)
111: *> The diagonal elements of the diagonal M-by-N sign matrix S,
112: *> D(i) = S(i,i), where 1 <= i <= min(M,N). The elements can
113: *> be only plus or minus one.
114: *> \endverbatim
115: *>
116: *> \param[out] INFO
117: *> \verbatim
118: *> INFO is INTEGER
119: *> = 0: successful exit
120: *> < 0: if INFO = -i, the i-th argument had an illegal value
121: *> \endverbatim
122: *>
123: * Authors:
124: * ========
125: *
126: *> \author Univ. of Tennessee
127: *> \author Univ. of California Berkeley
128: *> \author Univ. of Colorado Denver
129: *> \author NAG Ltd.
130: *
131: *> \date November 2019
132: *
133: *> \ingroup doubleGEcomputational
134: *
135: *> \par Contributors:
136: * ==================
137: *>
138: *> \verbatim
139: *>
140: *> November 2019, Igor Kozachenko,
141: *> Computer Science Division,
142: *> University of California, Berkeley
143: *>
144: *> \endverbatim
145: *
146: * =====================================================================
147: SUBROUTINE DLAORHR_COL_GETRFNP( M, N, A, LDA, D, INFO )
148: IMPLICIT NONE
149: *
150: * -- LAPACK computational routine (version 3.9.0) --
151: * -- LAPACK is a software package provided by Univ. of Tennessee, --
152: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
153: * November 2019
154: *
155: * .. Scalar Arguments ..
156: INTEGER INFO, LDA, M, N
157: * ..
158: * .. Array Arguments ..
159: DOUBLE PRECISION A( LDA, * ), D( * )
160: * ..
161: *
162: * =====================================================================
163: *
164: * .. Parameters ..
165: DOUBLE PRECISION ONE
166: PARAMETER ( ONE = 1.0D+0 )
167: * ..
168: * .. Local Scalars ..
169: INTEGER IINFO, J, JB, NB
170: * ..
171: * .. External Subroutines ..
172: EXTERNAL DGEMM, DLAORHR_COL_GETRFNP2, DTRSM, XERBLA
173: * ..
174: * .. External Functions ..
175: INTEGER ILAENV
176: EXTERNAL ILAENV
177: * ..
178: * .. Intrinsic Functions ..
179: INTRINSIC MAX, MIN
180: * ..
181: * .. Executable Statements ..
182: *
183: * Test the input parameters.
184: *
185: INFO = 0
186: IF( M.LT.0 ) THEN
187: INFO = -1
188: ELSE IF( N.LT.0 ) THEN
189: INFO = -2
190: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
191: INFO = -4
192: END IF
193: IF( INFO.NE.0 ) THEN
194: CALL XERBLA( 'DLAORHR_COL_GETRFNP', -INFO )
195: RETURN
196: END IF
197: *
198: * Quick return if possible
199: *
200: IF( MIN( M, N ).EQ.0 )
201: $ RETURN
202: *
203: * Determine the block size for this environment.
204: *
205:
206: NB = ILAENV( 1, 'DLAORHR_COL_GETRFNP', ' ', M, N, -1, -1 )
207:
208: IF( NB.LE.1 .OR. NB.GE.MIN( M, N ) ) THEN
209: *
210: * Use unblocked code.
211: *
212: CALL DLAORHR_COL_GETRFNP2( M, N, A, LDA, D, INFO )
213: ELSE
214: *
215: * Use blocked code.
216: *
217: DO J = 1, MIN( M, N ), NB
218: JB = MIN( MIN( M, N )-J+1, NB )
219: *
220: * Factor diagonal and subdiagonal blocks.
221: *
222: CALL DLAORHR_COL_GETRFNP2( M-J+1, JB, A( J, J ), LDA,
223: $ D( J ), IINFO )
224: *
225: IF( J+JB.LE.N ) THEN
226: *
227: * Compute block row of U.
228: *
229: CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Unit', JB,
230: $ N-J-JB+1, ONE, A( J, J ), LDA, A( J, J+JB ),
231: $ LDA )
232: IF( J+JB.LE.M ) THEN
233: *
234: * Update trailing submatrix.
235: *
236: CALL DGEMM( 'No transpose', 'No transpose', M-J-JB+1,
237: $ N-J-JB+1, JB, -ONE, A( J+JB, J ), LDA,
238: $ A( J, J+JB ), LDA, ONE, A( J+JB, J+JB ),
239: $ LDA )
240: END IF
241: END IF
242: END DO
243: END IF
244: RETURN
245: *
246: * End of DLAORHR_COL_GETRFNP
247: *
248: END
CVSweb interface <joel.bertrand@systella.fr>