LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
claunhr_col_getrfnp2.f
Go to the documentation of this file.
1 *> \brief \b CLAUNHR_COL_GETRFNP2
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CLAUNHR_COL_GETRFNP2 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/claunhr_col_getrfnp2.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/claunhr_col_getrfnp2.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/claunhr_col_getrfnp2.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * RECURSIVE SUBROUTINE CLAUNHR_COL_GETRFNP2( M, N, A, LDA, D, INFO )
22 *
23 * .. Scalar Arguments ..
24 * INTEGER INFO, LDA, M, N
25 * ..
26 * .. Array Arguments ..
27 * COMPLEX A( LDA, * ), D( * )
28 * ..
29 *
30 *
31 *> \par Purpose:
32 * =============
33 *>
34 *> \verbatim
35 *>
36 *> CLAUNHR_COL_GETRFNP2 computes the modified LU factorization without
37 *> pivoting of a complex 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 at
48 *> least one in absolute value (so that division-by-zero not
49 *> 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 CUNHR_COL. In CUNHR_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 recursive version of the LU factorization algorithm.
70 *> Denote A - S by B. The algorithm divides the matrix B into four
71 *> submatrices:
72 *>
73 *> [ B11 | B12 ] where B11 is n1 by n1,
74 *> B = [ -----|----- ] B21 is (m-n1) by n1,
75 *> [ B21 | B22 ] B12 is n1 by n2,
76 *> B22 is (m-n1) by n2,
77 *> with n1 = min(m,n)/2, n2 = n-n1.
78 *>
79 *>
80 *> The subroutine calls itself to factor B11, solves for B21,
81 *> solves for B12, updates B22, then calls itself to factor B22.
82 *>
83 *> For more details on the recursive LU algorithm, see [2].
84 *>
85 *> CLAUNHR_COL_GETRFNP2 is called to factorize a block by the blocked
86 *> routine CLAUNHR_COL_GETRFNP, which uses blocked code calling
87 *> Level 3 BLAS to update the submatrix. However, CLAUNHR_COL_GETRFNP2
88 *> is self-sufficient and can be used without CLAUNHR_COL_GETRFNP.
89 *>
90 *> [1] "Reconstructing Householder vectors from tall-skinny QR",
91 *> G. Ballard, J. Demmel, L. Grigori, M. Jacquelin, H.D. Nguyen,
92 *> E. Solomonik, J. Parallel Distrib. Comput.,
93 *> vol. 85, pp. 3-31, 2015.
94 *>
95 *> [2] "Recursion leads to automatic variable blocking for dense linear
96 *> algebra algorithms", F. Gustavson, IBM J. of Res. and Dev.,
97 *> vol. 41, no. 6, pp. 737-755, 1997.
98 *> \endverbatim
99 *
100 * Arguments:
101 * ==========
102 *
103 *> \param[in] M
104 *> \verbatim
105 *> M is INTEGER
106 *> The number of rows of the matrix A. M >= 0.
107 *> \endverbatim
108 *>
109 *> \param[in] N
110 *> \verbatim
111 *> N is INTEGER
112 *> The number of columns of the matrix A. N >= 0.
113 *> \endverbatim
114 *>
115 *> \param[in,out] A
116 *> \verbatim
117 *> A is COMPLEX array, dimension (LDA,N)
118 *> On entry, the M-by-N matrix to be factored.
119 *> On exit, the factors L and U from the factorization
120 *> A-S=L*U; the unit diagonal elements of L are not stored.
121 *> \endverbatim
122 *>
123 *> \param[in] LDA
124 *> \verbatim
125 *> LDA is INTEGER
126 *> The leading dimension of the array A. LDA >= max(1,M).
127 *> \endverbatim
128 *>
129 *> \param[out] D
130 *> \verbatim
131 *> D is COMPLEX array, dimension min(M,N)
132 *> The diagonal elements of the diagonal M-by-N sign matrix S,
133 *> D(i) = S(i,i), where 1 <= i <= min(M,N). The elements can be
134 *> only ( +1.0, 0.0 ) or (-1.0, 0.0 ).
135 *> \endverbatim
136 *>
137 *> \param[out] INFO
138 *> \verbatim
139 *> INFO is INTEGER
140 *> = 0: successful exit
141 *> < 0: if INFO = -i, the i-th argument had an illegal value
142 *> \endverbatim
143 *>
144 * Authors:
145 * ========
146 *
147 *> \author Univ. of Tennessee
148 *> \author Univ. of California Berkeley
149 *> \author Univ. of Colorado Denver
150 *> \author NAG Ltd.
151 *
152 *> \ingroup complexGEcomputational
153 *
154 *> \par Contributors:
155 * ==================
156 *>
157 *> \verbatim
158 *>
159 *> November 2019, Igor Kozachenko,
160 *> Computer Science Division,
161 *> University of California, Berkeley
162 *>
163 *> \endverbatim
164 *
165 * =====================================================================
166  RECURSIVE SUBROUTINE claunhr_col_getrfnp2( M, N, A, LDA, D, INFO )
167  IMPLICIT NONE
168 *
169 * -- LAPACK computational routine --
170 * -- LAPACK is a software package provided by Univ. of Tennessee, --
171 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
172 *
173 * .. Scalar Arguments ..
174  INTEGER info, lda, m, n
175 * ..
176 * .. Array Arguments ..
177  COMPLEX a( lda, * ), d( * )
178 * ..
179 *
180 * =====================================================================
181 *
182 * .. Parameters ..
183  REAL one
184  parameter( one = 1.0e+0 )
185  COMPLEX cone
186  parameter( cone = ( 1.0e+0, 0.0e+0 ) )
187 * ..
188 * .. Local Scalars ..
189  REAL sfmin
190  INTEGER i, iinfo, n1, n2
191  COMPLEX z
192 * ..
193 * .. External Functions ..
194  REAL slamch
195  EXTERNAL slamch
196 * ..
197 * .. External Subroutines ..
198  EXTERNAL cgemm, cscal, ctrsm, xerbla
199 * ..
200 * .. Intrinsic Functions ..
201  INTRINSIC abs, real, cmplx, aimag, sign, max, min
202 * ..
203 * .. Statement Functions ..
204  DOUBLE PRECISION cabs1
205 * ..
206 * .. Statement Function definitions ..
207  cabs1( z ) = abs( real( z ) ) + abs( aimag( z ) )
208 * ..
209 * .. Executable Statements ..
210 *
211 * Test the input parameters
212 *
213  info = 0
214  IF( m.LT.0 ) THEN
215  info = -1
216  ELSE IF( n.LT.0 ) THEN
217  info = -2
218  ELSE IF( lda.LT.max( 1, m ) ) THEN
219  info = -4
220  END IF
221  IF( info.NE.0 ) THEN
222  CALL xerbla( 'CLAUNHR_COL_GETRFNP2', -info )
223  RETURN
224  END IF
225 *
226 * Quick return if possible
227 *
228  IF( min( m, n ).EQ.0 )
229  $ RETURN
230 
231  IF ( m.EQ.1 ) THEN
232 *
233 * One row case, (also recursion termination case),
234 * use unblocked code
235 *
236 * Transfer the sign
237 *
238  d( 1 ) = cmplx( -sign( one, real( a( 1, 1 ) ) ) )
239 *
240 * Construct the row of U
241 *
242  a( 1, 1 ) = a( 1, 1 ) - d( 1 )
243 *
244  ELSE IF( n.EQ.1 ) THEN
245 *
246 * One column case, (also recursion termination case),
247 * use unblocked code
248 *
249 * Transfer the sign
250 *
251  d( 1 ) = cmplx( -sign( one, real( a( 1, 1 ) ) ) )
252 *
253 * Construct the row of U
254 *
255  a( 1, 1 ) = a( 1, 1 ) - d( 1 )
256 *
257 * Scale the elements 2:M of the column
258 *
259 * Determine machine safe minimum
260 *
261  sfmin = slamch('S')
262 *
263 * Construct the subdiagonal elements of L
264 *
265  IF( cabs1( a( 1, 1 ) ) .GE. sfmin ) THEN
266  CALL cscal( m-1, cone / a( 1, 1 ), a( 2, 1 ), 1 )
267  ELSE
268  DO i = 2, m
269  a( i, 1 ) = a( i, 1 ) / a( 1, 1 )
270  END DO
271  END IF
272 *
273  ELSE
274 *
275 * Divide the matrix B into four submatrices
276 *
277  n1 = min( m, n ) / 2
278  n2 = n-n1
279 
280 *
281 * Factor B11, recursive call
282 *
283  CALL claunhr_col_getrfnp2( n1, n1, a, lda, d, iinfo )
284 *
285 * Solve for B21
286 *
287  CALL ctrsm( 'R', 'U', 'N', 'N', m-n1, n1, cone, a, lda,
288  $ a( n1+1, 1 ), lda )
289 *
290 * Solve for B12
291 *
292  CALL ctrsm( 'L', 'L', 'N', 'U', n1, n2, cone, a, lda,
293  $ a( 1, n1+1 ), lda )
294 *
295 * Update B22, i.e. compute the Schur complement
296 * B22 := B22 - B21*B12
297 *
298  CALL cgemm( 'N', 'N', m-n1, n2, n1, -cone, a( n1+1, 1 ), lda,
299  $ a( 1, n1+1 ), lda, cone, a( n1+1, n1+1 ), lda )
300 *
301 * Factor B22, recursive call
302 *
303  CALL claunhr_col_getrfnp2( m-n1, n2, a( n1+1, n1+1 ), lda,
304  $ d( n1+1 ), iinfo )
305 *
306  END IF
307  RETURN
308 *
309 * End of CLAUNHR_COL_GETRFNP2
310 *
311  END
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:78
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:187
subroutine ctrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
CTRSM
Definition: ctrsm.f:180
recursive subroutine claunhr_col_getrfnp2(M, N, A, LDA, D, INFO)
CLAUNHR_COL_GETRFNP2
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68