 LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ cgsvj1()

 subroutine cgsvj1 ( character*1 JOBV, integer M, integer N, integer N1, 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 )

CGSVJ1 pre-processor for the routine cgesvj, applies Jacobi rotations targeting only particular pivots.

Purpose:
``` CGSVJ1 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 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
~~~~~~~~~~~~~~~
CGSVJ1 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 postmulyiplying the N-by-N array V. (See the description of V.) = 'A': the product of the Jacobi rotations is accumulated by postmulyiplying 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 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 * 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 COMPLEX 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 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 onexit*diag(D_onexit).``` [in] MV ``` MV is INTEGER If JOBV = 'A', then MV rows of V are post-multipled 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-multipled by a sequence of Jacobi rotations. If JOBV = 'A' then MV rows of V are post-multipled 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```
Contributor:
Zlatko Drmac (Zagreb, Croatia)

Definition at line 234 of file cgsvj1.f.

236*
237* -- LAPACK computational routine --
238* -- LAPACK is a software package provided by Univ. of Tennessee, --
239* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
240*
241* .. Scalar Arguments ..
242 REAL EPS, SFMIN, TOL
243 INTEGER INFO, LDA, LDV, LWORK, M, MV, N, N1, NSWEEP
244 CHARACTER*1 JOBV
245* ..
246* .. Array Arguments ..
247 COMPLEX A( LDA, * ), D( N ), V( LDV, * ), WORK( LWORK )
248 REAL SVA( N )
249* ..
250*
251* =====================================================================
252*
253* .. Local Parameters ..
254 REAL ZERO, HALF, ONE
255 parameter( zero = 0.0e0, half = 0.5e0, one = 1.0e0)
256* ..
257* .. Local Scalars ..
258 COMPLEX AAPQ, OMPQ
259 REAL AAPP, AAPP0, AAPQ1, AAQQ, APOAQ, AQOAP, BIG,
260 \$ BIGTHETA, CS, MXAAPQ, MXSINJ, ROOTBIG,
261 \$ ROOTEPS, ROOTSFMIN, ROOTTOL, SMALL, SN, T,
262 \$ TEMP1, THETA, THSIGN
263 INTEGER BLSKIP, EMPTSW, i, ibr, igl, IERR, IJBLSK,
264 \$ ISWROT, jbc, jgl, KBL, MVL, NOTROT, nblc, nblr,
265 \$ p, PSKIPPED, q, ROWSKIP, SWBAND
266 LOGICAL APPLV, ROTOK, RSVEC
267* ..
268* ..
269* .. Intrinsic Functions ..
270 INTRINSIC abs, max, conjg, real, min, sign, sqrt
271* ..
272* .. External Functions ..
273 REAL SCNRM2
274 COMPLEX CDOTC
275 INTEGER ISAMAX
276 LOGICAL LSAME
277 EXTERNAL isamax, lsame, cdotc, scnrm2
278* ..
279* .. External Subroutines ..
280* .. from BLAS
281 EXTERNAL ccopy, crot, cswap, caxpy
282* .. from LAPACK
283 EXTERNAL clascl, classq, xerbla
284* ..
285* .. Executable Statements ..
286*
287* Test the input parameters.
288*
289 applv = lsame( jobv, 'A' )
290 rsvec = lsame( jobv, 'V' )
291 IF( .NOT.( rsvec .OR. applv .OR. lsame( jobv, 'N' ) ) ) THEN
292 info = -1
293 ELSE IF( m.LT.0 ) THEN
294 info = -2
295 ELSE IF( ( n.LT.0 ) .OR. ( n.GT.m ) ) THEN
296 info = -3
297 ELSE IF( n1.LT.0 ) THEN
298 info = -4
299 ELSE IF( lda.LT.m ) THEN
300 info = -6
301 ELSE IF( ( rsvec.OR.applv ) .AND. ( mv.LT.0 ) ) THEN
302 info = -9
303 ELSE IF( ( rsvec.AND.( ldv.LT.n ) ).OR.
304 \$ ( applv.AND.( ldv.LT.mv ) ) ) THEN
305 info = -11
306 ELSE IF( tol.LE.eps ) THEN
307 info = -14
308 ELSE IF( nsweep.LT.0 ) THEN
309 info = -15
310 ELSE IF( lwork.LT.m ) THEN
311 info = -17
312 ELSE
313 info = 0
314 END IF
315*
316* #:(
317 IF( info.NE.0 ) THEN
318 CALL xerbla( 'CGSVJ1', -info )
319 RETURN
320 END IF
321*
322 IF( rsvec ) THEN
323 mvl = n
324 ELSE IF( applv ) THEN
325 mvl = mv
326 END IF
327 rsvec = rsvec .OR. applv
328
329 rooteps = sqrt( eps )
330 rootsfmin = sqrt( sfmin )
331 small = sfmin / eps
332 big = one / sfmin
333 rootbig = one / rootsfmin
334* LARGE = BIG / SQRT( REAL( M*N ) )
335 bigtheta = one / rooteps
336 roottol = sqrt( tol )
337*
338* .. Initialize the right singular vector matrix ..
339*
340* RSVEC = LSAME( JOBV, 'Y' )
341*
342 emptsw = n1*( n-n1 )
343 notrot = 0
344*
345* .. Row-cyclic pivot strategy with de Rijk's pivoting ..
346*
347 kbl = min( 8, n )
348 nblr = n1 / kbl
349 IF( ( nblr*kbl ).NE.n1 )nblr = nblr + 1
350
351* .. the tiling is nblr-by-nblc [tiles]
352
353 nblc = ( n-n1 ) / kbl
354 IF( ( nblc*kbl ).NE.( n-n1 ) )nblc = nblc + 1
355 blskip = ( kbl**2 ) + 1
356*[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
357
358 rowskip = min( 5, kbl )
359*[TP] ROWSKIP is a tuning parameter.
360 swband = 0
361*[TP] SWBAND is a tuning parameter. It is meaningful and effective
362* if CGESVJ is used as a computational routine in the preconditioned
363* Jacobi SVD algorithm CGEJSV.
364*
365*
366* | * * * [x] [x] [x]|
367* | * * * [x] [x] [x]| Row-cycling in the nblr-by-nblc [x] blocks.
368* | * * * [x] [x] [x]| Row-cyclic pivoting inside each [x] block.
369* |[x] [x] [x] * * * |
370* |[x] [x] [x] * * * |
371* |[x] [x] [x] * * * |
372*
373*
374 DO 1993 i = 1, nsweep
375*
376* .. go go go ...
377*
378 mxaapq = zero
379 mxsinj = zero
380 iswrot = 0
381*
382 notrot = 0
383 pskipped = 0
384*
385* Each sweep is unrolled using KBL-by-KBL tiles over the pivot pairs
386* 1 <= p < q <= N. This is the first step toward a blocked implementation
387* of the rotations. New implementation, based on block transformations,
388* is under development.
389*
390 DO 2000 ibr = 1, nblr
391*
392 igl = ( ibr-1 )*kbl + 1
393*
394
395*
396* ... go to the off diagonal blocks
397*
398 igl = ( ibr-1 )*kbl + 1
399*
400* DO 2010 jbc = ibr + 1, NBL
401 DO 2010 jbc = 1, nblc
402*
403 jgl = ( jbc-1 )*kbl + n1 + 1
404*
405* doing the block at ( ibr, jbc )
406*
407 ijblsk = 0
408 DO 2100 p = igl, min( igl+kbl-1, n1 )
409*
410 aapp = sva( p )
411 IF( aapp.GT.zero ) THEN
412*
413 pskipped = 0
414*
415 DO 2200 q = jgl, min( jgl+kbl-1, n )
416*
417 aaqq = sva( q )
418 IF( aaqq.GT.zero ) THEN
419 aapp0 = aapp
420*
421* .. M x 2 Jacobi SVD ..
422*
423* Safe Gram matrix computation
424*
425 IF( aaqq.GE.one ) THEN
426 IF( aapp.GE.aaqq ) THEN
427 rotok = ( small*aapp ).LE.aaqq
428 ELSE
429 rotok = ( small*aaqq ).LE.aapp
430 END IF
431 IF( aapp.LT.( big / aaqq ) ) THEN
432 aapq = ( cdotc( m, a( 1, p ), 1,
433 \$ a( 1, q ), 1 ) / aaqq ) / aapp
434 ELSE
435 CALL ccopy( m, a( 1, p ), 1,
436 \$ work, 1 )
437 CALL clascl( 'G', 0, 0, aapp,
438 \$ one, m, 1,
439 \$ work, lda, ierr )
440 aapq = cdotc( m, work, 1,
441 \$ a( 1, q ), 1 ) / aaqq
442 END IF
443 ELSE
444 IF( aapp.GE.aaqq ) THEN
445 rotok = aapp.LE.( aaqq / small )
446 ELSE
447 rotok = aaqq.LE.( aapp / small )
448 END IF
449 IF( aapp.GT.( small / aaqq ) ) THEN
450 aapq = ( cdotc( m, a( 1, p ), 1,
451 \$ a( 1, q ), 1 ) / max(aaqq,aapp) )
452 \$ / min(aaqq,aapp)
453 ELSE
454 CALL ccopy( m, a( 1, q ), 1,
455 \$ work, 1 )
456 CALL clascl( 'G', 0, 0, aaqq,
457 \$ one, m, 1,
458 \$ work, lda, ierr )
459 aapq = cdotc( m, a( 1, p ), 1,
460 \$ work, 1 ) / aapp
461 END IF
462 END IF
463*
464* AAPQ = AAPQ * CONJG(CWORK(p))*CWORK(q)
465 aapq1 = -abs(aapq)
466 mxaapq = max( mxaapq, -aapq1 )
467*
468* TO rotate or NOT to rotate, THAT is the question ...
469*
470 IF( abs( aapq1 ).GT.tol ) THEN
471 ompq = aapq / abs(aapq)
472 notrot = 0
473*[RTD] ROTATED = ROTATED + 1
474 pskipped = 0
475 iswrot = iswrot + 1
476*
477 IF( rotok ) THEN
478*
479 aqoap = aaqq / aapp
480 apoaq = aapp / aaqq
481 theta = -half*abs( aqoap-apoaq )/ aapq1
482 IF( aaqq.GT.aapp0 )theta = -theta
483*
484 IF( abs( theta ).GT.bigtheta ) THEN
485 t = half / theta
486 cs = one
487 CALL crot( m, a(1,p), 1, a(1,q), 1,
488 \$ cs, conjg(ompq)*t )
489 IF( rsvec ) THEN
490 CALL crot( mvl, v(1,p), 1,
491 \$ v(1,q), 1, cs, conjg(ompq)*t )
492 END IF
493 sva( q ) = aaqq*sqrt( max( zero,
494 \$ one+t*apoaq*aapq1 ) )
495 aapp = aapp*sqrt( max( zero,
496 \$ one-t*aqoap*aapq1 ) )
497 mxsinj = max( mxsinj, abs( t ) )
498 ELSE
499*
500* .. choose correct signum for THETA and rotate
501*
502 thsign = -sign( one, aapq1 )
503 IF( aaqq.GT.aapp0 )thsign = -thsign
504 t = one / ( theta+thsign*
505 \$ sqrt( one+theta*theta ) )
506 cs = sqrt( one / ( one+t*t ) )
507 sn = t*cs
508 mxsinj = max( mxsinj, abs( sn ) )
509 sva( q ) = aaqq*sqrt( max( zero,
510 \$ one+t*apoaq*aapq1 ) )
511 aapp = aapp*sqrt( max( zero,
512 \$ one-t*aqoap*aapq1 ) )
513*
514 CALL crot( m, a(1,p), 1, a(1,q), 1,
515 \$ cs, conjg(ompq)*sn )
516 IF( rsvec ) THEN
517 CALL crot( mvl, v(1,p), 1,
518 \$ v(1,q), 1, cs, conjg(ompq)*sn )
519 END IF
520 END IF
521 d(p) = -d(q) * ompq
522*
523 ELSE
524* .. have to use modified Gram-Schmidt like transformation
525 IF( aapp.GT.aaqq ) THEN
526 CALL ccopy( m, a( 1, p ), 1,
527 \$ work, 1 )
528 CALL clascl( 'G', 0, 0, aapp, one,
529 \$ m, 1, work,lda,
530 \$ ierr )
531 CALL clascl( 'G', 0, 0, aaqq, one,
532 \$ m, 1, a( 1, q ), lda,
533 \$ ierr )
534 CALL caxpy( m, -aapq, work,
535 \$ 1, a( 1, q ), 1 )
536 CALL clascl( 'G', 0, 0, one, aaqq,
537 \$ m, 1, a( 1, q ), lda,
538 \$ ierr )
539 sva( q ) = aaqq*sqrt( max( zero,
540 \$ one-aapq1*aapq1 ) )
541 mxsinj = max( mxsinj, sfmin )
542 ELSE
543 CALL ccopy( m, a( 1, q ), 1,
544 \$ work, 1 )
545 CALL clascl( 'G', 0, 0, aaqq, one,
546 \$ m, 1, work,lda,
547 \$ ierr )
548 CALL clascl( 'G', 0, 0, aapp, one,
549 \$ m, 1, a( 1, p ), lda,
550 \$ ierr )
551 CALL caxpy( m, -conjg(aapq),
552 \$ work, 1, a( 1, p ), 1 )
553 CALL clascl( 'G', 0, 0, one, aapp,
554 \$ m, 1, a( 1, p ), lda,
555 \$ ierr )
556 sva( p ) = aapp*sqrt( max( zero,
557 \$ one-aapq1*aapq1 ) )
558 mxsinj = max( mxsinj, sfmin )
559 END IF
560 END IF
561* END IF ROTOK THEN ... ELSE
562*
563* In the case of cancellation in updating SVA(q), SVA(p)
564* .. recompute SVA(q), SVA(p)
565 IF( ( sva( q ) / aaqq )**2.LE.rooteps )
566 \$ THEN
567 IF( ( aaqq.LT.rootbig ) .AND.
568 \$ ( aaqq.GT.rootsfmin ) ) THEN
569 sva( q ) = scnrm2( m, a( 1, q ), 1)
570 ELSE
571 t = zero
572 aaqq = one
573 CALL classq( m, a( 1, q ), 1, t,
574 \$ aaqq )
575 sva( q ) = t*sqrt( aaqq )
576 END IF
577 END IF
578 IF( ( aapp / aapp0 )**2.LE.rooteps ) THEN
579 IF( ( aapp.LT.rootbig ) .AND.
580 \$ ( aapp.GT.rootsfmin ) ) THEN
581 aapp = scnrm2( m, a( 1, p ), 1 )
582 ELSE
583 t = zero
584 aapp = one
585 CALL classq( m, a( 1, p ), 1, t,
586 \$ aapp )
587 aapp = t*sqrt( aapp )
588 END IF
589 sva( p ) = aapp
590 END IF
591* end of OK rotation
592 ELSE
593 notrot = notrot + 1
594*[RTD] SKIPPED = SKIPPED + 1
595 pskipped = pskipped + 1
596 ijblsk = ijblsk + 1
597 END IF
598 ELSE
599 notrot = notrot + 1
600 pskipped = pskipped + 1
601 ijblsk = ijblsk + 1
602 END IF
603*
604 IF( ( i.LE.swband ) .AND. ( ijblsk.GE.blskip ) )
605 \$ THEN
606 sva( p ) = aapp
607 notrot = 0
608 GO TO 2011
609 END IF
610 IF( ( i.LE.swband ) .AND.
611 \$ ( pskipped.GT.rowskip ) ) THEN
612 aapp = -aapp
613 notrot = 0
614 GO TO 2203
615 END IF
616*
617 2200 CONTINUE
618* end of the q-loop
619 2203 CONTINUE
620*
621 sva( p ) = aapp
622*
623 ELSE
624*
625 IF( aapp.EQ.zero )notrot = notrot +
626 \$ min( jgl+kbl-1, n ) - jgl + 1
627 IF( aapp.LT.zero )notrot = 0
628*
629 END IF
630*
631 2100 CONTINUE
632* end of the p-loop
633 2010 CONTINUE
634* end of the jbc-loop
635 2011 CONTINUE
636*2011 bailed out of the jbc-loop
637 DO 2012 p = igl, min( igl+kbl-1, n )
638 sva( p ) = abs( sva( p ) )
639 2012 CONTINUE
640***
641 2000 CONTINUE
642*2000 :: end of the ibr-loop
643*
644* .. update SVA(N)
645 IF( ( sva( n ).LT.rootbig ) .AND. ( sva( n ).GT.rootsfmin ) )
646 \$ THEN
647 sva( n ) = scnrm2( m, a( 1, n ), 1 )
648 ELSE
649 t = zero
650 aapp = one
651 CALL classq( m, a( 1, n ), 1, t, aapp )
652 sva( n ) = t*sqrt( aapp )
653 END IF
654*
656*
657 IF( ( i.LT.swband ) .AND. ( ( mxaapq.LE.roottol ) .OR.
658 \$ ( iswrot.LE.n ) ) )swband = i
659*
660 IF( ( i.GT.swband+1 ) .AND. ( mxaapq.LT.sqrt( real( n ) )*
661 \$ tol ) .AND. ( real( n )*mxaapq*mxsinj.LT.tol ) ) THEN
662 GO TO 1994
663 END IF
664*
665 IF( notrot.GE.emptsw )GO TO 1994
666*
667 1993 CONTINUE
668* end i=1:NSWEEP loop
669*
670* #:( Reaching this point means that the procedure has not converged.
671 info = nsweep - 1
672 GO TO 1995
673*
674 1994 CONTINUE
675* #:) Reaching this point means numerical convergence after the i-th
676* sweep.
677*
678 info = 0
679* #:) INFO = 0 confirms successful iterations.
680 1995 CONTINUE
681*
682* Sort the vector SVA() of column norms.
683 DO 5991 p = 1, n - 1
684 q = isamax( n-p+1, sva( p ), 1 ) + p - 1
685 IF( p.NE.q ) THEN
686 temp1 = sva( p )
687 sva( p ) = sva( q )
688 sva( q ) = temp1
689 aapq = d( p )
690 d( p ) = d( q )
691 d( q ) = aapq
692 CALL cswap( m, a( 1, p ), 1, a( 1, q ), 1 )
693 IF( rsvec )CALL cswap( mvl, v( 1, p ), 1, v( 1, q ), 1 )
694 END IF
695 5991 CONTINUE
696*
697*
698 RETURN
699* ..
700* .. END OF CGSVJ1
701* ..
subroutine classq(n, x, incx, scl, sumsq)
CLASSQ updates a sum of squares represented in scaled form.
Definition: classq.f90:137
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:71
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
complex function cdotc(N, CX, INCX, CY, INCY)
CDOTC
Definition: cdotc.f:83
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
Definition: ccopy.f:81
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:88
subroutine cswap(N, CX, INCX, CY, INCY)
CSWAP
Definition: cswap.f:81
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
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
real(wp) function scnrm2(n, x, incx)
SCNRM2
Definition: scnrm2.f90:90
Here is the call graph for this function:
Here is the caller graph for this function: