181
182
183
184
185
186
187 CHARACTER UPLO
188 INTEGER INFO, LDB, N, NRHS, RANK, SMLSIZ
189 DOUBLE PRECISION RCOND
190
191
192 INTEGER IWORK( * )
193 DOUBLE PRECISION D( * ), E( * ), RWORK( * )
194 COMPLEX*16 B( LDB, * ), WORK( * )
195
196
197
198
199
200 DOUBLE PRECISION ZERO, ONE, TWO
201 parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
202 COMPLEX*16 CZERO
203 parameter( czero = ( 0.0d0, 0.0d0 ) )
204
205
206 INTEGER BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM,
207 $ GIVPTR, I, ICMPQ1, ICMPQ2, IRWB, IRWIB, IRWRB,
208 $ IRWU, IRWVT, IRWWRK, IWK, J, JCOL, JIMAG,
209 $ JREAL, JROW, K, NLVL, NM1, NRWORK, NSIZE, NSUB,
210 $ PERM, POLES, S, SIZEI, SMLSZP, SQRE, ST, ST1,
211 $ U, VT, Z
212 DOUBLE PRECISION CS, EPS, ORGNRM, RCND, R, SN, TOL
213
214
215 INTEGER IDAMAX
216 DOUBLE PRECISION DLAMCH, DLANST
218
219
223
224
225 INTRINSIC abs, dble, dcmplx, dimag, int, log, sign
226
227
228
229
230
231 info = 0
232
233 IF( n.LT.0 ) THEN
234 info = -3
235 ELSE IF( nrhs.LT.1 ) THEN
236 info = -4
237 ELSE IF( ( ldb.LT.1 ) .OR. ( ldb.LT.n ) ) THEN
238 info = -8
239 END IF
240 IF( info.NE.0 ) THEN
241 CALL xerbla(
'ZLALSD', -info )
242 RETURN
243 END IF
244
246
247
248
249 IF( ( rcond.LE.zero ) .OR. ( rcond.GE.one ) ) THEN
250 rcnd = eps
251 ELSE
252 rcnd = rcond
253 END IF
254
255 rank = 0
256
257
258
259 IF( n.EQ.0 ) THEN
260 RETURN
261 ELSE IF( n.EQ.1 ) THEN
262 IF( d( 1 ).EQ.zero ) THEN
263 CALL zlaset(
'A', 1, nrhs, czero, czero, b, ldb )
264 ELSE
265 rank = 1
266 CALL zlascl(
'G', 0, 0, d( 1 ), one, 1, nrhs, b, ldb, info )
267 d( 1 ) = abs( d( 1 ) )
268 END IF
269 RETURN
270 END IF
271
272
273
274 IF( uplo.EQ.'L' ) THEN
275 DO 10 i = 1, n - 1
276 CALL dlartg( d( i ), e( i ), cs, sn, r )
277 d( i ) = r
278 e( i ) = sn*d( i+1 )
279 d( i+1 ) = cs*d( i+1 )
280 IF( nrhs.EQ.1 ) THEN
281 CALL zdrot( 1, b( i, 1 ), 1, b( i+1, 1 ), 1, cs, sn )
282 ELSE
283 rwork( i*2-1 ) = cs
284 rwork( i*2 ) = sn
285 END IF
286 10 CONTINUE
287 IF( nrhs.GT.1 ) THEN
288 DO 30 i = 1, nrhs
289 DO 20 j = 1, n - 1
290 cs = rwork( j*2-1 )
291 sn = rwork( j*2 )
292 CALL zdrot( 1, b( j, i ), 1, b( j+1, i ), 1, cs, sn )
293 20 CONTINUE
294 30 CONTINUE
295 END IF
296 END IF
297
298
299
300 nm1 = n - 1
301 orgnrm =
dlanst(
'M', n, d, e )
302 IF( orgnrm.EQ.zero ) THEN
303 CALL zlaset(
'A', n, nrhs, czero, czero, b, ldb )
304 RETURN
305 END IF
306
307 CALL dlascl(
'G', 0, 0, orgnrm, one, n, 1, d, n, info )
308 CALL dlascl(
'G', 0, 0, orgnrm, one, nm1, 1, e, nm1, info )
309
310
311
312
313 IF( n.LE.smlsiz ) THEN
314 irwu = 1
315 irwvt = irwu + n*n
316 irwwrk = irwvt + n*n
317 irwrb = irwwrk
318 irwib = irwrb + n*nrhs
319 irwb = irwib + n*nrhs
320 CALL dlaset(
'A', n, n, zero, one, rwork( irwu ), n )
321 CALL dlaset(
'A', n, n, zero, one, rwork( irwvt ), n )
322 CALL dlasdq(
'U', 0, n, n, n, 0, d, e, rwork( irwvt ), n,
323 $ rwork( irwu ), n, rwork( irwwrk ), 1,
324 $ rwork( irwwrk ), info )
325 IF( info.NE.0 ) THEN
326 RETURN
327 END IF
328
329
330
331
332
333 j = irwb - 1
334 DO 50 jcol = 1, nrhs
335 DO 40 jrow = 1, n
336 j = j + 1
337 rwork( j ) = dble( b( jrow, jcol ) )
338 40 CONTINUE
339 50 CONTINUE
340 CALL dgemm(
'T',
'N', n, nrhs, n, one, rwork( irwu ), n,
341 $ rwork( irwb ), n, zero, rwork( irwrb ), n )
342 j = irwb - 1
343 DO 70 jcol = 1, nrhs
344 DO 60 jrow = 1, n
345 j = j + 1
346 rwork( j ) = dimag( b( jrow, jcol ) )
347 60 CONTINUE
348 70 CONTINUE
349 CALL dgemm(
'T',
'N', n, nrhs, n, one, rwork( irwu ), n,
350 $ rwork( irwb ), n, zero, rwork( irwib ), n )
351 jreal = irwrb - 1
352 jimag = irwib - 1
353 DO 90 jcol = 1, nrhs
354 DO 80 jrow = 1, n
355 jreal = jreal + 1
356 jimag = jimag + 1
357 b( jrow, jcol ) = dcmplx( rwork( jreal ),
358 $ rwork( jimag ) )
359 80 CONTINUE
360 90 CONTINUE
361
362 tol = rcnd*abs( d(
idamax( n, d, 1 ) ) )
363 DO 100 i = 1, n
364 IF( d( i ).LE.tol ) THEN
365 CALL zlaset(
'A', 1, nrhs, czero, czero, b( i, 1 ), ldb )
366 ELSE
367 CALL zlascl(
'G', 0, 0, d( i ), one, 1, nrhs, b( i, 1 ),
368 $ ldb, info )
369 rank = rank + 1
370 END IF
371 100 CONTINUE
372
373
374
375
376
377
378
379
380 j = irwb - 1
381 DO 120 jcol = 1, nrhs
382 DO 110 jrow = 1, n
383 j = j + 1
384 rwork( j ) = dble( b( jrow, jcol ) )
385 110 CONTINUE
386 120 CONTINUE
387 CALL dgemm(
'T',
'N', n, nrhs, n, one, rwork( irwvt ), n,
388 $ rwork( irwb ), n, zero, rwork( irwrb ), n )
389 j = irwb - 1
390 DO 140 jcol = 1, nrhs
391 DO 130 jrow = 1, n
392 j = j + 1
393 rwork( j ) = dimag( b( jrow, jcol ) )
394 130 CONTINUE
395 140 CONTINUE
396 CALL dgemm(
'T',
'N', n, nrhs, n, one, rwork( irwvt ), n,
397 $ rwork( irwb ), n, zero, rwork( irwib ), n )
398 jreal = irwrb - 1
399 jimag = irwib - 1
400 DO 160 jcol = 1, nrhs
401 DO 150 jrow = 1, n
402 jreal = jreal + 1
403 jimag = jimag + 1
404 b( jrow, jcol ) = dcmplx( rwork( jreal ),
405 $ rwork( jimag ) )
406 150 CONTINUE
407 160 CONTINUE
408
409
410
411 CALL dlascl(
'G', 0, 0, one, orgnrm, n, 1, d, n, info )
412 CALL dlasrt(
'D', n, d, info )
413 CALL zlascl(
'G', 0, 0, orgnrm, one, n, nrhs, b, ldb, info )
414
415 RETURN
416 END IF
417
418
419
420 nlvl = int( log( dble( n ) / dble( smlsiz+1 ) ) / log( two ) ) + 1
421
422 smlszp = smlsiz + 1
423
424 u = 1
425 vt = 1 + smlsiz*n
426 difl = vt + smlszp*n
427 difr = difl + nlvl*n
428 z = difr + nlvl*n*2
429 c = z + nlvl*n
430 s = c + n
431 poles = s + n
432 givnum = poles + 2*nlvl*n
433 nrwork = givnum + 2*nlvl*n
434 bx = 1
435
436 irwrb = nrwork
437 irwib = irwrb + smlsiz*nrhs
438 irwb = irwib + smlsiz*nrhs
439
440 sizei = 1 + n
441 k = sizei + n
442 givptr = k + n
443 perm = givptr + n
444 givcol = perm + nlvl*n
445 iwk = givcol + nlvl*n*2
446
447 st = 1
448 sqre = 0
449 icmpq1 = 1
450 icmpq2 = 0
451 nsub = 0
452
453 DO 170 i = 1, n
454 IF( abs( d( i ) ).LT.eps ) THEN
455 d( i ) = sign( eps, d( i ) )
456 END IF
457 170 CONTINUE
458
459 DO 240 i = 1, nm1
460 IF( ( abs( e( i ) ).LT.eps ) .OR. ( i.EQ.nm1 ) ) THEN
461 nsub = nsub + 1
462 iwork( nsub ) = st
463
464
465
466
467 IF( i.LT.nm1 ) THEN
468
469
470
471 nsize = i - st + 1
472 iwork( sizei+nsub-1 ) = nsize
473 ELSE IF( abs( e( i ) ).GE.eps ) THEN
474
475
476
477 nsize = n - st + 1
478 iwork( sizei+nsub-1 ) = nsize
479 ELSE
480
481
482
483
484
485 nsize = i - st + 1
486 iwork( sizei+nsub-1 ) = nsize
487 nsub = nsub + 1
488 iwork( nsub ) = n
489 iwork( sizei+nsub-1 ) = 1
490 CALL zcopy( nrhs, b( n, 1 ), ldb, work( bx+nm1 ), n )
491 END IF
492 st1 = st - 1
493 IF( nsize.EQ.1 ) THEN
494
495
496
497
498 CALL zcopy( nrhs, b( st, 1 ), ldb, work( bx+st1 ), n )
499 ELSE IF( nsize.LE.smlsiz ) THEN
500
501
502
503 CALL dlaset(
'A', nsize, nsize, zero, one,
504 $ rwork( vt+st1 ), n )
505 CALL dlaset(
'A', nsize, nsize, zero, one,
506 $ rwork( u+st1 ), n )
507 CALL dlasdq(
'U', 0, nsize, nsize, nsize, 0, d( st ),
508 $ e( st ), rwork( vt+st1 ), n, rwork( u+st1 ),
509 $ n, rwork( nrwork ), 1, rwork( nrwork ),
510 $ info )
511 IF( info.NE.0 ) THEN
512 RETURN
513 END IF
514
515
516
517
518
519 j = irwb - 1
520 DO 190 jcol = 1, nrhs
521 DO 180 jrow = st, st + nsize - 1
522 j = j + 1
523 rwork( j ) = dble( b( jrow, jcol ) )
524 180 CONTINUE
525 190 CONTINUE
526 CALL dgemm(
'T',
'N', nsize, nrhs, nsize, one,
527 $ rwork( u+st1 ), n, rwork( irwb ), nsize,
528 $ zero, rwork( irwrb ), nsize )
529 j = irwb - 1
530 DO 210 jcol = 1, nrhs
531 DO 200 jrow = st, st + nsize - 1
532 j = j + 1
533 rwork( j ) = dimag( b( jrow, jcol ) )
534 200 CONTINUE
535 210 CONTINUE
536 CALL dgemm(
'T',
'N', nsize, nrhs, nsize, one,
537 $ rwork( u+st1 ), n, rwork( irwb ), nsize,
538 $ zero, rwork( irwib ), nsize )
539 jreal = irwrb - 1
540 jimag = irwib - 1
541 DO 230 jcol = 1, nrhs
542 DO 220 jrow = st, st + nsize - 1
543 jreal = jreal + 1
544 jimag = jimag + 1
545 b( jrow, jcol ) = dcmplx( rwork( jreal ),
546 $ rwork( jimag ) )
547 220 CONTINUE
548 230 CONTINUE
549
550 CALL zlacpy(
'A', nsize, nrhs, b( st, 1 ), ldb,
551 $ work( bx+st1 ), n )
552 ELSE
553
554
555
556 CALL dlasda( icmpq1, smlsiz, nsize, sqre, d( st ),
557 $ e( st ), rwork( u+st1 ), n, rwork( vt+st1 ),
558 $ iwork( k+st1 ), rwork( difl+st1 ),
559 $ rwork( difr+st1 ), rwork( z+st1 ),
560 $ rwork( poles+st1 ), iwork( givptr+st1 ),
561 $ iwork( givcol+st1 ), n, iwork( perm+st1 ),
562 $ rwork( givnum+st1 ), rwork( c+st1 ),
563 $ rwork( s+st1 ), rwork( nrwork ),
564 $ iwork( iwk ), info )
565 IF( info.NE.0 ) THEN
566 RETURN
567 END IF
568 bxst = bx + st1
569 CALL zlalsa( icmpq2, smlsiz, nsize, nrhs, b( st, 1 ),
570 $ ldb, work( bxst ), n, rwork( u+st1 ), n,
571 $ rwork( vt+st1 ), iwork( k+st1 ),
572 $ rwork( difl+st1 ), rwork( difr+st1 ),
573 $ rwork( z+st1 ), rwork( poles+st1 ),
574 $ iwork( givptr+st1 ), iwork( givcol+st1 ), n,
575 $ iwork( perm+st1 ), rwork( givnum+st1 ),
576 $ rwork( c+st1 ), rwork( s+st1 ),
577 $ rwork( nrwork ), iwork( iwk ), info )
578 IF( info.NE.0 ) THEN
579 RETURN
580 END IF
581 END IF
582 st = i + 1
583 END IF
584 240 CONTINUE
585
586
587
588 tol = rcnd*abs( d(
idamax( n, d, 1 ) ) )
589
590 DO 250 i = 1, n
591
592
593
594
595 IF( abs( d( i ) ).LE.tol ) THEN
596 CALL zlaset(
'A', 1, nrhs, czero, czero, work( bx+i-1 ), n )
597 ELSE
598 rank = rank + 1
599 CALL zlascl(
'G', 0, 0, d( i ), one, 1, nrhs,
600 $ work( bx+i-1 ), n, info )
601 END IF
602 d( i ) = abs( d( i ) )
603 250 CONTINUE
604
605
606
607 icmpq2 = 1
608 DO 320 i = 1, nsub
609 st = iwork( i )
610 st1 = st - 1
611 nsize = iwork( sizei+i-1 )
612 bxst = bx + st1
613 IF( nsize.EQ.1 ) THEN
614 CALL zcopy( nrhs, work( bxst ), n, b( st, 1 ), ldb )
615 ELSE IF( nsize.LE.smlsiz ) THEN
616
617
618
619
620
621
622
623
624 j = bxst - n - 1
625 jreal = irwb - 1
626 DO 270 jcol = 1, nrhs
627 j = j + n
628 DO 260 jrow = 1, nsize
629 jreal = jreal + 1
630 rwork( jreal ) = dble( work( j+jrow ) )
631 260 CONTINUE
632 270 CONTINUE
633 CALL dgemm(
'T',
'N', nsize, nrhs, nsize, one,
634 $ rwork( vt+st1 ), n, rwork( irwb ), nsize, zero,
635 $ rwork( irwrb ), nsize )
636 j = bxst - n - 1
637 jimag = irwb - 1
638 DO 290 jcol = 1, nrhs
639 j = j + n
640 DO 280 jrow = 1, nsize
641 jimag = jimag + 1
642 rwork( jimag ) = dimag( work( j+jrow ) )
643 280 CONTINUE
644 290 CONTINUE
645 CALL dgemm(
'T',
'N', nsize, nrhs, nsize, one,
646 $ rwork( vt+st1 ), n, rwork( irwb ), nsize, zero,
647 $ rwork( irwib ), nsize )
648 jreal = irwrb - 1
649 jimag = irwib - 1
650 DO 310 jcol = 1, nrhs
651 DO 300 jrow = st, st + nsize - 1
652 jreal = jreal + 1
653 jimag = jimag + 1
654 b( jrow, jcol ) = dcmplx( rwork( jreal ),
655 $ rwork( jimag ) )
656 300 CONTINUE
657 310 CONTINUE
658 ELSE
659 CALL zlalsa( icmpq2, smlsiz, nsize, nrhs, work( bxst ), n,
660 $ b( st, 1 ), ldb, rwork( u+st1 ), n,
661 $ rwork( vt+st1 ), iwork( k+st1 ),
662 $ rwork( difl+st1 ), rwork( difr+st1 ),
663 $ rwork( z+st1 ), rwork( poles+st1 ),
664 $ iwork( givptr+st1 ), iwork( givcol+st1 ), n,
665 $ iwork( perm+st1 ), rwork( givnum+st1 ),
666 $ rwork( c+st1 ), rwork( s+st1 ),
667 $ rwork( nrwork ), iwork( iwk ), info )
668 IF( info.NE.0 ) THEN
669 RETURN
670 END IF
671 END IF
672 320 CONTINUE
673
674
675
676 CALL dlascl(
'G', 0, 0, one, orgnrm, n, 1, d, n, info )
677 CALL dlasrt(
'D', n, d, info )
678 CALL zlascl(
'G', 0, 0, orgnrm, one, n, nrhs, b, ldb, info )
679
680 RETURN
681
682
683
subroutine xerbla(srname, info)
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
integer function idamax(n, dx, incx)
IDAMAX
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zlalsa(icompq, smlsiz, n, nrhs, b, ldb, bx, ldbx, u, ldu, vt, k, difl, difr, z, poles, givptr, givcol, ldgcol, perm, givnum, c, s, rwork, iwork, info)
ZLALSA computes the SVD of the coefficient matrix in compact form. Used by sgelsd.
double precision function dlamch(cmach)
DLAMCH
double precision function dlanst(norm, n, d, e)
DLANST returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
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 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 dlasda(icompq, smlsiz, n, sqre, d, e, u, ldu, vt, k, difl, difr, z, poles, givptr, givcol, ldgcol, perm, givnum, c, s, work, iwork, info)
DLASDA computes the singular value decomposition (SVD) of a real upper bidiagonal matrix with diagona...
subroutine dlasdq(uplo, sqre, n, ncvt, nru, ncc, d, e, vt, ldvt, u, ldu, c, ldc, work, info)
DLASDQ computes the SVD of a real bidiagonal matrix with diagonal d and off-diagonal e....
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 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 dlasrt(id, n, d, info)
DLASRT sorts numbers in increasing or decreasing order.
subroutine zdrot(n, zx, incx, zy, incy, c, s)
ZDROT