 LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ 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

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```

Definition at line 176 of file sgelsx.f.

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