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