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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
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: