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
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