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