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: *> \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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgelsx.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack