236
237
238
239
240
241
242 DOUBLE PRECISION EPS, SFMIN, TOL
243 INTEGER INFO, LDA, LDV, LWORK, M, MV, N, N1, NSWEEP
244 CHARACTER*1 JOBV
245
246
247 DOUBLE PRECISION A( LDA, * ), D( N ), SVA( N ), V( LDV, * ),
248 $ WORK( LWORK )
249
250
251
252
253
254 DOUBLE PRECISION ZERO, HALF, ONE
255 parameter( zero = 0.0d0, half = 0.5d0, one = 1.0d0 )
256
257
258 DOUBLE PRECISION AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG,
259 $ BIGTHETA, CS, LARGE, MXAAPQ, MXSINJ, ROOTBIG,
260 $ ROOTEPS, ROOTSFMIN, ROOTTOL, SMALL, SN, T,
261 $ TEMP1, THETA, THSIGN
262 INTEGER BLSKIP, EMPTSW, i, ibr, igl, IERR, IJBLSK,
263 $ ISWROT, jbc, jgl, KBL, MVL, NOTROT, nblc, nblr,
264 $ p, PSKIPPED, q, ROWSKIP, SWBAND
265 LOGICAL APPLV, ROTOK, RSVEC
266
267
268 DOUBLE PRECISION FASTR( 5 )
269
270
271 INTRINSIC dabs, max, dble, min, dsign, dsqrt
272
273
274 DOUBLE PRECISION DDOT, DNRM2
275 INTEGER IDAMAX
276 LOGICAL LSAME
278
279
282
283
284
285
286
287 applv =
lsame( jobv,
'A' )
288 rsvec =
lsame( jobv,
'V' )
289 IF( .NOT.( rsvec .OR. applv .OR.
lsame( jobv,
'N' ) ) )
THEN
290 info = -1
291 ELSE IF( m.LT.0 ) THEN
292 info = -2
293 ELSE IF( ( n.LT.0 ) .OR. ( n.GT.m ) ) THEN
294 info = -3
295 ELSE IF( n1.LT.0 ) THEN
296 info = -4
297 ELSE IF( lda.LT.m ) THEN
298 info = -6
299 ELSE IF( ( rsvec.OR.applv ) .AND. ( mv.LT.0 ) ) THEN
300 info = -9
301 ELSE IF( ( rsvec.AND.( ldv.LT.n ) ).OR.
302 $ ( applv.AND.( ldv.LT.mv ) ) ) THEN
303 info = -11
304 ELSE IF( tol.LE.eps ) THEN
305 info = -14
306 ELSE IF( nsweep.LT.0 ) THEN
307 info = -15
308 ELSE IF( lwork.LT.m ) THEN
309 info = -17
310 ELSE
311 info = 0
312 END IF
313
314
315 IF( info.NE.0 ) THEN
316 CALL xerbla(
'DGSVJ1', -info )
317 RETURN
318 END IF
319
320 IF( rsvec ) THEN
321 mvl = n
322 ELSE IF( applv ) THEN
323 mvl = mv
324 END IF
325 rsvec = rsvec .OR. applv
326
327 rooteps = dsqrt( eps )
328 rootsfmin = dsqrt( sfmin )
329 small = sfmin / eps
330 big = one / sfmin
331 rootbig = one / rootsfmin
332 large = big / dsqrt( dble( m*n ) )
333 bigtheta = one / rooteps
334 roottol = dsqrt( tol )
335
336
337
338
339
340 emptsw = n1*( n-n1 )
341 notrot = 0
342 fastr( 1 ) = zero
343
344
345
346 kbl = min( 8, n )
347 nblr = n1 / kbl
348 IF( ( nblr*kbl ).NE.n1 )nblr = nblr + 1
349
350
351
352 nblc = ( n-n1 ) / kbl
353 IF( ( nblc*kbl ).NE.( n-n1 ) )nblc = nblc + 1
354 blskip = ( kbl**2 ) + 1
355
356
357 rowskip = min( 5, kbl )
358
359 swband = 0
360
361
362
363
364
365
366
367
368
369
370
371
372
373 DO 1993 i = 1, nsweep
374
375
376 mxaapq = zero
377 mxsinj = zero
378 iswrot = 0
379
380 notrot = 0
381 pskipped = 0
382
383 DO 2000 ibr = 1, nblr
384
385 igl = ( ibr-1 )*kbl + 1
386
387
388
389
390
391 igl = ( ibr-1 )*kbl + 1
392
393 DO 2010 jbc = 1, nblc
394
395 jgl = n1 + ( jbc-1 )*kbl + 1
396
397
398
399 ijblsk = 0
400 DO 2100 p = igl, min( igl+kbl-1, n1 )
401
402 aapp = sva( p )
403
404 IF( aapp.GT.zero ) THEN
405
406 pskipped = 0
407
408 DO 2200 q = jgl, min( jgl+kbl-1, n )
409
410 aaqq = sva( q )
411
412 IF( aaqq.GT.zero ) THEN
413 aapp0 = aapp
414
415
416
417
418
419 IF( aaqq.GE.one ) THEN
420 IF( aapp.GE.aaqq ) THEN
421 rotok = ( small*aapp ).LE.aaqq
422 ELSE
423 rotok = ( small*aaqq ).LE.aapp
424 END IF
425 IF( aapp.LT.( big / aaqq ) ) THEN
426 aapq = (
ddot( m, a( 1, p ), 1, a( 1,
427 $ q ), 1 )*d( p )*d( q ) / aaqq )
428 $ / aapp
429 ELSE
430 CALL dcopy( m, a( 1, p ), 1, work, 1 )
431 CALL dlascl(
'G', 0, 0, aapp, d( p ),
432 $ m, 1, work, lda, ierr )
433 aapq =
ddot( m, work, 1, a( 1, q ),
434 $ 1 )*d( q ) / aaqq
435 END IF
436 ELSE
437 IF( aapp.GE.aaqq ) THEN
438 rotok = aapp.LE.( aaqq / small )
439 ELSE
440 rotok = aaqq.LE.( aapp / small )
441 END IF
442 IF( aapp.GT.( small / aaqq ) ) THEN
443 aapq = (
ddot( m, a( 1, p ), 1, a( 1,
444 $ q ), 1 )*d( p )*d( q ) / aaqq )
445 $ / aapp
446 ELSE
447 CALL dcopy( m, a( 1, q ), 1, work, 1 )
448 CALL dlascl(
'G', 0, 0, aaqq, d( q ),
449 $ m, 1, work, lda, ierr )
450 aapq =
ddot( m, work, 1, a( 1, p ),
451 $ 1 )*d( p ) / aapp
452 END IF
453 END IF
454
455 mxaapq = max( mxaapq, dabs( aapq ) )
456
457
458
459 IF( dabs( aapq ).GT.tol ) THEN
460 notrot = 0
461
462 pskipped = 0
463 iswrot = iswrot + 1
464
465 IF( rotok ) THEN
466
467 aqoap = aaqq / aapp
468 apoaq = aapp / aaqq
469 theta = -half*dabs(aqoap-apoaq) / aapq
470 IF( aaqq.GT.aapp0 )theta = -theta
471
472 IF( dabs( theta ).GT.bigtheta ) THEN
473 t = half / theta
474 fastr( 3 ) = t*d( p ) / d( q )
475 fastr( 4 ) = -t*d( q ) / d( p )
476 CALL drotm( m, a( 1, p ), 1,
477 $ a( 1, q ), 1, fastr )
478 IF( rsvec )
CALL drotm( mvl,
479 $ v( 1, p ), 1,
480 $ v( 1, q ), 1,
481 $ fastr )
482 sva( q ) = aaqq*dsqrt( max( zero,
483 $ one+t*apoaq*aapq ) )
484 aapp = aapp*dsqrt( max( zero,
485 $ one-t*aqoap*aapq ) )
486 mxsinj = max( mxsinj, dabs( t ) )
487 ELSE
488
489
490
491 thsign = -dsign( one, aapq )
492 IF( aaqq.GT.aapp0 )thsign = -thsign
493 t = one / ( theta+thsign*
494 $ dsqrt( one+theta*theta ) )
495 cs = dsqrt( one / ( one+t*t ) )
496 sn = t*cs
497 mxsinj = max( mxsinj, dabs( sn ) )
498 sva( q ) = aaqq*dsqrt( max( zero,
499 $ one+t*apoaq*aapq ) )
500 aapp = aapp*dsqrt( max( zero,
501 $ one-t*aqoap*aapq ) )
502
503 apoaq = d( p ) / d( q )
504 aqoap = d( q ) / d( p )
505 IF( d( p ).GE.one ) THEN
506
507 IF( d( q ).GE.one ) THEN
508 fastr( 3 ) = t*apoaq
509 fastr( 4 ) = -t*aqoap
510 d( p ) = d( p )*cs
511 d( q ) = d( q )*cs
512 CALL drotm( m, a( 1, p ), 1,
513 $ a( 1, q ), 1,
514 $ fastr )
515 IF( rsvec )
CALL drotm( mvl,
516 $ v( 1, p ), 1, v( 1, q ),
517 $ 1, fastr )
518 ELSE
519 CALL daxpy( m, -t*aqoap,
520 $ a( 1, q ), 1,
521 $ a( 1, p ), 1 )
522 CALL daxpy( m, cs*sn*apoaq,
523 $ a( 1, p ), 1,
524 $ a( 1, q ), 1 )
525 IF( rsvec ) THEN
526 CALL daxpy( mvl, -t*aqoap,
527 $ v( 1, q ), 1,
528 $ v( 1, p ), 1 )
530 $ cs*sn*apoaq,
531 $ v( 1, p ), 1,
532 $ v( 1, q ), 1 )
533 END IF
534 d( p ) = d( p )*cs
535 d( q ) = d( q ) / cs
536 END IF
537 ELSE
538 IF( d( q ).GE.one ) THEN
539 CALL daxpy( m, t*apoaq,
540 $ a( 1, p ), 1,
541 $ a( 1, q ), 1 )
542 CALL daxpy( m, -cs*sn*aqoap,
543 $ a( 1, q ), 1,
544 $ a( 1, p ), 1 )
545 IF( rsvec ) THEN
546 CALL daxpy( mvl, t*apoaq,
547 $ v( 1, p ), 1,
548 $ v( 1, q ), 1 )
550 $ -cs*sn*aqoap,
551 $ v( 1, q ), 1,
552 $ v( 1, p ), 1 )
553 END IF
554 d( p ) = d( p ) / cs
555 d( q ) = d( q )*cs
556 ELSE
557 IF( d( p ).GE.d( q ) ) THEN
558 CALL daxpy( m, -t*aqoap,
559 $ a( 1, q ), 1,
560 $ a( 1, p ), 1 )
561 CALL daxpy( m, cs*sn*apoaq,
562 $ a( 1, p ), 1,
563 $ a( 1, q ), 1 )
564 d( p ) = d( p )*cs
565 d( q ) = d( q ) / cs
566 IF( rsvec ) THEN
568 $ -t*aqoap,
569 $ v( 1, q ), 1,
570 $ v( 1, p ), 1 )
572 $ cs*sn*apoaq,
573 $ v( 1, p ), 1,
574 $ v( 1, q ), 1 )
575 END IF
576 ELSE
577 CALL daxpy( m, t*apoaq,
578 $ a( 1, p ), 1,
579 $ a( 1, q ), 1 )
581 $ -cs*sn*aqoap,
582 $ a( 1, q ), 1,
583 $ a( 1, p ), 1 )
584 d( p ) = d( p ) / cs
585 d( q ) = d( q )*cs
586 IF( rsvec ) THEN
588 $ t*apoaq, v( 1, p ),
589 $ 1, v( 1, q ), 1 )
591 $ -cs*sn*aqoap,
592 $ v( 1, q ), 1,
593 $ v( 1, p ), 1 )
594 END IF
595 END IF
596 END IF
597 END IF
598 END IF
599
600 ELSE
601 IF( aapp.GT.aaqq ) THEN
602 CALL dcopy( m, a( 1, p ), 1, work,
603 $ 1 )
604 CALL dlascl(
'G', 0, 0, aapp, one,
605 $ m, 1, work, lda, ierr )
606 CALL dlascl(
'G', 0, 0, aaqq, one,
607 $ m, 1, a( 1, q ), lda,
608 $ ierr )
609 temp1 = -aapq*d( p ) / d( q )
610 CALL daxpy( m, temp1, work, 1,
611 $ a( 1, q ), 1 )
612 CALL dlascl(
'G', 0, 0, one, aaqq,
613 $ m, 1, a( 1, q ), lda,
614 $ ierr )
615 sva( q ) = aaqq*dsqrt( max( zero,
616 $ one-aapq*aapq ) )
617 mxsinj = max( mxsinj, sfmin )
618 ELSE
619 CALL dcopy( m, a( 1, q ), 1, work,
620 $ 1 )
621 CALL dlascl(
'G', 0, 0, aaqq, one,
622 $ m, 1, work, lda, ierr )
623 CALL dlascl(
'G', 0, 0, aapp, one,
624 $ m, 1, a( 1, p ), lda,
625 $ ierr )
626 temp1 = -aapq*d( q ) / d( p )
627 CALL daxpy( m, temp1, work, 1,
628 $ a( 1, p ), 1 )
629 CALL dlascl(
'G', 0, 0, one, aapp,
630 $ m, 1, a( 1, p ), lda,
631 $ ierr )
632 sva( p ) = aapp*dsqrt( max( zero,
633 $ one-aapq*aapq ) )
634 mxsinj = max( mxsinj, sfmin )
635 END IF
636 END IF
637
638
639
640
641 IF( ( sva( q ) / aaqq )**2.LE.rooteps )
642 $ THEN
643 IF( ( aaqq.LT.rootbig ) .AND.
644 $ ( aaqq.GT.rootsfmin ) ) THEN
645 sva( q ) =
dnrm2( m, a( 1, q ), 1 )*
646 $ d( q )
647 ELSE
648 t = zero
649 aaqq = one
650 CALL dlassq( m, a( 1, q ), 1, t,
651 $ aaqq )
652 sva( q ) = t*dsqrt( aaqq )*d( q )
653 END IF
654 END IF
655 IF( ( aapp / aapp0 )**2.LE.rooteps ) THEN
656 IF( ( aapp.LT.rootbig ) .AND.
657 $ ( aapp.GT.rootsfmin ) ) THEN
658 aapp =
dnrm2( m, a( 1, p ), 1 )*
659 $ d( p )
660 ELSE
661 t = zero
662 aapp = one
663 CALL dlassq( m, a( 1, p ), 1, t,
664 $ aapp )
665 aapp = t*dsqrt( aapp )*d( p )
666 END IF
667 sva( p ) = aapp
668 END IF
669
670 ELSE
671 notrot = notrot + 1
672
673 pskipped = pskipped + 1
674 ijblsk = ijblsk + 1
675 END IF
676 ELSE
677 notrot = notrot + 1
678 pskipped = pskipped + 1
679 ijblsk = ijblsk + 1
680 END IF
681
682
683 IF( ( i.LE.swband ) .AND. ( ijblsk.GE.blskip ) )
684 $ THEN
685 sva( p ) = aapp
686 notrot = 0
687 GO TO 2011
688 END IF
689 IF( ( i.LE.swband ) .AND.
690 $ ( pskipped.GT.rowskip ) ) THEN
691 aapp = -aapp
692 notrot = 0
693 GO TO 2203
694 END IF
695
696
697 2200 CONTINUE
698
699 2203 CONTINUE
700
701 sva( p ) = aapp
702
703 ELSE
704 IF( aapp.EQ.zero )notrot = notrot +
705 $ min( jgl+kbl-1, n ) - jgl + 1
706 IF( aapp.LT.zero )notrot = 0
707
708 END IF
709
710 2100 CONTINUE
711
712 2010 CONTINUE
713
714 2011 CONTINUE
715
716 DO 2012 p = igl, min( igl+kbl-1, n )
717 sva( p ) = dabs( sva( p ) )
718 2012 CONTINUE
719
720 2000 CONTINUE
721
722
723
724 IF( ( sva( n ).LT.rootbig ) .AND. ( sva( n ).GT.rootsfmin ) )
725 $ THEN
726 sva( n ) =
dnrm2( m, a( 1, n ), 1 )*d( n )
727 ELSE
728 t = zero
729 aapp = one
730 CALL dlassq( m, a( 1, n ), 1, t, aapp )
731 sva( n ) = t*dsqrt( aapp )*d( n )
732 END IF
733
734
735
736 IF( ( i.LT.swband ) .AND. ( ( mxaapq.LE.roottol ) .OR.
737 $ ( iswrot.LE.n ) ) )swband = i
738
739 IF( ( i.GT.swband+1 ) .AND. ( mxaapq.LT.dble( n )*tol ) .AND.
740 $ ( dble( n )*mxaapq*mxsinj.LT.tol ) ) THEN
741 GO TO 1994
742 END IF
743
744
745 IF( notrot.GE.emptsw )GO TO 1994
746
747 1993 CONTINUE
748
749
750
751 info = nsweep - 1
752 GO TO 1995
753 1994 CONTINUE
754
755
756
757 info = 0
758
759 1995 CONTINUE
760
761
762
763 DO 5991 p = 1, n - 1
764 q =
idamax( n-p+1, sva( p ), 1 ) + p - 1
765 IF( p.NE.q ) THEN
766 temp1 = sva( p )
767 sva( p ) = sva( q )
768 sva( q ) = temp1
769 temp1 = d( p )
770 d( p ) = d( q )
771 d( q ) = temp1
772 CALL dswap( m, a( 1, p ), 1, a( 1, q ), 1 )
773 IF( rsvec )
CALL dswap( mvl, v( 1, p ), 1, v( 1, q ), 1 )
774 END IF
775 5991 CONTINUE
776
777 RETURN
778
779
780
subroutine xerbla(srname, info)
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
double precision function ddot(n, dx, incx, dy, incy)
DDOT
integer function idamax(n, dx, incx)
IDAMAX
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 dlassq(n, x, incx, scale, sumsq)
DLASSQ updates a sum of squares represented in scaled form.
logical function lsame(ca, cb)
LSAME
real(wp) function dnrm2(n, x, incx)
DNRM2
subroutine drotm(n, dx, incx, dy, incy, dparam)
DROTM
subroutine dswap(n, dx, incx, dy, incy)
DSWAP