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