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

◆ cgsvj0()

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

CGSVJ0 pre-processor for the routine cgesvj.

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

Purpose:
 CGSVJ0 is called from CGESVJ as a pre-processor and that is its main
 purpose. It applies Jacobi rotations in the same way as CGESVJ 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 COMPLEX array, dimension (LDA,N)
          On entry, M-by-N matrix A, such that A*diag(D) represents
          the input matrix.
          On exit,
          A_onexit * diag(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 COMPLEX array, dimension (N)
          The array D accumulates the scaling factors from the complex 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 REAL 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 A_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 COMPLEX 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 REAL
          EPS = SLAMCH('Epsilon')
[in]SFMIN
          SFMIN is REAL
          SFMIN = SLAMCH('Safe Minimum')
[in]TOL
          TOL is REAL
          TOL is the threshold for Jacobi rotations. For a pair
          A(:,p), A(:,q) of pivot columns, the Jacobi rotation is
          applied only if ABS(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 COMPLEX 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:
CGSVJ0 is used just to enable CGESVJ to call a simplified version of itself to work on a submatrix of the original matrix.
Contributor:
Zlatko Drmac (Zagreb, Croatia)
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 216 of file cgsvj0.f.

218*
219* -- LAPACK computational routine --
220* -- LAPACK is a software package provided by Univ. of Tennessee, --
221* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
222*
223 IMPLICIT NONE
224* .. Scalar Arguments ..
225 INTEGER INFO, LDA, LDV, LWORK, M, MV, N, NSWEEP
226 REAL EPS, SFMIN, TOL
227 CHARACTER*1 JOBV
228* ..
229* .. Array Arguments ..
230 COMPLEX A( LDA, * ), D( N ), V( LDV, * ), WORK( LWORK )
231 REAL SVA( N )
232* ..
233*
234* =====================================================================
235*
236* .. Local Parameters ..
237 REAL ZERO, HALF, ONE
238 parameter( zero = 0.0e0, half = 0.5e0, one = 1.0e0)
239 COMPLEX CZERO, CONE
240 parameter( czero = (0.0e0, 0.0e0), cone = (1.0e0, 0.0e0) )
241* ..
242* .. Local Scalars ..
243 COMPLEX AAPQ, OMPQ
244 REAL AAPP, AAPP0, AAPQ1, AAQQ, APOAQ, AQOAP, BIG,
245 $ BIGTHETA, CS, MXAAPQ, MXSINJ, ROOTBIG, ROOTEPS,
246 $ ROOTSFMIN, ROOTTOL, SMALL, SN, T, TEMP1, THETA,
247 $ THSIGN
248 INTEGER BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1,
249 $ ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, NBL,
250 $ NOTROT, p, PSKIPPED, q, ROWSKIP, SWBAND
251 LOGICAL APPLV, ROTOK, RSVEC
252* ..
253* ..
254* .. Intrinsic Functions ..
255 INTRINSIC abs, max, conjg, real, min, sign, sqrt
256* ..
257* .. External Functions ..
258 REAL SCNRM2
259 COMPLEX CDOTC
260 INTEGER ISAMAX
261 LOGICAL LSAME
262 EXTERNAL isamax, lsame, cdotc, scnrm2
263* ..
264* ..
265* .. External Subroutines ..
266* ..
267* from BLAS
268 EXTERNAL ccopy, crot, cswap, caxpy
269* from LAPACK
270 EXTERNAL clascl, classq, xerbla
271* ..
272* .. Executable Statements ..
273*
274* Test the input parameters.
275*
276 applv = lsame( jobv, 'A' )
277 rsvec = lsame( jobv, 'V' )
278 IF( .NOT.( rsvec .OR. applv .OR. lsame( jobv, 'N' ) ) ) THEN
279 info = -1
280 ELSE IF( m.LT.0 ) THEN
281 info = -2
282 ELSE IF( ( n.LT.0 ) .OR. ( n.GT.m ) ) THEN
283 info = -3
284 ELSE IF( lda.LT.m ) THEN
285 info = -5
286 ELSE IF( ( rsvec.OR.applv ) .AND. ( mv.LT.0 ) ) THEN
287 info = -8
288 ELSE IF( ( rsvec.AND.( ldv.LT.n ) ).OR.
289 $ ( applv.AND.( ldv.LT.mv ) ) ) THEN
290 info = -10
291 ELSE IF( tol.LE.eps ) THEN
292 info = -13
293 ELSE IF( nsweep.LT.0 ) THEN
294 info = -14
295 ELSE IF( lwork.LT.m ) THEN
296 info = -16
297 ELSE
298 info = 0
299 END IF
300*
301* #:(
302 IF( info.NE.0 ) THEN
303 CALL xerbla( 'CGSVJ0', -info )
304 RETURN
305 END IF
306*
307 IF( rsvec ) THEN
308 mvl = n
309 ELSE IF( applv ) THEN
310 mvl = mv
311 END IF
312 rsvec = rsvec .OR. applv
313
314 rooteps = sqrt( eps )
315 rootsfmin = sqrt( sfmin )
316 small = sfmin / eps
317 big = one / sfmin
318 rootbig = one / rootsfmin
319 bigtheta = one / rooteps
320 roottol = sqrt( tol )
321*
322* .. Row-cyclic Jacobi SVD algorithm with column pivoting ..
323*
324 emptsw = ( n*( n-1 ) ) / 2
325 notrot = 0
326*
327* .. Row-cyclic pivot strategy with de Rijk's pivoting ..
328*
329
330 swband = 0
331*[TP] SWBAND is a tuning parameter [TP]. It is meaningful and effective
332* if CGESVJ is used as a computational routine in the preconditioned
333* Jacobi SVD algorithm CGEJSV. For sweeps i=1:SWBAND the procedure
334* works on pivots inside a band-like region around the diagonal.
335* The boundaries are determined dynamically, based on the number of
336* pivots above a threshold.
337*
338 kbl = min( 8, n )
339*[TP] KBL is a tuning parameter that defines the tile size in the
340* tiling of the p-q loops of pivot pairs. In general, an optimal
341* value of KBL depends on the matrix dimensions and on the
342* parameters of the computer's memory.
343*
344 nbl = n / kbl
345 IF( ( nbl*kbl ).NE.n )nbl = nbl + 1
346*
347 blskip = kbl**2
348*[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
349*
350 rowskip = min( 5, kbl )
351*[TP] ROWSKIP is a tuning parameter.
352*
353 lkahead = 1
354*[TP] LKAHEAD is a tuning parameter.
355*
356* Quasi block transformations, using the lower (upper) triangular
357* structure of the input matrix. The quasi-block-cycling usually
358* invokes cubic convergence. Big part of this cycle is done inside
359* canonical subspaces of dimensions less than M.
360*
361*
362* .. Row-cyclic pivot strategy with de Rijk's pivoting ..
363*
364 DO 1993 i = 1, nsweep
365*
366* .. go go go ...
367*
368 mxaapq = zero
369 mxsinj = zero
370 iswrot = 0
371*
372 notrot = 0
373 pskipped = 0
374*
375* Each sweep is unrolled using KBL-by-KBL tiles over the pivot pairs
376* 1 <= p < q <= N. This is the first step toward a blocked implementation
377* of the rotations. New implementation, based on block transformations,
378* is under development.
379*
380 DO 2000 ibr = 1, nbl
381*
382 igl = ( ibr-1 )*kbl + 1
383*
384 DO 1002 ir1 = 0, min( lkahead, nbl-ibr )
385*
386 igl = igl + ir1*kbl
387*
388 DO 2001 p = igl, min( igl+kbl-1, n-1 )
389*
390* .. de Rijk's pivoting
391*
392 q = isamax( n-p+1, sva( p ), 1 ) + p - 1
393 IF( p.NE.q ) THEN
394 CALL cswap( m, a( 1, p ), 1, a( 1, q ), 1 )
395 IF( rsvec )CALL cswap( mvl, v( 1, p ), 1,
396 $ v( 1, q ), 1 )
397 temp1 = sva( p )
398 sva( p ) = sva( q )
399 sva( q ) = temp1
400 aapq = d(p)
401 d(p) = d(q)
402 d(q) = aapq
403 END IF
404*
405 IF( ir1.EQ.0 ) THEN
406*
407* Column norms are periodically updated by explicit
408* norm computation.
409* Caveat:
410* Unfortunately, some BLAS implementations compute SNCRM2(M,A(1,p),1)
411* as SQRT(S=CDOTC(M,A(1,p),1,A(1,p),1)), which may cause the result to
412* overflow for ||A(:,p)||_2 > SQRT(overflow_threshold), and to
413* underflow for ||A(:,p)||_2 < SQRT(underflow_threshold).
414* Hence, SCNRM2 cannot be trusted, not even in the case when
415* the true norm is far from the under(over)flow boundaries.
416* If properly implemented SCNRM2 is available, the IF-THEN-ELSE-END IF
417* below should be replaced with "AAPP = SCNRM2( M, A(1,p), 1 )".
418*
419 IF( ( sva( p ).LT.rootbig ) .AND.
420 $ ( sva( p ).GT.rootsfmin ) ) THEN
421 sva( p ) = scnrm2( m, a( 1, p ), 1 )
422 ELSE
423 temp1 = zero
424 aapp = one
425 CALL classq( m, a( 1, p ), 1, temp1, aapp )
426 sva( p ) = temp1*sqrt( aapp )
427 END IF
428 aapp = sva( p )
429 ELSE
430 aapp = sva( p )
431 END IF
432*
433 IF( aapp.GT.zero ) THEN
434*
435 pskipped = 0
436*
437 DO 2002 q = p + 1, min( igl+kbl-1, n )
438*
439 aaqq = sva( q )
440*
441 IF( aaqq.GT.zero ) THEN
442*
443 aapp0 = aapp
444 IF( aaqq.GE.one ) THEN
445 rotok = ( small*aapp ).LE.aaqq
446 IF( aapp.LT.( big / aaqq ) ) THEN
447 aapq = ( cdotc( m, a( 1, p ), 1,
448 $ a( 1, q ), 1 ) / aaqq ) / aapp
449 ELSE
450 CALL ccopy( m, a( 1, p ), 1,
451 $ work, 1 )
452 CALL clascl( 'G', 0, 0, aapp, one,
453 $ m, 1, work, lda, ierr )
454 aapq = cdotc( m, work, 1,
455 $ a( 1, q ), 1 ) / aaqq
456 END IF
457 ELSE
458 rotok = aapp.LE.( aaqq / small )
459 IF( aapp.GT.( small / aaqq ) ) THEN
460 aapq = ( cdotc( m, a( 1, p ), 1,
461 $ a( 1, q ), 1 ) / aapp ) / aaqq
462 ELSE
463 CALL ccopy( m, a( 1, q ), 1,
464 $ work, 1 )
465 CALL clascl( 'G', 0, 0, aaqq,
466 $ one, m, 1,
467 $ work, lda, ierr )
468 aapq = cdotc( m, a( 1, p ), 1,
469 $ work, 1 ) / aapp
470 END IF
471 END IF
472*
473* AAPQ = AAPQ * CONJG( CWORK(p) ) * CWORK(q)
474 aapq1 = -abs(aapq)
475 mxaapq = max( mxaapq, -aapq1 )
476*
477* TO rotate or NOT to rotate, THAT is the question ...
478*
479 IF( abs( aapq1 ).GT.tol ) THEN
480 ompq = aapq / abs(aapq)
481*
482* .. rotate
483*[RTD] ROTATED = ROTATED + ONE
484*
485 IF( ir1.EQ.0 ) THEN
486 notrot = 0
487 pskipped = 0
488 iswrot = iswrot + 1
489 END IF
490*
491 IF( rotok ) THEN
492*
493 aqoap = aaqq / aapp
494 apoaq = aapp / aaqq
495 theta = -half*abs( aqoap-apoaq )/aapq1
496*
497 IF( abs( theta ).GT.bigtheta ) THEN
498*
499 t = half / theta
500 cs = one
501
502 CALL crot( m, a(1,p), 1, a(1,q), 1,
503 $ cs, conjg(ompq)*t )
504 IF ( rsvec ) THEN
505 CALL crot( mvl, v(1,p), 1,
506 $ v(1,q), 1, cs, conjg(ompq)*t )
507 END IF
508
509 sva( q ) = aaqq*sqrt( max( zero,
510 $ one+t*apoaq*aapq1 ) )
511 aapp = aapp*sqrt( max( zero,
512 $ one-t*aqoap*aapq1 ) )
513 mxsinj = max( mxsinj, abs( t ) )
514*
515 ELSE
516*
517* .. choose correct signum for THETA and rotate
518*
519 thsign = -sign( one, aapq1 )
520 t = one / ( theta+thsign*
521 $ sqrt( one+theta*theta ) )
522 cs = sqrt( one / ( one+t*t ) )
523 sn = t*cs
524*
525 mxsinj = max( mxsinj, abs( sn ) )
526 sva( q ) = aaqq*sqrt( max( zero,
527 $ one+t*apoaq*aapq1 ) )
528 aapp = aapp*sqrt( max( zero,
529 $ one-t*aqoap*aapq1 ) )
530*
531 CALL crot( m, a(1,p), 1, a(1,q), 1,
532 $ cs, conjg(ompq)*sn )
533 IF ( rsvec ) THEN
534 CALL crot( mvl, v(1,p), 1,
535 $ v(1,q), 1, cs, conjg(ompq)*sn )
536 END IF
537 END IF
538 d(p) = -d(q) * ompq
539*
540 ELSE
541* .. have to use modified Gram-Schmidt like transformation
542 CALL ccopy( m, a( 1, p ), 1,
543 $ work, 1 )
544 CALL clascl( 'G', 0, 0, aapp, one, m,
545 $ 1, work, lda,
546 $ ierr )
547 CALL clascl( 'G', 0, 0, aaqq, one, m,
548 $ 1, a( 1, q ), lda, ierr )
549 CALL caxpy( m, -aapq, work, 1,
550 $ a( 1, q ), 1 )
551 CALL clascl( 'G', 0, 0, one, aaqq, m,
552 $ 1, a( 1, q ), lda, ierr )
553 sva( q ) = aaqq*sqrt( max( zero,
554 $ one-aapq1*aapq1 ) )
555 mxsinj = max( mxsinj, sfmin )
556 END IF
557* END IF ROTOK THEN ... ELSE
558*
559* In the case of cancellation in updating SVA(q), SVA(p)
560* recompute SVA(q), SVA(p).
561*
562 IF( ( sva( q ) / aaqq )**2.LE.rooteps )
563 $ THEN
564 IF( ( aaqq.LT.rootbig ) .AND.
565 $ ( aaqq.GT.rootsfmin ) ) THEN
566 sva( q ) = scnrm2( m, a( 1, q ), 1 )
567 ELSE
568 t = zero
569 aaqq = one
570 CALL classq( m, a( 1, q ), 1, t,
571 $ aaqq )
572 sva( q ) = t*sqrt( aaqq )
573 END IF
574 END IF
575 IF( ( aapp / aapp0 ).LE.rooteps ) THEN
576 IF( ( aapp.LT.rootbig ) .AND.
577 $ ( aapp.GT.rootsfmin ) ) THEN
578 aapp = scnrm2( m, a( 1, p ), 1 )
579 ELSE
580 t = zero
581 aapp = one
582 CALL classq( m, a( 1, p ), 1, t,
583 $ aapp )
584 aapp = t*sqrt( aapp )
585 END IF
586 sva( p ) = aapp
587 END IF
588*
589 ELSE
590* A(:,p) and A(:,q) already numerically orthogonal
591 IF( ir1.EQ.0 )notrot = notrot + 1
592*[RTD] SKIPPED = SKIPPED + 1
593 pskipped = pskipped + 1
594 END IF
595 ELSE
596* A(:,q) is zero column
597 IF( ir1.EQ.0 )notrot = notrot + 1
598 pskipped = pskipped + 1
599 END IF
600*
601 IF( ( i.LE.swband ) .AND.
602 $ ( pskipped.GT.rowskip ) ) THEN
603 IF( ir1.EQ.0 )aapp = -aapp
604 notrot = 0
605 GO TO 2103
606 END IF
607*
608 2002 CONTINUE
609* END q-LOOP
610*
611 2103 CONTINUE
612* bailed out of q-loop
613*
614 sva( p ) = aapp
615*
616 ELSE
617 sva( p ) = aapp
618 IF( ( ir1.EQ.0 ) .AND. ( aapp.EQ.zero ) )
619 $ notrot = notrot + min( igl+kbl-1, n ) - p
620 END IF
621*
622 2001 CONTINUE
623* end of the p-loop
624* end of doing the block ( ibr, ibr )
625 1002 CONTINUE
626* end of ir1-loop
627*
628* ... go to the off diagonal blocks
629*
630 igl = ( ibr-1 )*kbl + 1
631*
632 DO 2010 jbc = ibr + 1, nbl
633*
634 jgl = ( jbc-1 )*kbl + 1
635*
636* doing the block at ( ibr, jbc )
637*
638 ijblsk = 0
639 DO 2100 p = igl, min( igl+kbl-1, n )
640*
641 aapp = sva( p )
642 IF( aapp.GT.zero ) THEN
643*
644 pskipped = 0
645*
646 DO 2200 q = jgl, min( jgl+kbl-1, n )
647*
648 aaqq = sva( q )
649 IF( aaqq.GT.zero ) THEN
650 aapp0 = aapp
651*
652* .. M x 2 Jacobi SVD ..
653*
654* Safe Gram matrix computation
655*
656 IF( aaqq.GE.one ) THEN
657 IF( aapp.GE.aaqq ) THEN
658 rotok = ( small*aapp ).LE.aaqq
659 ELSE
660 rotok = ( small*aaqq ).LE.aapp
661 END IF
662 IF( aapp.LT.( big / aaqq ) ) THEN
663 aapq = ( cdotc( m, a( 1, p ), 1,
664 $ a( 1, q ), 1 ) / aaqq ) / aapp
665 ELSE
666 CALL ccopy( m, a( 1, p ), 1,
667 $ work, 1 )
668 CALL clascl( 'G', 0, 0, aapp,
669 $ one, m, 1,
670 $ work, lda, ierr )
671 aapq = cdotc( m, work, 1,
672 $ a( 1, q ), 1 ) / aaqq
673 END IF
674 ELSE
675 IF( aapp.GE.aaqq ) THEN
676 rotok = aapp.LE.( aaqq / small )
677 ELSE
678 rotok = aaqq.LE.( aapp / small )
679 END IF
680 IF( aapp.GT.( small / aaqq ) ) THEN
681 aapq = ( cdotc( m, a( 1, p ), 1,
682 $ a( 1, q ), 1 ) / max(aaqq,aapp) )
683 $ / min(aaqq,aapp)
684 ELSE
685 CALL ccopy( m, a( 1, q ), 1,
686 $ work, 1 )
687 CALL clascl( 'G', 0, 0, aaqq,
688 $ one, m, 1,
689 $ work, lda, ierr )
690 aapq = cdotc( m, a( 1, p ), 1,
691 $ work, 1 ) / aapp
692 END IF
693 END IF
694*
695* AAPQ = AAPQ * CONJG(CWORK(p))*CWORK(q)
696 aapq1 = -abs(aapq)
697 mxaapq = max( mxaapq, -aapq1 )
698*
699* TO rotate or NOT to rotate, THAT is the question ...
700*
701 IF( abs( aapq1 ).GT.tol ) THEN
702 ompq = aapq / abs(aapq)
703 notrot = 0
704*[RTD] ROTATED = ROTATED + 1
705 pskipped = 0
706 iswrot = iswrot + 1
707*
708 IF( rotok ) THEN
709*
710 aqoap = aaqq / aapp
711 apoaq = aapp / aaqq
712 theta = -half*abs( aqoap-apoaq )/ aapq1
713 IF( aaqq.GT.aapp0 )theta = -theta
714*
715 IF( abs( theta ).GT.bigtheta ) THEN
716 t = half / theta
717 cs = one
718 CALL crot( m, a(1,p), 1, a(1,q), 1,
719 $ cs, conjg(ompq)*t )
720 IF( rsvec ) THEN
721 CALL crot( mvl, v(1,p), 1,
722 $ v(1,q), 1, cs, conjg(ompq)*t )
723 END IF
724 sva( q ) = aaqq*sqrt( max( zero,
725 $ one+t*apoaq*aapq1 ) )
726 aapp = aapp*sqrt( max( zero,
727 $ one-t*aqoap*aapq1 ) )
728 mxsinj = max( mxsinj, abs( t ) )
729 ELSE
730*
731* .. choose correct signum for THETA and rotate
732*
733 thsign = -sign( one, aapq1 )
734 IF( aaqq.GT.aapp0 )thsign = -thsign
735 t = one / ( theta+thsign*
736 $ sqrt( one+theta*theta ) )
737 cs = sqrt( one / ( one+t*t ) )
738 sn = t*cs
739 mxsinj = max( mxsinj, abs( sn ) )
740 sva( q ) = aaqq*sqrt( max( zero,
741 $ one+t*apoaq*aapq1 ) )
742 aapp = aapp*sqrt( max( zero,
743 $ one-t*aqoap*aapq1 ) )
744*
745 CALL crot( m, a(1,p), 1, a(1,q), 1,
746 $ cs, conjg(ompq)*sn )
747 IF( rsvec ) THEN
748 CALL crot( mvl, v(1,p), 1,
749 $ v(1,q), 1, cs, conjg(ompq)*sn )
750 END IF
751 END IF
752 d(p) = -d(q) * ompq
753*
754 ELSE
755* .. have to use modified Gram-Schmidt like transformation
756 IF( aapp.GT.aaqq ) THEN
757 CALL ccopy( m, a( 1, p ), 1,
758 $ work, 1 )
759 CALL clascl( 'G', 0, 0, aapp, one,
760 $ m, 1, work,lda,
761 $ ierr )
762 CALL clascl( 'G', 0, 0, aaqq, one,
763 $ m, 1, a( 1, q ), lda,
764 $ ierr )
765 CALL caxpy( m, -aapq, work,
766 $ 1, a( 1, q ), 1 )
767 CALL clascl( 'G', 0, 0, one, aaqq,
768 $ m, 1, a( 1, q ), lda,
769 $ ierr )
770 sva( q ) = aaqq*sqrt( max( zero,
771 $ one-aapq1*aapq1 ) )
772 mxsinj = max( mxsinj, sfmin )
773 ELSE
774 CALL ccopy( m, a( 1, q ), 1,
775 $ work, 1 )
776 CALL clascl( 'G', 0, 0, aaqq, one,
777 $ m, 1, work,lda,
778 $ ierr )
779 CALL clascl( 'G', 0, 0, aapp, one,
780 $ m, 1, a( 1, p ), lda,
781 $ ierr )
782 CALL caxpy( m, -conjg(aapq),
783 $ work, 1, a( 1, p ), 1 )
784 CALL clascl( 'G', 0, 0, one, aapp,
785 $ m, 1, a( 1, p ), lda,
786 $ ierr )
787 sva( p ) = aapp*sqrt( max( zero,
788 $ one-aapq1*aapq1 ) )
789 mxsinj = max( mxsinj, sfmin )
790 END IF
791 END IF
792* END IF ROTOK THEN ... ELSE
793*
794* In the case of cancellation in updating SVA(q), SVA(p)
795* .. recompute SVA(q), SVA(p)
796 IF( ( sva( q ) / aaqq )**2.LE.rooteps )
797 $ THEN
798 IF( ( aaqq.LT.rootbig ) .AND.
799 $ ( aaqq.GT.rootsfmin ) ) THEN
800 sva( q ) = scnrm2( m, a( 1, q ), 1)
801 ELSE
802 t = zero
803 aaqq = one
804 CALL classq( m, a( 1, q ), 1, t,
805 $ aaqq )
806 sva( q ) = t*sqrt( aaqq )
807 END IF
808 END IF
809 IF( ( aapp / aapp0 )**2.LE.rooteps ) THEN
810 IF( ( aapp.LT.rootbig ) .AND.
811 $ ( aapp.GT.rootsfmin ) ) THEN
812 aapp = scnrm2( m, a( 1, p ), 1 )
813 ELSE
814 t = zero
815 aapp = one
816 CALL classq( m, a( 1, p ), 1, t,
817 $ aapp )
818 aapp = t*sqrt( aapp )
819 END IF
820 sva( p ) = aapp
821 END IF
822* end of OK rotation
823 ELSE
824 notrot = notrot + 1
825*[RTD] SKIPPED = SKIPPED + 1
826 pskipped = pskipped + 1
827 ijblsk = ijblsk + 1
828 END IF
829 ELSE
830 notrot = notrot + 1
831 pskipped = pskipped + 1
832 ijblsk = ijblsk + 1
833 END IF
834*
835 IF( ( i.LE.swband ) .AND. ( ijblsk.GE.blskip ) )
836 $ THEN
837 sva( p ) = aapp
838 notrot = 0
839 GO TO 2011
840 END IF
841 IF( ( i.LE.swband ) .AND.
842 $ ( pskipped.GT.rowskip ) ) THEN
843 aapp = -aapp
844 notrot = 0
845 GO TO 2203
846 END IF
847*
848 2200 CONTINUE
849* end of the q-loop
850 2203 CONTINUE
851*
852 sva( p ) = aapp
853*
854 ELSE
855*
856 IF( aapp.EQ.zero )notrot = notrot +
857 $ min( jgl+kbl-1, n ) - jgl + 1
858 IF( aapp.LT.zero )notrot = 0
859*
860 END IF
861*
862 2100 CONTINUE
863* end of the p-loop
864 2010 CONTINUE
865* end of the jbc-loop
866 2011 CONTINUE
867*2011 bailed out of the jbc-loop
868 DO 2012 p = igl, min( igl+kbl-1, n )
869 sva( p ) = abs( sva( p ) )
870 2012 CONTINUE
871***
872 2000 CONTINUE
873*2000 :: end of the ibr-loop
874*
875* .. update SVA(N)
876 IF( ( sva( n ).LT.rootbig ) .AND. ( sva( n ).GT.rootsfmin ) )
877 $ THEN
878 sva( n ) = scnrm2( m, a( 1, n ), 1 )
879 ELSE
880 t = zero
881 aapp = one
882 CALL classq( m, a( 1, n ), 1, t, aapp )
883 sva( n ) = t*sqrt( aapp )
884 END IF
885*
886* Additional steering devices
887*
888 IF( ( i.LT.swband ) .AND. ( ( mxaapq.LE.roottol ) .OR.
889 $ ( iswrot.LE.n ) ) )swband = i
890*
891 IF( ( i.GT.swband+1 ) .AND. ( mxaapq.LT.sqrt( real( n ) )*
892 $ tol ) .AND. ( real( n )*mxaapq*mxsinj.LT.tol ) ) THEN
893 GO TO 1994
894 END IF
895*
896 IF( notrot.GE.emptsw )GO TO 1994
897*
898 1993 CONTINUE
899* end i=1:NSWEEP loop
900*
901* #:( Reaching this point means that the procedure has not converged.
902 info = nsweep - 1
903 GO TO 1995
904*
905 1994 CONTINUE
906* #:) Reaching this point means numerical convergence after the i-th
907* sweep.
908*
909 info = 0
910* #:) INFO = 0 confirms successful iterations.
911 1995 CONTINUE
912*
913* Sort the vector SVA() of column norms.
914 DO 5991 p = 1, n - 1
915 q = isamax( n-p+1, sva( p ), 1 ) + p - 1
916 IF( p.NE.q ) THEN
917 temp1 = sva( p )
918 sva( p ) = sva( q )
919 sva( q ) = temp1
920 aapq = d( p )
921 d( p ) = d( q )
922 d( q ) = aapq
923 CALL cswap( m, a( 1, p ), 1, a( 1, q ), 1 )
924 IF( rsvec )CALL cswap( mvl, v( 1, p ), 1, v( 1, q ), 1 )
925 END IF
926 5991 CONTINUE
927*
928 RETURN
929* ..
930* .. END OF CGSVJ0
931* ..
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
caxpy
subroutine caxpy(n, ca, cx, incx, cy, incy)
CAXPY
Definition caxpy.f:88
ccopy
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
cdotc
complex function cdotc(n, cx, incx, cy, incy)
CDOTC
Definition cdotc.f:83
isamax
integer function isamax(n, sx, incx)
ISAMAX
Definition isamax.f:71
clascl
subroutine clascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition clascl.f:143
classq
subroutine classq(n, x, incx, scale, sumsq)
CLASSQ updates a sum of squares represented in scaled form.
Definition classq.f90:124
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
scnrm2
real(wp) function scnrm2(n, x, incx)
SCNRM2
Definition scnrm2.f90:90
crot
subroutine crot(n, cx, incx, cy, incy, c, s)
CROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors.
Definition crot.f:103
cswap
subroutine cswap(n, cx, incx, cy, incy)
CSWAP
Definition cswap.f:81
Here is the call graph for this function:
Here is the caller graph for this function: