 LAPACK  3.10.1 LAPACK: Linear Algebra PACKage

## ◆ slaorhr_col_getrfnp2()

 recursive subroutine slaorhr_col_getrfnp2 ( integer M, integer N, real, dimension( lda, * ) A, integer LDA, real, dimension( * ) D, integer INFO )

SLAORHR_COL_GETRFNP2

Purpose:
``` SLAORHR_COL_GETRFNP2 computes the modified LU factorization without
pivoting of a real general M-by-N matrix A. The factorization has
the form:

A - S = L * U,

where:
S is a m-by-n diagonal sign matrix with the diagonal D, so that
D(i) = S(i,i), 1 <= i <= min(M,N). The diagonal D is constructed
as D(i)=-SIGN(A(i,i)), where A(i,i) is the value after performing
i-1 steps of Gaussian elimination. This means that the diagonal
element at each step of "modified" Gaussian elimination is at
least one in absolute value (so that division-by-zero not
possible during the division by the diagonal element);

L is a M-by-N lower triangular matrix with unit diagonal elements
(lower trapezoidal if M > N);

and U is a M-by-N upper triangular matrix
(upper trapezoidal if M < N).

This routine is an auxiliary routine used in the Householder
reconstruction routine SORHR_COL. In SORHR_COL, this routine is
applied to an M-by-N matrix A with orthonormal columns, where each
element is bounded by one in absolute value. With the choice of
the matrix S above, one can show that the diagonal element at each
step of Gaussian elimination is the largest (in absolute value) in
the column on or below the diagonal, so that no pivoting is required
for numerical stability .

For more details on the Householder reconstruction algorithm,
including the modified LU factorization, see .

This is the recursive version of the LU factorization algorithm.
Denote A - S by B. The algorithm divides the matrix B into four
submatrices:

[  B11 | B12  ]  where B11 is n1 by n1,
B = [ -----|----- ]        B21 is (m-n1) by n1,
[  B21 | B22  ]        B12 is n1 by n2,
B22 is (m-n1) by n2,
with n1 = min(m,n)/2, n2 = n-n1.

The subroutine calls itself to factor B11, solves for B21,
solves for B12, updates B22, then calls itself to factor B22.

For more details on the recursive LU algorithm, see .

SLAORHR_COL_GETRFNP2 is called to factorize a block by the blocked
routine SLAORHR_COL_GETRFNP, which uses blocked code calling
Level 3 BLAS to update the submatrix. However, SLAORHR_COL_GETRFNP2
is self-sufficient and can be used without SLAORHR_COL_GETRFNP.

 "Reconstructing Householder vectors from tall-skinny QR",
G. Ballard, J. Demmel, L. Grigori, M. Jacquelin, H.D. Nguyen,
E. Solomonik, J. Parallel Distrib. Comput.,
vol. 85, pp. 3-31, 2015.

 "Recursion leads to automatic variable blocking for dense linear
algebra algorithms", F. Gustavson, IBM J. of Res. and Dev.,
vol. 41, no. 6, pp. 737-755, 1997.```
Parameters
 [in] M ``` M is INTEGER The number of rows of the matrix A. M >= 0.``` [in] N ``` N is INTEGER The number of columns of the matrix A. N >= 0.``` [in,out] A ``` A is REAL array, dimension (LDA,N) On entry, the M-by-N matrix to be factored. On exit, the factors L and U from the factorization A-S=L*U; the unit diagonal elements of L are not stored.``` [in] LDA ``` LDA is INTEGER The leading dimension of the array A. LDA >= max(1,M).``` [out] D ``` D is REAL array, dimension min(M,N) The diagonal elements of the diagonal M-by-N sign matrix S, D(i) = S(i,i), where 1 <= i <= min(M,N). The elements can be only plus or minus one.``` [out] INFO ``` INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value```
Contributors:
``` November 2019, Igor Kozachenko,
Computer Science Division,
University of California, Berkeley```

Definition at line 166 of file slaorhr_col_getrfnp2.f.

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  REAL A( LDA, * ), D( * )
178 * ..
179 *
180 * =====================================================================
181 *
182 * .. Parameters ..
183  REAL ONE
184  parameter( one = 1.0e+0 )
185 * ..
186 * .. Local Scalars ..
187  REAL SFMIN
188  INTEGER I, IINFO, N1, N2
189 * ..
190 * .. External Functions ..
191  REAL SLAMCH
192  EXTERNAL slamch
193 * ..
194 * .. External Subroutines ..
195  EXTERNAL sgemm, sscal, strsm, xerbla
196 * ..
197 * .. Intrinsic Functions ..
198  INTRINSIC abs, sign, max, min
199 * ..
200 * .. Executable Statements ..
201 *
202 * Test the input parameters
203 *
204  info = 0
205  IF( m.LT.0 ) THEN
206  info = -1
207  ELSE IF( n.LT.0 ) THEN
208  info = -2
209  ELSE IF( lda.LT.max( 1, m ) ) THEN
210  info = -4
211  END IF
212  IF( info.NE.0 ) THEN
213  CALL xerbla( 'SLAORHR_COL_GETRFNP2', -info )
214  RETURN
215  END IF
216 *
217 * Quick return if possible
218 *
219  IF( min( m, n ).EQ.0 )
220  \$ RETURN
221
222  IF ( m.EQ.1 ) THEN
223 *
224 * One row case, (also recursion termination case),
225 * use unblocked code
226 *
227 * Transfer the sign
228 *
229  d( 1 ) = -sign( one, a( 1, 1 ) )
230 *
231 * Construct the row of U
232 *
233  a( 1, 1 ) = a( 1, 1 ) - d( 1 )
234 *
235  ELSE IF( n.EQ.1 ) THEN
236 *
237 * One column case, (also recursion termination case),
238 * use unblocked code
239 *
240 * Transfer the sign
241 *
242  d( 1 ) = -sign( one, a( 1, 1 ) )
243 *
244 * Construct the row of U
245 *
246  a( 1, 1 ) = a( 1, 1 ) - d( 1 )
247 *
248 * Scale the elements 2:M of the column
249 *
250 * Determine machine safe minimum
251 *
252  sfmin = slamch('S')
253 *
254 * Construct the subdiagonal elements of L
255 *
256  IF( abs( a( 1, 1 ) ) .GE. sfmin ) THEN
257  CALL sscal( m-1, one / a( 1, 1 ), a( 2, 1 ), 1 )
258  ELSE
259  DO i = 2, m
260  a( i, 1 ) = a( i, 1 ) / a( 1, 1 )
261  END DO
262  END IF
263 *
264  ELSE
265 *
266 * Divide the matrix B into four submatrices
267 *
268  n1 = min( m, n ) / 2
269  n2 = n-n1
270
271 *
272 * Factor B11, recursive call
273 *
274  CALL slaorhr_col_getrfnp2( n1, n1, a, lda, d, iinfo )
275 *
276 * Solve for B21
277 *
278  CALL strsm( 'R', 'U', 'N', 'N', m-n1, n1, one, a, lda,
279  \$ a( n1+1, 1 ), lda )
280 *
281 * Solve for B12
282 *
283  CALL strsm( 'L', 'L', 'N', 'U', n1, n2, one, a, lda,
284  \$ a( 1, n1+1 ), lda )
285 *
286 * Update B22, i.e. compute the Schur complement
287 * B22 := B22 - B21*B12
288 *
289  CALL sgemm( 'N', 'N', m-n1, n2, n1, -one, a( n1+1, 1 ), lda,
290  \$ a( 1, n1+1 ), lda, one, a( n1+1, n1+1 ), lda )
291 *
292 * Factor B22, recursive call
293 *
294  CALL slaorhr_col_getrfnp2( m-n1, n2, a( n1+1, n1+1 ), lda,
295  \$ d( n1+1 ), iinfo )
296 *
297  END IF
298  RETURN
299 *
300 * End of SLAORHR_COL_GETRFNP2
301 *
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
recursive subroutine slaorhr_col_getrfnp2(M, N, A, LDA, D, INFO)
SLAORHR_COL_GETRFNP2
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
subroutine strsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
STRSM
Definition: strsm.f:181
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: