225
226
227
228
229
230
231 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
232 REAL RCOND
233
234
235 INTEGER IWORK( * )
236 REAL RWORK( * ), S( * )
237 COMPLEX A( LDA, * ), B( LDB, * ), WORK( * )
238
239
240
241
242
243 REAL ZERO, ONE, TWO
244 parameter( zero = 0.0e+0, one = 1.0e+0, two = 2.0e+0 )
245 COMPLEX CZERO
246 parameter( czero = ( 0.0e+0, 0.0e+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 REAL ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM
254
255
260
261
262 INTEGER ILAENV
263 REAL CLANGE, SLAMCH
265
266
267 INTRINSIC int, log, max, min, real
268
269
270
271
272
273 info = 0
274 minmn = min( m, n )
275 maxmn = max( m, n )
276 lquery = ( lwork.EQ.-1 )
277 IF( m.LT.0 ) THEN
278 info = -1
279 ELSE IF( n.LT.0 ) THEN
280 info = -2
281 ELSE IF( nrhs.LT.0 ) THEN
282 info = -3
283 ELSE IF( lda.LT.max( 1, m ) ) THEN
284 info = -5
285 ELSE IF( ldb.LT.max( 1, maxmn ) ) THEN
286 info = -7
287 END IF
288
289
290
291
292
293
294
295
296 IF( info.EQ.0 ) THEN
297 minwrk = 1
298 maxwrk = 1
299 liwork = 1
300 lrwork = 1
301 IF( minmn.GT.0 ) THEN
302 smlsiz =
ilaenv( 9,
'CGELSD',
' ', 0, 0, 0, 0 )
303 mnthr =
ilaenv( 6,
'CGELSD',
' ', m, n, nrhs, -1 )
304 nlvl = max( int( log( real( minmn ) / real( smlsiz + 1 ) ) /
305 $ log( two ) ) + 1, 0 )
306 liwork = 3*minmn*nlvl + 11*minmn
307 mm = m
308 IF( m.GE.n .AND. m.GE.mnthr ) THEN
309
310
311
312
313 mm = n
314 maxwrk = max( maxwrk, n*
ilaenv( 1,
'CGEQRF',
' ', m, n,
315 $ -1, -1 ) )
316 maxwrk = max( maxwrk, nrhs*
ilaenv( 1,
'CUNMQR',
'LC', m,
317 $ nrhs, n, -1 ) )
318 END IF
319 IF( m.GE.n ) THEN
320
321
322
323 lrwork = 10*n + 2*n*smlsiz + 8*n*nlvl + 3*smlsiz*nrhs +
324 $ max( (smlsiz+1)**2, n*(1+nrhs) + 2*nrhs )
325 maxwrk = max( maxwrk, 2*n + ( mm + n )*
ilaenv( 1,
326 $ 'CGEBRD', ' ', mm, n, -1, -1 ) )
327 maxwrk = max( maxwrk, 2*n + nrhs*
ilaenv( 1,
'CUNMBR',
328 $ 'QLC', mm, nrhs, n, -1 ) )
329 maxwrk = max( maxwrk, 2*n + ( n - 1 )*
ilaenv( 1,
330 $ 'CUNMBR', 'PLN', n, nrhs, n, -1 ) )
331 maxwrk = max( maxwrk, 2*n + n*nrhs )
332 minwrk = max( 2*n + mm, 2*n + n*nrhs )
333 END IF
334 IF( n.GT.m ) THEN
335 lrwork = 10*m + 2*m*smlsiz + 8*m*nlvl + 3*smlsiz*nrhs +
336 $ max( (smlsiz+1)**2, n*(1+nrhs) + 2*nrhs )
337 IF( n.GE.mnthr ) THEN
338
339
340
341
342 maxwrk = m + m*
ilaenv( 1,
'CGELQF',
' ', m, n, -1,
343 $ -1 )
344 maxwrk = max( maxwrk, m*m + 4*m + 2*m*
ilaenv( 1,
345 $ 'CGEBRD', ' ', m, m, -1, -1 ) )
346 maxwrk = max( maxwrk, m*m + 4*m + nrhs*
ilaenv( 1,
347 $ 'CUNMBR', 'QLC', m, nrhs, m, -1 ) )
348 maxwrk = max( maxwrk, m*m + 4*m + ( m - 1 )*
ilaenv( 1,
349 $ 'CUNMLQ', 'LC', n, nrhs, m, -1 ) )
350 IF( nrhs.GT.1 ) THEN
351 maxwrk = max( maxwrk, m*m + m + m*nrhs )
352 ELSE
353 maxwrk = max( maxwrk, m*m + 2*m )
354 END IF
355 maxwrk = max( maxwrk, m*m + 4*m + m*nrhs )
356
357
358 maxwrk = max( maxwrk,
359 $ 4*m+m*m+max( m, 2*m-4, nrhs, n-3*m ) )
360 ELSE
361
362
363
364 maxwrk = 2*m + ( n + m )*
ilaenv( 1,
'CGEBRD',
' ', m,
365 $ n, -1, -1 )
366 maxwrk = max( maxwrk, 2*m + nrhs*
ilaenv( 1,
'CUNMBR',
367 $ 'QLC', m, nrhs, m, -1 ) )
368 maxwrk = max( maxwrk, 2*m + m*
ilaenv( 1,
'CUNMBR',
369 $ 'PLN', n, nrhs, m, -1 ) )
370 maxwrk = max( maxwrk, 2*m + m*nrhs )
371 END IF
372 minwrk = max( 2*m + n, 2*m + m*nrhs )
373 END IF
374 END IF
375 minwrk = min( minwrk, maxwrk )
376 work( 1 ) = maxwrk
377 iwork( 1 ) = liwork
378 rwork( 1 ) = lrwork
379
380 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
381 info = -12
382 END IF
383 END IF
384
385 IF( info.NE.0 ) THEN
386 CALL xerbla(
'CGELSD', -info )
387 RETURN
388 ELSE IF( lquery ) THEN
389 RETURN
390 END IF
391
392
393
394 IF( m.EQ.0 .OR. n.EQ.0 ) THEN
395 rank = 0
396 RETURN
397 END IF
398
399
400
403 smlnum = sfmin / eps
404 bignum = one / smlnum
405 CALL slabad( smlnum, bignum )
406
407
408
409 anrm =
clange(
'M', m, n, a, lda, rwork )
410 iascl = 0
411 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
412
413
414
415 CALL clascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
416 iascl = 1
417 ELSE IF( anrm.GT.bignum ) THEN
418
419
420
421 CALL clascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, info )
422 iascl = 2
423 ELSE IF( anrm.EQ.zero ) THEN
424
425
426
427 CALL claset(
'F', max( m, n ), nrhs, czero, czero, b, ldb )
428 CALL slaset(
'F', minmn, 1, zero, zero, s, 1 )
429 rank = 0
430 GO TO 10
431 END IF
432
433
434
435 bnrm =
clange(
'M', m, nrhs, b, ldb, rwork )
436 ibscl = 0
437 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
438
439
440
441 CALL clascl(
'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
442 ibscl = 1
443 ELSE IF( bnrm.GT.bignum ) THEN
444
445
446
447 CALL clascl(
'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
448 ibscl = 2
449 END IF
450
451
452
453 IF( m.LT.n )
454 $
CALL claset(
'F', n-m, nrhs, czero, czero, b( m+1, 1 ), ldb )
455
456
457
458 IF( m.GE.n ) THEN
459
460
461
462 mm = m
463 IF( m.GE.mnthr ) THEN
464
465
466
467 mm = n
468 itau = 1
469 nwork = itau + n
470
471
472
473
474
475 CALL cgeqrf( m, n, a, lda, work( itau ), work( nwork ),
476 $ lwork-nwork+1, info )
477
478
479
480
481
482 CALL cunmqr(
'L',
'C', m, nrhs, n, a, lda, work( itau ), b,
483 $ ldb, work( nwork ), lwork-nwork+1, info )
484
485
486
487 IF( n.GT.1 ) THEN
488 CALL claset(
'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 cgebrd( 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 cunmbr(
'Q',
'L',
'C', mm, nrhs, n, a, lda, work( itauq ),
511 $ b, ldb, work( nwork ), lwork-nwork+1, info )
512
513
514
515 CALL clalsd(
'U', smlsiz, n, nrhs, s, rwork( ie ), b, ldb,
516 $ rcond, rank, work( nwork ), rwork( nrwork ),
517 $ iwork, info )
518 IF( info.NE.0 ) THEN
519 GO TO 10
520 END IF
521
522
523
524 CALL cunmbr(
'P',
'L',
'N', n, nrhs, n, a, lda, work( itaup ),
525 $ b, ldb, work( nwork ), lwork-nwork+1, info )
526
527 ELSE IF( n.GE.mnthr .AND. lwork.GE.4*m+m*m+
528 $ max( m, 2*m-4, nrhs, n-3*m ) ) THEN
529
530
531
532
533 ldwork = m
534 IF( lwork.GE.max( 4*m+m*lda+max( m, 2*m-4, nrhs, n-3*m ),
535 $ m*lda+m+m*nrhs ) )ldwork = lda
536 itau = 1
537 nwork = m + 1
538
539
540
541
542 CALL cgelqf( m, n, a, lda, work( itau ), work( nwork ),
543 $ lwork-nwork+1, info )
544 il = nwork
545
546
547
548 CALL clacpy(
'L', m, m, a, lda, work( il ), ldwork )
549 CALL claset(
'U', m-1, m-1, czero, czero, work( il+ldwork ),
550 $ ldwork )
551 itauq = il + ldwork*m
552 itaup = itauq + m
553 nwork = itaup + m
554 ie = 1
555 nrwork = ie + m
556
557
558
559
560
561 CALL cgebrd( m, m, work( il ), ldwork, s, rwork( ie ),
562 $ work( itauq ), work( itaup ), work( nwork ),
563 $ lwork-nwork+1, info )
564
565
566
567
568 CALL cunmbr(
'Q',
'L',
'C', m, nrhs, m, work( il ), ldwork,
569 $ work( itauq ), b, ldb, work( nwork ),
570 $ lwork-nwork+1, info )
571
572
573
574 CALL clalsd(
'U', smlsiz, m, nrhs, s, rwork( ie ), b, ldb,
575 $ rcond, rank, work( nwork ), rwork( nrwork ),
576 $ iwork, info )
577 IF( info.NE.0 ) THEN
578 GO TO 10
579 END IF
580
581
582
583 CALL cunmbr(
'P',
'L',
'N', m, nrhs, m, work( il ), ldwork,
584 $ work( itaup ), b, ldb, work( nwork ),
585 $ lwork-nwork+1, info )
586
587
588
589 CALL claset(
'F', n-m, nrhs, czero, czero, b( m+1, 1 ), ldb )
590 nwork = itau + m
591
592
593
594
595 CALL cunmlq(
'L',
'C', n, nrhs, m, a, lda, work( itau ), b,
596 $ ldb, work( nwork ), lwork-nwork+1, info )
597
598 ELSE
599
600
601
602 itauq = 1
603 itaup = itauq + m
604 nwork = itaup + m
605 ie = 1
606 nrwork = ie + m
607
608
609
610
611
612 CALL cgebrd( m, n, a, lda, s, rwork( ie ), work( itauq ),
613 $ work( itaup ), work( nwork ), lwork-nwork+1,
614 $ info )
615
616
617
618
619 CALL cunmbr(
'Q',
'L',
'C', m, nrhs, n, a, lda, work( itauq ),
620 $ b, ldb, work( nwork ), lwork-nwork+1, info )
621
622
623
624 CALL clalsd(
'L', smlsiz, m, nrhs, s, rwork( ie ), b, ldb,
625 $ rcond, rank, work( nwork ), rwork( nrwork ),
626 $ iwork, info )
627 IF( info.NE.0 ) THEN
628 GO TO 10
629 END IF
630
631
632
633 CALL cunmbr(
'P',
'L',
'N', n, nrhs, m, a, lda, work( itaup ),
634 $ b, ldb, work( nwork ), lwork-nwork+1, info )
635
636 END IF
637
638
639
640 IF( iascl.EQ.1 ) THEN
641 CALL clascl(
'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
642 CALL slascl(
'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
643 $ info )
644 ELSE IF( iascl.EQ.2 ) THEN
645 CALL clascl(
'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
646 CALL slascl(
'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
647 $ info )
648 END IF
649 IF( ibscl.EQ.1 ) THEN
650 CALL clascl(
'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
651 ELSE IF( ibscl.EQ.2 ) THEN
652 CALL clascl(
'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
653 END IF
654
655 10 CONTINUE
656 work( 1 ) = maxwrk
657 iwork( 1 ) = liwork
658 rwork( 1 ) = lrwork
659 RETURN
660
661
662
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET 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
real function clange(NORM, M, N, A, LDA, WORK)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine cgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGEQRF
subroutine cgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
CGEBRD
subroutine cgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGELQF
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
subroutine clascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
subroutine clalsd(UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND, RANK, WORK, RWORK, IWORK, INFO)
CLALSD uses the singular value decomposition of A to solve the least squares problem.
subroutine cunmbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMBR
subroutine cunmlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMLQ
subroutine cunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMQR
real function slamch(CMACH)
SLAMCH