LAPACK  3.10.1 LAPACK: Linear Algebra PACKage
dgelsx.f
Go to the documentation of this file.
1 *> \brief <b> DGELSX 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
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgelsx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgelsx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgelsx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DGELSX( 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 * 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 *> This routine is deprecated and has been replaced by routine DGELSY.
40 *>
41 *> DGELSX 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 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
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 doubleGEsolve
174 *
175 * =====================================================================
176  SUBROUTINE dgelsx( 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  DOUBLE PRECISION RCOND
186 * ..
187 * .. Array Arguments ..
188  INTEGER JPVT( * )
189  DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
190 * ..
191 *
192 * =====================================================================
193 *
194 * .. Parameters ..
195  INTEGER IMAX, IMIN
196  parameter( imax = 1, imin = 2 )
197  DOUBLE PRECISION ZERO, ONE, DONE, NTDONE
198  parameter( zero = 0.0d0, one = 1.0d0, done = zero,
199  \$ ntdone = one )
200 * ..
201 * .. Local Scalars ..
202  INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, K, MN
203  DOUBLE PRECISION ANRM, BIGNUM, BNRM, C1, C2, S1, S2, SMAX,
204  \$ smaxpr, smin, sminpr, smlnum, t1, t2
205 * ..
206 * .. External Functions ..
207  DOUBLE PRECISION DLAMCH, DLANGE
208  EXTERNAL dlamch, dlange
209 * ..
210 * .. External Subroutines ..
211  EXTERNAL dgeqpf, dlaic1, dlascl, dlaset, dlatzm, dorm2r,
212  \$ dtrsm, dtzrqf, xerbla
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( 'DGELSX', -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 = dlamch( 'S' ) / dlamch( 'P' )
253  bignum = one / smlnum
254  CALL dlabad( smlnum, bignum )
255 *
256 * Scale A, B if max elements outside range [SMLNUM,BIGNUM]
257 *
258  anrm = dlange( '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 dlascl( '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 dlascl( '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 dlaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
277  rank = 0
278  GO TO 100
279  END IF
280 *
281  bnrm = dlange( '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 dlascl( '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 dlascl( '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 dgeqpf( 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 dlaset( '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 dlaic1( imin, rank, work( ismin ), smin, a( 1, i ),
323  \$ a( i, i ), sminpr, s1, c1 )
324  CALL dlaic1( 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 dtzrqf( 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 dorm2r( '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 dtrsm( '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 dlatzm( '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 dlascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
413  CALL dlascl( 'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
414  \$ info )
415  ELSE IF( iascl.EQ.2 ) THEN
416  CALL dlascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
417  CALL dlascl( 'U', 0, 0, bignum, anrm, rank, rank, a, lda,
418  \$ info )
419  END IF
420  IF( ibscl.EQ.1 ) THEN
421  CALL dlascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
422  ELSE IF( ibscl.EQ.2 ) THEN
423  CALL dlascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
424  END IF
425 *
426  100 CONTINUE
427 *
428  RETURN
429 *
430 * End of DGELSX
431 *
432  END
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 dtrsm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
DTRSM
Definition: dtrsm.f:181
subroutine dgeqpf(M, N, A, LDA, JPVT, TAU, WORK, INFO)
DGEQPF
Definition: dgeqpf.f:142
subroutine dgelsx(M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK, WORK, INFO)
DGELSX solves overdetermined or underdetermined systems for GE matrices
Definition: dgelsx.f:178
subroutine dlaic1(JOB, J, X, SEST, W, GAMMA, SESTPR, S, C)
DLAIC1 applies one step of incremental condition estimation.
Definition: dlaic1.f:134
subroutine dorm2r(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
DORM2R multiplies a general matrix by the orthogonal matrix from a QR factorization determined by sge...
Definition: dorm2r.f:159
subroutine dtzrqf(M, N, A, LDA, TAU, INFO)
DTZRQF
Definition: dtzrqf.f:138
subroutine dlatzm(SIDE, M, N, V, INCV, TAU, C1, C2, LDC, WORK)
DLATZM
Definition: dlatzm.f:151