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