217
218
219
220
221
222
223 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
224 DOUBLE PRECISION RCOND
225
226
227 INTEGER IWORK( * )
228 DOUBLE PRECISION RWORK( * ), S( * )
229 COMPLEX*16 A( LDA, * ), B( LDB, * ), WORK( * )
230
231
232
233
234
235 DOUBLE PRECISION ZERO, ONE, TWO
236 parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0 )
237 COMPLEX*16 CZERO
238 parameter( czero = ( 0.0d+0, 0.0d+0 ) )
239
240
241 LOGICAL LQUERY
242 INTEGER IASCL, IBSCL, IE, IL, ITAU, ITAUP, ITAUQ,
243 $ LDWORK, LIWORK, LRWORK, MAXMN, MAXWRK, MINMN,
244 $ MINWRK, MM, MNTHR, NLVL, NRWORK, NWORK, SMLSIZ
245 DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM
246
247
252
253
254 INTEGER ILAENV
255 DOUBLE PRECISION DLAMCH, ZLANGE
257
258
259 INTRINSIC int, log, max, min, dble
260
261
262
263
264
265 info = 0
266 minmn = min( m, n )
267 maxmn = max( m, n )
268 lquery = ( lwork.EQ.-1 )
269 IF( m.LT.0 ) THEN
270 info = -1
271 ELSE IF( n.LT.0 ) THEN
272 info = -2
273 ELSE IF( nrhs.LT.0 ) THEN
274 info = -3
275 ELSE IF( lda.LT.max( 1, m ) ) THEN
276 info = -5
277 ELSE IF( ldb.LT.max( 1, maxmn ) ) THEN
278 info = -7
279 END IF
280
281
282
283
284
285
286
287
288 IF( info.EQ.0 ) THEN
289 minwrk = 1
290 maxwrk = 1
291 liwork = 1
292 lrwork = 1
293 IF( minmn.GT.0 ) THEN
294 smlsiz =
ilaenv( 9,
'ZGELSD',
' ', 0, 0, 0, 0 )
295 mnthr =
ilaenv( 6,
'ZGELSD',
' ', m, n, nrhs, -1 )
296 nlvl = max( int( log( dble( minmn ) / dble( smlsiz + 1 ) ) /
297 $ log( two ) ) + 1, 0 )
298 liwork = 3*minmn*nlvl + 11*minmn
299 mm = m
300 IF( m.GE.n .AND. m.GE.mnthr ) THEN
301
302
303
304
305 mm = n
306 maxwrk = max( maxwrk, n*
ilaenv( 1,
'ZGEQRF',
' ', m,
307 $ n,
308 $ -1, -1 ) )
309 maxwrk = max( maxwrk, nrhs*
ilaenv( 1,
'ZUNMQR',
'LC',
310 $ m,
311 $ nrhs, n, -1 ) )
312 END IF
313 IF( m.GE.n ) THEN
314
315
316
317 lrwork = 10*n + 2*n*smlsiz + 8*n*nlvl + 3*smlsiz*nrhs +
318 $ max( (smlsiz+1)**2, n*(1+nrhs) + 2*nrhs )
319 maxwrk = max( maxwrk, 2*n + ( mm + n )*
ilaenv( 1,
320 $ 'ZGEBRD', ' ', mm, n, -1, -1 ) )
321 maxwrk = max( maxwrk, 2*n + nrhs*
ilaenv( 1,
'ZUNMBR',
322 $ 'QLC', mm, nrhs, n, -1 ) )
323 maxwrk = max( maxwrk, 2*n + ( n - 1 )*
ilaenv( 1,
324 $ 'ZUNMBR', 'PLN', n, nrhs, n, -1 ) )
325 maxwrk = max( maxwrk, 2*n + n*nrhs )
326 minwrk = max( 2*n + mm, 2*n + n*nrhs )
327 END IF
328 IF( n.GT.m ) THEN
329 lrwork = 10*m + 2*m*smlsiz + 8*m*nlvl + 3*smlsiz*nrhs +
330 $ max( (smlsiz+1)**2, n*(1+nrhs) + 2*nrhs )
331 IF( n.GE.mnthr ) THEN
332
333
334
335
336 maxwrk = m + m*
ilaenv( 1,
'ZGELQF',
' ', m, n, -1,
337 $ -1 )
338 maxwrk = max( maxwrk, m*m + 4*m + 2*m*
ilaenv( 1,
339 $ 'ZGEBRD', ' ', m, m, -1, -1 ) )
340 maxwrk = max( maxwrk, m*m + 4*m + nrhs*
ilaenv( 1,
341 $ 'ZUNMBR', 'QLC', m, nrhs, m, -1 ) )
342 maxwrk = max( maxwrk,
343 $ m*m + 4*m + ( m - 1 )*
ilaenv( 1,
344 $ 'ZUNMLQ', 'LC', n, nrhs, m, -1 ) )
345 IF( nrhs.GT.1 ) THEN
346 maxwrk = max( maxwrk, m*m + m + m*nrhs )
347 ELSE
348 maxwrk = max( maxwrk, m*m + 2*m )
349 END IF
350 maxwrk = max( maxwrk, m*m + 4*m + m*nrhs )
351
352
353 maxwrk = max( maxwrk,
354 $ 4*m+m*m+max( m, 2*m-4, nrhs, n-3*m ) )
355 ELSE
356
357
358
359 maxwrk = 2*m + ( n + m )*
ilaenv( 1,
'ZGEBRD',
' ',
360 $ m,
361 $ n, -1, -1 )
362 maxwrk = max( maxwrk, 2*m + nrhs*
ilaenv( 1,
363 $ 'ZUNMBR',
364 $ 'QLC', m, nrhs, m, -1 ) )
365 maxwrk = max( maxwrk, 2*m + m*
ilaenv( 1,
'ZUNMBR',
366 $ 'PLN', n, nrhs, m, -1 ) )
367 maxwrk = max( maxwrk, 2*m + m*nrhs )
368 END IF
369 minwrk = max( 2*m + n, 2*m + m*nrhs )
370 END IF
371 END IF
372 minwrk = min( minwrk, maxwrk )
373 work( 1 ) = maxwrk
374 iwork( 1 ) = liwork
375 rwork( 1 ) = lrwork
376
377 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
378 info = -12
379 END IF
380 END IF
381
382 IF( info.NE.0 ) THEN
383 CALL xerbla(
'ZGELSD', -info )
384 RETURN
385 ELSE IF( lquery ) THEN
386 RETURN
387 END IF
388
389
390
391 IF( m.EQ.0 .OR. n.EQ.0 ) THEN
392 rank = 0
393 RETURN
394 END IF
395
396
397
400 smlnum = sfmin / eps
401 bignum = one / smlnum
402
403
404
405 anrm =
zlange(
'M', m, n, a, lda, rwork )
406 iascl = 0
407 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
408
409
410
411 CALL zlascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
412 iascl = 1
413 ELSE IF( anrm.GT.bignum ) THEN
414
415
416
417 CALL zlascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, info )
418 iascl = 2
419 ELSE IF( anrm.EQ.zero ) THEN
420
421
422
423 CALL zlaset(
'F', max( m, n ), nrhs, czero, czero, b, ldb )
424 CALL dlaset(
'F', minmn, 1, zero, zero, s, 1 )
425 rank = 0
426 GO TO 10
427 END IF
428
429
430
431 bnrm =
zlange(
'M', m, nrhs, b, ldb, rwork )
432 ibscl = 0
433 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
434
435
436
437 CALL zlascl(
'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb,
438 $ info )
439 ibscl = 1
440 ELSE IF( bnrm.GT.bignum ) THEN
441
442
443
444 CALL zlascl(
'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb,
445 $ info )
446 ibscl = 2
447 END IF
448
449
450
451 IF( m.LT.n )
452 $
CALL zlaset(
'F', n-m, nrhs, czero, czero, b( m+1, 1 ),
453 $ ldb )
454
455
456
457 IF( m.GE.n ) THEN
458
459
460
461 mm = m
462 IF( m.GE.mnthr ) THEN
463
464
465
466 mm = n
467 itau = 1
468 nwork = itau + n
469
470
471
472
473
474 CALL zgeqrf( m, n, a, lda, work( itau ), work( nwork ),
475 $ lwork-nwork+1, info )
476
477
478
479
480
481 CALL zunmqr(
'L',
'C', m, nrhs, n, a, lda, work( itau ),
482 $ b,
483 $ ldb, work( nwork ), lwork-nwork+1, info )
484
485
486
487 IF( n.GT.1 ) THEN
488 CALL zlaset(
'L', n-1, n-1, czero, czero, a( 2, 1 ),
489 $ lda )
490 END IF
491 END IF
492
493 itauq = 1
494 itaup = itauq + n
495 nwork = itaup + n
496 ie = 1
497 nrwork = ie + n
498
499
500
501
502
503 CALL zgebrd( mm, n, a, lda, s, rwork( ie ), work( itauq ),
504 $ work( itaup ), work( nwork ), lwork-nwork+1,
505 $ info )
506
507
508
509
510 CALL zunmbr(
'Q',
'L',
'C', mm, nrhs, n, a, lda,
511 $ work( itauq ),
512 $ b, ldb, work( nwork ), lwork-nwork+1, info )
513
514
515
516 CALL zlalsd(
'U', smlsiz, n, nrhs, s, rwork( ie ), b, ldb,
517 $ rcond, rank, work( nwork ), rwork( nrwork ),
518 $ iwork, info )
519 IF( info.NE.0 ) THEN
520 GO TO 10
521 END IF
522
523
524
525 CALL zunmbr(
'P',
'L',
'N', n, nrhs, n, a, lda,
526 $ work( itaup ),
527 $ b, ldb, work( nwork ), lwork-nwork+1, info )
528
529 ELSE IF( n.GE.mnthr .AND. lwork.GE.4*m+m*m+
530 $ max( m, 2*m-4, nrhs, n-3*m ) ) THEN
531
532
533
534
535 ldwork = m
536 IF( lwork.GE.max( 4*m+m*lda+max( m, 2*m-4, nrhs, n-3*m ),
537 $ m*lda+m+m*nrhs ) )ldwork = lda
538 itau = 1
539 nwork = m + 1
540
541
542
543
544 CALL zgelqf( m, n, a, lda, work( itau ), work( nwork ),
545 $ lwork-nwork+1, info )
546 il = nwork
547
548
549
550 CALL zlacpy(
'L', m, m, a, lda, work( il ), ldwork )
551 CALL zlaset(
'U', m-1, m-1, czero, czero, work( il+ldwork ),
552 $ ldwork )
553 itauq = il + ldwork*m
554 itaup = itauq + m
555 nwork = itaup + m
556 ie = 1
557 nrwork = ie + m
558
559
560
561
562
563 CALL zgebrd( m, m, work( il ), ldwork, s, rwork( ie ),
564 $ work( itauq ), work( itaup ), work( nwork ),
565 $ lwork-nwork+1, info )
566
567
568
569
570 CALL zunmbr(
'Q',
'L',
'C', m, nrhs, m, work( il ), ldwork,
571 $ work( itauq ), b, ldb, work( nwork ),
572 $ lwork-nwork+1, info )
573
574
575
576 CALL zlalsd(
'U', smlsiz, m, nrhs, s, rwork( ie ), b, ldb,
577 $ rcond, rank, work( nwork ), rwork( nrwork ),
578 $ iwork, info )
579 IF( info.NE.0 ) THEN
580 GO TO 10
581 END IF
582
583
584
585 CALL zunmbr(
'P',
'L',
'N', m, nrhs, m, work( il ), ldwork,
586 $ work( itaup ), b, ldb, work( nwork ),
587 $ lwork-nwork+1, info )
588
589
590
591 CALL zlaset(
'F', n-m, nrhs, czero, czero, b( m+1, 1 ),
592 $ ldb )
593 nwork = itau + m
594
595
596
597
598 CALL zunmlq(
'L',
'C', n, nrhs, m, a, lda, work( itau ), b,
599 $ ldb, work( nwork ), lwork-nwork+1, info )
600
601 ELSE
602
603
604
605 itauq = 1
606 itaup = itauq + m
607 nwork = itaup + m
608 ie = 1
609 nrwork = ie + m
610
611
612
613
614
615 CALL zgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
616 $ work( itaup ), work( nwork ), lwork-nwork+1,
617 $ info )
618
619
620
621
622 CALL zunmbr(
'Q',
'L',
'C', m, nrhs, n, a, lda,
623 $ work( itauq ),
624 $ b, ldb, work( nwork ), lwork-nwork+1, info )
625
626
627
628 CALL zlalsd(
'L', smlsiz, m, nrhs, s, rwork( ie ), b, ldb,
629 $ rcond, rank, work( nwork ), rwork( nrwork ),
630 $ iwork, info )
631 IF( info.NE.0 ) THEN
632 GO TO 10
633 END IF
634
635
636
637 CALL zunmbr(
'P',
'L',
'N', n, nrhs, m, a, lda,
638 $ work( itaup ),
639 $ b, ldb, work( nwork ), lwork-nwork+1, info )
640
641 END IF
642
643
644
645 IF( iascl.EQ.1 ) THEN
646 CALL zlascl(
'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb,
647 $ info )
648 CALL dlascl(
'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
649 $ info )
650 ELSE IF( iascl.EQ.2 ) THEN
651 CALL zlascl(
'G', 0, 0, anrm, bignum, n, nrhs, b, ldb,
652 $ info )
653 CALL dlascl(
'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
654 $ info )
655 END IF
656 IF( ibscl.EQ.1 ) THEN
657 CALL zlascl(
'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb,
658 $ info )
659 ELSE IF( ibscl.EQ.2 ) THEN
660 CALL zlascl(
'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb,
661 $ info )
662 END IF
663
664 10 CONTINUE
665 work( 1 ) = maxwrk
666 iwork( 1 ) = liwork
667 rwork( 1 ) = lrwork
668 RETURN
669
670
671
subroutine xerbla(srname, info)
subroutine zgebrd(m, n, a, lda, d, e, tauq, taup, work, lwork, info)
ZGEBRD
subroutine zgelqf(m, n, a, lda, tau, work, lwork, info)
ZGELQF
subroutine zgeqrf(m, n, a, lda, tau, work, lwork, info)
ZGEQRF
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zlalsd(uplo, smlsiz, n, nrhs, d, e, b, ldb, rcond, rank, work, rwork, iwork, info)
ZLALSD uses the singular value decomposition of A to solve the least squares problem.
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 dlascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
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.
subroutine zunmbr(vect, side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
ZUNMBR
subroutine zunmlq(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
ZUNMLQ
subroutine zunmqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
ZUNMQR