LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
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
9 *> Download DGELSX + dependencies
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 *> \date November 2011
174 *
175 *> \ingroup doubleGEsolve
176 *
177 * =====================================================================
178  SUBROUTINE dgelsx( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
179  $ work, info )
180 *
181 * -- LAPACK driver routine (version 3.4.0) --
182 * -- LAPACK is a software package provided by Univ. of Tennessee, --
183 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
184 * November 2011
185 *
186 * .. Scalar Arguments ..
187  INTEGER info, lda, ldb, m, n, nrhs, rank
188  DOUBLE PRECISION rcond
189 * ..
190 * .. Array Arguments ..
191  INTEGER jpvt( * )
192  DOUBLE PRECISION a( lda, * ), b( ldb, * ), work( * )
193 * ..
194 *
195 * =====================================================================
196 *
197 * .. Parameters ..
198  INTEGER imax, imin
199  parameter( imax = 1, imin = 2 )
200  DOUBLE PRECISION zero, one, done, ntdone
201  parameter( zero = 0.0d0, one = 1.0d0, done = zero,
202  $ ntdone = one )
203 * ..
204 * .. Local Scalars ..
205  INTEGER i, iascl, ibscl, ismax, ismin, j, k, mn
206  DOUBLE PRECISION anrm, bignum, bnrm, c1, c2, s1, s2, smax,
207  $ smaxpr, smin, sminpr, smlnum, t1, t2
208 * ..
209 * .. External Functions ..
210  DOUBLE PRECISION dlamch, dlange
211  EXTERNAL dlamch, dlange
212 * ..
213 * .. External Subroutines ..
214  EXTERNAL dgeqpf, dlaic1, dlascl, dlaset, dlatzm, dorm2r,
215  $ dtrsm, dtzrqf, xerbla
216 * ..
217 * .. Intrinsic Functions ..
218  INTRINSIC abs, max, min
219 * ..
220 * .. Executable Statements ..
221 *
222  mn = min( m, n )
223  ismin = mn + 1
224  ismax = 2*mn + 1
225 *
226 * Test the input arguments.
227 *
228  info = 0
229  IF( m.LT.0 ) THEN
230  info = -1
231  ELSE IF( n.LT.0 ) THEN
232  info = -2
233  ELSE IF( nrhs.LT.0 ) THEN
234  info = -3
235  ELSE IF( lda.LT.max( 1, m ) ) THEN
236  info = -5
237  ELSE IF( ldb.LT.max( 1, m, n ) ) THEN
238  info = -7
239  END IF
240 *
241  IF( info.NE.0 ) THEN
242  CALL xerbla( 'DGELSX', -info )
243  RETURN
244  END IF
245 *
246 * Quick return if possible
247 *
248  IF( min( m, n, nrhs ).EQ.0 ) THEN
249  rank = 0
250  RETURN
251  END IF
252 *
253 * Get machine parameters
254 *
255  smlnum = dlamch( 'S' ) / dlamch( 'P' )
256  bignum = one / smlnum
257  CALL dlabad( smlnum, bignum )
258 *
259 * Scale A, B if max elements outside range [SMLNUM,BIGNUM]
260 *
261  anrm = dlange( 'M', m, n, a, lda, work )
262  iascl = 0
263  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
264 *
265 * Scale matrix norm up to SMLNUM
266 *
267  CALL dlascl( 'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
268  iascl = 1
269  ELSE IF( anrm.GT.bignum ) THEN
270 *
271 * Scale matrix norm down to BIGNUM
272 *
273  CALL dlascl( 'G', 0, 0, anrm, bignum, m, n, a, lda, info )
274  iascl = 2
275  ELSE IF( anrm.EQ.zero ) THEN
276 *
277 * Matrix all zero. Return zero solution.
278 *
279  CALL dlaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
280  rank = 0
281  go to 100
282  END IF
283 *
284  bnrm = dlange( 'M', m, nrhs, b, ldb, work )
285  ibscl = 0
286  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
287 *
288 * Scale matrix norm up to SMLNUM
289 *
290  CALL dlascl( 'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
291  ibscl = 1
292  ELSE IF( bnrm.GT.bignum ) THEN
293 *
294 * Scale matrix norm down to BIGNUM
295 *
296  CALL dlascl( 'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
297  ibscl = 2
298  END IF
299 *
300 * Compute QR factorization with column pivoting of A:
301 * A * P = Q * R
302 *
303  CALL dgeqpf( m, n, a, lda, jpvt, work( 1 ), work( mn+1 ), info )
304 *
305 * workspace 3*N. Details of Householder rotations stored
306 * in WORK(1:MN).
307 *
308 * Determine RANK using incremental condition estimation
309 *
310  work( ismin ) = one
311  work( ismax ) = one
312  smax = abs( a( 1, 1 ) )
313  smin = smax
314  IF( abs( a( 1, 1 ) ).EQ.zero ) THEN
315  rank = 0
316  CALL dlaset( 'F', max( m, n ), nrhs, zero, zero, b, ldb )
317  go to 100
318  ELSE
319  rank = 1
320  END IF
321 *
322  10 CONTINUE
323  IF( rank.LT.mn ) THEN
324  i = rank + 1
325  CALL dlaic1( imin, rank, work( ismin ), smin, a( 1, i ),
326  $ a( i, i ), sminpr, s1, c1 )
327  CALL dlaic1( imax, rank, work( ismax ), smax, a( 1, i ),
328  $ a( i, i ), smaxpr, s2, c2 )
329 *
330  IF( smaxpr*rcond.LE.sminpr ) THEN
331  DO 20 i = 1, rank
332  work( ismin+i-1 ) = s1*work( ismin+i-1 )
333  work( ismax+i-1 ) = s2*work( ismax+i-1 )
334  20 CONTINUE
335  work( ismin+rank ) = c1
336  work( ismax+rank ) = c2
337  smin = sminpr
338  smax = smaxpr
339  rank = rank + 1
340  go to 10
341  END IF
342  END IF
343 *
344 * Logically partition R = [ R11 R12 ]
345 * [ 0 R22 ]
346 * where R11 = R(1:RANK,1:RANK)
347 *
348 * [R11,R12] = [ T11, 0 ] * Y
349 *
350  IF( rank.LT.n )
351  $ CALL dtzrqf( rank, n, a, lda, work( mn+1 ), info )
352 *
353 * Details of Householder rotations stored in WORK(MN+1:2*MN)
354 *
355 * B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
356 *
357  CALL dorm2r( 'Left', 'Transpose', m, nrhs, mn, a, lda, work( 1 ),
358  $ b, ldb, work( 2*mn+1 ), info )
359 *
360 * workspace NRHS
361 *
362 * B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
363 *
364  CALL dtrsm( 'Left', 'Upper', 'No transpose', 'Non-unit', rank,
365  $ nrhs, one, a, lda, b, ldb )
366 *
367  DO 40 i = rank + 1, n
368  DO 30 j = 1, nrhs
369  b( i, j ) = zero
370  30 CONTINUE
371  40 CONTINUE
372 *
373 * B(1:N,1:NRHS) := Y**T * B(1:N,1:NRHS)
374 *
375  IF( rank.LT.n ) THEN
376  DO 50 i = 1, rank
377  CALL dlatzm( 'Left', n-rank+1, nrhs, a( i, rank+1 ), lda,
378  $ work( mn+i ), b( i, 1 ), b( rank+1, 1 ), ldb,
379  $ work( 2*mn+1 ) )
380  50 CONTINUE
381  END IF
382 *
383 * workspace NRHS
384 *
385 * B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
386 *
387  DO 90 j = 1, nrhs
388  DO 60 i = 1, n
389  work( 2*mn+i ) = ntdone
390  60 CONTINUE
391  DO 80 i = 1, n
392  IF( work( 2*mn+i ).EQ.ntdone ) THEN
393  IF( jpvt( i ).NE.i ) THEN
394  k = i
395  t1 = b( k, j )
396  t2 = b( jpvt( k ), j )
397  70 CONTINUE
398  b( jpvt( k ), j ) = t1
399  work( 2*mn+k ) = done
400  t1 = t2
401  k = jpvt( k )
402  t2 = b( jpvt( k ), j )
403  IF( jpvt( k ).NE.i )
404  $ go to 70
405  b( i, j ) = t1
406  work( 2*mn+k ) = done
407  END IF
408  END IF
409  80 CONTINUE
410  90 CONTINUE
411 *
412 * Undo scaling
413 *
414  IF( iascl.EQ.1 ) THEN
415  CALL dlascl( 'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
416  CALL dlascl( 'U', 0, 0, smlnum, anrm, rank, rank, a, lda,
417  $ info )
418  ELSE IF( iascl.EQ.2 ) THEN
419  CALL dlascl( 'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
420  CALL dlascl( 'U', 0, 0, bignum, anrm, rank, rank, a, lda,
421  $ info )
422  END IF
423  IF( ibscl.EQ.1 ) THEN
424  CALL dlascl( 'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
425  ELSE IF( ibscl.EQ.2 ) THEN
426  CALL dlascl( 'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
427  END IF
428 *
429  100 CONTINUE
430 *
431  RETURN
432 *
433 * End of DGELSX
434 *
435  END