LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
sgsvj1.f
Go to the documentation of this file.
00001 *> \brief \b SGSVJ1
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download SGSVJ1 + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgsvj1.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgsvj1.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgsvj1.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE SGSVJ1( JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV,
00022 *                          EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       REAL               EPS, SFMIN, TOL
00026 *       INTEGER            INFO, LDA, LDV, LWORK, M, MV, N, N1, NSWEEP
00027 *       CHARACTER*1        JOBV
00028 *       ..
00029 *       .. Array Arguments ..
00030 *       REAL               A( LDA, * ), D( N ), SVA( N ), V( LDV, * ),
00031 *      $                   WORK( LWORK )
00032 *       ..
00033 *  
00034 *
00035 *> \par Purpose:
00036 *  =============
00037 *>
00038 *> \verbatim
00039 *>
00040 *> SGSVJ1 is called from SGESVJ as a pre-processor and that is its main
00041 *> purpose. It applies Jacobi rotations in the same way as SGESVJ does, but
00042 *> it targets only particular pivots and it does not check convergence
00043 *> (stopping criterion). Few tunning parameters (marked by [TP]) are
00044 *> available for the implementer.
00045 *>
00046 *> Further Details
00047 *> ~~~~~~~~~~~~~~~
00048 *> SGSVJ1 applies few sweeps of Jacobi rotations in the column space of
00049 *> the input M-by-N matrix A. The pivot pairs are taken from the (1,2)
00050 *> off-diagonal block in the corresponding N-by-N Gram matrix A^T * A. The
00051 *> block-entries (tiles) of the (1,2) off-diagonal block are marked by the
00052 *> [x]'s in the following scheme:
00053 *>
00054 *>    | *  *  * [x] [x] [x]|
00055 *>    | *  *  * [x] [x] [x]|    Row-cycling in the nblr-by-nblc [x] blocks.
00056 *>    | *  *  * [x] [x] [x]|    Row-cyclic pivoting inside each [x] block.
00057 *>    |[x] [x] [x] *  *  * |
00058 *>    |[x] [x] [x] *  *  * |
00059 *>    |[x] [x] [x] *  *  * |
00060 *>
00061 *> In terms of the columns of A, the first N1 columns are rotated 'against'
00062 *> the remaining N-N1 columns, trying to increase the angle between the
00063 *> corresponding subspaces. The off-diagonal block is N1-by(N-N1) and it is
00064 *> tiled using quadratic tiles of side KBL. Here, KBL is a tunning parmeter.
00065 *> The number of sweeps is given in NSWEEP and the orthogonality threshold
00066 *> is given in TOL.
00067 *> \endverbatim
00068 *
00069 *  Arguments:
00070 *  ==========
00071 *
00072 *> \param[in] JOBV
00073 *> \verbatim
00074 *>          JOBV is CHARACTER*1
00075 *>          Specifies whether the output from this procedure is used
00076 *>          to compute the matrix V:
00077 *>          = 'V': the product of the Jacobi rotations is accumulated
00078 *>                 by postmulyiplying the N-by-N array V.
00079 *>                (See the description of V.)
00080 *>          = 'A': the product of the Jacobi rotations is accumulated
00081 *>                 by postmulyiplying the MV-by-N array V.
00082 *>                (See the descriptions of MV and V.)
00083 *>          = 'N': the Jacobi rotations are not accumulated.
00084 *> \endverbatim
00085 *>
00086 *> \param[in] M
00087 *> \verbatim
00088 *>          M is INTEGER
00089 *>          The number of rows of the input matrix A.  M >= 0.
00090 *> \endverbatim
00091 *>
00092 *> \param[in] N
00093 *> \verbatim
00094 *>          N is INTEGER
00095 *>          The number of columns of the input matrix A.
00096 *>          M >= N >= 0.
00097 *> \endverbatim
00098 *>
00099 *> \param[in] N1
00100 *> \verbatim
00101 *>          N1 is INTEGER
00102 *>          N1 specifies the 2 x 2 block partition, the first N1 columns are
00103 *>          rotated 'against' the remaining N-N1 columns of A.
00104 *> \endverbatim
00105 *>
00106 *> \param[in,out] A
00107 *> \verbatim
00108 *>          A is REAL array, dimension (LDA,N)
00109 *>          On entry, M-by-N matrix A, such that A*diag(D) represents
00110 *>          the input matrix.
00111 *>          On exit,
00112 *>          A_onexit * D_onexit represents the input matrix A*diag(D)
00113 *>          post-multiplied by a sequence of Jacobi rotations, where the
00114 *>          rotation threshold and the total number of sweeps are given in
00115 *>          TOL and NSWEEP, respectively.
00116 *>          (See the descriptions of N1, D, TOL and NSWEEP.)
00117 *> \endverbatim
00118 *>
00119 *> \param[in] LDA
00120 *> \verbatim
00121 *>          LDA is INTEGER
00122 *>          The leading dimension of the array A.  LDA >= max(1,M).
00123 *> \endverbatim
00124 *>
00125 *> \param[in,out] D
00126 *> \verbatim
00127 *>          D is REAL array, dimension (N)
00128 *>          The array D accumulates the scaling factors from the fast scaled
00129 *>          Jacobi rotations.
00130 *>          On entry, A*diag(D) represents the input matrix.
00131 *>          On exit, A_onexit*diag(D_onexit) represents the input matrix
00132 *>          post-multiplied by a sequence of Jacobi rotations, where the
00133 *>          rotation threshold and the total number of sweeps are given in
00134 *>          TOL and NSWEEP, respectively.
00135 *>          (See the descriptions of N1, A, TOL and NSWEEP.)
00136 *> \endverbatim
00137 *>
00138 *> \param[in,out] SVA
00139 *> \verbatim
00140 *>          SVA is REAL array, dimension (N)
00141 *>          On entry, SVA contains the Euclidean norms of the columns of
00142 *>          the matrix A*diag(D).
00143 *>          On exit, SVA contains the Euclidean norms of the columns of
00144 *>          the matrix onexit*diag(D_onexit).
00145 *> \endverbatim
00146 *>
00147 *> \param[in] MV
00148 *> \verbatim
00149 *>          MV is INTEGER
00150 *>          If JOBV .EQ. 'A', then MV rows of V are post-multipled by a
00151 *>                           sequence of Jacobi rotations.
00152 *>          If JOBV = 'N',   then MV is not referenced.
00153 *> \endverbatim
00154 *>
00155 *> \param[in,out] V
00156 *> \verbatim
00157 *>          V is REAL array, dimension (LDV,N)
00158 *>          If JOBV .EQ. 'V' then N rows of V are post-multipled by a
00159 *>                           sequence of Jacobi rotations.
00160 *>          If JOBV .EQ. 'A' then MV rows of V are post-multipled by a
00161 *>                           sequence of Jacobi rotations.
00162 *>          If JOBV = 'N',   then V is not referenced.
00163 *> \endverbatim
00164 *>
00165 *> \param[in] LDV
00166 *> \verbatim
00167 *>          LDV is INTEGER
00168 *>          The leading dimension of the array V,  LDV >= 1.
00169 *>          If JOBV = 'V', LDV .GE. N.
00170 *>          If JOBV = 'A', LDV .GE. MV.
00171 *> \endverbatim
00172 *>
00173 *> \param[in] EPS
00174 *> \verbatim
00175 *>          EPS is INTEGER
00176 *>          EPS = SLAMCH('Epsilon')
00177 *> \endverbatim
00178 *>
00179 *> \param[in] SFMIN
00180 *> \verbatim
00181 *>          SFMIN is INTEGER
00182 *>          SFMIN = SLAMCH('Safe Minimum')
00183 *> \endverbatim
00184 *>
00185 *> \param[in] TOL
00186 *> \verbatim
00187 *>          TOL is REAL
00188 *>          TOL is the threshold for Jacobi rotations. For a pair
00189 *>          A(:,p), A(:,q) of pivot columns, the Jacobi rotation is
00190 *>          applied only if ABS(COS(angle(A(:,p),A(:,q)))) .GT. TOL.
00191 *> \endverbatim
00192 *>
00193 *> \param[in] NSWEEP
00194 *> \verbatim
00195 *>          NSWEEP is INTEGER
00196 *>          NSWEEP is the number of sweeps of Jacobi rotations to be
00197 *>          performed.
00198 *> \endverbatim
00199 *>
00200 *> \param[out] WORK
00201 *> \verbatim
00202 *>          WORK is REAL array, dimension LWORK.
00203 *> \endverbatim
00204 *>
00205 *> \param[in] LWORK
00206 *> \verbatim
00207 *>          LWORK is INTEGER
00208 *>          LWORK is the dimension of WORK. LWORK .GE. M.
00209 *> \endverbatim
00210 *>
00211 *> \param[out] INFO
00212 *> \verbatim
00213 *>          INFO is INTEGER
00214 *>          = 0 : successful exit.
00215 *>          < 0 : if INFO = -i, then the i-th argument had an illegal value
00216 *> \endverbatim
00217 *
00218 *  Authors:
00219 *  ========
00220 *
00221 *> \author Univ. of Tennessee 
00222 *> \author Univ. of California Berkeley 
00223 *> \author Univ. of Colorado Denver 
00224 *> \author NAG Ltd. 
00225 *
00226 *> \date November 2011
00227 *
00228 *> \ingroup realOTHERcomputational
00229 *
00230 *> \par Contributors:
00231 *  ==================
00232 *>
00233 *> Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany)
00234 *
00235 *  =====================================================================
00236       SUBROUTINE SGSVJ1( JOBV, M, N, N1, A, LDA, D, SVA, MV, V, LDV,
00237      $                   EPS, SFMIN, TOL, NSWEEP, WORK, LWORK, INFO )
00238 *
00239 *  -- LAPACK computational routine (version 3.4.0) --
00240 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00241 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00242 *     November 2011
00243 *
00244 *     .. Scalar Arguments ..
00245       REAL               EPS, SFMIN, TOL
00246       INTEGER            INFO, LDA, LDV, LWORK, M, MV, N, N1, NSWEEP
00247       CHARACTER*1        JOBV
00248 *     ..
00249 *     .. Array Arguments ..
00250       REAL               A( LDA, * ), D( N ), SVA( N ), V( LDV, * ),
00251      $                   WORK( LWORK )
00252 *     ..
00253 *
00254 *  =====================================================================
00255 *
00256 *     .. Local Parameters ..
00257       REAL               ZERO, HALF, ONE
00258       PARAMETER          ( ZERO = 0.0E0, HALF = 0.5E0, ONE = 1.0E0)
00259 *     ..
00260 *     .. Local Scalars ..
00261       REAL               AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG,
00262      $                   BIGTHETA, CS, LARGE, MXAAPQ, MXSINJ, ROOTBIG,
00263      $                   ROOTEPS, ROOTSFMIN, ROOTTOL, SMALL, SN, T,
00264      $                   TEMP1, THETA, THSIGN
00265       INTEGER            BLSKIP, EMPTSW, i, ibr, igl, IERR, IJBLSK,
00266      $                   ISWROT, jbc, jgl, KBL, MVL, NOTROT, nblc, nblr,
00267      $                   p, PSKIPPED, q, ROWSKIP, SWBAND
00268       LOGICAL            APPLV, ROTOK, RSVEC
00269 *     ..
00270 *     .. Local Arrays ..
00271       REAL               FASTR( 5 )
00272 *     ..
00273 *     .. Intrinsic Functions ..
00274       INTRINSIC          ABS, AMAX1, FLOAT, MIN0, SIGN, SQRT
00275 *     ..
00276 *     .. External Functions ..
00277       REAL               SDOT, SNRM2
00278       INTEGER            ISAMAX
00279       LOGICAL            LSAME
00280       EXTERNAL           ISAMAX, LSAME, SDOT, SNRM2
00281 *     ..
00282 *     .. External Subroutines ..
00283       EXTERNAL           SAXPY, SCOPY, SLASCL, SLASSQ, SROTM, SSWAP
00284 *     ..
00285 *     .. Executable Statements ..
00286 *
00287 *     Test the input parameters.
00288 *
00289       APPLV = LSAME( JOBV, 'A' )
00290       RSVEC = LSAME( JOBV, 'V' )
00291       IF( .NOT.( RSVEC .OR. APPLV .OR. LSAME( JOBV, 'N' ) ) ) THEN
00292          INFO = -1
00293       ELSE IF( M.LT.0 ) THEN
00294          INFO = -2
00295       ELSE IF( ( N.LT.0 ) .OR. ( N.GT.M ) ) THEN
00296          INFO = -3
00297       ELSE IF( N1.LT.0 ) THEN
00298          INFO = -4
00299       ELSE IF( LDA.LT.M ) THEN
00300          INFO = -6
00301       ELSE IF( ( RSVEC.OR.APPLV ) .AND. ( MV.LT.0 ) ) THEN
00302          INFO = -9
00303       ELSE IF( ( RSVEC.AND.( LDV.LT.N ) ).OR. 
00304      $         ( APPLV.AND.( LDV.LT.MV ) )  ) THEN
00305          INFO = -11
00306       ELSE IF( TOL.LE.EPS ) THEN
00307          INFO = -14
00308       ELSE IF( NSWEEP.LT.0 ) THEN
00309          INFO = -15
00310       ELSE IF( LWORK.LT.M ) THEN
00311          INFO = -17
00312       ELSE
00313          INFO = 0
00314       END IF
00315 *
00316 *     #:(
00317       IF( INFO.NE.0 ) THEN
00318          CALL XERBLA( 'SGSVJ1', -INFO )
00319          RETURN
00320       END IF
00321 *
00322       IF( RSVEC ) THEN
00323          MVL = N
00324       ELSE IF( APPLV ) THEN
00325          MVL = MV
00326       END IF
00327       RSVEC = RSVEC .OR. APPLV
00328 
00329       ROOTEPS = SQRT( EPS )
00330       ROOTSFMIN = SQRT( SFMIN )
00331       SMALL = SFMIN / EPS
00332       BIG = ONE / SFMIN
00333       ROOTBIG = ONE / ROOTSFMIN
00334       LARGE = BIG / SQRT( FLOAT( M*N ) )
00335       BIGTHETA = ONE / ROOTEPS
00336       ROOTTOL = SQRT( TOL )
00337 *
00338 *     .. Initialize the right singular vector matrix ..
00339 *
00340 *     RSVEC = LSAME( JOBV, 'Y' )
00341 *
00342       EMPTSW = N1*( N-N1 )
00343       NOTROT = 0
00344       FASTR( 1 ) = ZERO
00345 *
00346 *     .. Row-cyclic pivot strategy with de Rijk's pivoting ..
00347 *
00348       KBL = MIN0( 8, N )
00349       NBLR = N1 / KBL
00350       IF( ( NBLR*KBL ).NE.N1 )NBLR = NBLR + 1
00351 
00352 *     .. the tiling is nblr-by-nblc [tiles]
00353 
00354       NBLC = ( N-N1 ) / KBL
00355       IF( ( NBLC*KBL ).NE.( N-N1 ) )NBLC = NBLC + 1
00356       BLSKIP = ( KBL**2 ) + 1
00357 *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL.
00358 
00359       ROWSKIP = MIN0( 5, KBL )
00360 *[TP] ROWSKIP is a tuning parameter.
00361       SWBAND = 0
00362 *[TP] SWBAND is a tuning parameter. It is meaningful and effective
00363 *     if SGESVJ is used as a computational routine in the preconditioned
00364 *     Jacobi SVD algorithm SGESVJ.
00365 *
00366 *
00367 *     | *   *   * [x] [x] [x]|
00368 *     | *   *   * [x] [x] [x]|    Row-cycling in the nblr-by-nblc [x] blocks.
00369 *     | *   *   * [x] [x] [x]|    Row-cyclic pivoting inside each [x] block.
00370 *     |[x] [x] [x] *   *   * |
00371 *     |[x] [x] [x] *   *   * |
00372 *     |[x] [x] [x] *   *   * |
00373 *
00374 *
00375       DO 1993 i = 1, NSWEEP
00376 *     .. go go go ...
00377 *
00378          MXAAPQ = ZERO
00379          MXSINJ = ZERO
00380          ISWROT = 0
00381 *
00382          NOTROT = 0
00383          PSKIPPED = 0
00384 *
00385          DO 2000 ibr = 1, NBLR
00386 
00387             igl = ( ibr-1 )*KBL + 1
00388 *
00389 *
00390 *........................................................
00391 * ... go to the off diagonal blocks
00392 
00393             igl = ( ibr-1 )*KBL + 1
00394 
00395             DO 2010 jbc = 1, NBLC
00396 
00397                jgl = N1 + ( jbc-1 )*KBL + 1
00398 
00399 *        doing the block at ( ibr, jbc )
00400 
00401                IJBLSK = 0
00402                DO 2100 p = igl, MIN0( igl+KBL-1, N1 )
00403 
00404                   AAPP = SVA( p )
00405 
00406                   IF( AAPP.GT.ZERO ) THEN
00407 
00408                      PSKIPPED = 0
00409 
00410                      DO 2200 q = jgl, MIN0( jgl+KBL-1, N )
00411 *
00412                         AAQQ = SVA( q )
00413 
00414                         IF( AAQQ.GT.ZERO ) THEN
00415                            AAPP0 = AAPP
00416 *
00417 *     .. M x 2 Jacobi SVD ..
00418 *
00419 *        .. Safe Gram matrix computation ..
00420 *
00421                            IF( AAQQ.GE.ONE ) THEN
00422                               IF( AAPP.GE.AAQQ ) THEN
00423                                  ROTOK = ( SMALL*AAPP ).LE.AAQQ
00424                               ELSE
00425                                  ROTOK = ( SMALL*AAQQ ).LE.AAPP
00426                               END IF
00427                               IF( AAPP.LT.( BIG / AAQQ ) ) THEN
00428                                  AAPQ = ( SDOT( M, A( 1, p ), 1, A( 1,
00429      $                                  q ), 1 )*D( p )*D( q ) / AAQQ )
00430      $                                  / AAPP
00431                               ELSE
00432                                  CALL SCOPY( M, A( 1, p ), 1, WORK, 1 )
00433                                  CALL SLASCL( 'G', 0, 0, AAPP, D( p ),
00434      $                                        M, 1, WORK, LDA, IERR )
00435                                  AAPQ = SDOT( M, WORK, 1, A( 1, q ),
00436      $                                  1 )*D( q ) / AAQQ
00437                               END IF
00438                            ELSE
00439                               IF( AAPP.GE.AAQQ ) THEN
00440                                  ROTOK = AAPP.LE.( AAQQ / SMALL )
00441                               ELSE
00442                                  ROTOK = AAQQ.LE.( AAPP / SMALL )
00443                               END IF
00444                               IF( AAPP.GT.( SMALL / AAQQ ) ) THEN
00445                                  AAPQ = ( SDOT( M, A( 1, p ), 1, A( 1,
00446      $                                  q ), 1 )*D( p )*D( q ) / AAQQ )
00447      $                                  / AAPP
00448                               ELSE
00449                                  CALL SCOPY( M, A( 1, q ), 1, WORK, 1 )
00450                                  CALL SLASCL( 'G', 0, 0, AAQQ, D( q ),
00451      $                                        M, 1, WORK, LDA, IERR )
00452                                  AAPQ = SDOT( M, WORK, 1, A( 1, p ),
00453      $                                  1 )*D( p ) / AAPP
00454                               END IF
00455                            END IF
00456 
00457                            MXAAPQ = AMAX1( MXAAPQ, ABS( AAPQ ) )
00458 
00459 *        TO rotate or NOT to rotate, THAT is the question ...
00460 *
00461                            IF( ABS( AAPQ ).GT.TOL ) THEN
00462                               NOTROT = 0
00463 *           ROTATED  = ROTATED + 1
00464                               PSKIPPED = 0
00465                               ISWROT = ISWROT + 1
00466 *
00467                               IF( ROTOK ) THEN
00468 *
00469                                  AQOAP = AAQQ / AAPP
00470                                  APOAQ = AAPP / AAQQ
00471                                  THETA = -HALF*ABS( AQOAP-APOAQ ) / AAPQ
00472                                  IF( AAQQ.GT.AAPP0 )THETA = -THETA
00473 
00474                                  IF( ABS( THETA ).GT.BIGTHETA ) THEN
00475                                     T = HALF / THETA
00476                                     FASTR( 3 ) = T*D( p ) / D( q )
00477                                     FASTR( 4 ) = -T*D( q ) / D( p )
00478                                     CALL SROTM( M, A( 1, p ), 1,
00479      $                                          A( 1, q ), 1, FASTR )
00480                                     IF( RSVEC )CALL SROTM( MVL,
00481      $                                              V( 1, p ), 1,
00482      $                                              V( 1, q ), 1,
00483      $                                              FASTR )
00484                                     SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
00485      $                                         ONE+T*APOAQ*AAPQ ) )
00486                                     AAPP = AAPP*SQRT( AMAX1( ZERO,
00487      $                                     ONE-T*AQOAP*AAPQ ) )
00488                                     MXSINJ = AMAX1( MXSINJ, ABS( T ) )
00489                                  ELSE
00490 *
00491 *                 .. choose correct signum for THETA and rotate
00492 *
00493                                     THSIGN = -SIGN( ONE, AAPQ )
00494                                     IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN
00495                                     T = ONE / ( THETA+THSIGN*
00496      $                                  SQRT( ONE+THETA*THETA ) )
00497                                     CS = SQRT( ONE / ( ONE+T*T ) )
00498                                     SN = T*CS
00499                                     MXSINJ = AMAX1( MXSINJ, ABS( SN ) )
00500                                     SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
00501      $                                         ONE+T*APOAQ*AAPQ ) )
00502                                     AAPP = AAPP*SQRT( AMAX1( ZERO, 
00503      $                                         ONE-T*AQOAP*AAPQ ) )
00504 
00505                                     APOAQ = D( p ) / D( q )
00506                                     AQOAP = D( q ) / D( p )
00507                                     IF( D( p ).GE.ONE ) THEN
00508 *
00509                                        IF( D( q ).GE.ONE ) THEN
00510                                           FASTR( 3 ) = T*APOAQ
00511                                           FASTR( 4 ) = -T*AQOAP
00512                                           D( p ) = D( p )*CS
00513                                           D( q ) = D( q )*CS
00514                                           CALL SROTM( M, A( 1, p ), 1,
00515      $                                                A( 1, q ), 1,
00516      $                                                FASTR )
00517                                           IF( RSVEC )CALL SROTM( MVL,
00518      $                                        V( 1, p ), 1, V( 1, q ),
00519      $                                        1, FASTR )
00520                                        ELSE
00521                                           CALL SAXPY( M, -T*AQOAP,
00522      $                                                A( 1, q ), 1,
00523      $                                                A( 1, p ), 1 )
00524                                           CALL SAXPY( M, CS*SN*APOAQ,
00525      $                                                A( 1, p ), 1,
00526      $                                                A( 1, q ), 1 )
00527                                           IF( RSVEC ) THEN
00528                                              CALL SAXPY( MVL, -T*AQOAP,
00529      $                                                   V( 1, q ), 1,
00530      $                                                   V( 1, p ), 1 )
00531                                              CALL SAXPY( MVL,
00532      $                                                   CS*SN*APOAQ,
00533      $                                                   V( 1, p ), 1,
00534      $                                                   V( 1, q ), 1 )
00535                                           END IF
00536                                           D( p ) = D( p )*CS
00537                                           D( q ) = D( q ) / CS
00538                                        END IF
00539                                     ELSE
00540                                        IF( D( q ).GE.ONE ) THEN
00541                                           CALL SAXPY( M, T*APOAQ,
00542      $                                                A( 1, p ), 1,
00543      $                                                A( 1, q ), 1 )
00544                                           CALL SAXPY( M, -CS*SN*AQOAP,
00545      $                                                A( 1, q ), 1,
00546      $                                                A( 1, p ), 1 )
00547                                           IF( RSVEC ) THEN
00548                                              CALL SAXPY( MVL, T*APOAQ,
00549      $                                                   V( 1, p ), 1,
00550      $                                                   V( 1, q ), 1 )
00551                                              CALL SAXPY( MVL,
00552      $                                                   -CS*SN*AQOAP,
00553      $                                                   V( 1, q ), 1,
00554      $                                                   V( 1, p ), 1 )
00555                                           END IF
00556                                           D( p ) = D( p ) / CS
00557                                           D( q ) = D( q )*CS
00558                                        ELSE
00559                                           IF( D( p ).GE.D( q ) ) THEN
00560                                              CALL SAXPY( M, -T*AQOAP,
00561      $                                                   A( 1, q ), 1,
00562      $                                                   A( 1, p ), 1 )
00563                                              CALL SAXPY( M, CS*SN*APOAQ,
00564      $                                                   A( 1, p ), 1,
00565      $                                                   A( 1, q ), 1 )
00566                                              D( p ) = D( p )*CS
00567                                              D( q ) = D( q ) / CS
00568                                              IF( RSVEC ) THEN
00569                                                 CALL SAXPY( MVL,
00570      $                                               -T*AQOAP,
00571      $                                               V( 1, q ), 1,
00572      $                                               V( 1, p ), 1 )
00573                                                 CALL SAXPY( MVL,
00574      $                                               CS*SN*APOAQ,
00575      $                                               V( 1, p ), 1,
00576      $                                               V( 1, q ), 1 )
00577                                              END IF
00578                                           ELSE
00579                                              CALL SAXPY( M, T*APOAQ,
00580      $                                                   A( 1, p ), 1,
00581      $                                                   A( 1, q ), 1 )
00582                                              CALL SAXPY( M,
00583      $                                                   -CS*SN*AQOAP,
00584      $                                                   A( 1, q ), 1,
00585      $                                                   A( 1, p ), 1 )
00586                                              D( p ) = D( p ) / CS
00587                                              D( q ) = D( q )*CS
00588                                              IF( RSVEC ) THEN
00589                                                 CALL SAXPY( MVL,
00590      $                                               T*APOAQ, V( 1, p ),
00591      $                                               1, V( 1, q ), 1 )
00592                                                 CALL SAXPY( MVL,
00593      $                                               -CS*SN*AQOAP,
00594      $                                               V( 1, q ), 1,
00595      $                                               V( 1, p ), 1 )
00596                                              END IF
00597                                           END IF
00598                                        END IF
00599                                     END IF
00600                                  END IF
00601 
00602                               ELSE
00603                                  IF( AAPP.GT.AAQQ ) THEN
00604                                     CALL SCOPY( M, A( 1, p ), 1, WORK,
00605      $                                          1 )
00606                                     CALL SLASCL( 'G', 0, 0, AAPP, ONE,
00607      $                                           M, 1, WORK, LDA, IERR )
00608                                     CALL SLASCL( 'G', 0, 0, AAQQ, ONE,
00609      $                                           M, 1, A( 1, q ), LDA,
00610      $                                           IERR )
00611                                     TEMP1 = -AAPQ*D( p ) / D( q )
00612                                     CALL SAXPY( M, TEMP1, WORK, 1,
00613      $                                          A( 1, q ), 1 )
00614                                     CALL SLASCL( 'G', 0, 0, ONE, AAQQ,
00615      $                                           M, 1, A( 1, q ), LDA,
00616      $                                           IERR )
00617                                     SVA( q ) = AAQQ*SQRT( AMAX1( ZERO,
00618      $                                         ONE-AAPQ*AAPQ ) )
00619                                     MXSINJ = AMAX1( MXSINJ, SFMIN )
00620                                  ELSE
00621                                     CALL SCOPY( M, A( 1, q ), 1, WORK,
00622      $                                          1 )
00623                                     CALL SLASCL( 'G', 0, 0, AAQQ, ONE,
00624      $                                           M, 1, WORK, LDA, IERR )
00625                                     CALL SLASCL( 'G', 0, 0, AAPP, ONE,
00626      $                                           M, 1, A( 1, p ), LDA,
00627      $                                           IERR )
00628                                     TEMP1 = -AAPQ*D( q ) / D( p )
00629                                     CALL SAXPY( M, TEMP1, WORK, 1,
00630      $                                          A( 1, p ), 1 )
00631                                     CALL SLASCL( 'G', 0, 0, ONE, AAPP,
00632      $                                           M, 1, A( 1, p ), LDA,
00633      $                                           IERR )
00634                                     SVA( p ) = AAPP*SQRT( AMAX1( ZERO,
00635      $                                         ONE-AAPQ*AAPQ ) )
00636                                     MXSINJ = AMAX1( MXSINJ, SFMIN )
00637                                  END IF
00638                               END IF
00639 *           END IF ROTOK THEN ... ELSE
00640 *
00641 *           In the case of cancellation in updating SVA(q)
00642 *           .. recompute SVA(q)
00643                               IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS )
00644      $                            THEN
00645                                  IF( ( AAQQ.LT.ROOTBIG ) .AND.
00646      $                               ( AAQQ.GT.ROOTSFMIN ) ) THEN
00647                                     SVA( q ) = SNRM2( M, A( 1, q ), 1 )*
00648      $                                         D( q )
00649                                  ELSE
00650                                     T = ZERO
00651                                     AAQQ = ONE
00652                                     CALL SLASSQ( M, A( 1, q ), 1, T,
00653      $                                           AAQQ )
00654                                     SVA( q ) = T*SQRT( AAQQ )*D( q )
00655                                  END IF
00656                               END IF
00657                               IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN
00658                                  IF( ( AAPP.LT.ROOTBIG ) .AND.
00659      $                               ( AAPP.GT.ROOTSFMIN ) ) THEN
00660                                     AAPP = SNRM2( M, A( 1, p ), 1 )*
00661      $                                     D( p )
00662                                  ELSE
00663                                     T = ZERO
00664                                     AAPP = ONE
00665                                     CALL SLASSQ( M, A( 1, p ), 1, T,
00666      $                                           AAPP )
00667                                     AAPP = T*SQRT( AAPP )*D( p )
00668                                  END IF
00669                                  SVA( p ) = AAPP
00670                               END IF
00671 *              end of OK rotation
00672                            ELSE
00673                               NOTROT = NOTROT + 1
00674 *           SKIPPED  = SKIPPED  + 1
00675                               PSKIPPED = PSKIPPED + 1
00676                               IJBLSK = IJBLSK + 1
00677                            END IF
00678                         ELSE
00679                            NOTROT = NOTROT + 1
00680                            PSKIPPED = PSKIPPED + 1
00681                            IJBLSK = IJBLSK + 1
00682                         END IF
00683 
00684 *      IF ( NOTROT .GE. EMPTSW )  GO TO 2011
00685                         IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) )
00686      $                      THEN
00687                            SVA( p ) = AAPP
00688                            NOTROT = 0
00689                            GO TO 2011
00690                         END IF
00691                         IF( ( i.LE.SWBAND ) .AND.
00692      $                      ( PSKIPPED.GT.ROWSKIP ) ) THEN
00693                            AAPP = -AAPP
00694                            NOTROT = 0
00695                            GO TO 2203
00696                         END IF
00697 
00698 *
00699  2200                CONTINUE
00700 *        end of the q-loop
00701  2203                CONTINUE
00702 
00703                      SVA( p ) = AAPP
00704 *
00705                   ELSE
00706                      IF( AAPP.EQ.ZERO )NOTROT = NOTROT +
00707      $                   MIN0( jgl+KBL-1, N ) - jgl + 1
00708                      IF( AAPP.LT.ZERO )NOTROT = 0
00709 ***      IF ( NOTROT .GE. EMPTSW )  GO TO 2011
00710                   END IF
00711 
00712  2100          CONTINUE
00713 *     end of the p-loop
00714  2010       CONTINUE
00715 *     end of the jbc-loop
00716  2011       CONTINUE
00717 *2011 bailed out of the jbc-loop
00718             DO 2012 p = igl, MIN0( igl+KBL-1, N )
00719                SVA( p ) = ABS( SVA( p ) )
00720  2012       CONTINUE
00721 ***   IF ( NOTROT .GE. EMPTSW ) GO TO 1994
00722  2000    CONTINUE
00723 *2000 :: end of the ibr-loop
00724 *
00725 *     .. update SVA(N)
00726          IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) )
00727      $       THEN
00728             SVA( N ) = SNRM2( M, A( 1, N ), 1 )*D( N )
00729          ELSE
00730             T = ZERO
00731             AAPP = ONE
00732             CALL SLASSQ( M, A( 1, N ), 1, T, AAPP )
00733             SVA( N ) = T*SQRT( AAPP )*D( N )
00734          END IF
00735 *
00736 *     Additional steering devices
00737 *
00738          IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR.
00739      $       ( ISWROT.LE.N ) ) )SWBAND = i
00740 
00741          IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.FLOAT( N )*TOL ) .AND.
00742      $       ( FLOAT( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN
00743             GO TO 1994
00744          END IF
00745 
00746 *
00747          IF( NOTROT.GE.EMPTSW )GO TO 1994
00748 
00749  1993 CONTINUE
00750 *     end i=1:NSWEEP loop
00751 * #:) Reaching this point means that the procedure has completed the given
00752 *     number of sweeps.
00753       INFO = NSWEEP - 1
00754       GO TO 1995
00755  1994 CONTINUE
00756 * #:) Reaching this point means that during the i-th sweep all pivots were
00757 *     below the given threshold, causing early exit.
00758 
00759       INFO = 0
00760 * #:) INFO = 0 confirms successful iterations.
00761  1995 CONTINUE
00762 *
00763 *     Sort the vector D
00764 *
00765       DO 5991 p = 1, N - 1
00766          q = ISAMAX( N-p+1, SVA( p ), 1 ) + p - 1
00767          IF( p.NE.q ) THEN
00768             TEMP1 = SVA( p )
00769             SVA( p ) = SVA( q )
00770             SVA( q ) = TEMP1
00771             TEMP1 = D( p )
00772             D( p ) = D( q )
00773             D( q ) = TEMP1
00774             CALL SSWAP( M, A( 1, p ), 1, A( 1, q ), 1 )
00775             IF( RSVEC )CALL SSWAP( MVL, V( 1, p ), 1, V( 1, q ), 1 )
00776          END IF
00777  5991 CONTINUE
00778 *
00779       RETURN
00780 *     ..
00781 *     .. END OF SGSVJ1
00782 *     ..
00783       END
 All Files Functions