LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
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*> Download SGELSX + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgelsx.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgelsx.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgelsx.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE SGELSX( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
20* WORK, INFO )
21*
22* .. Scalar Arguments ..
23* INTEGER INFO, LDA, LDB, M, N, NRHS, RANK
24* REAL RCOND
25* ..
26* .. Array Arguments ..
27* INTEGER JPVT( * )
28* REAL A( LDA, * ), B( LDB, * ), WORK( * )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> This routine is deprecated and has been replaced by routine SGELSY.
38*>
39*> SGELSX computes the minimum-norm solution to a real linear least
40*> squares problem:
41*> minimize || A * X - B ||
42*> using a complete orthogonal factorization of A. A is an M-by-N
43*> matrix which may be rank-deficient.
44*>
45*> Several right hand side vectors b and solution vectors x can be
46*> handled in a single call; they are stored as the columns of the
47*> M-by-NRHS right hand side matrix B and the N-by-NRHS solution
48*> matrix X.
49*>
50*> The routine first computes a QR factorization with column pivoting:
51*> A * P = Q * [ R11 R12 ]
52*> [ 0 R22 ]
53*> with R11 defined as the largest leading submatrix whose estimated
54*> condition number is less than 1/RCOND. The order of R11, RANK,
55*> is the effective rank of A.
56*>
57*> Then, R22 is considered to be negligible, and R12 is annihilated
58*> by orthogonal transformations from the right, arriving at the
59*> complete orthogonal factorization:
60*> A * P = Q * [ T11 0 ] * Z
61*> [ 0 0 ]
62*> The minimum-norm solution is then
63*> X = P * Z**T [ inv(T11)*Q1**T*B ]
64*> [ 0 ]
65*> where Q1 consists of the first RANK columns of Q.
66*> \endverbatim
67*
68* Arguments:
69* ==========
70*
71*> \param[in] M
72*> \verbatim
73*> M is INTEGER
74*> The number of rows of the matrix A. M >= 0.
75*> \endverbatim
76*>
77*> \param[in] N
78*> \verbatim
79*> N is INTEGER
80*> The number of columns of the matrix A. N >= 0.
81*> \endverbatim
82*>
83*> \param[in] NRHS
84*> \verbatim
85*> NRHS is INTEGER
86*> The number of right hand sides, i.e., the number of
87*> columns of matrices B and X. NRHS >= 0.
88*> \endverbatim
89*>
90*> \param[in,out] A
91*> \verbatim
92*> A is REAL array, dimension (LDA,N)
93*> On entry, the M-by-N matrix A.
94*> On exit, A has been overwritten by details of its
95*> complete orthogonal factorization.
96*> \endverbatim
97*>
98*> \param[in] LDA
99*> \verbatim
100*> LDA is INTEGER
101*> The leading dimension of the array A. LDA >= max(1,M).
102*> \endverbatim
103*>
104*> \param[in,out] B
105*> \verbatim
106*> B is REAL array, dimension (LDB,NRHS)
107*> On entry, the M-by-NRHS right hand side matrix B.
108*> On exit, the N-by-NRHS solution matrix X.
109*> If m >= n and RANK = n, the residual sum-of-squares for
110*> the solution in the i-th column is given by the sum of
111*> squares of elements N+1:M in that column.
112*> \endverbatim
113*>
114*> \param[in] LDB
115*> \verbatim
116*> LDB is INTEGER
117*> The leading dimension of the array B. LDB >= max(1,M,N).
118*> \endverbatim
119*>
120*> \param[in,out] JPVT
121*> \verbatim
122*> JPVT is INTEGER array, dimension (N)
123*> On entry, if JPVT(i) .ne. 0, the i-th column of A is an
124*> initial column, otherwise it is a free column. Before
125*> the QR factorization of A, all initial columns are
126*> permuted to the leading positions; only the remaining
127*> free columns are moved as a result of column pivoting
128*> during the factorization.
129*> On exit, if JPVT(i) = k, then the i-th column of A*P
130*> was the k-th column of A.
131*> \endverbatim
132*>
133*> \param[in] RCOND
134*> \verbatim
135*> RCOND is REAL
136*> RCOND is used to determine the effective rank of A, which
137*> is defined as the order of the largest leading triangular
138*> submatrix R11 in the QR factorization with pivoting of A,
139*> whose estimated condition number < 1/RCOND.
140*> \endverbatim
141*>
142*> \param[out] RANK
143*> \verbatim
144*> RANK is INTEGER
145*> The effective rank of A, i.e., the order of the submatrix
146*> R11. This is the same as the order of the submatrix T11
147*> in the complete orthogonal factorization of A.
148*> \endverbatim
149*>
150*> \param[out] WORK
151*> \verbatim
152*> WORK is REAL array, dimension
153*> (max( min(M,N)+3*N, 2*min(M,N)+NRHS )),
154*> \endverbatim
155*>
156*> \param[out] INFO
157*> \verbatim
158*> INFO is INTEGER
159*> = 0: successful exit
160*> < 0: if INFO = -i, the i-th argument had an illegal value
161*> \endverbatim
162*
163* Authors:
164* ========
165*
166*> \author Univ. of Tennessee
167*> \author Univ. of California Berkeley
168*> \author Univ. of Colorado Denver
169*> \author NAG Ltd.
170*
171*> \ingroup realGEsolve
172*
173* =====================================================================
174 SUBROUTINE sgelsx( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND,
175 $ RANK, WORK, INFO )
176*
177* -- LAPACK driver routine --
178* -- LAPACK is a software package provided by Univ. of Tennessee, --
179* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
180*
181* .. Scalar Arguments ..
182 INTEGER INFO, LDA, LDB, M, N, NRHS, RANK
183 REAL RCOND
184* ..
185* .. Array Arguments ..
186 INTEGER JPVT( * )
187 REAL A( LDA, * ), B( LDB, * ), WORK( * )
188* ..
189*
190* =====================================================================
191*
192* .. Parameters ..
193 INTEGER IMAX, IMIN
194 parameter( imax = 1, imin = 2 )
195 REAL ZERO, ONE, DONE, NTDONE
196 parameter( zero = 0.0e0, one = 1.0e0, done = zero,
197 $ ntdone = one )
198* ..
199* .. Local Scalars ..
200 INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, K, MN
201 REAL ANRM, BIGNUM, BNRM, C1, C2, S1, S2, SMAX,
202 $ smaxpr, smin, sminpr, smlnum, t1, t2
203* ..
204* .. External Functions ..
205 REAL SLAMCH, SLANGE
206 EXTERNAL slamch, slange
207* ..
208* .. External Subroutines ..
209 EXTERNAL sgeqpf, slaic1, slascl, slaset, slatzm,
211* ..
212* .. Intrinsic Functions ..
213 INTRINSIC abs, max, min
214* ..
215* .. Executable Statements ..
216*
217 mn = min( m, n )
218 ismin = mn + 1
219 ismax = 2*mn + 1
220*
221* Test the input arguments.
222*
223 info = 0
224 IF( m.LT.0 ) THEN
225 info = -1
226 ELSE IF( n.LT.0 ) THEN
227 info = -2
228 ELSE IF( nrhs.LT.0 ) THEN
229 info = -3
230 ELSE IF( lda.LT.max( 1, m ) ) THEN
231 info = -5
232 ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
233 info = -7
234 END IF
235*
236 IF( info.NE.0 ) THEN
237 CALL xerbla( 'SGELSX', -info )
238 RETURN
239 END IF
240*
241* Quick return if possible
242*
243 IF( min( m, n, nrhs ).EQ.0 ) THEN
244 rank = 0
245 RETURN
246 END IF
247*
248* Get machine parameters
249*
250 smlnum = slamch( 'S' ) / slamch( 'P' )
251 bignum = one / smlnum
252*
253* Scale A, B if max elements outside range [SMLNUM,BIGNUM]
254*
255 anrm = slange( 'M', m, n, a, lda, work )
256 iascl = 0
257 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
258*
259* Scale matrix norm up to SMLNUM
260*
261 CALL slascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
262 iascl = 1
263 ELSE IF( anrm.GT.bignum ) THEN
264*
265* Scale matrix norm down to BIGNUM
266*
267 CALL slascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
268 iascl = 2
269 ELSE IF( anrm.EQ.zero ) THEN
270*
271* Matrix all zero. Return zero solution.
272*
273 CALL slaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
274 rank = 0
275 GO TO 100
276 END IF
277*
278 bnrm = slange( 'M', m, nrhs, b, ldb, work )
279 ibscl = 0
280 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
281*
282* Scale matrix norm up to SMLNUM
283*
284 CALL slascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb,
285 $ info )
286 ibscl = 1
287 ELSE IF( bnrm.GT.bignum ) THEN
288*
289* Scale matrix norm down to BIGNUM
290*
291 CALL slascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb,
292 $ info )
293 ibscl = 2
294 END IF
295*
296* Compute QR factorization with column pivoting of A:
297* A * P = Q * R
298*
299 CALL sgeqpf( m, n, a, lda, jpvt, work( 1 ), work( mn+1 ),
300 $ 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,
355 $ work( 1 ), 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',
362 $ rank, 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,
413 $ info )
414 CALL slascl( 'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
415 $ info )
416 ELSE IF( iascl.EQ.2 ) THEN
417 CALL slascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb,
418 $ info )
419 CALL slascl( 'U', 0, 0, bignum, anrm, rank, rank, a, lda,
420 $ info )
421 END IF
422 IF( ibscl.EQ.1 ) THEN
423 CALL slascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb,
424 $ info )
425 ELSE IF( ibscl.EQ.2 ) THEN
426 CALL slascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb,
427 $ info )
428 END IF
429*
430 100 CONTINUE
431*
432 RETURN
433*
434* End of SGELSX
435*
436 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine slaic1(job, j, x, sest, w, gamma, sestpr, s, c)
SLAIC1 applies one step of incremental condition estimation.
Definition slaic1.f:132
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:142
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:108
subroutine strsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
STRSM
Definition strsm.f:181
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:157
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:176
subroutine sgeqpf(m, n, a, lda, jpvt, tau, work, info)
SGEQPF
Definition sgeqpf.f:140
subroutine slatzm(side, m, n, v, incv, tau, c1, c2, ldc, work)
SLATZM
Definition slatzm.f:150
subroutine stzrqf(m, n, a, lda, tau, info)
STZRQF
Definition stzrqf.f:136