LAPACK  3.6.1 LAPACK: Linear Algebra PACKage
 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.``` [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.``` [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.```
Date
November 2011
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 206 of file sgelsy.f.

206 *
207 * -- LAPACK driver routine (version 3.4.0) --
208 * -- LAPACK is a software package provided by Univ. of Tennessee, --
209 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
210 * November 2011
211 *
212 * .. Scalar Arguments ..
213  INTEGER info, lda, ldb, lwork, m, n, nrhs, rank
214  REAL rcond
215 * ..
216 * .. Array Arguments ..
217  INTEGER jpvt( * )
218  REAL a( lda, * ), b( ldb, * ), work( * )
219 * ..
220 *
221 * =====================================================================
222 *
223 * .. Parameters ..
224  INTEGER imax, imin
225  parameter ( imax = 1, imin = 2 )
226  REAL zero, one
227  parameter ( zero = 0.0e+0, one = 1.0e+0 )
228 * ..
229 * .. Local Scalars ..
230  LOGICAL lquery
231  INTEGER i, iascl, ibscl, ismax, ismin, j, lwkmin,
232  \$ lwkopt, mn, nb, nb1, nb2, nb3, nb4
233  REAL anrm, bignum, bnrm, c1, c2, s1, s2, smax,
234  \$ smaxpr, smin, sminpr, smlnum, wsize
235 * ..
236 * .. External Functions ..
237  INTEGER ilaenv
238  REAL slamch, slange
239  EXTERNAL ilaenv, slamch, slange
240 * ..
241 * .. External Subroutines ..
242  EXTERNAL scopy, sgeqp3, slabad, slaic1, slascl, slaset,
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 ) = 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  CALL slabad( smlnum, bignum )
312 *
313 * Scale A, B if max entries outside range [SMLNUM,BIGNUM]
314 *
315  anrm = slange( 'M', m, n, a, lda, work )
316  iascl = 0
317  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
318 *
319 * Scale matrix norm up to SMLNUM
320 *
321  CALL slascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
322  iascl = 1
323  ELSE IF( anrm.GT.bignum ) THEN
324 *
325 * Scale matrix norm down to BIGNUM
326 *
327  CALL slascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
328  iascl = 2
329  ELSE IF( anrm.EQ.zero ) THEN
330 *
331 * Matrix all zero. Return zero solution.
332 *
333  CALL slaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
334  rank = 0
335  GO TO 70
336  END IF
337 *
338  bnrm = slange( 'M', m, nrhs, b, ldb, work )
339  ibscl = 0
340  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
341 *
342 * Scale matrix norm up to SMLNUM
343 *
344  CALL slascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, 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, info )
351  ibscl = 2
352  END IF
353 *
354 * Compute QR factorization with column pivoting of A:
355 * A * P = Q * R
356 *
357  CALL sgeqp3( m, n, a, lda, jpvt, work( 1 ), work( mn+1 ),
358  \$ lwork-mn, info )
359  wsize = mn + work( mn+1 )
360 *
361 * workspace: MN+2*N+NB*(N+1).
362 * Details of Householder rotations stored in WORK(1:MN).
363 *
364 * Determine RANK using incremental condition estimation
365 *
366  work( ismin ) = one
367  work( ismax ) = one
368  smax = abs( a( 1, 1 ) )
369  smin = smax
370  IF( abs( a( 1, 1 ) ).EQ.zero ) THEN
371  rank = 0
372  CALL slaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
373  GO TO 70
374  ELSE
375  rank = 1
376  END IF
377 *
378  10 CONTINUE
379  IF( rank.LT.mn ) THEN
380  i = rank + 1
381  CALL slaic1( imin, rank, work( ismin ), smin, a( 1, i ),
382  \$ a( i, i ), sminpr, s1, c1 )
383  CALL slaic1( imax, rank, work( ismax ), smax, a( 1, i ),
384  \$ a( i, i ), smaxpr, s2, c2 )
385 *
386  IF( smaxpr*rcond.LE.sminpr ) THEN
387  DO 20 i = 1, rank
388  work( ismin+i-1 ) = s1*work( ismin+i-1 )
389  work( ismax+i-1 ) = s2*work( ismax+i-1 )
390  20 CONTINUE
391  work( ismin+rank ) = c1
392  work( ismax+rank ) = c2
393  smin = sminpr
394  smax = smaxpr
395  rank = rank + 1
396  GO TO 10
397  END IF
398  END IF
399 *
400 * workspace: 3*MN.
401 *
402 * Logically partition R = [ R11 R12 ]
403 * [ 0 R22 ]
404 * where R11 = R(1:RANK,1:RANK)
405 *
406 * [R11,R12] = [ T11, 0 ] * Y
407 *
408  IF( rank.LT.n )
409  \$ CALL stzrzf( rank, n, a, lda, work( mn+1 ), work( 2*mn+1 ),
410  \$ lwork-2*mn, info )
411 *
412 * workspace: 2*MN.
413 * Details of Householder rotations stored in WORK(MN+1:2*MN)
414 *
415 * B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
416 *
417  CALL sormqr( 'Left', 'Transpose', m, nrhs, mn, a, lda, work( 1 ),
418  \$ b, ldb, work( 2*mn+1 ), lwork-2*mn, info )
419  wsize = max( wsize, 2*mn+work( 2*mn+1 ) )
420 *
421 * workspace: 2*MN+NB*NRHS.
422 *
423 * B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
424 *
425  CALL strsm( 'Left', 'Upper', 'No transpose', 'Non-unit', rank,
426  \$ nrhs, one, a, lda, b, ldb )
427 *
428  DO 40 j = 1, nrhs
429  DO 30 i = rank + 1, n
430  b( i, j ) = zero
431  30 CONTINUE
432  40 CONTINUE
433 *
434 * B(1:N,1:NRHS) := Y**T * B(1:N,1:NRHS)
435 *
436  IF( rank.LT.n ) THEN
437  CALL sormrz( 'Left', 'Transpose', n, nrhs, rank, n-rank, a,
438  \$ lda, work( mn+1 ), b, ldb, work( 2*mn+1 ),
439  \$ lwork-2*mn, info )
440  END IF
441 *
442 * workspace: 2*MN+NRHS.
443 *
444 * B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
445 *
446  DO 60 j = 1, nrhs
447  DO 50 i = 1, n
448  work( jpvt( i ) ) = b( i, j )
449  50 CONTINUE
450  CALL scopy( n, work( 1 ), 1, b( 1, j ), 1 )
451  60 CONTINUE
452 *
453 * workspace: N.
454 *
455 * Undo scaling
456 *
457  IF( iascl.EQ.1 ) THEN
458  CALL slascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
459  CALL slascl( 'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
460  \$ info )
461  ELSE IF( iascl.EQ.2 ) THEN
462  CALL slascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
463  CALL slascl( 'U', 0, 0, bignum, anrm, rank, rank, a, lda,
464  \$ info )
465  END IF
466  IF( ibscl.EQ.1 ) THEN
467  CALL slascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
468  ELSE IF( ibscl.EQ.2 ) THEN
469  CALL slascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
470  END IF
471 *
472  70 CONTINUE
473  work( 1 ) = lwkopt
474 *
475  RETURN
476 *
477 * End of SGELSY
478 *
subroutine stzrzf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
STZRZF
Definition: stzrzf.f:153
subroutine sormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMQR
Definition: sormqr.f:170
subroutine strsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
STRSM
Definition: strsm.f:183
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:145
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
subroutine slaic1(JOB, J, X, SEST, W, GAMMA, SESTPR, S, C)
SLAIC1 applies one step of incremental condition estimation.
Definition: slaic1.f:136
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
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:112
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:116
subroutine sormrz(SIDE, TRANS, M, N, K, L, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMRZ
Definition: sormrz.f:189
subroutine sgeqp3(M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO)
SGEQP3
Definition: sgeqp3.f:153
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53

Here is the call graph for this function:

Here is the caller graph for this function: