Annotation of rpl/lapack/lapack/zgetsls.f, revision 1.4
1.1 bertrand 1: * Definition:
2: * ===========
3: *
4: * SUBROUTINE ZGETSLS( 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: * COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
13: * ..
14: *
15: *
16: *> \par Purpose:
17: * =============
18: *>
19: *> \verbatim
20: *>
21: *> ZGETSLS solves overdetermined or underdetermined complex 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 = 'C' and m >= n: find the minimum norm solution of
37: *> an undetermined system A**T * X = B.
38: *>
39: *> 4. If TRANS = 'C' 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;
1.3 bertrand 56: *> = 'C': the linear system involves A**H.
1.1 bertrand 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 COMPLEX*16 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 ZGEQR or ZGELQ.
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 COMPLEX*16 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 = 'C'.
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 = 'C' and m >= n, rows 1 to M of B contain the
106: *> minimum norm solution vectors;
107: *> if TRANS = 'C' 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) COMPLEX*16 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 complex16GEsolve
158: *
159: * =====================================================================
160: SUBROUTINE ZGETSLS( 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: COMPLEX*16 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: COMPLEX*16 CZERO
183: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) )
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
1.3 bertrand 190: DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMLNUM, DUM( 1 )
191: COMPLEX*16 TQ( 5 ), WORKQ( 1 )
1.1 bertrand 192: * ..
193: * .. External Functions ..
194: LOGICAL LSAME
195: INTEGER ILAENV
196: DOUBLE PRECISION DLAMCH, ZLANGE
197: EXTERNAL LSAME, ILAENV, DLABAD, DLAMCH, ZLANGE
198: * ..
199: * .. External Subroutines ..
200: EXTERNAL ZGEQR, ZGEMQR, ZLASCL, ZLASET,
201: $ ZTRTRS, XERBLA, ZGELQ, ZGEMLQ
202: * ..
203: * .. Intrinsic Functions ..
204: INTRINSIC DBLE, MAX, MIN, INT
205: * ..
206: * .. Executable Statements ..
207: *
208: * Test the input arguments.
209: *
210: INFO = 0
211: MINMN = MIN( M, N )
212: MAXMN = MAX( M, N )
213: MNK = MAX( MINMN, NRHS )
214: TRAN = LSAME( TRANS, 'C' )
215: *
216: LQUERY = ( LWORK.EQ.-1 .OR. LWORK.EQ.-2 )
217: IF( .NOT.( LSAME( TRANS, 'N' ) .OR.
218: $ LSAME( TRANS, 'C' ) ) ) THEN
219: INFO = -1
220: ELSE IF( M.LT.0 ) THEN
221: INFO = -2
222: ELSE IF( N.LT.0 ) THEN
223: INFO = -3
224: ELSE IF( NRHS.LT.0 ) THEN
225: INFO = -4
226: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
227: INFO = -6
228: ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
229: INFO = -8
230: END IF
231: *
232: IF( INFO.EQ.0 ) THEN
233: *
234: * Determine the block size and minimum LWORK
235: *
236: IF( M.GE.N ) THEN
237: CALL ZGEQR( M, N, A, LDA, TQ, -1, WORKQ, -1, INFO2 )
238: TSZO = INT( TQ( 1 ) )
1.3 bertrand 239: LWO = INT( WORKQ( 1 ) )
1.1 bertrand 240: CALL ZGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, TQ,
241: $ TSZO, B, LDB, WORKQ, -1, INFO2 )
1.3 bertrand 242: LWO = MAX( LWO, INT( WORKQ( 1 ) ) )
1.1 bertrand 243: CALL ZGEQR( M, N, A, LDA, TQ, -2, WORKQ, -2, INFO2 )
244: TSZM = INT( TQ( 1 ) )
1.3 bertrand 245: LWM = INT( WORKQ( 1 ) )
1.1 bertrand 246: CALL ZGEMQR( 'L', TRANS, M, NRHS, N, A, LDA, TQ,
247: $ TSZM, B, LDB, WORKQ, -1, INFO2 )
1.3 bertrand 248: LWM = MAX( LWM, INT( WORKQ( 1 ) ) )
1.1 bertrand 249: WSIZEO = TSZO + LWO
250: WSIZEM = TSZM + LWM
251: ELSE
252: CALL ZGELQ( M, N, A, LDA, TQ, -1, WORKQ, -1, INFO2 )
253: TSZO = INT( TQ( 1 ) )
1.3 bertrand 254: LWO = INT( WORKQ( 1 ) )
1.1 bertrand 255: CALL ZGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, TQ,
256: $ TSZO, B, LDB, WORKQ, -1, INFO2 )
1.3 bertrand 257: LWO = MAX( LWO, INT( WORKQ( 1 ) ) )
1.1 bertrand 258: CALL ZGELQ( M, N, A, LDA, TQ, -2, WORKQ, -2, INFO2 )
259: TSZM = INT( TQ( 1 ) )
1.3 bertrand 260: LWM = INT( WORKQ( 1 ) )
1.1 bertrand 261: CALL ZGEMLQ( 'L', TRANS, N, NRHS, M, A, LDA, TQ,
262: $ TSZO, B, LDB, WORKQ, -1, INFO2 )
1.3 bertrand 263: LWM = MAX( LWM, INT( WORKQ( 1 ) ) )
1.1 bertrand 264: WSIZEO = TSZO + LWO
265: WSIZEM = TSZM + LWM
266: END IF
267: *
268: IF( ( LWORK.LT.WSIZEM ).AND.( .NOT.LQUERY ) ) THEN
269: INFO = -10
270: END IF
271: *
272: END IF
273: *
274: IF( INFO.NE.0 ) THEN
275: CALL XERBLA( 'ZGETSLS', -INFO )
276: WORK( 1 ) = DBLE( WSIZEO )
277: RETURN
278: END IF
279: IF( LQUERY ) THEN
280: IF( LWORK.EQ.-1 ) WORK( 1 ) = REAL( WSIZEO )
281: IF( LWORK.EQ.-2 ) WORK( 1 ) = REAL( WSIZEM )
282: RETURN
283: END IF
284: IF( LWORK.LT.WSIZEO ) THEN
285: LW1 = TSZM
286: LW2 = LWM
287: ELSE
288: LW1 = TSZO
289: LW2 = LWO
290: END IF
291: *
292: * Quick return if possible
293: *
294: IF( MIN( M, N, NRHS ).EQ.0 ) THEN
295: CALL ZLASET( 'FULL', MAX( M, N ), NRHS, CZERO, CZERO,
296: $ B, LDB )
297: RETURN
298: END IF
299: *
300: * Get machine parameters
301: *
302: SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
303: BIGNUM = ONE / SMLNUM
304: CALL DLABAD( SMLNUM, BIGNUM )
305: *
306: * Scale A, B if max element outside range [SMLNUM,BIGNUM]
307: *
1.3 bertrand 308: ANRM = ZLANGE( 'M', M, N, A, LDA, DUM )
1.1 bertrand 309: IASCL = 0
310: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
311: *
312: * Scale matrix norm up to SMLNUM
313: *
314: CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
315: IASCL = 1
316: ELSE IF( ANRM.GT.BIGNUM ) THEN
317: *
318: * Scale matrix norm down to BIGNUM
319: *
320: CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
321: IASCL = 2
322: ELSE IF( ANRM.EQ.ZERO ) THEN
323: *
324: * Matrix all zero. Return zero solution.
325: *
326: CALL ZLASET( 'F', MAXMN, NRHS, CZERO, CZERO, B, LDB )
327: GO TO 50
328: END IF
329: *
330: BROW = M
331: IF ( TRAN ) THEN
332: BROW = N
333: END IF
1.3 bertrand 334: BNRM = ZLANGE( 'M', BROW, NRHS, B, LDB, DUM )
1.1 bertrand 335: IBSCL = 0
336: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
337: *
338: * Scale matrix norm up to SMLNUM
339: *
340: CALL ZLASCL( 'G', 0, 0, BNRM, SMLNUM, BROW, NRHS, B, LDB,
341: $ INFO )
342: IBSCL = 1
343: ELSE IF( BNRM.GT.BIGNUM ) THEN
344: *
345: * Scale matrix norm down to BIGNUM
346: *
347: CALL ZLASCL( 'G', 0, 0, BNRM, BIGNUM, BROW, NRHS, B, LDB,
348: $ INFO )
349: IBSCL = 2
350: END IF
351: *
352: IF ( M.GE.N ) THEN
353: *
354: * compute QR factorization of A
355: *
356: CALL ZGEQR( M, N, A, LDA, WORK( LW2+1 ), LW1,
357: $ WORK( 1 ), LW2, INFO )
358: IF ( .NOT.TRAN ) THEN
359: *
360: * Least-Squares Problem min || A * X - B ||
361: *
362: * B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS)
363: *
364: CALL ZGEMQR( 'L' , 'C', M, NRHS, N, A, LDA,
365: $ WORK( LW2+1 ), LW1, B, LDB, WORK( 1 ), LW2,
366: $ INFO )
367: *
368: * B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
369: *
370: CALL ZTRTRS( 'U', 'N', 'N', N, NRHS,
371: $ A, LDA, B, LDB, INFO )
372: IF( INFO.GT.0 ) THEN
373: RETURN
374: END IF
375: SCLLEN = N
376: ELSE
377: *
378: * Overdetermined system of equations A**T * X = B
379: *
380: * B(1:N,1:NRHS) := inv(R**T) * B(1:N,1:NRHS)
381: *
382: CALL ZTRTRS( 'U', 'C', 'N', N, NRHS,
383: $ A, LDA, B, LDB, INFO )
384: *
385: IF( INFO.GT.0 ) THEN
386: RETURN
387: END IF
388: *
389: * B(N+1:M,1:NRHS) = CZERO
390: *
391: DO 20 J = 1, NRHS
392: DO 10 I = N + 1, M
393: B( I, J ) = CZERO
394: 10 CONTINUE
395: 20 CONTINUE
396: *
397: * B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS)
398: *
399: CALL ZGEMQR( 'L', 'N', M, NRHS, N, A, LDA,
400: $ WORK( LW2+1 ), LW1, B, LDB, WORK( 1 ), LW2,
401: $ INFO )
402: *
403: SCLLEN = M
404: *
405: END IF
406: *
407: ELSE
408: *
409: * Compute LQ factorization of A
410: *
411: CALL ZGELQ( M, N, A, LDA, WORK( LW2+1 ), LW1,
412: $ WORK( 1 ), LW2, INFO )
413: *
414: * workspace at least M, optimally M*NB.
415: *
416: IF( .NOT.TRAN ) THEN
417: *
418: * underdetermined system of equations A * X = B
419: *
420: * B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
421: *
422: CALL ZTRTRS( 'L', 'N', 'N', M, NRHS,
423: $ A, LDA, B, LDB, INFO )
424: *
425: IF( INFO.GT.0 ) THEN
426: RETURN
427: END IF
428: *
429: * B(M+1:N,1:NRHS) = 0
430: *
431: DO 40 J = 1, NRHS
432: DO 30 I = M + 1, N
433: B( I, J ) = CZERO
434: 30 CONTINUE
435: 40 CONTINUE
436: *
437: * B(1:N,1:NRHS) := Q(1:N,:)**T * B(1:M,1:NRHS)
438: *
439: CALL ZGEMLQ( 'L', 'C', N, NRHS, M, A, LDA,
440: $ WORK( LW2+1 ), LW1, B, LDB, WORK( 1 ), LW2,
441: $ INFO )
442: *
443: * workspace at least NRHS, optimally NRHS*NB
444: *
445: SCLLEN = N
446: *
447: ELSE
448: *
449: * overdetermined system min || A**T * X - B ||
450: *
451: * B(1:N,1:NRHS) := Q * B(1:N,1:NRHS)
452: *
453: CALL ZGEMLQ( 'L', 'N', N, NRHS, M, A, LDA,
454: $ WORK( LW2+1 ), LW1, B, LDB, WORK( 1 ), LW2,
455: $ INFO )
456: *
457: * workspace at least NRHS, optimally NRHS*NB
458: *
459: * B(1:M,1:NRHS) := inv(L**T) * B(1:M,1:NRHS)
460: *
461: CALL ZTRTRS( 'L', 'C', 'N', M, NRHS,
462: $ A, LDA, B, LDB, INFO )
463: *
464: IF( INFO.GT.0 ) THEN
465: RETURN
466: END IF
467: *
468: SCLLEN = M
469: *
470: END IF
471: *
472: END IF
473: *
474: * Undo scaling
475: *
476: IF( IASCL.EQ.1 ) THEN
477: CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
478: $ INFO )
479: ELSE IF( IASCL.EQ.2 ) THEN
480: CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
481: $ INFO )
482: END IF
483: IF( IBSCL.EQ.1 ) THEN
484: CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
485: $ INFO )
486: ELSE IF( IBSCL.EQ.2 ) THEN
487: CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
488: $ INFO )
489: END IF
490: *
491: 50 CONTINUE
492: WORK( 1 ) = DBLE( TSZO + LWO )
493: RETURN
494: *
495: * End of ZGETSLS
496: *
497: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgetsls.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack