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

◆ dgsvj1()

subroutine dgsvj1 ( character*1 jobv,
integer m,
integer n,
integer n1,
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 )

DGSVJ1 pre-processor for the routine dgesvj, applies Jacobi rotations targeting only particular pivots.

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

Purpose:
!>
!> DGSVJ1 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 targets only particular pivots and it does not check convergence
!> (stopping criterion). Few tuning parameters (marked by [TP]) are
!> available for the implementer.
!>
!> Further Details
!> ~~~~~~~~~~~~~~~
!> DGSVJ1 applies few sweeps of Jacobi rotations in the column space of
!> the input M-by-N matrix A. The pivot pairs are taken from the (1,2)
!> off-diagonal block in the corresponding N-by-N Gram matrix A^T * A. The
!> block-entries (tiles) of the (1,2) off-diagonal block are marked by the
!> [x]'s in the following scheme:
!>
!>    | *  *  * [x] [x] [x]|
!>    | *  *  * [x] [x] [x]|    Row-cycling in the nblr-by-nblc [x] blocks.
!>    | *  *  * [x] [x] [x]|    Row-cyclic pivoting inside each [x] block.
!>    |[x] [x] [x] *  *  * |
!>    |[x] [x] [x] *  *  * |
!>    |[x] [x] [x] *  *  * |
!>
!> In terms of the columns of A, the first N1 columns are rotated 'against'
!> the remaining N-N1 columns, trying to increase the angle between the
!> corresponding subspaces. The off-diagonal block is N1-by(N-N1) and it is
!> tiled using quadratic tiles of side KBL. Here, KBL is a tuning parameter.
!> The number of sweeps is given in NSWEEP and the orthogonality threshold
!> is given in TOL.
!>
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]N1
!>          N1 is INTEGER
!>          N1 specifies the 2 x 2 block partition, the first N1 columns are
!>          rotated 'against' the remaining N-N1 columns of A.
!>
[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 N1, 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 N1, 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.
Contributors:
Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)

Definition at line 232 of file dgsvj1.f.

234*
235* -- LAPACK computational routine --
236* -- LAPACK is a software package provided by Univ. of Tennessee, --
237* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
238*
239* .. Scalar Arguments ..
240 DOUBLE PRECISION EPS, SFMIN, TOL
241 INTEGER INFO, LDA, LDV, LWORK, M, MV, N, N1, NSWEEP
242 CHARACTER*1 JOBV
243* ..
244* .. Array Arguments ..
245 DOUBLE PRECISION A( LDA, * ), D( N ), SVA( N ), V( LDV, * ),
246 $ WORK( LWORK )
247* ..
248*
249* =====================================================================
250*
251* .. Local Parameters ..
252 DOUBLE PRECISION ZERO, HALF, ONE
253 parameter( zero = 0.0d0, half = 0.5d0, one = 1.0d0 )
254* ..
255* .. Local Scalars ..
256 DOUBLE PRECISION AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG,
257 $ BIGTHETA, CS, LARGE, MXAAPQ, MXSINJ, ROOTBIG,
258 $ ROOTEPS, ROOTSFMIN, ROOTTOL, SMALL, SN, T,
259 $ TEMP1, THETA, THSIGN
260 INTEGER BLSKIP, EMPTSW, i, ibr, igl, IERR, IJBLSK,
261 $ ISWROT, jbc, jgl, KBL, MVL, NOTROT, nblc, nblr,
262 $ p, PSKIPPED, q, ROWSKIP, SWBAND
263 LOGICAL APPLV, ROTOK, RSVEC
264* ..
265* .. Local Arrays ..
266 DOUBLE PRECISION FASTR( 5 )
267* ..
268* .. Intrinsic Functions ..
269 INTRINSIC dabs, max, dble, min, dsign, dsqrt
270* ..
271* .. External Functions ..
272 DOUBLE PRECISION DDOT, DNRM2
273 INTEGER IDAMAX
274 LOGICAL LSAME
275 EXTERNAL idamax, lsame, ddot, dnrm2
276* ..
277* .. External Subroutines ..
278 EXTERNAL daxpy, dcopy, dlascl, dlassq, drotm,
279 $ dswap,
280 $ xerbla
281* ..
282* .. Executable Statements ..
283*
284* Test the input parameters.
285*
286 applv = lsame( jobv, 'A' )
287 rsvec = lsame( jobv, 'V' )
288 IF( .NOT.( rsvec .OR. applv .OR. lsame( jobv, 'N' ) ) ) THEN
289 info = -1
290 ELSE IF( m.LT.0 ) THEN
291 info = -2
292 ELSE IF( ( n.LT.0 ) .OR. ( n.GT.m ) ) THEN
293 info = -3
294 ELSE IF( n1.LT.0 ) THEN
295 info = -4
296 ELSE IF( lda.LT.m ) THEN
297 info = -6
298 ELSE IF( ( rsvec.OR.applv ) .AND. ( mv.LT.0 ) ) THEN
299 info = -9
300 ELSE IF( ( rsvec.AND.( ldv.LT.n ) ).OR.
301 $ ( applv.AND.( ldv.LT.mv ) ) ) THEN
302 info = -11
303 ELSE IF( tol.LE.eps ) THEN
304 info = -14
305 ELSE IF( nsweep.LT.0 ) THEN
306 info = -15
307 ELSE IF( lwork.LT.m ) THEN
308 info = -17
309 ELSE
310 info = 0
311 END IF
312*
313* #:(
314 IF( info.NE.0 ) THEN
315 CALL xerbla( 'DGSVJ1', -info )
316 RETURN
317 END IF
318*
319 IF( rsvec ) THEN
320 mvl = n
321 ELSE IF( applv ) THEN
322 mvl = mv
323 END IF
324 rsvec = rsvec .OR. applv
325
326 rooteps = dsqrt( eps )
327 rootsfmin = dsqrt( sfmin )
328 small = sfmin / eps
329 big = one / sfmin
330 rootbig = one / rootsfmin
331 large = big / dsqrt( dble( m*n ) )
332 bigtheta = one / rooteps
333 roottol = dsqrt( tol )
334*
335* .. Initialize the right singular vector matrix ..
336*
337* RSVEC = LSAME( JOBV, 'Y' )
338*
339 emptsw = n1*( n-n1 )
340 notrot = 0
341 fastr( 1 ) = zero
342*
343* .. Row-cyclic pivot strategy with de Rijk's pivoting ..
344*
345 kbl = min( 8, n )
346 nblr = n1 / kbl
347 IF( ( nblr*kbl ).NE.n1 )nblr = nblr + 1
348
349* .. the tiling is nblr-by-nblc [tiles]
350
351 nblc = ( n-n1 ) / kbl
352 IF( ( nblc*kbl ).NE.( n-n1 ) )nblc = nblc + 1
353 blskip = ( kbl**2 ) + 1
354*[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
355
356 rowskip = min( 5, kbl )
357*[TP] ROWSKIP is a tuning parameter.
358 swband = 0
359*[TP] SWBAND is a tuning parameter. It is meaningful and effective
360* if SGESVJ is used as a computational routine in the preconditioned
361* Jacobi SVD algorithm SGESVJ.
362*
363*
364* | * * * [x] [x] [x]|
365* | * * * [x] [x] [x]| Row-cycling in the nblr-by-nblc [x] blocks.
366* | * * * [x] [x] [x]| Row-cyclic pivoting inside each [x] block.
367* |[x] [x] [x] * * * |
368* |[x] [x] [x] * * * |
369* |[x] [x] [x] * * * |
370*
371*
372 DO 1993 i = 1, nsweep
373* .. go go go ...
374*
375 mxaapq = zero
376 mxsinj = zero
377 iswrot = 0
378*
379 notrot = 0
380 pskipped = 0
381*
382 DO 2000 ibr = 1, nblr
383
384 igl = ( ibr-1 )*kbl + 1
385*
386*
387*........................................................
388* ... go to the off diagonal blocks
389
390 igl = ( ibr-1 )*kbl + 1
391
392 DO 2010 jbc = 1, nblc
393
394 jgl = n1 + ( jbc-1 )*kbl + 1
395
396* doing the block at ( ibr, jbc )
397
398 ijblsk = 0
399 DO 2100 p = igl, min( igl+kbl-1, n1 )
400
401 aapp = sva( p )
402
403 IF( aapp.GT.zero ) THEN
404
405 pskipped = 0
406
407 DO 2200 q = jgl, min( jgl+kbl-1, n )
408*
409 aaqq = sva( q )
410
411 IF( aaqq.GT.zero ) THEN
412 aapp0 = aapp
413*
414* .. M x 2 Jacobi SVD ..
415*
416* .. Safe Gram matrix computation ..
417*
418 IF( aaqq.GE.one ) THEN
419 IF( aapp.GE.aaqq ) THEN
420 rotok = ( small*aapp ).LE.aaqq
421 ELSE
422 rotok = ( small*aaqq ).LE.aapp
423 END IF
424 IF( aapp.LT.( big / aaqq ) ) THEN
425 aapq = ( ddot( m, a( 1, p ), 1,
426 $ a( 1,
427 $ q ), 1 )*d( p )*d( q ) / aaqq )
428 $ / aapp
429 ELSE
430 CALL dcopy( m, a( 1, p ), 1, work,
431 $ 1 )
432 CALL dlascl( 'G', 0, 0, aapp,
433 $ d( p ),
434 $ m, 1, work, lda, ierr )
435 aapq = ddot( m, work, 1, a( 1, q ),
436 $ 1 )*d( q ) / aaqq
437 END IF
438 ELSE
439 IF( aapp.GE.aaqq ) THEN
440 rotok = aapp.LE.( aaqq / small )
441 ELSE
442 rotok = aaqq.LE.( aapp / small )
443 END IF
444 IF( aapp.GT.( small / aaqq ) ) THEN
445 aapq = ( ddot( m, a( 1, p ), 1,
446 $ a( 1,
447 $ q ), 1 )*d( p )*d( q ) / aaqq )
448 $ / aapp
449 ELSE
450 CALL dcopy( m, a( 1, q ), 1, work,
451 $ 1 )
452 CALL dlascl( 'G', 0, 0, aaqq,
453 $ d( q ),
454 $ m, 1, work, lda, ierr )
455 aapq = ddot( m, work, 1, a( 1, p ),
456 $ 1 )*d( p ) / aapp
457 END IF
458 END IF
459
460 mxaapq = max( mxaapq, dabs( aapq ) )
461
462* TO rotate or NOT to rotate, THAT is the question ...
463*
464 IF( dabs( aapq ).GT.tol ) THEN
465 notrot = 0
466* ROTATED = ROTATED + 1
467 pskipped = 0
468 iswrot = iswrot + 1
469*
470 IF( rotok ) THEN
471*
472 aqoap = aaqq / aapp
473 apoaq = aapp / aaqq
474 theta = -half*dabs(aqoap-apoaq) / aapq
475 IF( aaqq.GT.aapp0 )theta = -theta
476
477 IF( dabs( theta ).GT.bigtheta ) THEN
478 t = half / theta
479 fastr( 3 ) = t*d( p ) / d( q )
480 fastr( 4 ) = -t*d( q ) / d( p )
481 CALL drotm( m, a( 1, p ), 1,
482 $ a( 1, q ), 1, fastr )
483 IF( rsvec )CALL drotm( mvl,
484 $ v( 1, p ), 1,
485 $ v( 1, q ), 1,
486 $ fastr )
487 sva( q ) = aaqq*dsqrt( max( zero,
488 $ one+t*apoaq*aapq ) )
489 aapp = aapp*dsqrt( max( zero,
490 $ one-t*aqoap*aapq ) )
491 mxsinj = max( mxsinj, dabs( t ) )
492 ELSE
493*
494* .. choose correct signum for THETA and rotate
495*
496 thsign = -dsign( one, aapq )
497 IF( aaqq.GT.aapp0 )thsign = -thsign
498 t = one / ( theta+thsign*
499 $ dsqrt( one+theta*theta ) )
500 cs = dsqrt( one / ( one+t*t ) )
501 sn = t*cs
502 mxsinj = max( mxsinj, dabs( sn ) )
503 sva( q ) = aaqq*dsqrt( max( zero,
504 $ one+t*apoaq*aapq ) )
505 aapp = aapp*dsqrt( max( zero,
506 $ one-t*aqoap*aapq ) )
507
508 apoaq = d( p ) / d( q )
509 aqoap = d( q ) / d( p )
510 IF( d( p ).GE.one ) THEN
511*
512 IF( d( q ).GE.one ) THEN
513 fastr( 3 ) = t*apoaq
514 fastr( 4 ) = -t*aqoap
515 d( p ) = d( p )*cs
516 d( q ) = d( q )*cs
517 CALL drotm( m, a( 1, p ),
518 $ 1,
519 $ a( 1, q ), 1,
520 $ fastr )
521 IF( rsvec )CALL drotm( mvl,
522 $ v( 1, p ), 1, v( 1, q ),
523 $ 1, fastr )
524 ELSE
525 CALL daxpy( m, -t*aqoap,
526 $ a( 1, q ), 1,
527 $ a( 1, p ), 1 )
528 CALL daxpy( m, cs*sn*apoaq,
529 $ a( 1, p ), 1,
530 $ a( 1, q ), 1 )
531 IF( rsvec ) THEN
532 CALL daxpy( mvl,
533 $ -t*aqoap,
534 $ v( 1, q ), 1,
535 $ v( 1, p ), 1 )
536 CALL daxpy( mvl,
537 $ cs*sn*apoaq,
538 $ v( 1, p ), 1,
539 $ v( 1, q ), 1 )
540 END IF
541 d( p ) = d( p )*cs
542 d( q ) = d( q ) / cs
543 END IF
544 ELSE
545 IF( d( q ).GE.one ) THEN
546 CALL daxpy( m, t*apoaq,
547 $ a( 1, p ), 1,
548 $ a( 1, q ), 1 )
549 CALL daxpy( m,
550 $ -cs*sn*aqoap,
551 $ a( 1, q ), 1,
552 $ a( 1, p ), 1 )
553 IF( rsvec ) THEN
554 CALL daxpy( mvl,
555 $ t*apoaq,
556 $ v( 1, p ), 1,
557 $ v( 1, q ), 1 )
558 CALL daxpy( mvl,
559 $ -cs*sn*aqoap,
560 $ v( 1, q ), 1,
561 $ v( 1, p ), 1 )
562 END IF
563 d( p ) = d( p ) / cs
564 d( q ) = d( q )*cs
565 ELSE
566 IF( d( p ).GE.d( q ) ) THEN
567 CALL daxpy( m, -t*aqoap,
568 $ a( 1, q ), 1,
569 $ a( 1, p ), 1 )
570 CALL daxpy( m,
571 $ cs*sn*apoaq,
572 $ a( 1, p ), 1,
573 $ a( 1, q ), 1 )
574 d( p ) = d( p )*cs
575 d( q ) = d( q ) / cs
576 IF( rsvec ) THEN
577 CALL daxpy( mvl,
578 $ -t*aqoap,
579 $ v( 1, q ), 1,
580 $ v( 1, p ), 1 )
581 CALL daxpy( mvl,
582 $ cs*sn*apoaq,
583 $ v( 1, p ), 1,
584 $ v( 1, q ), 1 )
585 END IF
586 ELSE
587 CALL daxpy( m, t*apoaq,
588 $ a( 1, p ), 1,
589 $ a( 1, q ), 1 )
590 CALL daxpy( m,
591 $ -cs*sn*aqoap,
592 $ a( 1, q ), 1,
593 $ a( 1, p ), 1 )
594 d( p ) = d( p ) / cs
595 d( q ) = d( q )*cs
596 IF( rsvec ) THEN
597 CALL daxpy( mvl,
598 $ t*apoaq, v( 1, p ),
599 $ 1, v( 1, q ), 1 )
600 CALL daxpy( mvl,
601 $ -cs*sn*aqoap,
602 $ v( 1, q ), 1,
603 $ v( 1, p ), 1 )
604 END IF
605 END IF
606 END IF
607 END IF
608 END IF
609
610 ELSE
611 IF( aapp.GT.aaqq ) THEN
612 CALL dcopy( m, a( 1, p ), 1,
613 $ work,
614 $ 1 )
615 CALL dlascl( 'G', 0, 0, aapp,
616 $ one,
617 $ m, 1, work, lda, ierr )
618 CALL dlascl( 'G', 0, 0, aaqq,
619 $ one,
620 $ m, 1, a( 1, q ), lda,
621 $ 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,
626 $ aaqq,
627 $ m, 1, a( 1, q ), lda,
628 $ ierr )
629 sva( q ) = aaqq*dsqrt( max( zero,
630 $ one-aapq*aapq ) )
631 mxsinj = max( mxsinj, sfmin )
632 ELSE
633 CALL dcopy( m, a( 1, q ), 1,
634 $ work,
635 $ 1 )
636 CALL dlascl( 'G', 0, 0, aaqq,
637 $ one,
638 $ m, 1, work, lda, ierr )
639 CALL dlascl( 'G', 0, 0, aapp,
640 $ one,
641 $ m, 1, a( 1, p ), lda,
642 $ ierr )
643 temp1 = -aapq*d( q ) / d( p )
644 CALL daxpy( m, temp1, work, 1,
645 $ a( 1, p ), 1 )
646 CALL dlascl( 'G', 0, 0, one,
647 $ aapp,
648 $ m, 1, a( 1, p ), lda,
649 $ ierr )
650 sva( p ) = aapp*dsqrt( max( zero,
651 $ one-aapq*aapq ) )
652 mxsinj = max( mxsinj, sfmin )
653 END IF
654 END IF
655* END IF ROTOK THEN ... ELSE
656*
657* In the case of cancellation in updating SVA(q)
658* .. recompute SVA(q)
659 IF( ( sva( q ) / aaqq )**2.LE.rooteps )
660 $ THEN
661 IF( ( aaqq.LT.rootbig ) .AND.
662 $ ( aaqq.GT.rootsfmin ) ) THEN
663 sva( q ) = dnrm2( m, a( 1, q ),
664 $ 1 )*
665 $ d( q )
666 ELSE
667 t = zero
668 aaqq = one
669 CALL dlassq( m, a( 1, q ), 1, t,
670 $ aaqq )
671 sva( q ) = t*dsqrt( aaqq )*d( q )
672 END IF
673 END IF
674 IF( ( aapp / aapp0 )**2.LE.rooteps ) THEN
675 IF( ( aapp.LT.rootbig ) .AND.
676 $ ( aapp.GT.rootsfmin ) ) THEN
677 aapp = dnrm2( m, a( 1, p ), 1 )*
678 $ d( p )
679 ELSE
680 t = zero
681 aapp = one
682 CALL dlassq( m, a( 1, p ), 1, t,
683 $ aapp )
684 aapp = t*dsqrt( aapp )*d( p )
685 END IF
686 sva( p ) = aapp
687 END IF
688* end of OK rotation
689 ELSE
690 notrot = notrot + 1
691* SKIPPED = SKIPPED + 1
692 pskipped = pskipped + 1
693 ijblsk = ijblsk + 1
694 END IF
695 ELSE
696 notrot = notrot + 1
697 pskipped = pskipped + 1
698 ijblsk = ijblsk + 1
699 END IF
700
701* IF ( NOTROT .GE. EMPTSW ) GO TO 2011
702 IF( ( i.LE.swband ) .AND. ( ijblsk.GE.blskip ) )
703 $ THEN
704 sva( p ) = aapp
705 notrot = 0
706 GO TO 2011
707 END IF
708 IF( ( i.LE.swband ) .AND.
709 $ ( pskipped.GT.rowskip ) ) THEN
710 aapp = -aapp
711 notrot = 0
712 GO TO 2203
713 END IF
714
715*
716 2200 CONTINUE
717* end of the q-loop
718 2203 CONTINUE
719
720 sva( p ) = aapp
721*
722 ELSE
723 IF( aapp.EQ.zero )notrot = notrot +
724 $ min( jgl+kbl-1, n ) - jgl + 1
725 IF( aapp.LT.zero )notrot = 0
726*** IF ( NOTROT .GE. EMPTSW ) GO TO 2011
727 END IF
728
729 2100 CONTINUE
730* end of the p-loop
731 2010 CONTINUE
732* end of the jbc-loop
733 2011 CONTINUE
734*2011 bailed out of the jbc-loop
735 DO 2012 p = igl, min( igl+kbl-1, n )
736 sva( p ) = dabs( sva( p ) )
737 2012 CONTINUE
738*** IF ( NOTROT .GE. EMPTSW ) GO TO 1994
739 2000 CONTINUE
740*2000 :: end of the ibr-loop
741*
742* .. update SVA(N)
743 IF( ( sva( n ).LT.rootbig ) .AND. ( sva( n ).GT.rootsfmin ) )
744 $ THEN
745 sva( n ) = dnrm2( m, a( 1, n ), 1 )*d( n )
746 ELSE
747 t = zero
748 aapp = one
749 CALL dlassq( m, a( 1, n ), 1, t, aapp )
750 sva( n ) = t*dsqrt( aapp )*d( n )
751 END IF
752*
753* Additional steering devices
754*
755 IF( ( i.LT.swband ) .AND. ( ( mxaapq.LE.roottol ) .OR.
756 $ ( iswrot.LE.n ) ) )swband = i
757
758 IF( ( i.GT.swband+1 ) .AND. ( mxaapq.LT.dble( n )*tol ) .AND.
759 $ ( dble( n )*mxaapq*mxsinj.LT.tol ) ) THEN
760 GO TO 1994
761 END IF
762
763*
764 IF( notrot.GE.emptsw )GO TO 1994
765
766 1993 CONTINUE
767* end i=1:NSWEEP loop
768* #:) Reaching this point means that the procedure has completed the given
769* number of sweeps.
770 info = nsweep - 1
771 GO TO 1995
772 1994 CONTINUE
773* #:) Reaching this point means that during the i-th sweep all pivots were
774* below the given threshold, causing early exit.
775
776 info = 0
777* #:) INFO = 0 confirms successful iterations.
778 1995 CONTINUE
779*
780* Sort the vector D
781*
782 DO 5991 p = 1, n - 1
783 q = idamax( n-p+1, sva( p ), 1 ) + p - 1
784 IF( p.NE.q ) THEN
785 temp1 = sva( p )
786 sva( p ) = sva( q )
787 sva( q ) = temp1
788 temp1 = d( p )
789 d( p ) = d( q )
790 d( q ) = temp1
791 CALL dswap( m, a( 1, p ), 1, a( 1, q ), 1 )
792 IF( rsvec )CALL dswap( mvl, v( 1, p ), 1, v( 1, q ), 1 )
793 END IF
794 5991 CONTINUE
795*
796 RETURN
797* ..
798* .. END OF DGSVJ1
799* ..
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: