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 underdetermined 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: *> \ingroup doubleGEsolve
179: *
180: * =====================================================================
181: SUBROUTINE DGELS( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK,
182: $ INFO )
183: *
184: * -- LAPACK driver routine --
185: * -- LAPACK is a software package provided by Univ. of Tennessee, --
186: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
187: *
188: * .. Scalar Arguments ..
189: CHARACTER TRANS
190: INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS
191: * ..
192: * .. Array Arguments ..
193: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), WORK( * )
194: * ..
195: *
196: * =====================================================================
197: *
198: * .. Parameters ..
199: DOUBLE PRECISION ZERO, ONE
200: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
201: * ..
202: * .. Local Scalars ..
203: LOGICAL LQUERY, TPSD
204: INTEGER BROW, I, IASCL, IBSCL, J, MN, NB, SCLLEN, WSIZE
205: DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMLNUM
206: * ..
207: * .. Local Arrays ..
208: DOUBLE PRECISION RWORK( 1 )
209: * ..
210: * .. External Functions ..
211: LOGICAL LSAME
212: INTEGER ILAENV
213: DOUBLE PRECISION DLAMCH, DLANGE
214: EXTERNAL LSAME, ILAENV, DLABAD, DLAMCH, DLANGE
215: * ..
216: * .. External Subroutines ..
217: EXTERNAL DGELQF, DGEQRF, DLASCL, DLASET, DORMLQ, DORMQR,
218: $ DTRTRS, XERBLA
219: * ..
220: * .. Intrinsic Functions ..
221: INTRINSIC DBLE, MAX, MIN
222: * ..
223: * .. Executable Statements ..
224: *
225: * Test the input arguments.
226: *
227: INFO = 0
228: MN = MIN( M, N )
229: LQUERY = ( LWORK.EQ.-1 )
230: IF( .NOT.( LSAME( TRANS, 'N' ) .OR. LSAME( TRANS, 'T' ) ) ) THEN
231: INFO = -1
232: ELSE IF( M.LT.0 ) THEN
233: INFO = -2
234: ELSE IF( N.LT.0 ) THEN
235: INFO = -3
236: ELSE IF( NRHS.LT.0 ) THEN
237: INFO = -4
238: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
239: INFO = -6
240: ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
241: INFO = -8
242: ELSE IF( LWORK.LT.MAX( 1, MN+MAX( MN, NRHS ) ) .AND. .NOT.LQUERY )
243: $ THEN
244: INFO = -10
245: END IF
246: *
247: * Figure out optimal block size
248: *
249: IF( INFO.EQ.0 .OR. INFO.EQ.-10 ) THEN
250: *
251: TPSD = .TRUE.
252: IF( LSAME( TRANS, 'N' ) )
253: $ TPSD = .FALSE.
254: *
255: IF( M.GE.N ) THEN
256: NB = ILAENV( 1, 'DGEQRF', ' ', M, N, -1, -1 )
257: IF( TPSD ) THEN
258: NB = MAX( NB, ILAENV( 1, 'DORMQR', 'LN', M, NRHS, N,
259: $ -1 ) )
260: ELSE
261: NB = MAX( NB, ILAENV( 1, 'DORMQR', 'LT', M, NRHS, N,
262: $ -1 ) )
263: END IF
264: ELSE
265: NB = ILAENV( 1, 'DGELQF', ' ', M, N, -1, -1 )
266: IF( TPSD ) THEN
267: NB = MAX( NB, ILAENV( 1, 'DORMLQ', 'LT', N, NRHS, M,
268: $ -1 ) )
269: ELSE
270: NB = MAX( NB, ILAENV( 1, 'DORMLQ', 'LN', N, NRHS, M,
271: $ -1 ) )
272: END IF
273: END IF
274: *
275: WSIZE = MAX( 1, MN+MAX( MN, NRHS )*NB )
276: WORK( 1 ) = DBLE( WSIZE )
277: *
278: END IF
279: *
280: IF( INFO.NE.0 ) THEN
281: CALL XERBLA( 'DGELS ', -INFO )
282: RETURN
283: ELSE IF( LQUERY ) THEN
284: RETURN
285: END IF
286: *
287: * Quick return if possible
288: *
289: IF( MIN( M, N, NRHS ).EQ.0 ) THEN
290: CALL DLASET( 'Full', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
291: RETURN
292: END IF
293: *
294: * Get machine parameters
295: *
296: SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
297: BIGNUM = ONE / SMLNUM
298: CALL DLABAD( SMLNUM, BIGNUM )
299: *
300: * Scale A, B if max element outside range [SMLNUM,BIGNUM]
301: *
302: ANRM = DLANGE( 'M', M, N, A, LDA, RWORK )
303: IASCL = 0
304: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
305: *
306: * Scale matrix norm up to SMLNUM
307: *
308: CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
309: IASCL = 1
310: ELSE IF( ANRM.GT.BIGNUM ) THEN
311: *
312: * Scale matrix norm down to BIGNUM
313: *
314: CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
315: IASCL = 2
316: ELSE IF( ANRM.EQ.ZERO ) THEN
317: *
318: * Matrix all zero. Return zero solution.
319: *
320: CALL DLASET( 'F', MAX( M, N ), NRHS, ZERO, ZERO, B, LDB )
321: GO TO 50
322: END IF
323: *
324: BROW = M
325: IF( TPSD )
326: $ BROW = N
327: BNRM = DLANGE( 'M', BROW, NRHS, B, LDB, RWORK )
328: IBSCL = 0
329: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
330: *
331: * Scale matrix norm up to SMLNUM
332: *
333: CALL DLASCL( 'G', 0, 0, BNRM, SMLNUM, BROW, NRHS, B, LDB,
334: $ INFO )
335: IBSCL = 1
336: ELSE IF( BNRM.GT.BIGNUM ) THEN
337: *
338: * Scale matrix norm down to BIGNUM
339: *
340: CALL DLASCL( 'G', 0, 0, BNRM, BIGNUM, BROW, NRHS, B, LDB,
341: $ INFO )
342: IBSCL = 2
343: END IF
344: *
345: IF( M.GE.N ) THEN
346: *
347: * compute QR factorization of A
348: *
349: CALL DGEQRF( M, N, A, LDA, WORK( 1 ), WORK( MN+1 ), LWORK-MN,
350: $ INFO )
351: *
352: * workspace at least N, optimally N*NB
353: *
354: IF( .NOT.TPSD ) THEN
355: *
356: * Least-Squares Problem min || A * X - B ||
357: *
358: * B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
359: *
360: CALL DORMQR( 'Left', 'Transpose', M, NRHS, N, A, LDA,
361: $ WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
362: $ INFO )
363: *
364: * workspace at least NRHS, optimally NRHS*NB
365: *
366: * B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
367: *
368: CALL DTRTRS( 'Upper', 'No transpose', 'Non-unit', N, NRHS,
369: $ A, LDA, B, LDB, INFO )
370: *
371: IF( INFO.GT.0 ) THEN
372: RETURN
373: END IF
374: *
375: SCLLEN = N
376: *
377: ELSE
378: *
379: * Underdetermined system of equations A**T * X = B
380: *
381: * B(1:N,1:NRHS) := inv(R**T) * B(1:N,1:NRHS)
382: *
383: CALL DTRTRS( 'Upper', 'Transpose', 'Non-unit', N, NRHS,
384: $ A, LDA, B, LDB, INFO )
385: *
386: IF( INFO.GT.0 ) THEN
387: RETURN
388: END IF
389: *
390: * B(N+1:M,1:NRHS) = ZERO
391: *
392: DO 20 J = 1, NRHS
393: DO 10 I = N + 1, M
394: B( I, J ) = ZERO
395: 10 CONTINUE
396: 20 CONTINUE
397: *
398: * B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS)
399: *
400: CALL DORMQR( 'Left', 'No transpose', M, NRHS, N, A, LDA,
401: $ WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
402: $ INFO )
403: *
404: * workspace at least NRHS, optimally NRHS*NB
405: *
406: SCLLEN = M
407: *
408: END IF
409: *
410: ELSE
411: *
412: * Compute LQ factorization of A
413: *
414: CALL DGELQF( M, N, A, LDA, WORK( 1 ), WORK( MN+1 ), LWORK-MN,
415: $ INFO )
416: *
417: * workspace at least M, optimally M*NB.
418: *
419: IF( .NOT.TPSD ) THEN
420: *
421: * underdetermined system of equations A * X = B
422: *
423: * B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
424: *
425: CALL DTRTRS( 'Lower', 'No transpose', 'Non-unit', M, NRHS,
426: $ A, LDA, B, LDB, INFO )
427: *
428: IF( INFO.GT.0 ) THEN
429: RETURN
430: END IF
431: *
432: * B(M+1:N,1:NRHS) = 0
433: *
434: DO 40 J = 1, NRHS
435: DO 30 I = M + 1, N
436: B( I, J ) = ZERO
437: 30 CONTINUE
438: 40 CONTINUE
439: *
440: * B(1:N,1:NRHS) := Q(1:N,:)**T * B(1:M,1:NRHS)
441: *
442: CALL DORMLQ( 'Left', 'Transpose', N, NRHS, M, A, LDA,
443: $ WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
444: $ INFO )
445: *
446: * workspace at least NRHS, optimally NRHS*NB
447: *
448: SCLLEN = N
449: *
450: ELSE
451: *
452: * overdetermined system min || A**T * X - B ||
453: *
454: * B(1:N,1:NRHS) := Q * B(1:N,1:NRHS)
455: *
456: CALL DORMLQ( 'Left', 'No transpose', N, NRHS, M, A, LDA,
457: $ WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
458: $ INFO )
459: *
460: * workspace at least NRHS, optimally NRHS*NB
461: *
462: * B(1:M,1:NRHS) := inv(L**T) * B(1:M,1:NRHS)
463: *
464: CALL DTRTRS( 'Lower', 'Transpose', 'Non-unit', M, NRHS,
465: $ A, LDA, B, LDB, INFO )
466: *
467: IF( INFO.GT.0 ) THEN
468: RETURN
469: END IF
470: *
471: SCLLEN = M
472: *
473: END IF
474: *
475: END IF
476: *
477: * Undo scaling
478: *
479: IF( IASCL.EQ.1 ) THEN
480: CALL DLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
481: $ INFO )
482: ELSE IF( IASCL.EQ.2 ) THEN
483: CALL DLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
484: $ INFO )
485: END IF
486: IF( IBSCL.EQ.1 ) THEN
487: CALL DLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
488: $ INFO )
489: ELSE IF( IBSCL.EQ.2 ) THEN
490: CALL DLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
491: $ INFO )
492: END IF
493: *
494: 50 CONTINUE
495: WORK( 1 ) = DBLE( WSIZE )
496: *
497: RETURN
498: *
499: * End of DGELS
500: *
501: 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