LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dgelsy.f
Go to the documentation of this file.
1*> \brief <b> DGELSY 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 DGELSY + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgelsy.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgelsy.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgelsy.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
22* WORK, LWORK, INFO )
23*
24* .. Scalar Arguments ..
25* INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
26* DOUBLE PRECISION RCOND
27* ..
28* .. Array Arguments ..
29* INTEGER JPVT( * )
30* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> DGELSY 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*>
67*> This routine is basically identical to the original xGELSX except
68*> three differences:
69*> o The call to the subroutine xGEQPF has been substituted by the
70*> the call to the subroutine xGEQP3. This subroutine is a Blas-3
71*> version of the QR factorization with column pivoting.
72*> o Matrix B (the right hand side) is updated with Blas-3.
73*> o The permutation of matrix B (the right hand side) is faster and
74*> more simple.
75*> \endverbatim
76*
77* Arguments:
78* ==========
79*
80*> \param[in] M
81*> \verbatim
82*> M is INTEGER
83*> The number of rows of the matrix A. M >= 0.
84*> \endverbatim
85*>
86*> \param[in] N
87*> \verbatim
88*> N is INTEGER
89*> The number of columns of the matrix A. N >= 0.
90*> \endverbatim
91*>
92*> \param[in] NRHS
93*> \verbatim
94*> NRHS is INTEGER
95*> The number of right hand sides, i.e., the number of
96*> columns of matrices B and X. NRHS >= 0.
97*> \endverbatim
98*>
99*> \param[in,out] A
100*> \verbatim
101*> A is DOUBLE PRECISION array, dimension (LDA,N)
102*> On entry, the M-by-N matrix A.
103*> On exit, A has been overwritten by details of its
104*> complete orthogonal factorization.
105*> \endverbatim
106*>
107*> \param[in] LDA
108*> \verbatim
109*> LDA is INTEGER
110*> The leading dimension of the array A. LDA >= max(1,M).
111*> \endverbatim
112*>
113*> \param[in,out] B
114*> \verbatim
115*> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
116*> On entry, the M-by-NRHS right hand side matrix B.
117*> On exit, the N-by-NRHS solution matrix X.
118*> If M = 0 or N = 0, B is not referenced.
119*> \endverbatim
120*>
121*> \param[in] LDB
122*> \verbatim
123*> LDB is INTEGER
124*> The leading dimension of the array B. LDB >= max(1,M,N).
125*> \endverbatim
126*>
127*> \param[in,out] JPVT
128*> \verbatim
129*> JPVT is INTEGER array, dimension (N)
130*> On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
131*> to the front of AP, otherwise column i is a free column.
132*> On exit, if JPVT(i) = k, then the i-th column of AP
133*> was the k-th column of A.
134*> \endverbatim
135*>
136*> \param[in] RCOND
137*> \verbatim
138*> RCOND is DOUBLE PRECISION
139*> RCOND is used to determine the effective rank of A, which
140*> is defined as the order of the largest leading triangular
141*> submatrix R11 in the QR factorization with pivoting of A,
142*> whose estimated condition number < 1/RCOND.
143*> \endverbatim
144*>
145*> \param[out] RANK
146*> \verbatim
147*> RANK is INTEGER
148*> The effective rank of A, i.e., the order of the submatrix
149*> R11. This is the same as the order of the submatrix T11
150*> in the complete orthogonal factorization of A.
151*> If NRHS = 0, RANK = 0 on output.
152*> \endverbatim
153*>
154*> \param[out] WORK
155*> \verbatim
156*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
157*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
158*> \endverbatim
159*>
160*> \param[in] LWORK
161*> \verbatim
162*> LWORK is INTEGER
163*> The dimension of the array WORK.
164*> The unblocked strategy requires that:
165*> LWORK >= MAX( MN+3*N+1, 2*MN+NRHS ),
166*> where MN = min( M, N ).
167*> The block algorithm requires that:
168*> LWORK >= MAX( MN+2*N+NB*(N+1), 2*MN+NB*NRHS ),
169*> where NB is an upper bound on the blocksize returned
170*> by ILAENV for the routines DGEQP3, DTZRZF, STZRQF, DORMQR,
171*> and DORMRZ.
172*>
173*> If LWORK = -1, then a workspace query is assumed; the routine
174*> only calculates the optimal size of the WORK array, returns
175*> this value as the first entry of the WORK array, and no error
176*> message related to LWORK is issued by XERBLA.
177*> \endverbatim
178*>
179*> \param[out] INFO
180*> \verbatim
181*> INFO is INTEGER
182*> = 0: successful exit
183*> < 0: If INFO = -i, the i-th argument had an illegal value.
184*> \endverbatim
185*
186* Authors:
187* ========
188*
189*> \author Univ. of Tennessee
190*> \author Univ. of California Berkeley
191*> \author Univ. of Colorado Denver
192*> \author NAG Ltd.
193*
194*> \ingroup gelsy
195*
196*> \par Contributors:
197* ==================
198*>
199*> A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA \n
200*> E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain \n
201*> G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain \n
202*>
203* =====================================================================
204 SUBROUTINE dgelsy( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
205 $ WORK, LWORK, INFO )
206*
207* -- LAPACK driver routine --
208* -- LAPACK is a software package provided by Univ. of Tennessee, --
209* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
210*
211* .. Scalar Arguments ..
212 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
213 DOUBLE PRECISION RCOND
214* ..
215* .. Array Arguments ..
216 INTEGER JPVT( * )
217 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
218* ..
219*
220* =====================================================================
221*
222* .. Parameters ..
223 INTEGER IMAX, IMIN
224 parameter( imax = 1, imin = 2 )
225 DOUBLE PRECISION ZERO, ONE
226 parameter( zero = 0.0d+0, one = 1.0d+0 )
227* ..
228* .. Local Scalars ..
229 LOGICAL LQUERY
230 INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, LWKMIN,
231 $ lwkopt, mn, nb, nb1, nb2, nb3, nb4
232 DOUBLE PRECISION ANRM, BIGNUM, BNRM, C1, C2, S1, S2, SMAX,
233 $ smaxpr, smin, sminpr, smlnum, wsize
234* ..
235* .. External Functions ..
236 INTEGER ILAENV
237 DOUBLE PRECISION DLAMCH, DLANGE
238 EXTERNAL ilaenv, dlamch, dlange
239* ..
240* .. External Subroutines ..
241 EXTERNAL dcopy, dgeqp3, dlaic1, dlascl, dlaset,
243* ..
244* .. Intrinsic Functions ..
245 INTRINSIC abs, max, min
246* ..
247* .. Executable Statements ..
248*
249 mn = min( m, n )
250 ismin = mn + 1
251 ismax = 2*mn + 1
252*
253* Test the input arguments.
254*
255 info = 0
256 lquery = ( lwork.EQ.-1 )
257 IF( m.LT.0 ) THEN
258 info = -1
259 ELSE IF( n.LT.0 ) THEN
260 info = -2
261 ELSE IF( nrhs.LT.0 ) THEN
262 info = -3
263 ELSE IF( lda.LT.max( 1, m ) ) THEN
264 info = -5
265 ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
266 info = -7
267 END IF
268*
269* Figure out optimal block size
270*
271 IF( info.EQ.0 ) THEN
272 IF( mn.EQ.0 .OR. nrhs.EQ.0 ) THEN
273 lwkmin = 1
274 lwkopt = 1
275 ELSE
276 nb1 = ilaenv( 1, 'DGEQRF', ' ', m, n, -1, -1 )
277 nb2 = ilaenv( 1, 'DGERQF', ' ', m, n, -1, -1 )
278 nb3 = ilaenv( 1, 'DORMQR', ' ', m, n, nrhs, -1 )
279 nb4 = ilaenv( 1, 'DORMRQ', ' ', m, n, nrhs, -1 )
280 nb = max( nb1, nb2, nb3, nb4 )
281 lwkmin = mn + max( 2*mn, n + 1, mn + nrhs )
282 lwkopt = max( lwkmin,
283 $ mn + 2*n + nb*( n + 1 ), 2*mn + nb*nrhs )
284 END IF
285 work( 1 ) = lwkopt
286*
287 IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
288 info = -12
289 END IF
290 END IF
291*
292 IF( info.NE.0 ) THEN
293 CALL xerbla( 'DGELSY', -info )
294 RETURN
295 ELSE IF( lquery ) THEN
296 RETURN
297 END IF
298*
299* Quick return if possible
300*
301 IF( mn.EQ.0 .OR. nrhs.EQ.0 ) THEN
302 rank = 0
303 RETURN
304 END IF
305*
306* Get machine parameters
307*
308 smlnum = dlamch( 'S' ) / dlamch( 'P' )
309 bignum = one / smlnum
310*
311* Scale A, B if max entries outside range [SMLNUM,BIGNUM]
312*
313 anrm = dlange( 'M', m, n, a, lda, work )
314 iascl = 0
315 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
316*
317* Scale matrix norm up to SMLNUM
318*
319 CALL dlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
320 iascl = 1
321 ELSE IF( anrm.GT.bignum ) THEN
322*
323* Scale matrix norm down to BIGNUM
324*
325 CALL dlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
326 iascl = 2
327 ELSE IF( anrm.EQ.zero ) THEN
328*
329* Matrix all zero. Return zero solution.
330*
331 CALL dlaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
332 rank = 0
333 GO TO 70
334 END IF
335*
336 bnrm = dlange( 'M', m, nrhs, b, ldb, work )
337 ibscl = 0
338 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
339*
340* Scale matrix norm up to SMLNUM
341*
342 CALL dlascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
343 ibscl = 1
344 ELSE IF( bnrm.GT.bignum ) THEN
345*
346* Scale matrix norm down to BIGNUM
347*
348 CALL dlascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
349 ibscl = 2
350 END IF
351*
352* Compute QR factorization with column pivoting of A:
353* A * P = Q * R
354*
355 CALL dgeqp3( m, n, a, lda, jpvt, work( 1 ), work( mn+1 ),
356 $ lwork-mn, info )
357 wsize = mn + work( mn+1 )
358*
359* workspace: MN+2*N+NB*(N+1).
360* Details of Householder rotations stored in WORK(1:MN).
361*
362* Determine RANK using incremental condition estimation
363*
364 work( ismin ) = one
365 work( ismax ) = one
366 smax = abs( a( 1, 1 ) )
367 smin = smax
368 IF( abs( a( 1, 1 ) ).EQ.zero ) THEN
369 rank = 0
370 CALL dlaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
371 GO TO 70
372 ELSE
373 rank = 1
374 END IF
375*
376 10 CONTINUE
377 IF( rank.LT.mn ) THEN
378 i = rank + 1
379 CALL dlaic1( imin, rank, work( ismin ), smin, a( 1, i ),
380 $ a( i, i ), sminpr, s1, c1 )
381 CALL dlaic1( imax, rank, work( ismax ), smax, a( 1, i ),
382 $ a( i, i ), smaxpr, s2, c2 )
383*
384 IF( smaxpr*rcond.LE.sminpr ) THEN
385 DO 20 i = 1, rank
386 work( ismin+i-1 ) = s1*work( ismin+i-1 )
387 work( ismax+i-1 ) = s2*work( ismax+i-1 )
388 20 CONTINUE
389 work( ismin+rank ) = c1
390 work( ismax+rank ) = c2
391 smin = sminpr
392 smax = smaxpr
393 rank = rank + 1
394 GO TO 10
395 END IF
396 END IF
397*
398* workspace: 3*MN.
399*
400* Logically partition R = [ R11 R12 ]
401* [ 0 R22 ]
402* where R11 = R(1:RANK,1:RANK)
403*
404* [R11,R12] = [ T11, 0 ] * Y
405*
406 IF( rank.LT.n )
407 $ CALL dtzrzf( rank, n, a, lda, work( mn+1 ), work( 2*mn+1 ),
408 $ lwork-2*mn, info )
409*
410* workspace: 2*MN.
411* Details of Householder rotations stored in WORK(MN+1:2*MN)
412*
413* B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
414*
415 CALL dormqr( 'Left', 'Transpose', m, nrhs, mn, a, lda, work( 1 ),
416 $ b, ldb, work( 2*mn+1 ), lwork-2*mn, info )
417 wsize = max( wsize, 2*mn+work( 2*mn+1 ) )
418*
419* workspace: 2*MN+NB*NRHS.
420*
421* B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
422*
423 CALL dtrsm( 'Left', 'Upper', 'No transpose', 'Non-unit', rank,
424 $ nrhs, one, a, lda, b, ldb )
425*
426 DO 40 j = 1, nrhs
427 DO 30 i = rank + 1, n
428 b( i, j ) = zero
429 30 CONTINUE
430 40 CONTINUE
431*
432* B(1:N,1:NRHS) := Y**T * B(1:N,1:NRHS)
433*
434 IF( rank.LT.n ) THEN
435 CALL dormrz( 'Left', 'Transpose', n, nrhs, rank, n-rank, a,
436 $ lda, work( mn+1 ), b, ldb, work( 2*mn+1 ),
437 $ lwork-2*mn, info )
438 END IF
439*
440* workspace: 2*MN+NRHS.
441*
442* B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
443*
444 DO 60 j = 1, nrhs
445 DO 50 i = 1, n
446 work( jpvt( i ) ) = b( i, j )
447 50 CONTINUE
448 CALL dcopy( n, work( 1 ), 1, b( 1, j ), 1 )
449 60 CONTINUE
450*
451* workspace: N.
452*
453* Undo scaling
454*
455 IF( iascl.EQ.1 ) THEN
456 CALL dlascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
457 CALL dlascl( 'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
458 $ info )
459 ELSE IF( iascl.EQ.2 ) THEN
460 CALL dlascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
461 CALL dlascl( 'U', 0, 0, bignum, anrm, rank, rank, a, lda,
462 $ info )
463 END IF
464 IF( ibscl.EQ.1 ) THEN
465 CALL dlascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
466 ELSE IF( ibscl.EQ.2 ) THEN
467 CALL dlascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
468 END IF
469*
470 70 CONTINUE
471 work( 1 ) = lwkopt
472*
473 RETURN
474*
475* End of DGELSY
476*
477 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
subroutine dgelsy(m, n, nrhs, a, lda, b, ldb, jpvt, rcond, rank, work, lwork, info)
DGELSY solves overdetermined or underdetermined systems for GE matrices
Definition dgelsy.f:206
subroutine dgeqp3(m, n, a, lda, jpvt, tau, work, lwork, info)
DGEQP3
Definition dgeqp3.f:151
subroutine dlaic1(job, j, x, sest, w, gamma, sestpr, s, c)
DLAIC1 applies one step of incremental condition estimation.
Definition dlaic1.f:134
subroutine dlascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition dlascl.f:143
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:110
subroutine dtrsm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
DTRSM
Definition dtrsm.f:181
subroutine dtzrzf(m, n, a, lda, tau, work, lwork, info)
DTZRZF
Definition dtzrzf.f:151
subroutine dormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
DORMQR
Definition dormqr.f:167
subroutine dormrz(side, trans, m, n, k, l, a, lda, tau, c, ldc, work, lwork, info)
DORMRZ
Definition dormrz.f:187