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