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