LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
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 *> \endverbatim
119 *>
120 *> \param[in] LDB
121 *> \verbatim
122 *> LDB is INTEGER
123 *> The leading dimension of the array B. LDB >= max(1,M,N).
124 *> \endverbatim
125 *>
126 *> \param[in,out] JPVT
127 *> \verbatim
128 *> JPVT is INTEGER array, dimension (N)
129 *> On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
130 *> to the front of AP, otherwise column i is a free column.
131 *> On exit, if JPVT(i) = k, then the i-th column of AP
132 *> was the k-th column of A.
133 *> \endverbatim
134 *>
135 *> \param[in] RCOND
136 *> \verbatim
137 *> RCOND is DOUBLE PRECISION
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 DOUBLE PRECISION array, dimension (MAX(1,LWORK))
155 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
156 *> \endverbatim
157 *>
158 *> \param[in] LWORK
159 *> \verbatim
160 *> LWORK is INTEGER
161 *> The dimension of the array WORK.
162 *> The unblocked strategy requires that:
163 *> LWORK >= MAX( MN+3*N+1, 2*MN+NRHS ),
164 *> where MN = min( M, N ).
165 *> The block algorithm requires that:
166 *> LWORK >= MAX( MN+2*N+NB*(N+1), 2*MN+NB*NRHS ),
167 *> where NB is an upper bound on the blocksize returned
168 *> by ILAENV for the routines DGEQP3, DTZRZF, STZRQF, DORMQR,
169 *> and DORMRZ.
170 *>
171 *> If LWORK = -1, then a workspace query is assumed; the routine
172 *> only calculates the optimal size of the WORK array, returns
173 *> this value as the first entry of the WORK array, and no error
174 *> message related to LWORK is issued by XERBLA.
175 *> \endverbatim
176 *>
177 *> \param[out] INFO
178 *> \verbatim
179 *> INFO is INTEGER
180 *> = 0: successful exit
181 *> < 0: If INFO = -i, the i-th argument had an illegal value.
182 *> \endverbatim
183 *
184 * Authors:
185 * ========
186 *
187 *> \author Univ. of Tennessee
188 *> \author Univ. of California Berkeley
189 *> \author Univ. of Colorado Denver
190 *> \author NAG Ltd.
191 *
192 *> \ingroup doubleGEsolve
193 *
194 *> \par Contributors:
195 * ==================
196 *>
197 *> A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA \n
198 *> E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain \n
199 *> G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain \n
200 *>
201 * =====================================================================
202  SUBROUTINE dgelsy( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
203  $ WORK, LWORK, INFO )
204 *
205 * -- LAPACK driver routine --
206 * -- LAPACK is a software package provided by Univ. of Tennessee, --
207 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
208 *
209 * .. Scalar Arguments ..
210  INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
211  DOUBLE PRECISION RCOND
212 * ..
213 * .. Array Arguments ..
214  INTEGER JPVT( * )
215  DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
216 * ..
217 *
218 * =====================================================================
219 *
220 * .. Parameters ..
221  INTEGER IMAX, IMIN
222  parameter( imax = 1, imin = 2 )
223  DOUBLE PRECISION ZERO, ONE
224  parameter( zero = 0.0d+0, one = 1.0d+0 )
225 * ..
226 * .. Local Scalars ..
227  LOGICAL LQUERY
228  INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, LWKMIN,
229  $ lwkopt, mn, nb, nb1, nb2, nb3, nb4
230  DOUBLE PRECISION ANRM, BIGNUM, BNRM, C1, C2, S1, S2, SMAX,
231  $ smaxpr, smin, sminpr, smlnum, wsize
232 * ..
233 * .. External Functions ..
234  INTEGER ILAENV
235  DOUBLE PRECISION DLAMCH, DLANGE
236  EXTERNAL ilaenv, dlamch, dlange
237 * ..
238 * .. External Subroutines ..
239  EXTERNAL dcopy, dgeqp3, dlabad, dlaic1, dlascl, dlaset,
241 * ..
242 * .. Intrinsic Functions ..
243  INTRINSIC abs, max, min
244 * ..
245 * .. Executable Statements ..
246 *
247  mn = min( m, n )
248  ismin = mn + 1
249  ismax = 2*mn + 1
250 *
251 * Test the input arguments.
252 *
253  info = 0
254  lquery = ( lwork.EQ.-1 )
255  IF( m.LT.0 ) THEN
256  info = -1
257  ELSE IF( n.LT.0 ) THEN
258  info = -2
259  ELSE IF( nrhs.LT.0 ) THEN
260  info = -3
261  ELSE IF( lda.LT.max( 1, m ) ) THEN
262  info = -5
263  ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
264  info = -7
265  END IF
266 *
267 * Figure out optimal block size
268 *
269  IF( info.EQ.0 ) THEN
270  IF( mn.EQ.0 .OR. nrhs.EQ.0 ) THEN
271  lwkmin = 1
272  lwkopt = 1
273  ELSE
274  nb1 = ilaenv( 1, 'DGEQRF', ' ', m, n, -1, -1 )
275  nb2 = ilaenv( 1, 'DGERQF', ' ', m, n, -1, -1 )
276  nb3 = ilaenv( 1, 'DORMQR', ' ', m, n, nrhs, -1 )
277  nb4 = ilaenv( 1, 'DORMRQ', ' ', m, n, nrhs, -1 )
278  nb = max( nb1, nb2, nb3, nb4 )
279  lwkmin = mn + max( 2*mn, n + 1, mn + nrhs )
280  lwkopt = max( lwkmin,
281  $ mn + 2*n + nb*( n + 1 ), 2*mn + nb*nrhs )
282  END IF
283  work( 1 ) = lwkopt
284 *
285  IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
286  info = -12
287  END IF
288  END IF
289 *
290  IF( info.NE.0 ) THEN
291  CALL xerbla( 'DGELSY', -info )
292  RETURN
293  ELSE IF( lquery ) THEN
294  RETURN
295  END IF
296 *
297 * Quick return if possible
298 *
299  IF( mn.EQ.0 .OR. nrhs.EQ.0 ) THEN
300  rank = 0
301  RETURN
302  END IF
303 *
304 * Get machine parameters
305 *
306  smlnum = dlamch( 'S' ) / dlamch( 'P' )
307  bignum = one / smlnum
308  CALL dlabad( smlnum, bignum )
309 *
310 * Scale A, B if max entries outside range [SMLNUM,BIGNUM]
311 *
312  anrm = dlange( 'M', m, n, a, lda, work )
313  iascl = 0
314  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
315 *
316 * Scale matrix norm up to SMLNUM
317 *
318  CALL dlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
319  iascl = 1
320  ELSE IF( anrm.GT.bignum ) THEN
321 *
322 * Scale matrix norm down to BIGNUM
323 *
324  CALL dlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
325  iascl = 2
326  ELSE IF( anrm.EQ.zero ) THEN
327 *
328 * Matrix all zero. Return zero solution.
329 *
330  CALL dlaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
331  rank = 0
332  GO TO 70
333  END IF
334 *
335  bnrm = dlange( 'M', m, nrhs, b, ldb, work )
336  ibscl = 0
337  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
338 *
339 * Scale matrix norm up to SMLNUM
340 *
341  CALL dlascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
342  ibscl = 1
343  ELSE IF( bnrm.GT.bignum ) THEN
344 *
345 * Scale matrix norm down to BIGNUM
346 *
347  CALL dlascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
348  ibscl = 2
349  END IF
350 *
351 * Compute QR factorization with column pivoting of A:
352 * A * P = Q * R
353 *
354  CALL dgeqp3( m, n, a, lda, jpvt, work( 1 ), work( mn+1 ),
355  $ lwork-mn, info )
356  wsize = mn + work( mn+1 )
357 *
358 * workspace: MN+2*N+NB*(N+1).
359 * Details of Householder rotations stored in WORK(1:MN).
360 *
361 * Determine RANK using incremental condition estimation
362 *
363  work( ismin ) = one
364  work( ismax ) = one
365  smax = abs( a( 1, 1 ) )
366  smin = smax
367  IF( abs( a( 1, 1 ) ).EQ.zero ) THEN
368  rank = 0
369  CALL dlaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
370  GO TO 70
371  ELSE
372  rank = 1
373  END IF
374 *
375  10 CONTINUE
376  IF( rank.LT.mn ) THEN
377  i = rank + 1
378  CALL dlaic1( imin, rank, work( ismin ), smin, a( 1, i ),
379  $ a( i, i ), sminpr, s1, c1 )
380  CALL dlaic1( imax, rank, work( ismax ), smax, a( 1, i ),
381  $ a( i, i ), smaxpr, s2, c2 )
382 *
383  IF( smaxpr*rcond.LE.sminpr ) THEN
384  DO 20 i = 1, rank
385  work( ismin+i-1 ) = s1*work( ismin+i-1 )
386  work( ismax+i-1 ) = s2*work( ismax+i-1 )
387  20 CONTINUE
388  work( ismin+rank ) = c1
389  work( ismax+rank ) = c2
390  smin = sminpr
391  smax = smaxpr
392  rank = rank + 1
393  GO TO 10
394  END IF
395  END IF
396 *
397 * workspace: 3*MN.
398 *
399 * Logically partition R = [ R11 R12 ]
400 * [ 0 R22 ]
401 * where R11 = R(1:RANK,1:RANK)
402 *
403 * [R11,R12] = [ T11, 0 ] * Y
404 *
405  IF( rank.LT.n )
406  $ CALL dtzrzf( rank, n, a, lda, work( mn+1 ), work( 2*mn+1 ),
407  $ lwork-2*mn, info )
408 *
409 * workspace: 2*MN.
410 * Details of Householder rotations stored in WORK(MN+1:2*MN)
411 *
412 * B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
413 *
414  CALL dormqr( 'Left', 'Transpose', m, nrhs, mn, a, lda, work( 1 ),
415  $ b, ldb, work( 2*mn+1 ), lwork-2*mn, info )
416  wsize = max( wsize, 2*mn+work( 2*mn+1 ) )
417 *
418 * workspace: 2*MN+NB*NRHS.
419 *
420 * B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
421 *
422  CALL dtrsm( 'Left', 'Upper', 'No transpose', 'Non-unit', rank,
423  $ nrhs, one, a, lda, b, ldb )
424 *
425  DO 40 j = 1, nrhs
426  DO 30 i = rank + 1, n
427  b( i, j ) = zero
428  30 CONTINUE
429  40 CONTINUE
430 *
431 * B(1:N,1:NRHS) := Y**T * B(1:N,1:NRHS)
432 *
433  IF( rank.LT.n ) THEN
434  CALL dormrz( 'Left', 'Transpose', n, nrhs, rank, n-rank, a,
435  $ lda, work( mn+1 ), b, ldb, work( 2*mn+1 ),
436  $ lwork-2*mn, info )
437  END IF
438 *
439 * workspace: 2*MN+NRHS.
440 *
441 * B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
442 *
443  DO 60 j = 1, nrhs
444  DO 50 i = 1, n
445  work( jpvt( i ) ) = b( i, j )
446  50 CONTINUE
447  CALL dcopy( n, work( 1 ), 1, b( 1, j ), 1 )
448  60 CONTINUE
449 *
450 * workspace: N.
451 *
452 * Undo scaling
453 *
454  IF( iascl.EQ.1 ) THEN
455  CALL dlascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
456  CALL dlascl( 'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
457  $ info )
458  ELSE IF( iascl.EQ.2 ) THEN
459  CALL dlascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
460  CALL dlascl( 'U', 0, 0, bignum, anrm, rank, rank, a, lda,
461  $ info )
462  END IF
463  IF( ibscl.EQ.1 ) THEN
464  CALL dlascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
465  ELSE IF( ibscl.EQ.2 ) THEN
466  CALL dlascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
467  END IF
468 *
469  70 CONTINUE
470  work( 1 ) = lwkopt
471 *
472  RETURN
473 *
474 * End of DGELSY
475 *
476  END
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
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 xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine dtrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
DTRSM
Definition: dtrsm.f:181
subroutine dgeqp3(M, N, A, LDA, JPVT, TAU, WORK, LWORK, INFO)
DGEQP3
Definition: dgeqp3.f:151
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:204
subroutine dlaic1(JOB, J, X, SEST, W, GAMMA, SESTPR, S, C)
DLAIC1 applies one step of incremental condition estimation.
Definition: dlaic1.f:134
subroutine dormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMQR
Definition: dormqr.f:167
subroutine dtzrzf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DTZRZF
Definition: dtzrzf.f:151
subroutine dormrz(SIDE, TRANS, M, N, K, L, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMRZ
Definition: dormrz.f:187