209
210
211
212
213
214
215 INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK
216 DOUBLE PRECISION RCOND
217
218
219 INTEGER IWORK( * )
220 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), S( * ), WORK( * )
221
222
223
224
225
226 DOUBLE PRECISION ZERO, ONE, TWO
227 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
228
229
230 LOGICAL LQUERY
231 INTEGER IASCL, IBSCL, IE, IL, ITAU, ITAUP, ITAUQ,
232 $ LDWORK, LIWORK, MAXMN, MAXWRK, MINMN, MINWRK,
233 $ MM, MNTHR, NLVL, NWORK, SMLSIZ, WLALSD
234 DOUBLE PRECISION ANRM, BIGNUM, BNRM, EPS, SFMIN, SMLNUM
235
236
239
240
241 INTEGER ILAENV
242 DOUBLE PRECISION DLAMCH, DLANGE
244
245
246 INTRINSIC dble, int, log, max, min
247
248
249
250
251
252 info = 0
253 minmn = min( m, n )
254 maxmn = max( m, n )
255 mnthr =
ilaenv( 6,
'DGELSD',
' ', m, n, nrhs, -1 )
256 lquery = ( lwork.EQ.-1 )
257 IF( m.LT.0 ) THEN
258 info = -1
259 ELSE IF( n.LT.0 ) THEN
260 info = -2
261 ELSE IF( nrhs.LT.0 ) THEN
262 info = -3
263 ELSE IF( lda.LT.max( 1, m ) ) THEN
264 info = -5
265 ELSE IF( ldb.LT.max( 1, maxmn ) ) THEN
266 info = -7
267 END IF
268
269 smlsiz =
ilaenv( 9,
'DGELSD',
' ', 0, 0, 0, 0 )
270
271
272
273
274
275
276
277
278 minwrk = 1
279 liwork = 1
280 minmn = max( 1, minmn )
281 nlvl = max( int( log( dble( minmn ) / dble( smlsiz+1 ) ) /
282 $ log( two ) ) + 1, 0 )
283
284 IF( info.EQ.0 ) THEN
285 maxwrk = 0
286 liwork = 3*minmn*nlvl + 11*minmn
287 mm = m
288 IF( m.GE.n .AND. m.GE.mnthr ) THEN
289
290
291
292 mm = n
293 maxwrk = max( maxwrk, n+n*
ilaenv( 1,
'DGEQRF',
' ', m, n,
294 $ -1, -1 ) )
295 maxwrk = max( maxwrk, n+nrhs*
296 $
ilaenv( 1,
'DORMQR',
'LT', m, nrhs, n, -1 ) )
297 END IF
298 IF( m.GE.n ) THEN
299
300
301
302 maxwrk = max( maxwrk, 3*n+( mm+n )*
303 $
ilaenv( 1,
'DGEBRD',
' ', mm, n, -1, -1 ) )
304 maxwrk = max( maxwrk, 3*n+nrhs*
305 $
ilaenv( 1,
'DORMBR',
'QLT', mm, nrhs, n, -1 ) )
306 maxwrk = max( maxwrk, 3*n+( n-1 )*
307 $
ilaenv( 1,
'DORMBR',
'PLN', n, nrhs, n, -1 ) )
308 wlalsd = 9*n+2*n*smlsiz+8*n*nlvl+n*nrhs+(smlsiz+1)**2
309 maxwrk = max( maxwrk, 3*n+wlalsd )
310 minwrk = max( 3*n+mm, 3*n+nrhs, 3*n+wlalsd )
311 END IF
312 IF( n.GT.m ) THEN
313 wlalsd = 9*m+2*m*smlsiz+8*m*nlvl+m*nrhs+(smlsiz+1)**2
314 IF( n.GE.mnthr ) THEN
315
316
317
318
319 maxwrk = m + m*
ilaenv( 1,
'DGELQF',
' ', m, n, -1, -1 )
320 maxwrk = max( maxwrk, m*m+4*m+2*m*
321 $
ilaenv( 1,
'DGEBRD',
' ', m, m, -1, -1 ) )
322 maxwrk = max( maxwrk, m*m+4*m+nrhs*
323 $
ilaenv( 1,
'DORMBR',
'QLT', m, nrhs, m, -1 ) )
324 maxwrk = max( maxwrk, m*m+4*m+( m-1 )*
325 $
ilaenv( 1,
'DORMBR',
'PLN', m, nrhs, m, -1 ) )
326 IF( nrhs.GT.1 ) THEN
327 maxwrk = max( maxwrk, m*m+m+m*nrhs )
328 ELSE
329 maxwrk = max( maxwrk, m*m+2*m )
330 END IF
331 maxwrk = max( maxwrk, m+nrhs*
332 $
ilaenv( 1,
'DORMLQ',
'LT', n, nrhs, m, -1 ) )
333 maxwrk = max( maxwrk, m*m+4*m+wlalsd )
334
335
336 maxwrk = max( maxwrk,
337 $ 4*m+m*m+max( m, 2*m-4, nrhs, n-3*m ) )
338 ELSE
339
340
341
342 maxwrk = 3*m + ( n+m )*
ilaenv( 1,
'DGEBRD',
' ', m, n,
343 $ -1, -1 )
344 maxwrk = max( maxwrk, 3*m+nrhs*
345 $
ilaenv( 1,
'DORMBR',
'QLT', m, nrhs, n, -1 ) )
346 maxwrk = max( maxwrk, 3*m+m*
347 $
ilaenv( 1,
'DORMBR',
'PLN', n, nrhs, m, -1 ) )
348 maxwrk = max( maxwrk, 3*m+wlalsd )
349 END IF
350 minwrk = max( 3*m+nrhs, 3*m+m, 3*m+wlalsd )
351 END IF
352 minwrk = min( minwrk, maxwrk )
353 work( 1 ) = maxwrk
354 iwork( 1 ) = liwork
355
356 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
357 info = -12
358 END IF
359 END IF
360
361 IF( info.NE.0 ) THEN
362 CALL xerbla(
'DGELSD', -info )
363 RETURN
364 ELSE IF( lquery ) THEN
365 GO TO 10
366 END IF
367
368
369
370 IF( m.EQ.0 .OR. n.EQ.0 ) THEN
371 rank = 0
372 RETURN
373 END IF
374
375
376
379 smlnum = sfmin / eps
380 bignum = one / smlnum
381 CALL dlabad( smlnum, bignum )
382
383
384
385 anrm =
dlange(
'M', m, n, a, lda, work )
386 iascl = 0
387 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
388
389
390
391 CALL dlascl(
'G', 0, 0, anrm, smlnum, m, n, a, lda, info )
392 iascl = 1
393 ELSE IF( anrm.GT.bignum ) THEN
394
395
396
397 CALL dlascl(
'G', 0, 0, anrm, bignum, m, n, a, lda, info )
398 iascl = 2
399 ELSE IF( anrm.EQ.zero ) THEN
400
401
402
403 CALL dlaset(
'F', max( m, n ), nrhs, zero, zero, b, ldb )
404 CALL dlaset(
'F', minmn, 1, zero, zero, s, 1 )
405 rank = 0
406 GO TO 10
407 END IF
408
409
410
411 bnrm =
dlange(
'M', m, nrhs, b, ldb, work )
412 ibscl = 0
413 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
414
415
416
417 CALL dlascl(
'G', 0, 0, bnrm, smlnum, m, nrhs, b, ldb, info )
418 ibscl = 1
419 ELSE IF( bnrm.GT.bignum ) THEN
420
421
422
423 CALL dlascl(
'G', 0, 0, bnrm, bignum, m, nrhs, b, ldb, info )
424 ibscl = 2
425 END IF
426
427
428
429 IF( m.LT.n )
430 $
CALL dlaset(
'F', n-m, nrhs, zero, zero, b( m+1, 1 ), ldb )
431
432
433
434 IF( m.GE.n ) THEN
435
436
437
438 mm = m
439 IF( m.GE.mnthr ) THEN
440
441
442
443 mm = n
444 itau = 1
445 nwork = itau + n
446
447
448
449
450 CALL dgeqrf( m, n, a, lda, work( itau ), work( nwork ),
451 $ lwork-nwork+1, info )
452
453
454
455
456 CALL dormqr(
'L',
'T', m, nrhs, n, a, lda, work( itau ), b,
457 $ ldb, work( nwork ), lwork-nwork+1, info )
458
459
460
461 IF( n.GT.1 ) THEN
462 CALL dlaset(
'L', n-1, n-1, zero, zero, a( 2, 1 ), lda )
463 END IF
464 END IF
465
466 ie = 1
467 itauq = ie + n
468 itaup = itauq + n
469 nwork = itaup + n
470
471
472
473
474 CALL dgebrd( mm, n, a, lda, s, work( ie ), work( itauq ),
475 $ work( itaup ), work( nwork ), lwork-nwork+1,
476 $ info )
477
478
479
480
481 CALL dormbr(
'Q',
'L',
'T', mm, nrhs, n, a, lda, work( itauq ),
482 $ b, ldb, work( nwork ), lwork-nwork+1, info )
483
484
485
486 CALL dlalsd(
'U', smlsiz, n, nrhs, s, work( ie ), b, ldb,
487 $ rcond, rank, work( nwork ), iwork, info )
488 IF( info.NE.0 ) THEN
489 GO TO 10
490 END IF
491
492
493
494 CALL dormbr(
'P',
'L',
'N', n, nrhs, n, a, lda, work( itaup ),
495 $ b, ldb, work( nwork ), lwork-nwork+1, info )
496
497 ELSE IF( n.GE.mnthr .AND. lwork.GE.4*m+m*m+
498 $ max( m, 2*m-4, nrhs, n-3*m, wlalsd ) ) THEN
499
500
501
502
503 ldwork = m
504 IF( lwork.GE.max( 4*m+m*lda+max( m, 2*m-4, nrhs, n-3*m ),
505 $ m*lda+m+m*nrhs, 4*m+m*lda+wlalsd ) )ldwork = lda
506 itau = 1
507 nwork = m + 1
508
509
510
511
512 CALL dgelqf( m, n, a, lda, work( itau ), work( nwork ),
513 $ lwork-nwork+1, info )
514 il = nwork
515
516
517
518 CALL dlacpy(
'L', m, m, a, lda, work( il ), ldwork )
519 CALL dlaset(
'U', m-1, m-1, zero, zero, work( il+ldwork ),
520 $ ldwork )
521 ie = il + ldwork*m
522 itauq = ie + m
523 itaup = itauq + m
524 nwork = itaup + m
525
526
527
528
529 CALL dgebrd( m, m, work( il ), ldwork, s, work( ie ),
530 $ work( itauq ), work( itaup ), work( nwork ),
531 $ lwork-nwork+1, info )
532
533
534
535
536 CALL dormbr(
'Q',
'L',
'T', m, nrhs, m, work( il ), ldwork,
537 $ work( itauq ), b, ldb, work( nwork ),
538 $ lwork-nwork+1, info )
539
540
541
542 CALL dlalsd(
'U', smlsiz, m, nrhs, s, work( ie ), b, ldb,
543 $ rcond, rank, work( nwork ), iwork, info )
544 IF( info.NE.0 ) THEN
545 GO TO 10
546 END IF
547
548
549
550 CALL dormbr(
'P',
'L',
'N', m, nrhs, m, work( il ), ldwork,
551 $ work( itaup ), b, ldb, work( nwork ),
552 $ lwork-nwork+1, info )
553
554
555
556 CALL dlaset(
'F', n-m, nrhs, zero, zero, b( m+1, 1 ), ldb )
557 nwork = itau + m
558
559
560
561
562 CALL dormlq(
'L',
'T', n, nrhs, m, a, lda, work( itau ), b,
563 $ ldb, work( nwork ), lwork-nwork+1, info )
564
565 ELSE
566
567
568
569 ie = 1
570 itauq = ie + m
571 itaup = itauq + m
572 nwork = itaup + m
573
574
575
576
577 CALL dgebrd( m, n, a, lda, s, work( ie ), work( itauq ),
578 $ work( itaup ), work( nwork ), lwork-nwork+1,
579 $ info )
580
581
582
583
584 CALL dormbr(
'Q',
'L',
'T', m, nrhs, n, a, lda, work( itauq ),
585 $ b, ldb, work( nwork ), lwork-nwork+1, info )
586
587
588
589 CALL dlalsd(
'L', smlsiz, m, nrhs, s, work( ie ), b, ldb,
590 $ rcond, rank, work( nwork ), iwork, info )
591 IF( info.NE.0 ) THEN
592 GO TO 10
593 END IF
594
595
596
597 CALL dormbr(
'P',
'L',
'N', n, nrhs, m, a, lda, work( itaup ),
598 $ b, ldb, work( nwork ), lwork-nwork+1, info )
599
600 END IF
601
602
603
604 IF( iascl.EQ.1 ) THEN
605 CALL dlascl(
'G', 0, 0, anrm, smlnum, n, nrhs, b, ldb, info )
606 CALL dlascl(
'G', 0, 0, smlnum, anrm, minmn, 1, s, minmn,
607 $ info )
608 ELSE IF( iascl.EQ.2 ) THEN
609 CALL dlascl(
'G', 0, 0, anrm, bignum, n, nrhs, b, ldb, info )
610 CALL dlascl(
'G', 0, 0, bignum, anrm, minmn, 1, s, minmn,
611 $ info )
612 END IF
613 IF( ibscl.EQ.1 ) THEN
614 CALL dlascl(
'G', 0, 0, smlnum, bnrm, n, nrhs, b, ldb, info )
615 ELSE IF( ibscl.EQ.2 ) THEN
616 CALL dlascl(
'G', 0, 0, bignum, bnrm, n, nrhs, b, ldb, info )
617 END IF
618
619 10 CONTINUE
620 work( 1 ) = maxwrk
621 iwork( 1 ) = liwork
622 RETURN
623
624
625
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 dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
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 dlange(NORM, M, N, A, LDA, WORK)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine dgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGEQRF
subroutine dgelqf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGELQF
subroutine dgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
DGEBRD
subroutine dormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMQR
subroutine dormlq(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMLQ
subroutine dlalsd(UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND, RANK, WORK, IWORK, INFO)
DLALSD uses the singular value decomposition of A to solve the least squares problem.
subroutine dormbr(VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMBR