LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
sgelsx.f
Go to the documentation of this file.
1*> \brief <b> SGELSX solves overdetermined or underdetermined systems for GE matrices</b>
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download SGELSX + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgelsx.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgelsx.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgelsx.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SGELSX( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
22* WORK, INFO )
23*
24* .. Scalar Arguments ..
25* INTEGER INFO, LDA, LDB, M, N, NRHS, RANK
26* REAL RCOND
27* ..
28* .. Array Arguments ..
29* INTEGER JPVT( * )
30* REAL A( LDA, * ), B( LDB, * ), WORK( * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> This routine is deprecated and has been replaced by routine SGELSY.
40*>
41*> SGELSX computes the minimum-norm solution to a real linear least
42*> squares problem:
43*> minimize || A * X - B ||
44*> using a complete orthogonal factorization of A. A is an M-by-N
45*> matrix which may be rank-deficient.
46*>
47*> Several right hand side vectors b and solution vectors x can be
48*> handled in a single call; they are stored as the columns of the
49*> M-by-NRHS right hand side matrix B and the N-by-NRHS solution
50*> matrix X.
51*>
52*> The routine first computes a QR factorization with column pivoting:
53*> A * P = Q * [ R11 R12 ]
54*> [ 0 R22 ]
55*> with R11 defined as the largest leading submatrix whose estimated
56*> condition number is less than 1/RCOND. The order of R11, RANK,
57*> is the effective rank of A.
58*>
59*> Then, R22 is considered to be negligible, and R12 is annihilated
60*> by orthogonal transformations from the right, arriving at the
61*> complete orthogonal factorization:
62*> A * P = Q * [ T11 0 ] * Z
63*> [ 0 0 ]
64*> The minimum-norm solution is then
65*> X = P * Z**T [ inv(T11)*Q1**T*B ]
66*> [ 0 ]
67*> where Q1 consists of the first RANK columns of Q.
68*> \endverbatim
69*
70* Arguments:
71* ==========
72*
73*> \param[in] M
74*> \verbatim
75*> M is INTEGER
76*> The number of rows of the matrix A. M >= 0.
77*> \endverbatim
78*>
79*> \param[in] N
80*> \verbatim
81*> N is INTEGER
82*> The number of columns of the matrix A. N >= 0.
83*> \endverbatim
84*>
85*> \param[in] NRHS
86*> \verbatim
87*> NRHS is INTEGER
88*> The number of right hand sides, i.e., the number of
89*> columns of matrices B and X. NRHS >= 0.
90*> \endverbatim
91*>
92*> \param[in,out] A
93*> \verbatim
94*> A is REAL array, dimension (LDA,N)
95*> On entry, the M-by-N matrix A.
96*> On exit, A has been overwritten by details of its
97*> complete orthogonal factorization.
98*> \endverbatim
99*>
100*> \param[in] LDA
101*> \verbatim
102*> LDA is INTEGER
103*> The leading dimension of the array A. LDA >= max(1,M).
104*> \endverbatim
105*>
106*> \param[in,out] B
107*> \verbatim
108*> B is REAL array, dimension (LDB,NRHS)
109*> On entry, the M-by-NRHS right hand side matrix B.
110*> On exit, the N-by-NRHS solution matrix X.
111*> If m >= n and RANK = n, the residual sum-of-squares for
112*> the solution in the i-th column is given by the sum of
113*> squares of elements N+1:M in that column.
114*> \endverbatim
115*>
116*> \param[in] LDB
117*> \verbatim
118*> LDB is INTEGER
119*> The leading dimension of the array B. LDB >= max(1,M,N).
120*> \endverbatim
121*>
122*> \param[in,out] JPVT
123*> \verbatim
124*> JPVT is INTEGER array, dimension (N)
125*> On entry, if JPVT(i) .ne. 0, the i-th column of A is an
126*> initial column, otherwise it is a free column. Before
127*> the QR factorization of A, all initial columns are
128*> permuted to the leading positions; only the remaining
129*> free columns are moved as a result of column pivoting
130*> during the factorization.
131*> On exit, if JPVT(i) = k, then the i-th column of A*P
132*> was the k-th column of A.
133*> \endverbatim
134*>
135*> \param[in] RCOND
136*> \verbatim
137*> RCOND is REAL
138*> RCOND is used to determine the effective rank of A, which
139*> is defined as the order of the largest leading triangular
140*> submatrix R11 in the QR factorization with pivoting of A,
141*> whose estimated condition number < 1/RCOND.
142*> \endverbatim
143*>
144*> \param[out] RANK
145*> \verbatim
146*> RANK is INTEGER
147*> The effective rank of A, i.e., the order of the submatrix
148*> R11. This is the same as the order of the submatrix T11
149*> in the complete orthogonal factorization of A.
150*> \endverbatim
151*>
152*> \param[out] WORK
153*> \verbatim
154*> WORK is REAL array, dimension
155*> (max( min(M,N)+3*N, 2*min(M,N)+NRHS )),
156*> \endverbatim
157*>
158*> \param[out] INFO
159*> \verbatim
160*> INFO is INTEGER
161*> = 0: successful exit
162*> < 0: if INFO = -i, the i-th argument had an illegal value
163*> \endverbatim
164*
165* Authors:
166* ========
167*
168*> \author Univ. of Tennessee
169*> \author Univ. of California Berkeley
170*> \author Univ. of Colorado Denver
171*> \author NAG Ltd.
172*
173*> \ingroup realGEsolve
174*
175* =====================================================================
176 SUBROUTINE sgelsx( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
177 $ WORK, INFO )
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*
432 END
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
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
subroutine sgeqpf(M, N, A, LDA, JPVT, TAU, WORK, INFO)
SGEQPF
Definition: sgeqpf.f:142
subroutine sgelsx(M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, WORK, INFO)
SGELSX solves overdetermined or underdetermined systems for GE matrices
Definition: sgelsx.f:178
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