Annotation of rpl/lapack/lapack/dgels.f, revision 1.9
1.9 ! bertrand 1: *> \brief <b> DGELS 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 DGELS + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgels.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgels.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgels.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DGELS( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK,
! 22: * INFO )
! 23: *
! 24: * .. Scalar Arguments ..
! 25: * CHARACTER TRANS
! 26: * INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS
! 27: * ..
! 28: * .. Array Arguments ..
! 29: * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
! 30: * ..
! 31: *
! 32: *
! 33: *> \par Purpose:
! 34: * =============
! 35: *>
! 36: *> \verbatim
! 37: *>
! 38: *> DGELS solves overdetermined or underdetermined real linear systems
! 39: *> involving an M-by-N matrix A, or its transpose, using a QR or LQ
! 40: *> factorization of A. It is assumed that A has full rank.
! 41: *>
! 42: *> The following options are provided:
! 43: *>
! 44: *> 1. If TRANS = 'N' and m >= n: find the least squares solution of
! 45: *> an overdetermined system, i.e., solve the least squares problem
! 46: *> minimize || B - A*X ||.
! 47: *>
! 48: *> 2. If TRANS = 'N' and m < n: find the minimum norm solution of
! 49: *> an underdetermined system A * X = B.
! 50: *>
! 51: *> 3. If TRANS = 'T' and m >= n: find the minimum norm solution of
! 52: *> an undetermined system A**T * X = B.
! 53: *>
! 54: *> 4. If TRANS = 'T' and m < n: find the least squares solution of
! 55: *> an overdetermined system, i.e., solve the least squares problem
! 56: *> minimize || B - A**T * X ||.
! 57: *>
! 58: *> Several right hand side vectors b and solution vectors x can be
! 59: *> handled in a single call; they are stored as the columns of the
! 60: *> M-by-NRHS right hand side matrix B and the N-by-NRHS solution
! 61: *> matrix X.
! 62: *> \endverbatim
! 63: *
! 64: * Arguments:
! 65: * ==========
! 66: *
! 67: *> \param[in] TRANS
! 68: *> \verbatim
! 69: *> TRANS is CHARACTER*1
! 70: *> = 'N': the linear system involves A;
! 71: *> = 'T': the linear system involves A**T.
! 72: *> \endverbatim
! 73: *>
! 74: *> \param[in] M
! 75: *> \verbatim
! 76: *> M is INTEGER
! 77: *> The number of rows of the matrix A. M >= 0.
! 78: *> \endverbatim
! 79: *>
! 80: *> \param[in] N
! 81: *> \verbatim
! 82: *> N is INTEGER
! 83: *> The number of columns of the matrix A. N >= 0.
! 84: *> \endverbatim
! 85: *>
! 86: *> \param[in] NRHS
! 87: *> \verbatim
! 88: *> NRHS is INTEGER
! 89: *> The number of right hand sides, i.e., the number of
! 90: *> columns of the matrices B and X. NRHS >=0.
! 91: *> \endverbatim
! 92: *>
! 93: *> \param[in,out] A
! 94: *> \verbatim
! 95: *> A is DOUBLE PRECISION array, dimension (LDA,N)
! 96: *> On entry, the M-by-N matrix A.
! 97: *> On exit,
! 98: *> if M >= N, A is overwritten by details of its QR
! 99: *> factorization as returned by DGEQRF;
! 100: *> if M < N, A is overwritten by details of its LQ
! 101: *> factorization as returned by DGELQF.
! 102: *> \endverbatim
! 103: *>
! 104: *> \param[in] LDA
! 105: *> \verbatim
! 106: *> LDA is INTEGER
! 107: *> The leading dimension of the array A. LDA >= max(1,M).
! 108: *> \endverbatim
! 109: *>
! 110: *> \param[in,out] B
! 111: *> \verbatim
! 112: *> B is DOUBLE PRECISION array, dimension (LDB,NRHS)
! 113: *> On entry, the matrix B of right hand side vectors, stored
! 114: *> columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS
! 115: *> if TRANS = 'T'.
! 116: *> On exit, if INFO = 0, B is overwritten by the solution
! 117: *> vectors, stored columnwise:
! 118: *> if TRANS = 'N' and m >= n, rows 1 to n of B contain the least
! 119: *> squares solution vectors; the residual sum of squares for the
! 120: *> solution in each column is given by the sum of squares of
! 121: *> elements N+1 to M in that column;
! 122: *> if TRANS = 'N' and m < n, rows 1 to N of B contain the
! 123: *> minimum norm solution vectors;
! 124: *> if TRANS = 'T' and m >= n, rows 1 to M of B contain the
! 125: *> minimum norm solution vectors;
! 126: *> if TRANS = 'T' and m < n, rows 1 to M of B contain the
! 127: *> least squares solution vectors; the residual sum of squares
! 128: *> for the solution in each column is given by the sum of
! 129: *> squares of elements M+1 to N in that column.
! 130: *> \endverbatim
! 131: *>
! 132: *> \param[in] LDB
! 133: *> \verbatim
! 134: *> LDB is INTEGER
! 135: *> The leading dimension of the array B. LDB >= MAX(1,M,N).
! 136: *> \endverbatim
! 137: *>
! 138: *> \param[out] WORK
! 139: *> \verbatim
! 140: *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
! 141: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
! 142: *> \endverbatim
! 143: *>
! 144: *> \param[in] LWORK
! 145: *> \verbatim
! 146: *> LWORK is INTEGER
! 147: *> The dimension of the array WORK.
! 148: *> LWORK >= max( 1, MN + max( MN, NRHS ) ).
! 149: *> For optimal performance,
! 150: *> LWORK >= max( 1, MN + max( MN, NRHS )*NB ).
! 151: *> where MN = min(M,N) and NB is the optimum block size.
! 152: *>
! 153: *> If LWORK = -1, then a workspace query is assumed; the routine
! 154: *> only calculates the optimal size of the WORK array, returns
! 155: *> this value as the first entry of the WORK array, and no error
! 156: *> message related to LWORK is issued by XERBLA.
! 157: *> \endverbatim
! 158: *>
! 159: *> \param[out] INFO
! 160: *> \verbatim
! 161: *> INFO is INTEGER
! 162: *> = 0: successful exit
! 163: *> < 0: if INFO = -i, the i-th argument had an illegal value
! 164: *> > 0: if INFO = i, the i-th diagonal element of the
! 165: *> triangular factor of A is zero, so that A does not have
! 166: *> full rank; the least squares solution could not be
! 167: *> computed.
! 168: *> \endverbatim
! 169: *
! 170: * Authors:
! 171: * ========
! 172: *
! 173: *> \author Univ. of Tennessee
! 174: *> \author Univ. of California Berkeley
! 175: *> \author Univ. of Colorado Denver
! 176: *> \author NAG Ltd.
! 177: *
! 178: *> \date November 2011
! 179: *
! 180: *> \ingroup doubleGEsolve
! 181: *
! 182: * =====================================================================
1.1 bertrand 183: SUBROUTINE DGELS( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK,
184: $ INFO )
185: *
1.9 ! bertrand 186: * -- LAPACK driver routine (version 3.4.0) --
1.1 bertrand 187: * -- LAPACK is a software package provided by Univ. of Tennessee, --
188: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 189: * November 2011
1.1 bertrand 190: *
191: * .. Scalar Arguments ..
192: CHARACTER TRANS
193: INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS
194: * ..
195: * .. Array Arguments ..
196: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
197: * ..
198: *
199: * =====================================================================
200: *
201: * .. Parameters ..
202: DOUBLE PRECISION ZERO, ONE
203: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
204: * ..
205: * .. Local Scalars ..
206: LOGICAL LQUERY, TPSD
207: INTEGER BROW, I, IASCL, IBSCL, J, MN, NB, SCLLEN, WSIZE
208: DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMLNUM
209: * ..
210: * .. Local Arrays ..
211: DOUBLE PRECISION RWORK( 1 )
212: * ..
213: * .. External Functions ..
214: LOGICAL LSAME
215: INTEGER ILAENV
216: DOUBLE PRECISION DLAMCH, DLANGE
217: EXTERNAL LSAME, ILAENV, DLABAD, DLAMCH, DLANGE
218: * ..
219: * .. External Subroutines ..
220: EXTERNAL DGELQF, DGEQRF, DLASCL, DLASET, DORMLQ, DORMQR,
221: $ DTRTRS, XERBLA
222: * ..
223: * .. Intrinsic Functions ..
224: INTRINSIC DBLE, MAX, MIN
225: * ..
226: * .. Executable Statements ..
227: *
228: * Test the input arguments.
229: *
230: INFO = 0
231: MN = MIN( M, N )
232: LQUERY = ( LWORK.EQ.-1 )
233: IF( .NOT.( LSAME( TRANS, 'N' ) .OR. LSAME( TRANS, 'T' ) ) ) THEN
234: INFO = -1
235: ELSE IF( M.LT.0 ) THEN
236: INFO = -2
237: ELSE IF( N.LT.0 ) THEN
238: INFO = -3
239: ELSE IF( NRHS.LT.0 ) THEN
240: INFO = -4
241: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
242: INFO = -6
243: ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
244: INFO = -8
245: ELSE IF( LWORK.LT.MAX( 1, MN+MAX( MN, NRHS ) ) .AND. .NOT.LQUERY )
246: $ THEN
247: INFO = -10
248: END IF
249: *
250: * Figure out optimal block size
251: *
252: IF( INFO.EQ.0 .OR. INFO.EQ.-10 ) THEN
253: *
254: TPSD = .TRUE.
255: IF( LSAME( TRANS, 'N' ) )
256: $ TPSD = .FALSE.
257: *
258: IF( M.GE.N ) THEN
259: NB = ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
260: IF( TPSD ) THEN
261: NB = MAX( NB, ILAENV( 1, 'DORMQR', 'LN', M, NRHS, N,
262: $ -1 ) )
263: ELSE
264: NB = MAX( NB, ILAENV( 1, 'DORMQR', 'LT', M, NRHS, N,
265: $ -1 ) )
266: END IF
267: ELSE
268: NB = ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
269: IF( TPSD ) THEN
270: NB = MAX( NB, ILAENV( 1, 'DORMLQ', 'LT', N, NRHS, M,
271: $ -1 ) )
272: ELSE
273: NB = MAX( NB, ILAENV( 1, 'DORMLQ', 'LN', N, NRHS, M,
274: $ -1 ) )
275: END IF
276: END IF
277: *
278: WSIZE = MAX( 1, MN+MAX( MN, NRHS )*NB )
279: WORK( 1 ) = DBLE( WSIZE )
280: *
281: END IF
282: *
283: IF( INFO.NE.0 ) THEN
284: CALL XERBLA( 'DGELS ', -INFO )
285: RETURN
286: ELSE IF( LQUERY ) THEN
287: RETURN
288: END IF
289: *
290: * Quick return if possible
291: *
292: IF( MIN( M, N, NRHS ).EQ.0 ) THEN
293: CALL DLASET( 'Full', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
294: RETURN
295: END IF
296: *
297: * Get machine parameters
298: *
299: SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
300: BIGNUM = ONE / SMLNUM
301: CALL DLABAD( SMLNUM, BIGNUM )
302: *
303: * Scale A, B if max element outside range [SMLNUM,BIGNUM]
304: *
305: ANRM = DLANGE( 'M', M, N, A, LDA, RWORK )
306: IASCL = 0
307: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
308: *
309: * Scale matrix norm up to SMLNUM
310: *
311: CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
312: IASCL = 1
313: ELSE IF( ANRM.GT.BIGNUM ) THEN
314: *
315: * Scale matrix norm down to BIGNUM
316: *
317: CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
318: IASCL = 2
319: ELSE IF( ANRM.EQ.ZERO ) THEN
320: *
321: * Matrix all zero. Return zero solution.
322: *
323: CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
324: GO TO 50
325: END IF
326: *
327: BROW = M
328: IF( TPSD )
329: $ BROW = N
330: BNRM = DLANGE( 'M', BROW, NRHS, B, LDB, RWORK )
331: IBSCL = 0
332: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
333: *
334: * Scale matrix norm up to SMLNUM
335: *
336: CALL DLASCL( 'G', 0, 0, BNRM, SMLNUM, BROW, NRHS, B, LDB,
337: $ INFO )
338: IBSCL = 1
339: ELSE IF( BNRM.GT.BIGNUM ) THEN
340: *
341: * Scale matrix norm down to BIGNUM
342: *
343: CALL DLASCL( 'G', 0, 0, BNRM, BIGNUM, BROW, NRHS, B, LDB,
344: $ INFO )
345: IBSCL = 2
346: END IF
347: *
348: IF( M.GE.N ) THEN
349: *
350: * compute QR factorization of A
351: *
352: CALL DGEQRF( M, N, A, LDA, WORK( 1 ), WORK( MN+1 ), LWORK-MN,
353: $ INFO )
354: *
355: * workspace at least N, optimally N*NB
356: *
357: IF( .NOT.TPSD ) THEN
358: *
359: * Least-Squares Problem min || A * X - B ||
360: *
1.8 bertrand 361: * B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
1.1 bertrand 362: *
363: CALL DORMQR( 'Left', 'Transpose', M, NRHS, N, A, LDA,
364: $ WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
365: $ INFO )
366: *
367: * workspace at least NRHS, optimally NRHS*NB
368: *
369: * B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
370: *
371: CALL DTRTRS( 'Upper', 'No transpose', 'Non-unit', N, NRHS,
372: $ A, LDA, B, LDB, INFO )
373: *
374: IF( INFO.GT.0 ) THEN
375: RETURN
376: END IF
377: *
378: SCLLEN = N
379: *
380: ELSE
381: *
1.8 bertrand 382: * Overdetermined system of equations A**T * X = B
1.1 bertrand 383: *
1.8 bertrand 384: * B(1:N,1:NRHS) := inv(R**T) * B(1:N,1:NRHS)
1.1 bertrand 385: *
386: CALL DTRTRS( 'Upper', 'Transpose', 'Non-unit', N, NRHS,
387: $ A, LDA, B, LDB, INFO )
388: *
389: IF( INFO.GT.0 ) THEN
390: RETURN
391: END IF
392: *
393: * B(N+1:M,1:NRHS) = ZERO
394: *
395: DO 20 J = 1, NRHS
396: DO 10 I = N + 1, M
397: B( I, J ) = ZERO
398: 10 CONTINUE
399: 20 CONTINUE
400: *
401: * B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS)
402: *
403: CALL DORMQR( 'Left', 'No transpose', M, NRHS, N, A, LDA,
404: $ WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
405: $ INFO )
406: *
407: * workspace at least NRHS, optimally NRHS*NB
408: *
409: SCLLEN = M
410: *
411: END IF
412: *
413: ELSE
414: *
415: * Compute LQ factorization of A
416: *
417: CALL DGELQF( M, N, A, LDA, WORK( 1 ), WORK( MN+1 ), LWORK-MN,
418: $ INFO )
419: *
420: * workspace at least M, optimally M*NB.
421: *
422: IF( .NOT.TPSD ) THEN
423: *
424: * underdetermined system of equations A * X = B
425: *
426: * B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
427: *
428: CALL DTRTRS( 'Lower', 'No transpose', 'Non-unit', M, NRHS,
429: $ A, LDA, B, LDB, INFO )
430: *
431: IF( INFO.GT.0 ) THEN
432: RETURN
433: END IF
434: *
435: * B(M+1:N,1:NRHS) = 0
436: *
437: DO 40 J = 1, NRHS
438: DO 30 I = M + 1, N
439: B( I, J ) = ZERO
440: 30 CONTINUE
441: 40 CONTINUE
442: *
1.8 bertrand 443: * B(1:N,1:NRHS) := Q(1:N,:)**T * B(1:M,1:NRHS)
1.1 bertrand 444: *
445: CALL DORMLQ( 'Left', 'Transpose', N, NRHS, M, A, LDA,
446: $ WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
447: $ INFO )
448: *
449: * workspace at least NRHS, optimally NRHS*NB
450: *
451: SCLLEN = N
452: *
453: ELSE
454: *
1.8 bertrand 455: * overdetermined system min || A**T * X - B ||
1.1 bertrand 456: *
457: * B(1:N,1:NRHS) := Q * B(1:N,1:NRHS)
458: *
459: CALL DORMLQ( 'Left', 'No transpose', N, NRHS, M, A, LDA,
460: $ WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
461: $ INFO )
462: *
463: * workspace at least NRHS, optimally NRHS*NB
464: *
1.8 bertrand 465: * B(1:M,1:NRHS) := inv(L**T) * B(1:M,1:NRHS)
1.1 bertrand 466: *
467: CALL DTRTRS( 'Lower', 'Transpose', 'Non-unit', M, NRHS,
468: $ A, LDA, B, LDB, INFO )
469: *
470: IF( INFO.GT.0 ) THEN
471: RETURN
472: END IF
473: *
474: SCLLEN = M
475: *
476: END IF
477: *
478: END IF
479: *
480: * Undo scaling
481: *
482: IF( IASCL.EQ.1 ) THEN
483: CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
484: $ INFO )
485: ELSE IF( IASCL.EQ.2 ) THEN
486: CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
487: $ INFO )
488: END IF
489: IF( IBSCL.EQ.1 ) THEN
490: CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
491: $ INFO )
492: ELSE IF( IBSCL.EQ.2 ) THEN
493: CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
494: $ INFO )
495: END IF
496: *
497: 50 CONTINUE
498: WORK( 1 ) = DBLE( WSIZE )
499: *
500: RETURN
501: *
502: * End of DGELS
503: *
504: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgels.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack