LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine sgelsx ( 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  INFO 
)

SGELSX solves overdetermined or underdetermined systems for GE matrices

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

Purpose: 
 This routine is deprecated and has been replaced by routine SGELSY.

 SGELSX 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.
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 >= n and RANK = n, the residual sum-of-squares for
          the solution in the i-th column is given by the sum of
          squares of elements N+1:M in that column.
[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 an
          initial column, otherwise it is a free column.  Before
          the QR factorization of A, all initial columns are
          permuted to the leading positions; only the remaining
          free columns are moved as a result of column pivoting
          during the factorization.
          On exit, if JPVT(i) = k, then the i-th column of A*P
          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( min(M,N)+3*N, 2*min(M,N)+NRHS )),
[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

Definition at line 180 of file sgelsx.f.

180 *
181 * -- LAPACK driver routine (version 3.4.0) --
182 * -- LAPACK is a software package provided by Univ. of Tennessee, --
183 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
184 * November 2011
185 *
186 * .. Scalar Arguments ..
187  INTEGER info, lda, ldb, m, n, nrhs, rank
188  REAL rcond
189 * ..
190 * .. Array Arguments ..
191  INTEGER jpvt( * )
192  REAL a( lda, * ), b( ldb, * ), work( * )
193 * ..
194 *
195 * =====================================================================
196 *
197 * .. Parameters ..
198  INTEGER imax, imin
199  parameter ( imax = 1, imin = 2 )
200  REAL zero, one, done, ntdone
201  parameter ( zero = 0.0e0, one = 1.0e0, done = zero,
202  $ ntdone = one )
203 * ..
204 * .. Local Scalars ..
205  INTEGER i, iascl, ibscl, ismax, ismin, j, k, mn
206  REAL anrm, bignum, bnrm, c1, c2, s1, s2, smax,
207  $ smaxpr, smin, sminpr, smlnum, t1, t2
208 * ..
209 * .. External Functions ..
210  REAL slamch, slange
211  EXTERNAL slamch, slange
212 * ..
213 * .. External Subroutines ..
214  EXTERNAL sgeqpf, slabad, slaic1, slascl, slaset, slatzm,
215  $ sorm2r, strsm, stzrqf, xerbla
216 * ..
217 * .. Intrinsic Functions ..
218  INTRINSIC abs, max, min
219 * ..
220 * .. Executable Statements ..
221 *
222  mn = min( m, n )
223  ismin = mn + 1
224  ismax = 2*mn + 1
225 *
226 * Test the input arguments.
227 *
228  info = 0
229  IF( m.LT.0 ) THEN
230  info = -1
231  ELSE IF( n.LT.0 ) THEN
232  info = -2
233  ELSE IF( nrhs.LT.0 ) THEN
234  info = -3
235  ELSE IF( lda.LT.max( 1, m ) ) THEN
236  info = -5
237  ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
238  info = -7
239  END IF
240 *
241  IF( info.NE.0 ) THEN
242  CALL xerbla( 'SGELSX', -info )
243  RETURN
244  END IF
245 *
246 * Quick return if possible
247 *
248  IF( min( m, n, nrhs ).EQ.0 ) THEN
249  rank = 0
250  RETURN
251  END IF
252 *
253 * Get machine parameters
254 *
255  smlnum = slamch( 'S' ) / slamch( 'P' )
256  bignum = one / smlnum
257  CALL slabad( smlnum, bignum )
258 *
259 * Scale A, B if max elements outside range [SMLNUM,BIGNUM]
260 *
261  anrm = slange( 'M', m, n, a, lda, work )
262  iascl = 0
263  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
264 *
265 * Scale matrix norm up to SMLNUM
266 *
267  CALL slascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
268  iascl = 1
269  ELSE IF( anrm.GT.bignum ) THEN
270 *
271 * Scale matrix norm down to BIGNUM
272 *
273  CALL slascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
274  iascl = 2
275  ELSE IF( anrm.EQ.zero ) THEN
276 *
277 * Matrix all zero. Return zero solution.
278 *
279  CALL slaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
280  rank = 0
281  GO TO 100
282  END IF
283 *
284  bnrm = slange( 'M', m, nrhs, b, ldb, work )
285  ibscl = 0
286  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
287 *
288 * Scale matrix norm up to SMLNUM
289 *
290  CALL slascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
291  ibscl = 1
292  ELSE IF( bnrm.GT.bignum ) THEN
293 *
294 * Scale matrix norm down to BIGNUM
295 *
296  CALL slascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
297  ibscl = 2
298  END IF
299 *
300 * Compute QR factorization with column pivoting of A:
301 * A * P = Q * R
302 *
303  CALL sgeqpf( m, n, a, lda, jpvt, work( 1 ), work( mn+1 ), info )
304 *
305 * workspace 3*N. Details of Householder rotations stored
306 * in WORK(1:MN).
307 *
308 * Determine RANK using incremental condition estimation
309 *
310  work( ismin ) = one
311  work( ismax ) = one
312  smax = abs( a( 1, 1 ) )
313  smin = smax
314  IF( abs( a( 1, 1 ) ).EQ.zero ) THEN
315  rank = 0
316  CALL slaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
317  GO TO 100
318  ELSE
319  rank = 1
320  END IF
321 *
322  10 CONTINUE
323  IF( rank.LT.mn ) THEN
324  i = rank + 1
325  CALL slaic1( imin, rank, work( ismin ), smin, a( 1, i ),
326  $ a( i, i ), sminpr, s1, c1 )
327  CALL slaic1( imax, rank, work( ismax ), smax, a( 1, i ),
328  $ a( i, i ), smaxpr, s2, c2 )
329 *
330  IF( smaxpr*rcond.LE.sminpr ) THEN
331  DO 20 i = 1, rank
332  work( ismin+i-1 ) = s1*work( ismin+i-1 )
333  work( ismax+i-1 ) = s2*work( ismax+i-1 )
334  20 CONTINUE
335  work( ismin+rank ) = c1
336  work( ismax+rank ) = c2
337  smin = sminpr
338  smax = smaxpr
339  rank = rank + 1
340  GO TO 10
341  END IF
342  END IF
343 *
344 * Logically partition R = [ R11 R12 ]
345 * [ 0 R22 ]
346 * where R11 = R(1:RANK,1:RANK)
347 *
348 * [R11,R12] = [ T11, 0 ] * Y
349 *
350  IF( rank.LT.n )
351  $ CALL stzrqf( rank, n, a, lda, work( mn+1 ), info )
352 *
353 * Details of Householder rotations stored in WORK(MN+1:2*MN)
354 *
355 * B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
356 *
357  CALL sorm2r( 'Left', 'Transpose', m, nrhs, mn, a, lda, work( 1 ),
358  $ b, ldb, work( 2*mn+1 ), info )
359 *
360 * workspace NRHS
361 *
362 * B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
363 *
364  CALL strsm( 'Left', 'Upper', 'No transpose', 'Non-unit', rank,
365  $ nrhs, one, a, lda, b, ldb )
366 *
367  DO 40 i = rank + 1, n
368  DO 30 j = 1, nrhs
369  b( i, j ) = zero
370  30 CONTINUE
371  40 CONTINUE
372 *
373 * B(1:N,1:NRHS) := Y**T * B(1:N,1:NRHS)
374 *
375  IF( rank.LT.n ) THEN
376  DO 50 i = 1, rank
377  CALL slatzm( 'Left', n-rank+1, nrhs, a( i, rank+1 ), lda,
378  $ work( mn+i ), b( i, 1 ), b( rank+1, 1 ), ldb,
379  $ work( 2*mn+1 ) )
380  50 CONTINUE
381  END IF
382 *
383 * workspace NRHS
384 *
385 * B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
386 *
387  DO 90 j = 1, nrhs
388  DO 60 i = 1, n
389  work( 2*mn+i ) = ntdone
390  60 CONTINUE
391  DO 80 i = 1, n
392  IF( work( 2*mn+i ).EQ.ntdone ) THEN
393  IF( jpvt( i ).NE.i ) THEN
394  k = i
395  t1 = b( k, j )
396  t2 = b( jpvt( k ), j )
397  70 CONTINUE
398  b( jpvt( k ), j ) = t1
399  work( 2*mn+k ) = done
400  t1 = t2
401  k = jpvt( k )
402  t2 = b( jpvt( k ), j )
403  IF( jpvt( k ).NE.i )
404  $ GO TO 70
405  b( i, j ) = t1
406  work( 2*mn+k ) = done
407  END IF
408  END IF
409  80 CONTINUE
410  90 CONTINUE
411 *
412 * Undo scaling
413 *
414  IF( iascl.EQ.1 ) THEN
415  CALL slascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
416  CALL slascl( 'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
417  $ info )
418  ELSE IF( iascl.EQ.2 ) THEN
419  CALL slascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
420  CALL slascl( 'U', 0, 0, bignum, anrm, rank, rank, a, lda,
421  $ info )
422  END IF
423  IF( ibscl.EQ.1 ) THEN
424  CALL slascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
425  ELSE IF( ibscl.EQ.2 ) THEN
426  CALL slascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
427  END IF
428 *
429  100 CONTINUE
430 *
431  RETURN
432 *
433 * End of SGELSX
434 *
strsm
subroutine strsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
STRSM
Definition: strsm.f:183
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:145
slabad
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
sgeqpf
subroutine sgeqpf(M, N, A, LDA, JPVT, TAU, WORK, INFO)
SGEQPF
Definition: sgeqpf.f:144
slaic1
subroutine slaic1(JOB, J, X, SEST, W, GAMMA, SESTPR, S, C)
SLAIC1 applies one step of incremental condition estimation.
Definition: slaic1.f:136
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
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:112
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:116
sorm2r
subroutine sorm2r(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
SORM2R multiplies a general matrix by the orthogonal matrix from a QR factorization determined by sge...
Definition: sorm2r.f:161
slatzm
subroutine slatzm(SIDE, M, N, V, INCV, TAU, C1, C2, LDC, WORK)
SLATZM
Definition: slatzm.f:153
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
stzrqf
subroutine stzrqf(M, N, A, LDA, TAU, INFO)
STZRQF
Definition: stzrqf.f:140

Here is the call graph for this function: