1: *> \brief <b> ZGELST solves overdetermined or underdetermined systems for GE matrices using QR or LQ factorization with compact WY representation of Q.</b>
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZGELST + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgelst.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgelst.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgelst.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZGELST( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK,
22: * INFO )
23: *
24: * .. Scalar Arguments ..
25: * CHARACTER TRANS
26: * INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS
27: * ..
28: * .. Array Arguments ..
29: * COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
30: * ..
31: *
32: *
33: *> \par Purpose:
34: * =============
35: *>
36: *> \verbatim
37: *>
38: *> ZGELST solves overdetermined or underdetermined real linear systems
39: *> involving an M-by-N matrix A, or its conjugate-transpose, using a QR
40: *> or LQ factorization of A with compact WY representation of Q.
41: *> It is assumed that A has full rank.
42: *>
43: *> The following options are provided:
44: *>
45: *> 1. If TRANS = 'N' and m >= n: find the least squares solution of
46: *> an overdetermined system, i.e., solve the least squares problem
47: *> minimize || B - A*X ||.
48: *>
49: *> 2. If TRANS = 'N' and m < n: find the minimum norm solution of
50: *> an underdetermined system A * X = B.
51: *>
52: *> 3. If TRANS = 'C' and m >= n: find the minimum norm solution of
53: *> an underdetermined system A**T * X = B.
54: *>
55: *> 4. If TRANS = 'C' and m < n: find the least squares solution of
56: *> an overdetermined system, i.e., solve the least squares problem
57: *> minimize || B - A**T * X ||.
58: *>
59: *> Several right hand side vectors b and solution vectors x can be
60: *> handled in a single call; they are stored as the columns of the
61: *> M-by-NRHS right hand side matrix B and the N-by-NRHS solution
62: *> matrix X.
63: *> \endverbatim
64: *
65: * Arguments:
66: * ==========
67: *
68: *> \param[in] TRANS
69: *> \verbatim
70: *> TRANS is CHARACTER*1
71: *> = 'N': the linear system involves A;
72: *> = 'C': the linear system involves A**H.
73: *> \endverbatim
74: *>
75: *> \param[in] M
76: *> \verbatim
77: *> M is INTEGER
78: *> The number of rows of the matrix A. M >= 0.
79: *> \endverbatim
80: *>
81: *> \param[in] N
82: *> \verbatim
83: *> N is INTEGER
84: *> The number of columns of the matrix A. N >= 0.
85: *> \endverbatim
86: *>
87: *> \param[in] NRHS
88: *> \verbatim
89: *> NRHS is INTEGER
90: *> The number of right hand sides, i.e., the number of
91: *> columns of the matrices B and X. NRHS >=0.
92: *> \endverbatim
93: *>
94: *> \param[in,out] A
95: *> \verbatim
96: *> A is COMPLEX*16 array, dimension (LDA,N)
97: *> On entry, the M-by-N matrix A.
98: *> On exit,
99: *> if M >= N, A is overwritten by details of its QR
100: *> factorization as returned by ZGEQRT;
101: *> if M < N, A is overwritten by details of its LQ
102: *> factorization as returned by ZGELQT.
103: *> \endverbatim
104: *>
105: *> \param[in] LDA
106: *> \verbatim
107: *> LDA is INTEGER
108: *> The leading dimension of the array A. LDA >= max(1,M).
109: *> \endverbatim
110: *>
111: *> \param[in,out] B
112: *> \verbatim
113: *> B is COMPLEX*16 array, dimension (LDB,NRHS)
114: *> On entry, the matrix B of right hand side vectors, stored
115: *> columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS
116: *> if TRANS = 'C'.
117: *> On exit, if INFO = 0, B is overwritten by the solution
118: *> vectors, stored columnwise:
119: *> if TRANS = 'N' and m >= n, rows 1 to n of B contain the least
120: *> squares solution vectors; the residual sum of squares for the
121: *> solution in each column is given by the sum of squares of
122: *> modulus of elements N+1 to M in that column;
123: *> if TRANS = 'N' and m < n, rows 1 to N of B contain the
124: *> minimum norm solution vectors;
125: *> if TRANS = 'C' and m >= n, rows 1 to M of B contain the
126: *> minimum norm solution vectors;
127: *> if TRANS = 'C' and m < n, rows 1 to M of B contain the
128: *> least squares solution vectors; the residual sum of squares
129: *> for the solution in each column is given by the sum of
130: *> squares of the modulus of elements M+1 to N in that column.
131: *> \endverbatim
132: *>
133: *> \param[in] LDB
134: *> \verbatim
135: *> LDB is INTEGER
136: *> The leading dimension of the array B. LDB >= MAX(1,M,N).
137: *> \endverbatim
138: *>
139: *> \param[out] WORK
140: *> \verbatim
141: *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
142: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
143: *> \endverbatim
144: *>
145: *> \param[in] LWORK
146: *> \verbatim
147: *> LWORK is INTEGER
148: *> The dimension of the array WORK.
149: *> LWORK >= max( 1, MN + max( MN, NRHS ) ).
150: *> For optimal performance,
151: *> LWORK >= max( 1, (MN + max( MN, NRHS ))*NB ).
152: *> where MN = min(M,N) and NB is the optimum block size.
153: *>
154: *> If LWORK = -1, then a workspace query is assumed; the routine
155: *> only calculates the optimal size of the WORK array, returns
156: *> this value as the first entry of the WORK array, and no error
157: *> message related to LWORK is issued by XERBLA.
158: *> \endverbatim
159: *>
160: *> \param[out] INFO
161: *> \verbatim
162: *> INFO is INTEGER
163: *> = 0: successful exit
164: *> < 0: if INFO = -i, the i-th argument had an illegal value
165: *> > 0: if INFO = i, the i-th diagonal element of the
166: *> triangular factor of A is zero, so that A does not have
167: *> full rank; the least squares solution could not be
168: *> computed.
169: *> \endverbatim
170: *
171: * Authors:
172: * ========
173: *
174: *> \author Univ. of Tennessee
175: *> \author Univ. of California Berkeley
176: *> \author Univ. of Colorado Denver
177: *> \author NAG Ltd.
178: *
179: *> \ingroup complex16GEsolve
180: *
181: *> \par Contributors:
182: * ==================
183: *>
184: *> \verbatim
185: *>
186: *> November 2022, Igor Kozachenko,
187: *> Computer Science Division,
188: *> University of California, Berkeley
189: *> \endverbatim
190: *
191: * =====================================================================
192: SUBROUTINE ZGELST( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK,
193: $ INFO )
194: *
195: * -- LAPACK driver routine --
196: * -- LAPACK is a software package provided by Univ. of Tennessee, --
197: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
198: *
199: * .. Scalar Arguments ..
200: CHARACTER TRANS
201: INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS
202: * ..
203: * .. Array Arguments ..
204: COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
205: * ..
206: *
207: * =====================================================================
208: *
209: * .. Parameters ..
210: DOUBLE PRECISION ZERO, ONE
211: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
212: COMPLEX*16 CZERO
213: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) )
214: * ..
215: * .. Local Scalars ..
216: LOGICAL LQUERY, TPSD
217: INTEGER BROW, I, IASCL, IBSCL, J, LWOPT, MN, MNNRHS,
218: $ NB, NBMIN, SCLLEN
219: DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMLNUM
220: * ..
221: * .. Local Arrays ..
222: DOUBLE PRECISION RWORK( 1 )
223: * ..
224: * .. External Functions ..
225: LOGICAL LSAME
226: INTEGER ILAENV
227: DOUBLE PRECISION DLAMCH, ZLANGE
228: EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE
229: * ..
230: * .. External Subroutines ..
231: EXTERNAL ZGELQT, ZGEQRT, ZGEMLQT, ZGEMQRT, DLABAD,
232: $ ZLASCL, ZLASET, ZTRTRS, XERBLA
233: * ..
234: * .. Intrinsic Functions ..
235: INTRINSIC DBLE, MAX, MIN
236: * ..
237: * .. Executable Statements ..
238: *
239: * Test the input arguments.
240: *
241: INFO = 0
242: MN = MIN( M, N )
243: LQUERY = ( LWORK.EQ.-1 )
244: IF( .NOT.( LSAME( TRANS, 'N' ) .OR. LSAME( TRANS, 'C' ) ) ) THEN
245: INFO = -1
246: ELSE IF( M.LT.0 ) THEN
247: INFO = -2
248: ELSE IF( N.LT.0 ) THEN
249: INFO = -3
250: ELSE IF( NRHS.LT.0 ) THEN
251: INFO = -4
252: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
253: INFO = -6
254: ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
255: INFO = -8
256: ELSE IF( LWORK.LT.MAX( 1, MN+MAX( MN, NRHS ) ) .AND. .NOT.LQUERY )
257: $ THEN
258: INFO = -10
259: END IF
260: *
261: * Figure out optimal block size and optimal workspace size
262: *
263: IF( INFO.EQ.0 .OR. INFO.EQ.-10 ) THEN
264: *
265: TPSD = .TRUE.
266: IF( LSAME( TRANS, 'N' ) )
267: $ TPSD = .FALSE.
268: *
269: NB = ILAENV( 1, 'ZGELST', ' ', M, N, -1, -1 )
270: *
271: MNNRHS = MAX( MN, NRHS )
272: LWOPT = MAX( 1, (MN+MNNRHS)*NB )
273: WORK( 1 ) = DBLE( LWOPT )
274: *
275: END IF
276: *
277: IF( INFO.NE.0 ) THEN
278: CALL XERBLA( 'ZGELST ', -INFO )
279: RETURN
280: ELSE IF( LQUERY ) THEN
281: RETURN
282: END IF
283: *
284: * Quick return if possible
285: *
286: IF( MIN( M, N, NRHS ).EQ.0 ) THEN
287: CALL ZLASET( 'Full', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
288: WORK( 1 ) = DBLE( LWOPT )
289: RETURN
290: END IF
291: *
292: * *GEQRT and *GELQT routines cannot accept NB larger than min(M,N)
293: *
294: IF( NB.GT.MN ) NB = MN
295: *
296: * Determine the block size from the supplied LWORK
297: * ( at this stage we know that LWORK >= (minimum required workspace,
298: * but it may be less than optimal)
299: *
300: NB = MIN( NB, LWORK/( MN + MNNRHS ) )
301: *
302: * The minimum value of NB, when blocked code is used
303: *
304: NBMIN = MAX( 2, ILAENV( 2, 'ZGELST', ' ', M, N, -1, -1 ) )
305: *
306: IF( NB.LT.NBMIN ) THEN
307: NB = 1
308: END IF
309: *
310: * Get machine parameters
311: *
312: SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
313: BIGNUM = ONE / SMLNUM
314: CALL DLABAD( SMLNUM, BIGNUM )
315: *
316: * Scale A, B if max element outside range [SMLNUM,BIGNUM]
317: *
318: ANRM = ZLANGE( 'M', M, N, A, LDA, RWORK )
319: IASCL = 0
320: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
321: *
322: * Scale matrix norm up to SMLNUM
323: *
324: CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
325: IASCL = 1
326: ELSE IF( ANRM.GT.BIGNUM ) THEN
327: *
328: * Scale matrix norm down to BIGNUM
329: *
330: CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
331: IASCL = 2
332: ELSE IF( ANRM.EQ.ZERO ) THEN
333: *
334: * Matrix all zero. Return zero solution.
335: *
336: CALL ZLASET( 'Full', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
337: WORK( 1 ) = DBLE( LWOPT )
338: RETURN
339: END IF
340: *
341: BROW = M
342: IF( TPSD )
343: $ BROW = N
344: BNRM = ZLANGE( 'M', BROW, NRHS, B, LDB, RWORK )
345: IBSCL = 0
346: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
347: *
348: * Scale matrix norm up to SMLNUM
349: *
350: CALL ZLASCL( 'G', 0, 0, BNRM, SMLNUM, BROW, NRHS, B, LDB,
351: $ INFO )
352: IBSCL = 1
353: ELSE IF( BNRM.GT.BIGNUM ) THEN
354: *
355: * Scale matrix norm down to BIGNUM
356: *
357: CALL ZLASCL( 'G', 0, 0, BNRM, BIGNUM, BROW, NRHS, B, LDB,
358: $ INFO )
359: IBSCL = 2
360: END IF
361: *
362: IF( M.GE.N ) THEN
363: *
364: * M > N:
365: * Compute the blocked QR factorization of A,
366: * using the compact WY representation of Q,
367: * workspace at least N, optimally N*NB.
368: *
369: CALL ZGEQRT( M, N, NB, A, LDA, WORK( 1 ), NB,
370: $ WORK( MN*NB+1 ), INFO )
371: *
372: IF( .NOT.TPSD ) THEN
373: *
374: * M > N, A is not transposed:
375: * Overdetermined system of equations,
376: * least-squares problem, min || A * X - B ||.
377: *
378: * Compute B(1:M,1:NRHS) := Q**T * B(1:M,1:NRHS),
379: * using the compact WY representation of Q,
380: * workspace at least NRHS, optimally NRHS*NB.
381: *
382: CALL ZGEMQRT( 'Left', 'Conjugate transpose', M, NRHS, N, NB,
383: $ A, LDA, WORK( 1 ), NB, B, LDB,
384: $ WORK( MN*NB+1 ), INFO )
385: *
386: * Compute B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
387: *
388: CALL ZTRTRS( 'Upper', 'No transpose', 'Non-unit', N, NRHS,
389: $ A, LDA, B, LDB, INFO )
390: *
391: IF( INFO.GT.0 ) THEN
392: RETURN
393: END IF
394: *
395: SCLLEN = N
396: *
397: ELSE
398: *
399: * M > N, A is transposed:
400: * Underdetermined system of equations,
401: * minimum norm solution of A**T * X = B.
402: *
403: * Compute B := inv(R**T) * B in two row blocks of B.
404: *
405: * Block 1: B(1:N,1:NRHS) := inv(R**T) * B(1:N,1:NRHS)
406: *
407: CALL ZTRTRS( 'Upper', 'Conjugate transpose', 'Non-unit',
408: $ N, NRHS, A, LDA, B, LDB, INFO )
409: *
410: IF( INFO.GT.0 ) THEN
411: RETURN
412: END IF
413: *
414: * Block 2: Zero out all rows below the N-th row in B:
415: * B(N+1:M,1:NRHS) = ZERO
416: *
417: DO J = 1, NRHS
418: DO I = N + 1, M
419: B( I, J ) = ZERO
420: END DO
421: END DO
422: *
423: * Compute B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS),
424: * using the compact WY representation of Q,
425: * workspace at least NRHS, optimally NRHS*NB.
426: *
427: CALL ZGEMQRT( 'Left', 'No transpose', M, NRHS, N, NB,
428: $ A, LDA, WORK( 1 ), NB, B, LDB,
429: $ WORK( MN*NB+1 ), INFO )
430: *
431: SCLLEN = M
432: *
433: END IF
434: *
435: ELSE
436: *
437: * M < N:
438: * Compute the blocked LQ factorization of A,
439: * using the compact WY representation of Q,
440: * workspace at least M, optimally M*NB.
441: *
442: CALL ZGELQT( M, N, NB, A, LDA, WORK( 1 ), NB,
443: $ WORK( MN*NB+1 ), INFO )
444: *
445: IF( .NOT.TPSD ) THEN
446: *
447: * M < N, A is not transposed:
448: * Underdetermined system of equations,
449: * minimum norm solution of A * X = B.
450: *
451: * Compute B := inv(L) * B in two row blocks of B.
452: *
453: * Block 1: B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
454: *
455: CALL ZTRTRS( 'Lower', 'No transpose', 'Non-unit', M, NRHS,
456: $ A, LDA, B, LDB, INFO )
457: *
458: IF( INFO.GT.0 ) THEN
459: RETURN
460: END IF
461: *
462: * Block 2: Zero out all rows below the M-th row in B:
463: * B(M+1:N,1:NRHS) = ZERO
464: *
465: DO J = 1, NRHS
466: DO I = M + 1, N
467: B( I, J ) = ZERO
468: END DO
469: END DO
470: *
471: * Compute B(1:N,1:NRHS) := Q(1:N,:)**T * B(1:M,1:NRHS),
472: * using the compact WY representation of Q,
473: * workspace at least NRHS, optimally NRHS*NB.
474: *
475: CALL ZGEMLQT( 'Left', 'Conjugate transpose', N, NRHS, M, NB,
476: $ A, LDA, WORK( 1 ), NB, B, LDB,
477: $ WORK( MN*NB+1 ), INFO )
478: *
479: SCLLEN = N
480: *
481: ELSE
482: *
483: * M < N, A is transposed:
484: * Overdetermined system of equations,
485: * least-squares problem, min || A**T * X - B ||.
486: *
487: * Compute B(1:N,1:NRHS) := Q * B(1:N,1:NRHS),
488: * using the compact WY representation of Q,
489: * workspace at least NRHS, optimally NRHS*NB.
490: *
491: CALL ZGEMLQT( 'Left', 'No transpose', N, NRHS, M, NB,
492: $ A, LDA, WORK( 1 ), NB, B, LDB,
493: $ WORK( MN*NB+1), INFO )
494: *
495: * Compute B(1:M,1:NRHS) := inv(L**T) * B(1:M,1:NRHS)
496: *
497: CALL ZTRTRS( 'Lower', 'Conjugate transpose', 'Non-unit',
498: $ M, NRHS, A, LDA, B, LDB, INFO )
499: *
500: IF( INFO.GT.0 ) THEN
501: RETURN
502: END IF
503: *
504: SCLLEN = M
505: *
506: END IF
507: *
508: END IF
509: *
510: * Undo scaling
511: *
512: IF( IASCL.EQ.1 ) THEN
513: CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
514: $ INFO )
515: ELSE IF( IASCL.EQ.2 ) THEN
516: CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
517: $ INFO )
518: END IF
519: IF( IBSCL.EQ.1 ) THEN
520: CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
521: $ INFO )
522: ELSE IF( IBSCL.EQ.2 ) THEN
523: CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
524: $ INFO )
525: END IF
526: *
527: WORK( 1 ) = DBLE( LWOPT )
528: *
529: RETURN
530: *
531: * End of ZGELST
532: *
533: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgelst.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack