 LAPACK  3.8.0 LAPACK: Linear Algebra PACKage

## ◆ sgelsx()

 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```
Date
December 2016

Definition at line 180 of file sgelsx.f.

180 *
181 * -- LAPACK driver routine (version 3.7.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 * December 2016
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,
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 *
subroutine strsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
STRSM
Definition: strsm.f:183
subroutine slaic1(JOB, J, X, SEST, W, GAMMA, SESTPR, S, C)
SLAIC1 applies one step of incremental condition estimation.
Definition: slaic1.f:136
subroutine sgeqpf(M, N, A, LDA, JPVT, TAU, WORK, INFO)
SGEQPF
Definition: sgeqpf.f:144
subroutine stzrqf(M, N, A, LDA, TAU, INFO)
STZRQF
Definition: stzrqf.f:140
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 xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slatzm(SIDE, M, N, V, INCV, TAU, C1, C2, LDC, WORK)
SLATZM
Definition: slatzm.f:153
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
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
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
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
Here is the call graph for this function: