162
163
164
165
166
167
168 CHARACTER TRANS
169 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS
170
171
172 COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
173
174
175
176
177
178
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
185 LOGICAL LQUERY, TRAN
186 INTEGER I, IASCL, IBSCL, J, MAXMN, BROW,
187 $ SCLLEN, TSZO, TSZM, LWO, LWM, LW1, LW2,
188 $ WSIZEO, WSIZEM, INFO2
189 DOUBLE PRECISION ANRM, BIGNUM, BNRM, SMLNUM, DUM( 1 )
190 COMPLEX*16 TQ( 5 ), WORKQ( 1 )
191
192
193 LOGICAL LSAME
194 DOUBLE PRECISION DLAMCH, ZLANGE
196
197
200
201
202 INTRINSIC dble, max, min, int
203
204
205
206
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
230
231
232 IF( m.GE.n ) THEN
233 CALL zgeqr( m, n, a, lda, tq, -1, workq, -1, info2 )
234 tszo = int( tq( 1 ) )
235 lwo = int( workq( 1 ) )
236 CALL zgemqr(
'L', trans, m, nrhs, n, a, lda, tq,
237 $ tszo, b, ldb, workq, -1, info2 )
238 lwo = max( lwo, int( workq( 1 ) ) )
239 CALL zgeqr( m, n, a, lda, tq, -2, workq, -2, info2 )
240 tszm = int( tq( 1 ) )
241 lwm = int( workq( 1 ) )
242 CALL zgemqr(
'L', trans, m, nrhs, n, a, lda, tq,
243 $ tszm, b, ldb, workq, -1, info2 )
244 lwm = max( lwm, int( workq( 1 ) ) )
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 ) )
250 lwo = int( workq( 1 ) )
251 CALL zgemlq(
'L', trans, n, nrhs, m, a, lda, tq,
252 $ tszo, b, ldb, workq, -1, info2 )
253 lwo = max( lwo, int( workq( 1 ) ) )
254 CALL zgelq( m, n, a, lda, tq, -2, workq, -2, info2 )
255 tszm = int( tq( 1 ) )
256 lwm = int( workq( 1 ) )
257 CALL zgemlq(
'L', trans, n, nrhs, m, a, lda, tq,
258 $ tszm, b, ldb, workq, -1, info2 )
259 lwm = max( lwm, int( workq( 1 ) ) )
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
268 work( 1 ) = dble( wsizeo )
269
270 END IF
271
272 IF( info.NE.0 ) THEN
273 CALL xerbla(
'ZGETSLS', -info )
274 RETURN
275 END IF
276 IF( lquery ) THEN
277 IF( lwork.EQ.-2 ) work( 1 ) = dble( wsizem )
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
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
297
299 bignum = one / smlnum
300
301
302
303 anrm =
zlange(
'M', m, n, a, lda, dum )
304 iascl = 0
305 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
306
307
308
309 CALL zlascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
310 iascl = 1
311 ELSE IF( anrm.GT.bignum ) THEN
312
313
314
315 CALL zlascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, info )
316 iascl = 2
317 ELSE IF( anrm.EQ.zero ) THEN
318
319
320
321 CALL zlaset(
'F', maxmn, nrhs, czero, czero, b, ldb )
322 GO TO 50
323 END IF
324
325 brow = m
326 IF ( tran ) THEN
327 brow = n
328 END IF
329 bnrm =
zlange(
'M', brow, nrhs, b, ldb, dum )
330 ibscl = 0
331 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
332
333
334
335 CALL zlascl(
'G', 0, 0, bnrm, smlnum, brow, nrhs, b, ldb,
336 $ info )
337 ibscl = 1
338 ELSE IF( bnrm.GT.bignum ) THEN
339
340
341
342 CALL zlascl(
'G', 0, 0, bnrm, bignum, brow, nrhs, b, ldb,
343 $ info )
344 ibscl = 2
345 END IF
346
347 IF ( m.GE.n ) THEN
348
349
350
351 CALL zgeqr( m, n, a, lda, work( lw2+1 ), lw1,
352 $ work( 1 ), lw2, info )
353 IF ( .NOT.tran ) THEN
354
355
356
357
358
359 CALL zgemqr(
'L' ,
'C', m, nrhs, n, a, lda,
360 $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
361 $ info )
362
363
364
365 CALL ztrtrs(
'U',
'N',
'N', n, nrhs,
366 $ a, lda, b, ldb, info )
367 IF( info.GT.0 ) THEN
368 RETURN
369 END IF
370 scllen = n
371 ELSE
372
373
374
375
376
377 CALL ztrtrs(
'U',
'C',
'N', n, nrhs,
378 $ a, lda, b, ldb, info )
379
380 IF( info.GT.0 ) THEN
381 RETURN
382 END IF
383
384
385
386 DO 20 j = 1, nrhs
387 DO 10 i = n + 1, m
388 b( i, j ) = czero
389 10 CONTINUE
390 20 CONTINUE
391
392
393
394 CALL zgemqr(
'L',
'N', m, nrhs, n, a, lda,
395 $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
396 $ info )
397
398 scllen = m
399
400 END IF
401
402 ELSE
403
404
405
406 CALL zgelq( m, n, a, lda, work( lw2+1 ), lw1,
407 $ work( 1 ), lw2, info )
408
409
410
411 IF( .NOT.tran ) THEN
412
413
414
415
416
417 CALL ztrtrs(
'L',
'N',
'N', m, nrhs,
418 $ a, lda, b, ldb, info )
419
420 IF( info.GT.0 ) THEN
421 RETURN
422 END IF
423
424
425
426 DO 40 j = 1, nrhs
427 DO 30 i = m + 1, n
428 b( i, j ) = czero
429 30 CONTINUE
430 40 CONTINUE
431
432
433
434 CALL zgemlq(
'L',
'C', n, nrhs, m, a, lda,
435 $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
436 $ info )
437
438
439
440 scllen = n
441
442 ELSE
443
444
445
446
447
448 CALL zgemlq(
'L',
'N', n, nrhs, m, a, lda,
449 $ work( lw2+1 ), lw1, b, ldb, work( 1 ), lw2,
450 $ info )
451
452
453
454
455
456 CALL ztrtrs(
'L',
'C',
'N', m, nrhs,
457 $ a, lda, b, ldb, info )
458
459 IF( info.GT.0 ) THEN
460 RETURN
461 END IF
462
463 scllen = m
464
465 END IF
466
467 END IF
468
469
470
471 IF( iascl.EQ.1 ) THEN
472 CALL zlascl(
'G', 0, 0, anrm, smlnum, scllen, nrhs, b, ldb,
473 $ info )
474 ELSE IF( iascl.EQ.2 ) THEN
475 CALL zlascl(
'G', 0, 0, anrm, bignum, scllen, nrhs, b, ldb,
476 $ info )
477 END IF
478 IF( ibscl.EQ.1 ) THEN
479 CALL zlascl(
'G', 0, 0, smlnum, bnrm, scllen, nrhs, b, ldb,
480 $ info )
481 ELSE IF( ibscl.EQ.2 ) THEN
482 CALL zlascl(
'G', 0, 0, bignum, bnrm, scllen, nrhs, b, ldb,
483 $ info )
484 END IF
485
486 50 CONTINUE
487 work( 1 ) = dble( tszo + lwo )
488 RETURN
489
490
491
subroutine xerbla(srname, info)
subroutine zgelq(m, n, a, lda, t, tsize, work, lwork, info)
ZGELQ
subroutine zgemlq(side, trans, m, n, k, a, lda, t, tsize, c, ldc, work, lwork, info)
ZGEMLQ
subroutine zgemqr(side, trans, m, n, k, a, lda, t, tsize, c, ldc, work, lwork, info)
ZGEMQR
subroutine zgeqr(m, n, a, lda, t, tsize, work, lwork, info)
ZGEQR
double precision function dlamch(cmach)
DLAMCH
double precision function zlange(norm, m, n, a, lda, work)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine zlascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine zlaset(uplo, m, n, alpha, beta, a, lda)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
logical function lsame(ca, cb)
LSAME
subroutine ztrtrs(uplo, trans, diag, n, nrhs, a, lda, b, ldb, info)
ZTRTRS