LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dgsvj0()

subroutine dgsvj0 ( character*1 jobv,
integer m,
integer n,
double precision, dimension( lda, * ) a,
integer lda,
double precision, dimension( n ) d,
double precision, dimension( n ) sva,
integer mv,
double precision, dimension( ldv, * ) v,
integer ldv,
double precision eps,
double precision sfmin,
double precision tol,
integer nsweep,
double precision, dimension( lwork ) work,
integer lwork,
integer info )

DGSVJ0 pre-processor for the routine dgesvj.

Download DGSVJ0 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> DGSVJ0 is called from DGESVJ as a pre-processor and that is its main
!> purpose. It applies Jacobi rotations in the same way as DGESVJ does, but
!> it does not check convergence (stopping criterion). Few tuning
!> parameters (marked by [TP]) are available for the implementer.
!>
Parameters
[in]JOBV
!>          JOBV is CHARACTER*1
!>          Specifies whether the output from this procedure is used
!>          to compute the matrix V:
!>          = 'V': the product of the Jacobi rotations is accumulated
!>                 by postmultiplying the N-by-N array V.
!>                (See the description of V.)
!>          = 'A': the product of the Jacobi rotations is accumulated
!>                 by postmultiplying the MV-by-N array V.
!>                (See the descriptions of MV and V.)
!>          = 'N': the Jacobi rotations are not accumulated.
!>
[in]M
!>          M is INTEGER
!>          The number of rows of the input matrix A.  M >= 0.
!>
[in]N
!>          N is INTEGER
!>          The number of columns of the input matrix A.
!>          M >= N >= 0.
!>
[in,out]A
!>          A is DOUBLE PRECISION array, dimension (LDA,N)
!>          On entry, M-by-N matrix A, such that A*diag(D) represents
!>          the input matrix.
!>          On exit,
!>          A_onexit * D_onexit represents the input matrix A*diag(D)
!>          post-multiplied by a sequence of Jacobi rotations, where the
!>          rotation threshold and the total number of sweeps are given in
!>          TOL and NSWEEP, respectively.
!>          (See the descriptions of D, TOL and NSWEEP.)
!>
[in]LDA
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,M).
!>
[in,out]D
!>          D is DOUBLE PRECISION array, dimension (N)
!>          The array D accumulates the scaling factors from the fast scaled
!>          Jacobi rotations.
!>          On entry, A*diag(D) represents the input matrix.
!>          On exit, A_onexit*diag(D_onexit) represents the input matrix
!>          post-multiplied by a sequence of Jacobi rotations, where the
!>          rotation threshold and the total number of sweeps are given in
!>          TOL and NSWEEP, respectively.
!>          (See the descriptions of A, TOL and NSWEEP.)
!>
[in,out]SVA
!>          SVA is DOUBLE PRECISION array, dimension (N)
!>          On entry, SVA contains the Euclidean norms of the columns of
!>          the matrix A*diag(D).
!>          On exit, SVA contains the Euclidean norms of the columns of
!>          the matrix onexit*diag(D_onexit).
!>
[in]MV
!>          MV is INTEGER
!>          If JOBV = 'A', then MV rows of V are post-multiplied by a
!>                           sequence of Jacobi rotations.
!>          If JOBV = 'N',   then MV is not referenced.
!>
[in,out]V
!>          V is DOUBLE PRECISION array, dimension (LDV,N)
!>          If JOBV = 'V' then N rows of V are post-multiplied by a
!>                           sequence of Jacobi rotations.
!>          If JOBV = 'A' then MV rows of V are post-multiplied by a
!>                           sequence of Jacobi rotations.
!>          If JOBV = 'N',   then V is not referenced.
!>
[in]LDV
!>          LDV is INTEGER
!>          The leading dimension of the array V,  LDV >= 1.
!>          If JOBV = 'V', LDV >= N.
!>          If JOBV = 'A', LDV >= MV.
!>
[in]EPS
!>          EPS is DOUBLE PRECISION
!>          EPS = DLAMCH('Epsilon')
!>
[in]SFMIN
!>          SFMIN is DOUBLE PRECISION
!>          SFMIN = DLAMCH('Safe Minimum')
!>
[in]TOL
!>          TOL is DOUBLE PRECISION
!>          TOL is the threshold for Jacobi rotations. For a pair
!>          A(:,p), A(:,q) of pivot columns, the Jacobi rotation is
!>          applied only if DABS(COS(angle(A(:,p),A(:,q)))) > TOL.
!>
[in]NSWEEP
!>          NSWEEP is INTEGER
!>          NSWEEP is the number of sweeps of Jacobi rotations to be
!>          performed.
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array, dimension (LWORK)
!>
[in]LWORK
!>          LWORK is INTEGER
!>          LWORK is the dimension of WORK. LWORK >= M.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit.
!>          < 0:  if INFO = -i, then the i-th argument had an illegal value
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Further Details:
DGSVJ0 is used just to enable DGESVJ to call a simplified version of itself to work on a submatrix of the original matrix.
Contributors:
Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
Bugs, Examples and Comments:
Please report all bugs and send interesting test examples and comments to drmac.nosp@m.@mat.nosp@m.h.hr. Thank you.

Definition at line 214 of file dgsvj0.f.

216*
217* -- LAPACK computational routine --
218* -- LAPACK is a software package provided by Univ. of Tennessee, --
219* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
220*
221* .. Scalar Arguments ..
222 INTEGER INFO, LDA, LDV, LWORK, M, MV, N, NSWEEP
223 DOUBLE PRECISION EPS, SFMIN, TOL
224 CHARACTER*1 JOBV
225* ..
226* .. Array Arguments ..
227 DOUBLE PRECISION A( LDA, * ), SVA( N ), D( N ), V( LDV, * ),
228 $ WORK( LWORK )
229* ..
230*
231* =====================================================================
232*
233* .. Local Parameters ..
234 DOUBLE PRECISION ZERO, HALF, ONE
235 parameter( zero = 0.0d0, half = 0.5d0, one = 1.0d0)
236* ..
237* .. Local Scalars ..
238 DOUBLE PRECISION AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG,
239 $ BIGTHETA, CS, MXAAPQ, MXSINJ, ROOTBIG, ROOTEPS,
240 $ ROOTSFMIN, ROOTTOL, SMALL, SN, T, TEMP1, THETA,
241 $ THSIGN
242 INTEGER BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1,
243 $ ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, NBL,
244 $ NOTROT, p, PSKIPPED, q, ROWSKIP, SWBAND
245 LOGICAL APPLV, ROTOK, RSVEC
246* ..
247* .. Local Arrays ..
248 DOUBLE PRECISION FASTR( 5 )
249* ..
250* .. Intrinsic Functions ..
251 INTRINSIC dabs, max, dble, min, dsign, dsqrt
252* ..
253* .. External Functions ..
254 DOUBLE PRECISION DDOT, DNRM2
255 INTEGER IDAMAX
256 LOGICAL LSAME
257 EXTERNAL idamax, lsame, ddot, dnrm2
258* ..
259* .. External Subroutines ..
260 EXTERNAL daxpy, dcopy, dlascl, dlassq, drotm,
261 $ dswap,
262 $ xerbla
263* ..
264* .. Executable Statements ..
265*
266* Test the input parameters.
267*
268 applv = lsame( jobv, 'A' )
269 rsvec = lsame( jobv, 'V' )
270 IF( .NOT.( rsvec .OR. applv .OR. lsame( jobv, 'N' ) ) ) THEN
271 info = -1
272 ELSE IF( m.LT.0 ) THEN
273 info = -2
274 ELSE IF( ( n.LT.0 ) .OR. ( n.GT.m ) ) THEN
275 info = -3
276 ELSE IF( lda.LT.m ) THEN
277 info = -5
278 ELSE IF( ( rsvec.OR.applv ) .AND. ( mv.LT.0 ) ) THEN
279 info = -8
280 ELSE IF( ( rsvec.AND.( ldv.LT.n ) ).OR.
281 $ ( applv.AND.( ldv.LT.mv ) ) ) THEN
282 info = -10
283 ELSE IF( tol.LE.eps ) THEN
284 info = -13
285 ELSE IF( nsweep.LT.0 ) THEN
286 info = -14
287 ELSE IF( lwork.LT.m ) THEN
288 info = -16
289 ELSE
290 info = 0
291 END IF
292*
293* #:(
294 IF( info.NE.0 ) THEN
295 CALL xerbla( 'DGSVJ0', -info )
296 RETURN
297 END IF
298*
299 IF( rsvec ) THEN
300 mvl = n
301 ELSE IF( applv ) THEN
302 mvl = mv
303 END IF
304 rsvec = rsvec .OR. applv
305
306 rooteps = dsqrt( eps )
307 rootsfmin = dsqrt( sfmin )
308 small = sfmin / eps
309 big = one / sfmin
310 rootbig = one / rootsfmin
311 bigtheta = one / rooteps
312 roottol = dsqrt( tol )
313*
314* -#- Row-cyclic Jacobi SVD algorithm with column pivoting -#-
315*
316 emptsw = ( n*( n-1 ) ) / 2
317 notrot = 0
318 fastr( 1 ) = zero
319*
320* -#- Row-cyclic pivot strategy with de Rijk's pivoting -#-
321*
322
323 swband = 0
324*[TP] SWBAND is a tuning parameter. It is meaningful and effective
325* if SGESVJ is used as a computational routine in the preconditioned
326* Jacobi SVD algorithm SGESVJ. For sweeps i=1:SWBAND the procedure
327* ......
328
329 kbl = min( 8, n )
330*[TP] KBL is a tuning parameter that defines the tile size in the
331* tiling of the p-q loops of pivot pairs. In general, an optimal
332* value of KBL depends on the matrix dimensions and on the
333* parameters of the computer's memory.
334*
335 nbl = n / kbl
336 IF( ( nbl*kbl ).NE.n )nbl = nbl + 1
337
338 blskip = ( kbl**2 ) + 1
339*[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
340
341 rowskip = min( 5, kbl )
342*[TP] ROWSKIP is a tuning parameter.
343
344 lkahead = 1
345*[TP] LKAHEAD is a tuning parameter.
346 swband = 0
347 pskipped = 0
348*
349 DO 1993 i = 1, nsweep
350* .. go go go ...
351*
352 mxaapq = zero
353 mxsinj = zero
354 iswrot = 0
355*
356 notrot = 0
357 pskipped = 0
358*
359 DO 2000 ibr = 1, nbl
360
361 igl = ( ibr-1 )*kbl + 1
362*
363 DO 1002 ir1 = 0, min( lkahead, nbl-ibr )
364*
365 igl = igl + ir1*kbl
366*
367 DO 2001 p = igl, min( igl+kbl-1, n-1 )
368
369* .. de Rijk's pivoting
370 q = idamax( n-p+1, sva( p ), 1 ) + p - 1
371 IF( p.NE.q ) THEN
372 CALL dswap( m, a( 1, p ), 1, a( 1, q ), 1 )
373 IF( rsvec )CALL dswap( mvl, v( 1, p ), 1,
374 $ v( 1, q ), 1 )
375 temp1 = sva( p )
376 sva( p ) = sva( q )
377 sva( q ) = temp1
378 temp1 = d( p )
379 d( p ) = d( q )
380 d( q ) = temp1
381 END IF
382*
383 IF( ir1.EQ.0 ) THEN
384*
385* Column norms are periodically updated by explicit
386* norm computation.
387* Caveat:
388* Some BLAS implementations compute DNRM2(M,A(1,p),1)
389* as DSQRT(DDOT(M,A(1,p),1,A(1,p),1)), which may result in
390* overflow for ||A(:,p)||_2 > DSQRT(overflow_threshold), and
391* underflow for ||A(:,p)||_2 < DSQRT(underflow_threshold).
392* Hence, DNRM2 cannot be trusted, not even in the case when
393* the true norm is far from the under(over)flow boundaries.
394* If properly implemented DNRM2 is available, the IF-THEN-ELSE
395* below should read "AAPP = DNRM2( M, A(1,p), 1 ) * D(p)".
396*
397 IF( ( sva( p ).LT.rootbig ) .AND.
398 $ ( sva( p ).GT.rootsfmin ) ) THEN
399 sva( p ) = dnrm2( m, a( 1, p ), 1 )*d( p )
400 ELSE
401 temp1 = zero
402 aapp = one
403 CALL dlassq( m, a( 1, p ), 1, temp1, aapp )
404 sva( p ) = temp1*dsqrt( aapp )*d( p )
405 END IF
406 aapp = sva( p )
407 ELSE
408 aapp = sva( p )
409 END IF
410
411*
412 IF( aapp.GT.zero ) THEN
413*
414 pskipped = 0
415*
416 DO 2002 q = p + 1, min( igl+kbl-1, n )
417*
418 aaqq = sva( q )
419
420 IF( aaqq.GT.zero ) THEN
421*
422 aapp0 = aapp
423 IF( aaqq.GE.one ) THEN
424 rotok = ( small*aapp ).LE.aaqq
425 IF( aapp.LT.( big / aaqq ) ) THEN
426 aapq = ( ddot( m, a( 1, p ), 1,
427 $ a( 1,
428 $ q ), 1 )*d( p )*d( q ) / aaqq )
429 $ / aapp
430 ELSE
431 CALL dcopy( m, a( 1, p ), 1, work,
432 $ 1 )
433 CALL dlascl( 'G', 0, 0, aapp,
434 $ d( p ),
435 $ m, 1, work, lda, ierr )
436 aapq = ddot( m, work, 1, a( 1, q ),
437 $ 1 )*d( q ) / aaqq
438 END IF
439 ELSE
440 rotok = aapp.LE.( aaqq / small )
441 IF( aapp.GT.( small / aaqq ) ) THEN
442 aapq = ( ddot( m, a( 1, p ), 1,
443 $ a( 1,
444 $ q ), 1 )*d( p )*d( q ) / aaqq )
445 $ / aapp
446 ELSE
447 CALL dcopy( m, a( 1, q ), 1, work,
448 $ 1 )
449 CALL dlascl( 'G', 0, 0, aaqq,
450 $ d( q ),
451 $ m, 1, work, lda, ierr )
452 aapq = ddot( m, work, 1, a( 1, p ),
453 $ 1 )*d( p ) / aapp
454 END IF
455 END IF
456*
457 mxaapq = max( mxaapq, dabs( aapq ) )
458*
459* TO rotate or NOT to rotate, THAT is the question ...
460*
461 IF( dabs( aapq ).GT.tol ) THEN
462*
463* .. rotate
464* ROTATED = ROTATED + ONE
465*
466 IF( ir1.EQ.0 ) THEN
467 notrot = 0
468 pskipped = 0
469 iswrot = iswrot + 1
470 END IF
471*
472 IF( rotok ) THEN
473*
474 aqoap = aaqq / aapp
475 apoaq = aapp / aaqq
476 theta = -half*dabs( aqoap-apoaq )/aapq
477*
478 IF( dabs( theta ).GT.bigtheta ) THEN
479*
480 t = half / theta
481 fastr( 3 ) = t*d( p ) / d( q )
482 fastr( 4 ) = -t*d( q ) / d( p )
483 CALL drotm( m, a( 1, p ), 1,
484 $ a( 1, q ), 1, fastr )
485 IF( rsvec )CALL drotm( mvl,
486 $ v( 1, p ), 1,
487 $ v( 1, q ), 1,
488 $ fastr )
489 sva( q ) = aaqq*dsqrt( max( zero,
490 $ one+t*apoaq*aapq ) )
491 aapp = aapp*dsqrt( max( zero,
492 $ one-t*aqoap*aapq ) )
493 mxsinj = max( mxsinj, dabs( t ) )
494*
495 ELSE
496*
497* .. choose correct signum for THETA and rotate
498*
499 thsign = -dsign( one, aapq )
500 t = one / ( theta+thsign*
501 $ dsqrt( one+theta*theta ) )
502 cs = dsqrt( one / ( one+t*t ) )
503 sn = t*cs
504*
505 mxsinj = max( mxsinj, dabs( sn ) )
506 sva( q ) = aaqq*dsqrt( max( zero,
507 $ one+t*apoaq*aapq ) )
508 aapp = aapp*dsqrt( max( zero,
509 $ one-t*aqoap*aapq ) )
510*
511 apoaq = d( p ) / d( q )
512 aqoap = d( q ) / d( p )
513 IF( d( p ).GE.one ) THEN
514 IF( d( q ).GE.one ) THEN
515 fastr( 3 ) = t*apoaq
516 fastr( 4 ) = -t*aqoap
517 d( p ) = d( p )*cs
518 d( q ) = d( q )*cs
519 CALL drotm( m, a( 1, p ),
520 $ 1,
521 $ a( 1, q ), 1,
522 $ fastr )
523 IF( rsvec )CALL drotm( mvl,
524 $ v( 1, p ), 1, v( 1, q ),
525 $ 1, fastr )
526 ELSE
527 CALL daxpy( m, -t*aqoap,
528 $ a( 1, q ), 1,
529 $ a( 1, p ), 1 )
530 CALL daxpy( m, cs*sn*apoaq,
531 $ a( 1, p ), 1,
532 $ a( 1, q ), 1 )
533 d( p ) = d( p )*cs
534 d( q ) = d( q ) / cs
535 IF( rsvec ) THEN
536 CALL daxpy( mvl,
537 $ -t*aqoap,
538 $ v( 1, q ), 1,
539 $ v( 1, p ), 1 )
540 CALL daxpy( mvl,
541 $ cs*sn*apoaq,
542 $ v( 1, p ), 1,
543 $ v( 1, q ), 1 )
544 END IF
545 END IF
546 ELSE
547 IF( d( q ).GE.one ) THEN
548 CALL daxpy( m, t*apoaq,
549 $ a( 1, p ), 1,
550 $ a( 1, q ), 1 )
551 CALL daxpy( m,
552 $ -cs*sn*aqoap,
553 $ a( 1, q ), 1,
554 $ a( 1, p ), 1 )
555 d( p ) = d( p ) / cs
556 d( q ) = d( q )*cs
557 IF( rsvec ) THEN
558 CALL daxpy( mvl,
559 $ t*apoaq,
560 $ v( 1, p ), 1,
561 $ v( 1, q ), 1 )
562 CALL daxpy( mvl,
563 $ -cs*sn*aqoap,
564 $ v( 1, q ), 1,
565 $ v( 1, p ), 1 )
566 END IF
567 ELSE
568 IF( d( p ).GE.d( q ) ) THEN
569 CALL daxpy( m, -t*aqoap,
570 $ a( 1, q ), 1,
571 $ a( 1, p ), 1 )
572 CALL daxpy( m,
573 $ cs*sn*apoaq,
574 $ a( 1, p ), 1,
575 $ a( 1, q ), 1 )
576 d( p ) = d( p )*cs
577 d( q ) = d( q ) / cs
578 IF( rsvec ) THEN
579 CALL daxpy( mvl,
580 $ -t*aqoap,
581 $ v( 1, q ), 1,
582 $ v( 1, p ), 1 )
583 CALL daxpy( mvl,
584 $ cs*sn*apoaq,
585 $ v( 1, p ), 1,
586 $ v( 1, q ), 1 )
587 END IF
588 ELSE
589 CALL daxpy( m, t*apoaq,
590 $ a( 1, p ), 1,
591 $ a( 1, q ), 1 )
592 CALL daxpy( m,
593 $ -cs*sn*aqoap,
594 $ a( 1, q ), 1,
595 $ a( 1, p ), 1 )
596 d( p ) = d( p ) / cs
597 d( q ) = d( q )*cs
598 IF( rsvec ) THEN
599 CALL daxpy( mvl,
600 $ t*apoaq, v( 1, p ),
601 $ 1, v( 1, q ), 1 )
602 CALL daxpy( mvl,
603 $ -cs*sn*aqoap,
604 $ v( 1, q ), 1,
605 $ v( 1, p ), 1 )
606 END IF
607 END IF
608 END IF
609 END IF
610 END IF
611*
612 ELSE
613* .. have to use modified Gram-Schmidt like transformation
614 CALL dcopy( m, a( 1, p ), 1, work,
615 $ 1 )
616 CALL dlascl( 'G', 0, 0, aapp, one,
617 $ m,
618 $ 1, work, lda, ierr )
619 CALL dlascl( 'G', 0, 0, aaqq, one,
620 $ m,
621 $ 1, a( 1, q ), lda, ierr )
622 temp1 = -aapq*d( p ) / d( q )
623 CALL daxpy( m, temp1, work, 1,
624 $ a( 1, q ), 1 )
625 CALL dlascl( 'G', 0, 0, one, aaqq,
626 $ m,
627 $ 1, a( 1, q ), lda, ierr )
628 sva( q ) = aaqq*dsqrt( max( zero,
629 $ one-aapq*aapq ) )
630 mxsinj = max( mxsinj, sfmin )
631 END IF
632* END IF ROTOK THEN ... ELSE
633*
634* In the case of cancellation in updating SVA(q), SVA(p)
635* recompute SVA(q), SVA(p).
636 IF( ( sva( q ) / aaqq )**2.LE.rooteps )
637 $ THEN
638 IF( ( aaqq.LT.rootbig ) .AND.
639 $ ( aaqq.GT.rootsfmin ) ) THEN
640 sva( q ) = dnrm2( m, a( 1, q ),
641 $ 1 )*
642 $ d( q )
643 ELSE
644 t = zero
645 aaqq = one
646 CALL dlassq( m, a( 1, q ), 1, t,
647 $ aaqq )
648 sva( q ) = t*dsqrt( aaqq )*d( q )
649 END IF
650 END IF
651 IF( ( aapp / aapp0 ).LE.rooteps ) THEN
652 IF( ( aapp.LT.rootbig ) .AND.
653 $ ( aapp.GT.rootsfmin ) ) THEN
654 aapp = dnrm2( m, a( 1, p ), 1 )*
655 $ d( p )
656 ELSE
657 t = zero
658 aapp = one
659 CALL dlassq( m, a( 1, p ), 1, t,
660 $ aapp )
661 aapp = t*dsqrt( aapp )*d( p )
662 END IF
663 sva( p ) = aapp
664 END IF
665*
666 ELSE
667* A(:,p) and A(:,q) already numerically orthogonal
668 IF( ir1.EQ.0 )notrot = notrot + 1
669 pskipped = pskipped + 1
670 END IF
671 ELSE
672* A(:,q) is zero column
673 IF( ir1.EQ.0 )notrot = notrot + 1
674 pskipped = pskipped + 1
675 END IF
676*
677 IF( ( i.LE.swband ) .AND.
678 $ ( pskipped.GT.rowskip ) ) THEN
679 IF( ir1.EQ.0 )aapp = -aapp
680 notrot = 0
681 GO TO 2103
682 END IF
683*
684 2002 CONTINUE
685* END q-LOOP
686*
687 2103 CONTINUE
688* bailed out of q-loop
689
690 sva( p ) = aapp
691
692 ELSE
693 sva( p ) = aapp
694 IF( ( ir1.EQ.0 ) .AND. ( aapp.EQ.zero ) )
695 $ notrot = notrot + min( igl+kbl-1, n ) - p
696 END IF
697*
698 2001 CONTINUE
699* end of the p-loop
700* end of doing the block ( ibr, ibr )
701 1002 CONTINUE
702* end of ir1-loop
703*
704*........................................................
705* ... go to the off diagonal blocks
706*
707 igl = ( ibr-1 )*kbl + 1
708*
709 DO 2010 jbc = ibr + 1, nbl
710*
711 jgl = ( jbc-1 )*kbl + 1
712*
713* doing the block at ( ibr, jbc )
714*
715 ijblsk = 0
716 DO 2100 p = igl, min( igl+kbl-1, n )
717*
718 aapp = sva( p )
719*
720 IF( aapp.GT.zero ) THEN
721*
722 pskipped = 0
723*
724 DO 2200 q = jgl, min( jgl+kbl-1, n )
725*
726 aaqq = sva( q )
727*
728 IF( aaqq.GT.zero ) THEN
729 aapp0 = aapp
730*
731* -#- M x 2 Jacobi SVD -#-
732*
733* -#- Safe Gram matrix computation -#-
734*
735 IF( aaqq.GE.one ) THEN
736 IF( aapp.GE.aaqq ) THEN
737 rotok = ( small*aapp ).LE.aaqq
738 ELSE
739 rotok = ( small*aaqq ).LE.aapp
740 END IF
741 IF( aapp.LT.( big / aaqq ) ) THEN
742 aapq = ( ddot( m, a( 1, p ), 1,
743 $ a( 1,
744 $ q ), 1 )*d( p )*d( q ) / aaqq )
745 $ / aapp
746 ELSE
747 CALL dcopy( m, a( 1, p ), 1, work,
748 $ 1 )
749 CALL dlascl( 'G', 0, 0, aapp,
750 $ d( p ),
751 $ m, 1, work, lda, ierr )
752 aapq = ddot( m, work, 1, a( 1, q ),
753 $ 1 )*d( q ) / aaqq
754 END IF
755 ELSE
756 IF( aapp.GE.aaqq ) THEN
757 rotok = aapp.LE.( aaqq / small )
758 ELSE
759 rotok = aaqq.LE.( aapp / small )
760 END IF
761 IF( aapp.GT.( small / aaqq ) ) THEN
762 aapq = ( ddot( m, a( 1, p ), 1,
763 $ a( 1,
764 $ q ), 1 )*d( p )*d( q ) / aaqq )
765 $ / aapp
766 ELSE
767 CALL dcopy( m, a( 1, q ), 1, work,
768 $ 1 )
769 CALL dlascl( 'G', 0, 0, aaqq,
770 $ d( q ),
771 $ m, 1, work, lda, ierr )
772 aapq = ddot( m, work, 1, a( 1, p ),
773 $ 1 )*d( p ) / aapp
774 END IF
775 END IF
776*
777 mxaapq = max( mxaapq, dabs( aapq ) )
778*
779* TO rotate or NOT to rotate, THAT is the question ...
780*
781 IF( dabs( aapq ).GT.tol ) THEN
782 notrot = 0
783* ROTATED = ROTATED + 1
784 pskipped = 0
785 iswrot = iswrot + 1
786*
787 IF( rotok ) THEN
788*
789 aqoap = aaqq / aapp
790 apoaq = aapp / aaqq
791 theta = -half*dabs( aqoap-apoaq )/aapq
792 IF( aaqq.GT.aapp0 )theta = -theta
793*
794 IF( dabs( theta ).GT.bigtheta ) THEN
795 t = half / theta
796 fastr( 3 ) = t*d( p ) / d( q )
797 fastr( 4 ) = -t*d( q ) / d( p )
798 CALL drotm( m, a( 1, p ), 1,
799 $ a( 1, q ), 1, fastr )
800 IF( rsvec )CALL drotm( mvl,
801 $ v( 1, p ), 1,
802 $ v( 1, q ), 1,
803 $ fastr )
804 sva( q ) = aaqq*dsqrt( max( zero,
805 $ one+t*apoaq*aapq ) )
806 aapp = aapp*dsqrt( max( zero,
807 $ one-t*aqoap*aapq ) )
808 mxsinj = max( mxsinj, dabs( t ) )
809 ELSE
810*
811* .. choose correct signum for THETA and rotate
812*
813 thsign = -dsign( one, aapq )
814 IF( aaqq.GT.aapp0 )thsign = -thsign
815 t = one / ( theta+thsign*
816 $ dsqrt( one+theta*theta ) )
817 cs = dsqrt( one / ( one+t*t ) )
818 sn = t*cs
819 mxsinj = max( mxsinj, dabs( sn ) )
820 sva( q ) = aaqq*dsqrt( max( zero,
821 $ one+t*apoaq*aapq ) )
822 aapp = aapp*dsqrt( max( zero,
823 $ one-t*aqoap*aapq ) )
824*
825 apoaq = d( p ) / d( q )
826 aqoap = d( q ) / d( p )
827 IF( d( p ).GE.one ) THEN
828*
829 IF( d( q ).GE.one ) THEN
830 fastr( 3 ) = t*apoaq
831 fastr( 4 ) = -t*aqoap
832 d( p ) = d( p )*cs
833 d( q ) = d( q )*cs
834 CALL drotm( m, a( 1, p ),
835 $ 1,
836 $ a( 1, q ), 1,
837 $ fastr )
838 IF( rsvec )CALL drotm( mvl,
839 $ v( 1, p ), 1, v( 1, q ),
840 $ 1, fastr )
841 ELSE
842 CALL daxpy( m, -t*aqoap,
843 $ a( 1, q ), 1,
844 $ a( 1, p ), 1 )
845 CALL daxpy( m, cs*sn*apoaq,
846 $ a( 1, p ), 1,
847 $ a( 1, q ), 1 )
848 IF( rsvec ) THEN
849 CALL daxpy( mvl,
850 $ -t*aqoap,
851 $ v( 1, q ), 1,
852 $ v( 1, p ), 1 )
853 CALL daxpy( mvl,
854 $ cs*sn*apoaq,
855 $ v( 1, p ), 1,
856 $ v( 1, q ), 1 )
857 END IF
858 d( p ) = d( p )*cs
859 d( q ) = d( q ) / cs
860 END IF
861 ELSE
862 IF( d( q ).GE.one ) THEN
863 CALL daxpy( m, t*apoaq,
864 $ a( 1, p ), 1,
865 $ a( 1, q ), 1 )
866 CALL daxpy( m,
867 $ -cs*sn*aqoap,
868 $ a( 1, q ), 1,
869 $ a( 1, p ), 1 )
870 IF( rsvec ) THEN
871 CALL daxpy( mvl,
872 $ t*apoaq,
873 $ v( 1, p ), 1,
874 $ v( 1, q ), 1 )
875 CALL daxpy( mvl,
876 $ -cs*sn*aqoap,
877 $ v( 1, q ), 1,
878 $ v( 1, p ), 1 )
879 END IF
880 d( p ) = d( p ) / cs
881 d( q ) = d( q )*cs
882 ELSE
883 IF( d( p ).GE.d( q ) ) THEN
884 CALL daxpy( m, -t*aqoap,
885 $ a( 1, q ), 1,
886 $ a( 1, p ), 1 )
887 CALL daxpy( m,
888 $ cs*sn*apoaq,
889 $ a( 1, p ), 1,
890 $ a( 1, q ), 1 )
891 d( p ) = d( p )*cs
892 d( q ) = d( q ) / cs
893 IF( rsvec ) THEN
894 CALL daxpy( mvl,
895 $ -t*aqoap,
896 $ v( 1, q ), 1,
897 $ v( 1, p ), 1 )
898 CALL daxpy( mvl,
899 $ cs*sn*apoaq,
900 $ v( 1, p ), 1,
901 $ v( 1, q ), 1 )
902 END IF
903 ELSE
904 CALL daxpy( m, t*apoaq,
905 $ a( 1, p ), 1,
906 $ a( 1, q ), 1 )
907 CALL daxpy( m,
908 $ -cs*sn*aqoap,
909 $ a( 1, q ), 1,
910 $ a( 1, p ), 1 )
911 d( p ) = d( p ) / cs
912 d( q ) = d( q )*cs
913 IF( rsvec ) THEN
914 CALL daxpy( mvl,
915 $ t*apoaq, v( 1, p ),
916 $ 1, v( 1, q ), 1 )
917 CALL daxpy( mvl,
918 $ -cs*sn*aqoap,
919 $ v( 1, q ), 1,
920 $ v( 1, p ), 1 )
921 END IF
922 END IF
923 END IF
924 END IF
925 END IF
926*
927 ELSE
928 IF( aapp.GT.aaqq ) THEN
929 CALL dcopy( m, a( 1, p ), 1,
930 $ work,
931 $ 1 )
932 CALL dlascl( 'G', 0, 0, aapp,
933 $ one,
934 $ m, 1, work, lda, ierr )
935 CALL dlascl( 'G', 0, 0, aaqq,
936 $ one,
937 $ m, 1, a( 1, q ), lda,
938 $ ierr )
939 temp1 = -aapq*d( p ) / d( q )
940 CALL daxpy( m, temp1, work, 1,
941 $ a( 1, q ), 1 )
942 CALL dlascl( 'G', 0, 0, one,
943 $ aaqq,
944 $ m, 1, a( 1, q ), lda,
945 $ ierr )
946 sva( q ) = aaqq*dsqrt( max( zero,
947 $ one-aapq*aapq ) )
948 mxsinj = max( mxsinj, sfmin )
949 ELSE
950 CALL dcopy( m, a( 1, q ), 1,
951 $ work,
952 $ 1 )
953 CALL dlascl( 'G', 0, 0, aaqq,
954 $ one,
955 $ m, 1, work, lda, ierr )
956 CALL dlascl( 'G', 0, 0, aapp,
957 $ one,
958 $ m, 1, a( 1, p ), lda,
959 $ ierr )
960 temp1 = -aapq*d( q ) / d( p )
961 CALL daxpy( m, temp1, work, 1,
962 $ a( 1, p ), 1 )
963 CALL dlascl( 'G', 0, 0, one,
964 $ aapp,
965 $ m, 1, a( 1, p ), lda,
966 $ ierr )
967 sva( p ) = aapp*dsqrt( max( zero,
968 $ one-aapq*aapq ) )
969 mxsinj = max( mxsinj, sfmin )
970 END IF
971 END IF
972* END IF ROTOK THEN ... ELSE
973*
974* In the case of cancellation in updating SVA(q)
975* .. recompute SVA(q)
976 IF( ( sva( q ) / aaqq )**2.LE.rooteps )
977 $ THEN
978 IF( ( aaqq.LT.rootbig ) .AND.
979 $ ( aaqq.GT.rootsfmin ) ) THEN
980 sva( q ) = dnrm2( m, a( 1, q ),
981 $ 1 )*
982 $ d( q )
983 ELSE
984 t = zero
985 aaqq = one
986 CALL dlassq( m, a( 1, q ), 1, t,
987 $ aaqq )
988 sva( q ) = t*dsqrt( aaqq )*d( q )
989 END IF
990 END IF
991 IF( ( aapp / aapp0 )**2.LE.rooteps ) THEN
992 IF( ( aapp.LT.rootbig ) .AND.
993 $ ( aapp.GT.rootsfmin ) ) THEN
994 aapp = dnrm2( m, a( 1, p ), 1 )*
995 $ d( p )
996 ELSE
997 t = zero
998 aapp = one
999 CALL dlassq( m, a( 1, p ), 1, t,
1000 $ aapp )
1001 aapp = t*dsqrt( aapp )*d( p )
1002 END IF
1003 sva( p ) = aapp
1004 END IF
1005* end of OK rotation
1006 ELSE
1007 notrot = notrot + 1
1008 pskipped = pskipped + 1
1009 ijblsk = ijblsk + 1
1010 END IF
1011 ELSE
1012 notrot = notrot + 1
1013 pskipped = pskipped + 1
1014 ijblsk = ijblsk + 1
1015 END IF
1016*
1017 IF( ( i.LE.swband ) .AND. ( ijblsk.GE.blskip ) )
1018 $ THEN
1019 sva( p ) = aapp
1020 notrot = 0
1021 GO TO 2011
1022 END IF
1023 IF( ( i.LE.swband ) .AND.
1024 $ ( pskipped.GT.rowskip ) ) THEN
1025 aapp = -aapp
1026 notrot = 0
1027 GO TO 2203
1028 END IF
1029*
1030 2200 CONTINUE
1031* end of the q-loop
1032 2203 CONTINUE
1033*
1034 sva( p ) = aapp
1035*
1036 ELSE
1037 IF( aapp.EQ.zero )notrot = notrot +
1038 $ min( jgl+kbl-1, n ) - jgl + 1
1039 IF( aapp.LT.zero )notrot = 0
1040 END IF
1041
1042 2100 CONTINUE
1043* end of the p-loop
1044 2010 CONTINUE
1045* end of the jbc-loop
1046 2011 CONTINUE
1047*2011 bailed out of the jbc-loop
1048 DO 2012 p = igl, min( igl+kbl-1, n )
1049 sva( p ) = dabs( sva( p ) )
1050 2012 CONTINUE
1051*
1052 2000 CONTINUE
1053*2000 :: end of the ibr-loop
1054*
1055* .. update SVA(N)
1056 IF( ( sva( n ).LT.rootbig ) .AND. ( sva( n ).GT.rootsfmin ) )
1057 $ THEN
1058 sva( n ) = dnrm2( m, a( 1, n ), 1 )*d( n )
1059 ELSE
1060 t = zero
1061 aapp = one
1062 CALL dlassq( m, a( 1, n ), 1, t, aapp )
1063 sva( n ) = t*dsqrt( aapp )*d( n )
1064 END IF
1065*
1066* Additional steering devices
1067*
1068 IF( ( i.LT.swband ) .AND. ( ( mxaapq.LE.roottol ) .OR.
1069 $ ( iswrot.LE.n ) ) )swband = i
1070*
1071 IF( ( i.GT.swband+1 ) .AND. ( mxaapq.LT.dble( n )*tol ) .AND.
1072 $ ( dble( n )*mxaapq*mxsinj.LT.tol ) ) THEN
1073 GO TO 1994
1074 END IF
1075*
1076 IF( notrot.GE.emptsw )GO TO 1994
1077
1078 1993 CONTINUE
1079* end i=1:NSWEEP loop
1080* #:) Reaching this point means that the procedure has completed the given
1081* number of iterations.
1082 info = nsweep - 1
1083 GO TO 1995
1084 1994 CONTINUE
1085* #:) Reaching this point means that during the i-th sweep all pivots were
1086* below the given tolerance, causing early exit.
1087*
1088 info = 0
1089* #:) INFO = 0 confirms successful iterations.
1090 1995 CONTINUE
1091*
1092* Sort the vector D.
1093 DO 5991 p = 1, n - 1
1094 q = idamax( n-p+1, sva( p ), 1 ) + p - 1
1095 IF( p.NE.q ) THEN
1096 temp1 = sva( p )
1097 sva( p ) = sva( q )
1098 sva( q ) = temp1
1099 temp1 = d( p )
1100 d( p ) = d( q )
1101 d( q ) = temp1
1102 CALL dswap( m, a( 1, p ), 1, a( 1, q ), 1 )
1103 IF( rsvec )CALL dswap( mvl, v( 1, p ), 1, v( 1, q ), 1 )
1104 END IF
1105 5991 CONTINUE
1106*
1107 RETURN
1108* ..
1109* .. END OF DGSVJ0
1110* ..
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
daxpy
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
ddot
double precision function ddot(n, dx, incx, dy, incy)
DDOT
Definition ddot.f:82
idamax
integer function idamax(n, dx, incx)
IDAMAX
Definition idamax.f:71
dlascl
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.
Definition dlascl.f:142
dlassq
subroutine dlassq(n, x, incx, scale, sumsq)
DLASSQ updates a sum of squares represented in scaled form.
Definition dlassq.f90:122
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
dnrm2
real(wp) function dnrm2(n, x, incx)
DNRM2
Definition dnrm2.f90:89
drotm
subroutine drotm(n, dx, incx, dy, incy, dparam)
DROTM
Definition drotm.f:96
dswap
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82
Here is the call graph for this function:
Here is the caller graph for this function: