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