LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ sgelsy()

subroutine sgelsy ( integer m,
integer n,
integer nrhs,
real, dimension( lda, * ) a,
integer lda,
real, dimension( ldb, * ) b,
integer ldb,
integer, dimension( * ) jpvt,
real rcond,
integer rank,
real, dimension( * ) work,
integer lwork,
integer info )

SGELSY solves overdetermined or underdetermined systems for GE matrices

Download SGELSY + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> SGELSY computes the minimum-norm solution to a real linear least
!> squares problem:
!>     minimize || A * X - B ||
!> using a complete orthogonal factorization of A.  A is an M-by-N
!> matrix which may be rank-deficient.
!>
!> Several right hand side vectors b and solution vectors x can be
!> handled in a single call; they are stored as the columns of the
!> M-by-NRHS right hand side matrix B and the N-by-NRHS solution
!> matrix X.
!>
!> The routine first computes a QR factorization with column pivoting:
!>     A * P = Q * [ R11 R12 ]
!>                 [  0  R22 ]
!> with R11 defined as the largest leading submatrix whose estimated
!> condition number is less than 1/RCOND.  The order of R11, RANK,
!> is the effective rank of A.
!>
!> Then, R22 is considered to be negligible, and R12 is annihilated
!> by orthogonal transformations from the right, arriving at the
!> complete orthogonal factorization:
!>    A * P = Q * [ T11 0 ] * Z
!>                [  0  0 ]
!> The minimum-norm solution is then
!>    X = P * Z**T [ inv(T11)*Q1**T*B ]
!>                 [        0         ]
!> where Q1 consists of the first RANK columns of Q.
!>
!> This routine is basically identical to the original xGELSX except
!> three differences:
!>   o The call to the subroutine xGEQPF has been substituted by the
!>     the call to the subroutine xGEQP3. This subroutine is a Blas-3
!>     version of the QR factorization with column pivoting.
!>   o Matrix B (the right hand side) is updated with Blas-3.
!>   o The permutation of matrix B (the right hand side) is faster and
!>     more simple.
!>
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]NRHS
!>          NRHS is INTEGER
!>          The number of right hand sides, i.e., the number of
!>          columns of matrices B and X. NRHS >= 0.
!>
[in,out]A
!>          A is REAL array, dimension (LDA,N)
!>          On entry, the M-by-N matrix A.
!>          On exit, A has been overwritten by details of its
!>          complete orthogonal factorization.
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,M).
!>
[in,out]B
!>          B is REAL array, dimension (LDB,NRHS)
!>          On entry, the M-by-NRHS right hand side matrix B.
!>          On exit, the N-by-NRHS solution matrix X.
!>          If M = 0 or N = 0, B is not referenced.
!>
[in]LDB
!>          LDB is INTEGER
!>          The leading dimension of the array B. LDB >= max(1,M,N).
!>
[in,out]JPVT
!>          JPVT is INTEGER array, dimension (N)
!>          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
!>          to the front of AP, otherwise column i is a free column.
!>          On exit, if JPVT(i) = k, then the i-th column of AP
!>          was the k-th column of A.
!>
[in]RCOND
!>          RCOND is REAL
!>          RCOND is used to determine the effective rank of A, which
!>          is defined as the order of the largest leading triangular
!>          submatrix R11 in the QR factorization with pivoting of A,
!>          whose estimated condition number < 1/RCOND.
!>
[out]RANK
!>          RANK is INTEGER
!>          The effective rank of A, i.e., the order of the submatrix
!>          R11.  This is the same as the order of the submatrix T11
!>          in the complete orthogonal factorization of A.
!>          If NRHS = 0, RANK = 0 on output.
!>
[out]WORK
!>          WORK is REAL array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.
!>          The unblocked strategy requires that:
!>             LWORK >= MAX( MN+3*N+1, 2*MN+NRHS ),
!>          where MN = min( M, N ).
!>          The block algorithm requires that:
!>             LWORK >= MAX( MN+2*N+NB*(N+1), 2*MN+NB*NRHS ),
!>          where NB is an upper bound on the blocksize returned
!>          by ILAENV for the routines SGEQP3, STZRZF, STZRQF, SORMQR,
!>          and SORMRZ.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal size of the WORK array, returns
!>          this value as the first entry of the WORK array, and no error
!>          message related to LWORK is issued by XERBLA.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0: successful exit
!>          < 0: If INFO = -i, the i-th argument had an illegal value.
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain
G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain

Definition at line 202 of file sgelsy.f.

205*
206* -- LAPACK driver routine --
207* -- LAPACK is a software package provided by Univ. of Tennessee, --
208* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
209*
210* .. Scalar Arguments ..
211 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
212 REAL RCOND
213* ..
214* .. Array Arguments ..
215 INTEGER JPVT( * )
216 REAL A( LDA, * ), B( LDB, * ), WORK( * )
217* ..
218*
219* =====================================================================
220*
221* .. Parameters ..
222 INTEGER IMAX, IMIN
223 parameter( imax = 1, imin = 2 )
224 REAL ZERO, ONE
225 parameter( zero = 0.0e+0, one = 1.0e+0 )
226* ..
227* .. Local Scalars ..
228 LOGICAL LQUERY
229 INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, LWKMIN,
230 $ LWKOPT, MN, NB, NB1, NB2, NB3, NB4
231 REAL ANRM, BIGNUM, BNRM, C1, C2, S1, S2, SMAX,
232 $ SMAXPR, SMIN, SMINPR, SMLNUM, WSIZE
233* ..
234* .. External Functions ..
235 INTEGER ILAENV
236 REAL SLAMCH, SLANGE, SROUNDUP_LWORK
237 EXTERNAL ilaenv, slamch, slange,
238 $ sroundup_lwork
239* ..
240* .. External Subroutines ..
241 EXTERNAL scopy, sgeqp3, slaic1, slascl,
242 $ slaset,
243 $ sormqr, sormrz, strsm, stzrzf, xerbla
244* ..
245* .. Intrinsic Functions ..
246 INTRINSIC abs, max, min
247* ..
248* .. Executable Statements ..
249*
250 mn = min( m, n )
251 ismin = mn + 1
252 ismax = 2*mn + 1
253*
254* Test the input arguments.
255*
256 info = 0
257 lquery = ( lwork.EQ.-1 )
258 IF( m.LT.0 ) THEN
259 info = -1
260 ELSE IF( n.LT.0 ) THEN
261 info = -2
262 ELSE IF( nrhs.LT.0 ) THEN
263 info = -3
264 ELSE IF( lda.LT.max( 1, m ) ) THEN
265 info = -5
266 ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
267 info = -7
268 END IF
269*
270* Figure out optimal block size
271*
272 IF( info.EQ.0 ) THEN
273 IF( mn.EQ.0 .OR. nrhs.EQ.0 ) THEN
274 lwkmin = 1
275 lwkopt = 1
276 ELSE
277 nb1 = ilaenv( 1, 'SGEQRF', ' ', m, n, -1, -1 )
278 nb2 = ilaenv( 1, 'SGERQF', ' ', m, n, -1, -1 )
279 nb3 = ilaenv( 1, 'SORMQR', ' ', m, n, nrhs, -1 )
280 nb4 = ilaenv( 1, 'SORMRQ', ' ', m, n, nrhs, -1 )
281 nb = max( nb1, nb2, nb3, nb4 )
282 lwkmin = mn + max( 2*mn, n + 1, mn + nrhs )
283 lwkopt = max( lwkmin,
284 $ mn + 2*n + nb*( n + 1 ), 2*mn + nb*nrhs )
285 END IF
286 work( 1 ) = sroundup_lwork(lwkopt)
287*
288 IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
289 info = -12
290 END IF
291 END IF
292*
293 IF( info.NE.0 ) THEN
294 CALL xerbla( 'SGELSY', -info )
295 RETURN
296 ELSE IF( lquery ) THEN
297 RETURN
298 END IF
299*
300* Quick return if possible
301*
302 IF( mn.EQ.0 .OR. nrhs.EQ.0 ) THEN
303 rank = 0
304 RETURN
305 END IF
306*
307* Get machine parameters
308*
309 smlnum = slamch( 'S' ) / slamch( 'P' )
310 bignum = one / smlnum
311*
312* Scale A, B if max entries outside range [SMLNUM,BIGNUM]
313*
314 anrm = slange( 'M', m, n, a, lda, work )
315 iascl = 0
316 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
317*
318* Scale matrix norm up to SMLNUM
319*
320 CALL slascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
321 iascl = 1
322 ELSE IF( anrm.GT.bignum ) THEN
323*
324* Scale matrix norm down to BIGNUM
325*
326 CALL slascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
327 iascl = 2
328 ELSE IF( anrm.EQ.zero ) THEN
329*
330* Matrix all zero. Return zero solution.
331*
332 CALL slaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
333 rank = 0
334 GO TO 70
335 END IF
336*
337 bnrm = slange( 'M', m, nrhs, b, ldb, work )
338 ibscl = 0
339 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
340*
341* Scale matrix norm up to SMLNUM
342*
343 CALL slascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb,
344 $ info )
345 ibscl = 1
346 ELSE IF( bnrm.GT.bignum ) THEN
347*
348* Scale matrix norm down to BIGNUM
349*
350 CALL slascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb,
351 $ info )
352 ibscl = 2
353 END IF
354*
355* Compute QR factorization with column pivoting of A:
356* A * P = Q * R
357*
358 CALL sgeqp3( m, n, a, lda, jpvt, work( 1 ), work( mn+1 ),
359 $ lwork-mn, info )
360 wsize = real( mn ) + work( mn+1 )
361*
362* workspace: MN+2*N+NB*(N+1).
363* Details of Householder rotations stored in WORK(1:MN).
364*
365* Determine RANK using incremental condition estimation
366*
367 work( ismin ) = one
368 work( ismax ) = one
369 smax = abs( a( 1, 1 ) )
370 smin = smax
371 IF( abs( a( 1, 1 ) ).EQ.zero ) THEN
372 rank = 0
373 CALL slaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
374 GO TO 70
375 ELSE
376 rank = 1
377 END IF
378*
379 10 CONTINUE
380 IF( rank.LT.mn ) THEN
381 i = rank + 1
382 CALL slaic1( imin, rank, work( ismin ), smin, a( 1, i ),
383 $ a( i, i ), sminpr, s1, c1 )
384 CALL slaic1( imax, rank, work( ismax ), smax, a( 1, i ),
385 $ a( i, i ), smaxpr, s2, c2 )
386*
387 IF( smaxpr*rcond.LE.sminpr ) THEN
388 DO 20 i = 1, rank
389 work( ismin+i-1 ) = s1*work( ismin+i-1 )
390 work( ismax+i-1 ) = s2*work( ismax+i-1 )
391 20 CONTINUE
392 work( ismin+rank ) = c1
393 work( ismax+rank ) = c2
394 smin = sminpr
395 smax = smaxpr
396 rank = rank + 1
397 GO TO 10
398 END IF
399 END IF
400*
401* workspace: 3*MN.
402*
403* Logically partition R = [ R11 R12 ]
404* [ 0 R22 ]
405* where R11 = R(1:RANK,1:RANK)
406*
407* [R11,R12] = [ T11, 0 ] * Y
408*
409 IF( rank.LT.n )
410 $ CALL stzrzf( rank, n, a, lda, work( mn+1 ), work( 2*mn+1 ),
411 $ lwork-2*mn, info )
412*
413* workspace: 2*MN.
414* Details of Householder rotations stored in WORK(MN+1:2*MN)
415*
416* B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
417*
418 CALL sormqr( 'Left', 'Transpose', m, nrhs, mn, a, lda,
419 $ work( 1 ),
420 $ b, ldb, work( 2*mn+1 ), lwork-2*mn, info )
421 wsize = max( wsize, real( 2*mn )+work( 2*mn+1 ) )
422*
423* workspace: 2*MN+NB*NRHS.
424*
425* B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
426*
427 CALL strsm( 'Left', 'Upper', 'No transpose', 'Non-unit', rank,
428 $ nrhs, one, a, lda, b, ldb )
429*
430 DO 40 j = 1, nrhs
431 DO 30 i = rank + 1, n
432 b( i, j ) = zero
433 30 CONTINUE
434 40 CONTINUE
435*
436* B(1:N,1:NRHS) := Y**T * B(1:N,1:NRHS)
437*
438 IF( rank.LT.n ) THEN
439 CALL sormrz( 'Left', 'Transpose', n, nrhs, rank, n-rank, a,
440 $ lda, work( mn+1 ), b, ldb, work( 2*mn+1 ),
441 $ lwork-2*mn, info )
442 END IF
443*
444* workspace: 2*MN+NRHS.
445*
446* B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
447*
448 DO 60 j = 1, nrhs
449 DO 50 i = 1, n
450 work( jpvt( i ) ) = b( i, j )
451 50 CONTINUE
452 CALL scopy( n, work( 1 ), 1, b( 1, j ), 1 )
453 60 CONTINUE
454*
455* workspace: N.
456*
457* Undo scaling
458*
459 IF( iascl.EQ.1 ) THEN
460 CALL slascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb,
461 $ info )
462 CALL slascl( 'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
463 $ info )
464 ELSE IF( iascl.EQ.2 ) THEN
465 CALL slascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb,
466 $ info )
467 CALL slascl( 'U', 0, 0, bignum, anrm, rank, rank, a, lda,
468 $ info )
469 END IF
470 IF( ibscl.EQ.1 ) THEN
471 CALL slascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb,
472 $ info )
473 ELSE IF( ibscl.EQ.2 ) THEN
474 CALL slascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb,
475 $ info )
476 END IF
477*
478 70 CONTINUE
479 work( 1 ) = sroundup_lwork(lwkopt)
480*
481 RETURN
482*
483* End of SGELSY
484*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
sgeqp3
subroutine sgeqp3(m, n, a, lda, jpvt, tau, work, lwork, info)
SGEQP3
Definition sgeqp3.f:149
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:160
slaic1
subroutine slaic1(job, j, x, sest, w, gamma, sestpr, s, c)
SLAIC1 applies one step of incremental condition estimation.
Definition slaic1.f:132
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
slange
real function slange(norm, m, n, a, lda, work)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition slange.f:112
slascl
subroutine slascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition slascl.f:142
slaset
subroutine slaset(uplo, m, n, alpha, beta, a, lda)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition slaset.f:108
sroundup_lwork
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
Definition sroundup_lwork.f:59
strsm
subroutine strsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
STRSM
Definition strsm.f:181
stzrzf
subroutine stzrzf(m, n, a, lda, tau, work, lwork, info)
STZRZF
Definition stzrzf.f:149
sormqr
subroutine sormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
SORMQR
Definition sormqr.f:166
sormrz
subroutine sormrz(side, trans, m, n, k, l, a, lda, tau, c, ldc, work, lwork, info)
SORMRZ
Definition sormrz.f:186
Here is the call graph for this function:
Here is the caller graph for this function: