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