Annotation of rpl/lapack/lapack/dgelsy.f, revision 1.9
1.9 ! bertrand 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: *> \date November 2011
! 193: *
! 194: *> \ingroup doubleGEsolve
! 195: *
! 196: *> \par Contributors:
! 197: * ==================
! 198: *>
! 199: *> A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA \n
! 200: *> E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain \n
! 201: *> G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain \n
! 202: *>
! 203: * =====================================================================
1.1 bertrand 204: SUBROUTINE DGELSY( M, N, NRHS, A, LDA, B, LDB, JPVT, RCOND, RANK,
205: $ WORK, LWORK, INFO )
206: *
1.9 ! bertrand 207: * -- LAPACK driver routine (version 3.4.0) --
1.1 bertrand 208: * -- LAPACK is a software package provided by Univ. of Tennessee, --
209: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 210: * November 2011
1.1 bertrand 211: *
212: * .. Scalar Arguments ..
213: INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
214: DOUBLE PRECISION RCOND
215: * ..
216: * .. Array Arguments ..
217: INTEGER JPVT( * )
218: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
219: * ..
220: *
221: * =====================================================================
222: *
223: * .. Parameters ..
224: INTEGER IMAX, IMIN
225: PARAMETER ( IMAX = 1, IMIN = 2 )
226: DOUBLE PRECISION ZERO, ONE
227: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
228: * ..
229: * .. Local Scalars ..
230: LOGICAL LQUERY
231: INTEGER I, IASCL, IBSCL, ISMAX, ISMIN, J, LWKMIN,
232: $ LWKOPT, MN, NB, NB1, NB2, NB3, NB4
233: DOUBLE PRECISION ANRM, BIGNUM, BNRM, C1, C2, S1, S2, SMAX,
234: $ SMAXPR, SMIN, SMINPR, SMLNUM, WSIZE
235: * ..
236: * .. External Functions ..
237: INTEGER ILAENV
238: DOUBLE PRECISION DLAMCH, DLANGE
239: EXTERNAL ILAENV, DLAMCH, DLANGE
240: * ..
241: * .. External Subroutines ..
242: EXTERNAL DCOPY, DGEQP3, DLABAD, DLAIC1, DLASCL, DLASET,
243: $ DORMQR, DORMRZ, DTRSM, DTZRZF, XERBLA
244: * ..
245: * .. Intrinsic Functions ..
246: INTRINSIC ABS, MAX, MIN
247: * ..
248: * .. Executable Statements ..
249: *
250: MN = MIN( M, N )
251: ISMIN = MN + 1
252: ISMAX = 2*MN + 1
253: *
254: * Test the input arguments.
255: *
256: INFO = 0
257: LQUERY = ( LWORK.EQ.-1 )
258: IF( M.LT.0 ) THEN
259: INFO = -1
260: ELSE IF( N.LT.0 ) THEN
261: INFO = -2
262: ELSE IF( NRHS.LT.0 ) THEN
263: INFO = -3
264: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
265: INFO = -5
266: ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
267: INFO = -7
268: END IF
269: *
270: * Figure out optimal block size
271: *
272: IF( INFO.EQ.0 ) THEN
273: IF( MN.EQ.0 .OR. NRHS.EQ.0 ) THEN
274: LWKMIN = 1
275: LWKOPT = 1
276: ELSE
277: NB1 = ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
278: NB2 = ILAENV( 1, 'DGERQF', ' ', M, N, -1, -1 )
279: NB3 = ILAENV( 1, 'DORMQR', ' ', M, N, NRHS, -1 )
280: NB4 = ILAENV( 1, 'DORMRQ', ' ', M, N, NRHS, -1 )
281: NB = MAX( NB1, NB2, NB3, NB4 )
282: LWKMIN = MN + MAX( 2*MN, N + 1, MN + NRHS )
283: LWKOPT = MAX( LWKMIN,
284: $ MN + 2*N + NB*( N + 1 ), 2*MN + NB*NRHS )
285: END IF
286: WORK( 1 ) = LWKOPT
287: *
288: IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
289: INFO = -12
290: END IF
291: END IF
292: *
293: IF( INFO.NE.0 ) THEN
294: CALL XERBLA( 'DGELSY', -INFO )
295: RETURN
296: ELSE IF( LQUERY ) THEN
297: RETURN
298: END IF
299: *
300: * Quick return if possible
301: *
302: IF( MN.EQ.0 .OR. NRHS.EQ.0 ) THEN
303: RANK = 0
304: RETURN
305: END IF
306: *
307: * Get machine parameters
308: *
309: SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
310: BIGNUM = ONE / SMLNUM
311: CALL DLABAD( SMLNUM, BIGNUM )
312: *
313: * Scale A, B if max entries outside range [SMLNUM,BIGNUM]
314: *
315: ANRM = DLANGE( 'M', M, N, A, LDA, WORK )
316: IASCL = 0
317: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
318: *
319: * Scale matrix norm up to SMLNUM
320: *
321: CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
322: IASCL = 1
323: ELSE IF( ANRM.GT.BIGNUM ) THEN
324: *
325: * Scale matrix norm down to BIGNUM
326: *
327: CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
328: IASCL = 2
329: ELSE IF( ANRM.EQ.ZERO ) THEN
330: *
331: * Matrix all zero. Return zero solution.
332: *
333: CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
334: RANK = 0
335: GO TO 70
336: END IF
337: *
338: BNRM = DLANGE( 'M', M, NRHS, B, LDB, WORK )
339: IBSCL = 0
340: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
341: *
342: * Scale matrix norm up to SMLNUM
343: *
344: CALL DLASCL( 'G', 0, 0, BNRM, SMLNUM, M, NRHS, B, LDB, INFO )
345: IBSCL = 1
346: ELSE IF( BNRM.GT.BIGNUM ) THEN
347: *
348: * Scale matrix norm down to BIGNUM
349: *
350: CALL DLASCL( 'G', 0, 0, BNRM, BIGNUM, M, NRHS, B, LDB, INFO )
351: IBSCL = 2
352: END IF
353: *
354: * Compute QR factorization with column pivoting of A:
355: * A * P = Q * R
356: *
357: CALL DGEQP3( M, N, A, LDA, JPVT, WORK( 1 ), WORK( MN+1 ),
358: $ LWORK-MN, INFO )
359: WSIZE = MN + WORK( MN+1 )
360: *
361: * workspace: MN+2*N+NB*(N+1).
362: * Details of Householder rotations stored in WORK(1:MN).
363: *
364: * Determine RANK using incremental condition estimation
365: *
366: WORK( ISMIN ) = ONE
367: WORK( ISMAX ) = ONE
368: SMAX = ABS( A( 1, 1 ) )
369: SMIN = SMAX
370: IF( ABS( A( 1, 1 ) ).EQ.ZERO ) THEN
371: RANK = 0
372: CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
373: GO TO 70
374: ELSE
375: RANK = 1
376: END IF
377: *
378: 10 CONTINUE
379: IF( RANK.LT.MN ) THEN
380: I = RANK + 1
381: CALL DLAIC1( IMIN, RANK, WORK( ISMIN ), SMIN, A( 1, I ),
382: $ A( I, I ), SMINPR, S1, C1 )
383: CALL DLAIC1( IMAX, RANK, WORK( ISMAX ), SMAX, A( 1, I ),
384: $ A( I, I ), SMAXPR, S2, C2 )
385: *
386: IF( SMAXPR*RCOND.LE.SMINPR ) THEN
387: DO 20 I = 1, RANK
388: WORK( ISMIN+I-1 ) = S1*WORK( ISMIN+I-1 )
389: WORK( ISMAX+I-1 ) = S2*WORK( ISMAX+I-1 )
390: 20 CONTINUE
391: WORK( ISMIN+RANK ) = C1
392: WORK( ISMAX+RANK ) = C2
393: SMIN = SMINPR
394: SMAX = SMAXPR
395: RANK = RANK + 1
396: GO TO 10
397: END IF
398: END IF
399: *
400: * workspace: 3*MN.
401: *
402: * Logically partition R = [ R11 R12 ]
403: * [ 0 R22 ]
404: * where R11 = R(1:RANK,1:RANK)
405: *
406: * [R11,R12] = [ T11, 0 ] * Y
407: *
408: IF( RANK.LT.N )
409: $ CALL DTZRZF( RANK, N, A, LDA, WORK( MN+1 ), WORK( 2*MN+1 ),
410: $ LWORK-2*MN, INFO )
411: *
412: * workspace: 2*MN.
413: * Details of Householder rotations stored in WORK(MN+1:2*MN)
414: *
1.8 bertrand 415: * B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
1.1 bertrand 416: *
417: CALL DORMQR( 'Left', 'Transpose', M, NRHS, MN, A, LDA, WORK( 1 ),
418: $ B, LDB, WORK( 2*MN+1 ), LWORK-2*MN, INFO )
419: WSIZE = MAX( WSIZE, 2*MN+WORK( 2*MN+1 ) )
420: *
421: * workspace: 2*MN+NB*NRHS.
422: *
423: * B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS)
424: *
425: CALL DTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', RANK,
426: $ NRHS, ONE, A, LDA, B, LDB )
427: *
428: DO 40 J = 1, NRHS
429: DO 30 I = RANK + 1, N
430: B( I, J ) = ZERO
431: 30 CONTINUE
432: 40 CONTINUE
433: *
1.8 bertrand 434: * B(1:N,1:NRHS) := Y**T * B(1:N,1:NRHS)
1.1 bertrand 435: *
436: IF( RANK.LT.N ) THEN
437: CALL DORMRZ( 'Left', 'Transpose', N, NRHS, RANK, N-RANK, A,
438: $ LDA, WORK( MN+1 ), B, LDB, WORK( 2*MN+1 ),
439: $ LWORK-2*MN, INFO )
440: END IF
441: *
442: * workspace: 2*MN+NRHS.
443: *
444: * B(1:N,1:NRHS) := P * B(1:N,1:NRHS)
445: *
446: DO 60 J = 1, NRHS
447: DO 50 I = 1, N
448: WORK( JPVT( I ) ) = B( I, J )
449: 50 CONTINUE
450: CALL DCOPY( N, WORK( 1 ), 1, B( 1, J ), 1 )
451: 60 CONTINUE
452: *
453: * workspace: N.
454: *
455: * Undo scaling
456: *
457: IF( IASCL.EQ.1 ) THEN
458: CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, N, NRHS, B, LDB, INFO )
459: CALL DLASCL( 'U', 0, 0, SMLNUM, ANRM, RANK, RANK, A, LDA,
460: $ INFO )
461: ELSE IF( IASCL.EQ.2 ) THEN
462: CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, N, NRHS, B, LDB, INFO )
463: CALL DLASCL( 'U', 0, 0, BIGNUM, ANRM, RANK, RANK, A, LDA,
464: $ INFO )
465: END IF
466: IF( IBSCL.EQ.1 ) THEN
467: CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, N, NRHS, B, LDB, INFO )
468: ELSE IF( IBSCL.EQ.2 ) THEN
469: CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, N, NRHS, B, LDB, INFO )
470: END IF
471: *
472: 70 CONTINUE
473: WORK( 1 ) = LWKOPT
474: *
475: RETURN
476: *
477: * End of DGELSY
478: *
479: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgelsy.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack