1: SUBROUTINE ZGELS( TRANS, M, N, NRHS, A, LDA, B, LDB, WORK, LWORK,
2: $ INFO )
3: *
4: * -- LAPACK driver routine (version 3.3.1) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * -- April 2011 --
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: * Purpose
18: * =======
19: *
20: * ZGELS solves overdetermined or underdetermined complex linear systems
21: * involving an M-by-N matrix A, or its conjugate-transpose, using a QR
22: * or LQ factorization of A. It is assumed that A has full rank.
23: *
24: * The following options are provided:
25: *
26: * 1. If TRANS = 'N' and m >= n: find the least squares solution of
27: * an overdetermined system, i.e., solve the least squares problem
28: * minimize || B - A*X ||.
29: *
30: * 2. If TRANS = 'N' and m < n: find the minimum norm solution of
31: * an underdetermined system A * X = B.
32: *
33: * 3. If TRANS = 'C' and m >= n: find the minimum norm solution of
34: * an undetermined system A**H * X = B.
35: *
36: * 4. If TRANS = 'C' and m < n: find the least squares solution of
37: * an overdetermined system, i.e., solve the least squares problem
38: * minimize || B - A**H * X ||.
39: *
40: * Several right hand side vectors b and solution vectors x can be
41: * handled in a single call; they are stored as the columns of the
42: * M-by-NRHS right hand side matrix B and the N-by-NRHS solution
43: * matrix X.
44: *
45: * Arguments
46: * =========
47: *
48: * TRANS (input) CHARACTER*1
49: * = 'N': the linear system involves A;
50: * = 'C': the linear system involves A**H.
51: *
52: * M (input) INTEGER
53: * The number of rows of the matrix A. M >= 0.
54: *
55: * N (input) INTEGER
56: * The number of columns of the matrix A. N >= 0.
57: *
58: * NRHS (input) INTEGER
59: * The number of right hand sides, i.e., the number of
60: * columns of the matrices B and X. NRHS >= 0.
61: *
62: * A (input/output) COMPLEX*16 array, dimension (LDA,N)
63: * On entry, the M-by-N matrix A.
64: * if M >= N, A is overwritten by details of its QR
65: * factorization as returned by ZGEQRF;
66: * if M < N, A is overwritten by details of its LQ
67: * factorization as returned by ZGELQF.
68: *
69: * LDA (input) INTEGER
70: * The leading dimension of the array A. LDA >= max(1,M).
71: *
72: * B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
73: * On entry, the matrix B of right hand side vectors, stored
74: * columnwise; B is M-by-NRHS if TRANS = 'N', or N-by-NRHS
75: * if TRANS = 'C'.
76: * On exit, if INFO = 0, B is overwritten by the solution
77: * vectors, stored columnwise:
78: * if TRANS = 'N' and m >= n, rows 1 to n of B contain the least
79: * squares solution vectors; the residual sum of squares for the
80: * solution in each column is given by the sum of squares of the
81: * modulus of elements N+1 to M in that column;
82: * if TRANS = 'N' and m < n, rows 1 to N of B contain the
83: * minimum norm solution vectors;
84: * if TRANS = 'C' and m >= n, rows 1 to M of B contain the
85: * minimum norm solution vectors;
86: * if TRANS = 'C' and m < n, rows 1 to M of B contain the
87: * least squares solution vectors; the residual sum of squares
88: * for the solution in each column is given by the sum of
89: * squares of the modulus of elements M+1 to N in that column.
90: *
91: * LDB (input) INTEGER
92: * The leading dimension of the array B. LDB >= MAX(1,M,N).
93: *
94: * WORK (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
95: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
96: *
97: * LWORK (input) INTEGER
98: * The dimension of the array WORK.
99: * LWORK >= max( 1, MN + max( MN, NRHS ) ).
100: * For optimal performance,
101: * LWORK >= max( 1, MN + max( MN, NRHS )*NB ).
102: * where MN = min(M,N) and NB is the optimum block size.
103: *
104: * If LWORK = -1, then a workspace query is assumed; the routine
105: * only calculates the optimal size of the WORK array, returns
106: * this value as the first entry of the WORK array, and no error
107: * message related to LWORK is issued by XERBLA.
108: *
109: * INFO (output) INTEGER
110: * = 0: successful exit
111: * < 0: if INFO = -i, the i-th argument had an illegal value
112: * > 0: if INFO = i, the i-th diagonal element of the
113: * triangular factor of A is zero, so that A does not have
114: * full rank; the least squares solution could not be
115: * computed.
116: *
117: * =====================================================================
118: *
119: * .. Parameters ..
120: DOUBLE PRECISION ZERO, ONE
121: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
122: COMPLEX*16 CZERO
123: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) )
124: * ..
125: * .. Local Scalars ..
126: LOGICAL LQUERY, TPSD
127: INTEGER BROW, I, IASCL, IBSCL, J, MN, NB, SCLLEN, WSIZE
128: DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMLNUM
129: * ..
130: * .. Local Arrays ..
131: DOUBLE PRECISION RWORK( 1 )
132: * ..
133: * .. External Functions ..
134: LOGICAL LSAME
135: INTEGER ILAENV
136: DOUBLE PRECISION DLAMCH, ZLANGE
137: EXTERNAL LSAME, ILAENV, DLAMCH, ZLANGE
138: * ..
139: * .. External Subroutines ..
140: EXTERNAL DLABAD, XERBLA, ZGELQF, ZGEQRF, ZLASCL, ZLASET,
141: $ ZTRTRS, ZUNMLQ, ZUNMQR
142: * ..
143: * .. Intrinsic Functions ..
144: INTRINSIC DBLE, MAX, MIN
145: * ..
146: * .. Executable Statements ..
147: *
148: * Test the input arguments.
149: *
150: INFO = 0
151: MN = MIN( M, N )
152: LQUERY = ( LWORK.EQ.-1 )
153: IF( .NOT.( LSAME( TRANS, 'N' ) .OR. LSAME( TRANS, 'C' ) ) ) THEN
154: INFO = -1
155: ELSE IF( M.LT.0 ) THEN
156: INFO = -2
157: ELSE IF( N.LT.0 ) THEN
158: INFO = -3
159: ELSE IF( NRHS.LT.0 ) THEN
160: INFO = -4
161: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
162: INFO = -6
163: ELSE IF( LDB.LT.MAX( 1, M, N ) ) THEN
164: INFO = -8
165: ELSE IF( LWORK.LT.MAX( 1, MN+MAX( MN, NRHS ) ) .AND. .NOT.LQUERY )
166: $ THEN
167: INFO = -10
168: END IF
169: *
170: * Figure out optimal block size
171: *
172: IF( INFO.EQ.0 .OR. INFO.EQ.-10 ) THEN
173: *
174: TPSD = .TRUE.
175: IF( LSAME( TRANS, 'N' ) )
176: $ TPSD = .FALSE.
177: *
178: IF( M.GE.N ) THEN
179: NB = ILAENV( 1, 'ZGEQRF', ' ', M, N, -1, -1 )
180: IF( TPSD ) THEN
181: NB = MAX( NB, ILAENV( 1, 'ZUNMQR', 'LN', M, NRHS, N,
182: $ -1 ) )
183: ELSE
184: NB = MAX( NB, ILAENV( 1, 'ZUNMQR', 'LC', M, NRHS, N,
185: $ -1 ) )
186: END IF
187: ELSE
188: NB = ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 )
189: IF( TPSD ) THEN
190: NB = MAX( NB, ILAENV( 1, 'ZUNMLQ', 'LC', N, NRHS, M,
191: $ -1 ) )
192: ELSE
193: NB = MAX( NB, ILAENV( 1, 'ZUNMLQ', 'LN', N, NRHS, M,
194: $ -1 ) )
195: END IF
196: END IF
197: *
198: WSIZE = MAX( 1, MN+MAX( MN, NRHS )*NB )
199: WORK( 1 ) = DBLE( WSIZE )
200: *
201: END IF
202: *
203: IF( INFO.NE.0 ) THEN
204: CALL XERBLA( 'ZGELS ', -INFO )
205: RETURN
206: ELSE IF( LQUERY ) THEN
207: RETURN
208: END IF
209: *
210: * Quick return if possible
211: *
212: IF( MIN( M, N, NRHS ).EQ.0 ) THEN
213: CALL ZLASET( 'Full', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
214: RETURN
215: END IF
216: *
217: * Get machine parameters
218: *
219: SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
220: BIGNUM = ONE / SMLNUM
221: CALL DLABAD( SMLNUM, BIGNUM )
222: *
223: * Scale A, B if max element outside range [SMLNUM,BIGNUM]
224: *
225: ANRM = ZLANGE( 'M', M, N, A, LDA, RWORK )
226: IASCL = 0
227: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
228: *
229: * Scale matrix norm up to SMLNUM
230: *
231: CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, M, N, A, LDA, INFO )
232: IASCL = 1
233: ELSE IF( ANRM.GT.BIGNUM ) THEN
234: *
235: * Scale matrix norm down to BIGNUM
236: *
237: CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, M, N, A, LDA, INFO )
238: IASCL = 2
239: ELSE IF( ANRM.EQ.ZERO ) THEN
240: *
241: * Matrix all zero. Return zero solution.
242: *
243: CALL ZLASET( 'F', MAX( M, N ), NRHS, CZERO, CZERO, B, LDB )
244: GO TO 50
245: END IF
246: *
247: BROW = M
248: IF( TPSD )
249: $ BROW = N
250: BNRM = ZLANGE( 'M', BROW, NRHS, B, LDB, RWORK )
251: IBSCL = 0
252: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
253: *
254: * Scale matrix norm up to SMLNUM
255: *
256: CALL ZLASCL( 'G', 0, 0, BNRM, SMLNUM, BROW, NRHS, B, LDB,
257: $ INFO )
258: IBSCL = 1
259: ELSE IF( BNRM.GT.BIGNUM ) THEN
260: *
261: * Scale matrix norm down to BIGNUM
262: *
263: CALL ZLASCL( 'G', 0, 0, BNRM, BIGNUM, BROW, NRHS, B, LDB,
264: $ INFO )
265: IBSCL = 2
266: END IF
267: *
268: IF( M.GE.N ) THEN
269: *
270: * compute QR factorization of A
271: *
272: CALL ZGEQRF( M, N, A, LDA, WORK( 1 ), WORK( MN+1 ), LWORK-MN,
273: $ INFO )
274: *
275: * workspace at least N, optimally N*NB
276: *
277: IF( .NOT.TPSD ) THEN
278: *
279: * Least-Squares Problem min || A * X - B ||
280: *
281: * B(1:M,1:NRHS) := Q**H * B(1:M,1:NRHS)
282: *
283: CALL ZUNMQR( 'Left', 'Conjugate transpose', M, NRHS, N, A,
284: $ LDA, WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
285: $ INFO )
286: *
287: * workspace at least NRHS, optimally NRHS*NB
288: *
289: * B(1:N,1:NRHS) := inv(R) * B(1:N,1:NRHS)
290: *
291: CALL ZTRTRS( 'Upper', 'No transpose', 'Non-unit', N, NRHS,
292: $ A, LDA, B, LDB, INFO )
293: *
294: IF( INFO.GT.0 ) THEN
295: RETURN
296: END IF
297: *
298: SCLLEN = N
299: *
300: ELSE
301: *
302: * Overdetermined system of equations A**H * X = B
303: *
304: * B(1:N,1:NRHS) := inv(R**H) * B(1:N,1:NRHS)
305: *
306: CALL ZTRTRS( 'Upper', 'Conjugate transpose','Non-unit',
307: $ N, NRHS, A, LDA, B, LDB, INFO )
308: *
309: IF( INFO.GT.0 ) THEN
310: RETURN
311: END IF
312: *
313: * B(N+1:M,1:NRHS) = ZERO
314: *
315: DO 20 J = 1, NRHS
316: DO 10 I = N + 1, M
317: B( I, J ) = CZERO
318: 10 CONTINUE
319: 20 CONTINUE
320: *
321: * B(1:M,1:NRHS) := Q(1:N,:) * B(1:N,1:NRHS)
322: *
323: CALL ZUNMQR( 'Left', 'No transpose', M, NRHS, N, A, LDA,
324: $ WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
325: $ INFO )
326: *
327: * workspace at least NRHS, optimally NRHS*NB
328: *
329: SCLLEN = M
330: *
331: END IF
332: *
333: ELSE
334: *
335: * Compute LQ factorization of A
336: *
337: CALL ZGELQF( M, N, A, LDA, WORK( 1 ), WORK( MN+1 ), LWORK-MN,
338: $ INFO )
339: *
340: * workspace at least M, optimally M*NB.
341: *
342: IF( .NOT.TPSD ) THEN
343: *
344: * underdetermined system of equations A * X = B
345: *
346: * B(1:M,1:NRHS) := inv(L) * B(1:M,1:NRHS)
347: *
348: CALL ZTRTRS( 'Lower', 'No transpose', 'Non-unit', M, NRHS,
349: $ A, LDA, B, LDB, INFO )
350: *
351: IF( INFO.GT.0 ) THEN
352: RETURN
353: END IF
354: *
355: * B(M+1:N,1:NRHS) = 0
356: *
357: DO 40 J = 1, NRHS
358: DO 30 I = M + 1, N
359: B( I, J ) = CZERO
360: 30 CONTINUE
361: 40 CONTINUE
362: *
363: * B(1:N,1:NRHS) := Q(1:N,:)**H * B(1:M,1:NRHS)
364: *
365: CALL ZUNMLQ( 'Left', 'Conjugate transpose', N, NRHS, M, A,
366: $ LDA, WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
367: $ INFO )
368: *
369: * workspace at least NRHS, optimally NRHS*NB
370: *
371: SCLLEN = N
372: *
373: ELSE
374: *
375: * overdetermined system min || A**H * X - B ||
376: *
377: * B(1:N,1:NRHS) := Q * B(1:N,1:NRHS)
378: *
379: CALL ZUNMLQ( 'Left', 'No transpose', N, NRHS, M, A, LDA,
380: $ WORK( 1 ), B, LDB, WORK( MN+1 ), LWORK-MN,
381: $ INFO )
382: *
383: * workspace at least NRHS, optimally NRHS*NB
384: *
385: * B(1:M,1:NRHS) := inv(L**H) * B(1:M,1:NRHS)
386: *
387: CALL ZTRTRS( 'Lower', 'Conjugate transpose', 'Non-unit',
388: $ M, NRHS, A, LDA, B, LDB, INFO )
389: *
390: IF( INFO.GT.0 ) THEN
391: RETURN
392: END IF
393: *
394: SCLLEN = M
395: *
396: END IF
397: *
398: END IF
399: *
400: * Undo scaling
401: *
402: IF( IASCL.EQ.1 ) THEN
403: CALL ZLASCL( 'G', 0, 0, ANRM, SMLNUM, SCLLEN, NRHS, B, LDB,
404: $ INFO )
405: ELSE IF( IASCL.EQ.2 ) THEN
406: CALL ZLASCL( 'G', 0, 0, ANRM, BIGNUM, SCLLEN, NRHS, B, LDB,
407: $ INFO )
408: END IF
409: IF( IBSCL.EQ.1 ) THEN
410: CALL ZLASCL( 'G', 0, 0, SMLNUM, BNRM, SCLLEN, NRHS, B, LDB,
411: $ INFO )
412: ELSE IF( IBSCL.EQ.2 ) THEN
413: CALL ZLASCL( 'G', 0, 0, BIGNUM, BNRM, SCLLEN, NRHS, B, LDB,
414: $ INFO )
415: END IF
416: *
417: 50 CONTINUE
418: WORK( 1 ) = DBLE( WSIZE )
419: *
420: RETURN
421: *
422: * End of ZGELS
423: *
424: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zgels.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack