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