Annotation of rpl/lapack/lapack/dgelsx.f, revision 1.9
1.9 ! bertrand 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: * =====================================================================
1.1 bertrand 178: SUBROUTINE DGELSX( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
179: $ WORK, INFO )
180: *
1.9 ! bertrand 181: * -- LAPACK driver routine (version 3.4.0) --
1.1 bertrand 182: * -- LAPACK is a software package provided by Univ. of Tennessee, --
183: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 184: * November 2011
1.1 bertrand 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: *
1.8 bertrand 355: * B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
1.1 bertrand 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: *
1.8 bertrand 373: * B(1:N,1:NRHS) := Y**T * B(1:N,1:NRHS)
1.1 bertrand 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
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