 |
LAPACK
3.4.0
LAPACK: Linear Algebra PACKage
|
|
LAPACK
3.4.0
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b SGSVJ0 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download SGSVJ0 + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgsvj0.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgsvj0.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgsvj0.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, 00022 * SFMIN, TOL, NSWEEP, WORK, LWORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * INTEGER INFO, LDA, LDV, LWORK, M, MV, N, NSWEEP 00026 * REAL EPS, SFMIN, TOL 00027 * CHARACTER*1 JOBV 00028 * .. 00029 * .. Array Arguments .. 00030 * REAL A( LDA, * ), SVA( N ), D( N ), V( LDV, * ), 00031 * $ WORK( LWORK ) 00032 * .. 00033 * 00034 * 00035 *> \par Purpose: 00036 * ============= 00037 *> 00038 *> \verbatim 00039 *> 00040 *> SGSVJ0 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 does not check convergence (stopping criterion). Few tuning 00043 *> parameters (marked by [TP]) are available for the implementer. 00044 *> \endverbatim 00045 * 00046 * Arguments: 00047 * ========== 00048 * 00049 *> \param[in] JOBV 00050 *> \verbatim 00051 *> JOBV is CHARACTER*1 00052 *> Specifies whether the output from this procedure is used 00053 *> to compute the matrix V: 00054 *> = 'V': the product of the Jacobi rotations is accumulated 00055 *> by postmulyiplying the N-by-N array V. 00056 *> (See the description of V.) 00057 *> = 'A': the product of the Jacobi rotations is accumulated 00058 *> by postmulyiplying the MV-by-N array V. 00059 *> (See the descriptions of MV and V.) 00060 *> = 'N': the Jacobi rotations are not accumulated. 00061 *> \endverbatim 00062 *> 00063 *> \param[in] M 00064 *> \verbatim 00065 *> M is INTEGER 00066 *> The number of rows of the input matrix A. M >= 0. 00067 *> \endverbatim 00068 *> 00069 *> \param[in] N 00070 *> \verbatim 00071 *> N is INTEGER 00072 *> The number of columns of the input matrix A. 00073 *> M >= N >= 0. 00074 *> \endverbatim 00075 *> 00076 *> \param[in,out] A 00077 *> \verbatim 00078 *> A is REAL array, dimension (LDA,N) 00079 *> On entry, M-by-N matrix A, such that A*diag(D) represents 00080 *> the input matrix. 00081 *> On exit, 00082 *> A_onexit * D_onexit represents the input matrix A*diag(D) 00083 *> post-multiplied by a sequence of Jacobi rotations, where the 00084 *> rotation threshold and the total number of sweeps are given in 00085 *> TOL and NSWEEP, respectively. 00086 *> (See the descriptions of D, TOL and NSWEEP.) 00087 *> \endverbatim 00088 *> 00089 *> \param[in] LDA 00090 *> \verbatim 00091 *> LDA is INTEGER 00092 *> The leading dimension of the array A. LDA >= max(1,M). 00093 *> \endverbatim 00094 *> 00095 *> \param[in,out] D 00096 *> \verbatim 00097 *> D is REAL array, dimension (N) 00098 *> The array D accumulates the scaling factors from the fast scaled 00099 *> Jacobi rotations. 00100 *> On entry, A*diag(D) represents the input matrix. 00101 *> On exit, A_onexit*diag(D_onexit) represents the input matrix 00102 *> post-multiplied by a sequence of Jacobi rotations, where the 00103 *> rotation threshold and the total number of sweeps are given in 00104 *> TOL and NSWEEP, respectively. 00105 *> (See the descriptions of A, TOL and NSWEEP.) 00106 *> \endverbatim 00107 *> 00108 *> \param[in,out] SVA 00109 *> \verbatim 00110 *> SVA is REAL array, dimension (N) 00111 *> On entry, SVA contains the Euclidean norms of the columns of 00112 *> the matrix A*diag(D). 00113 *> On exit, SVA contains the Euclidean norms of the columns of 00114 *> the matrix onexit*diag(D_onexit). 00115 *> \endverbatim 00116 *> 00117 *> \param[in] MV 00118 *> \verbatim 00119 *> MV is INTEGER 00120 *> If JOBV .EQ. 'A', then MV rows of V are post-multipled by a 00121 *> sequence of Jacobi rotations. 00122 *> If JOBV = 'N', then MV is not referenced. 00123 *> \endverbatim 00124 *> 00125 *> \param[in,out] V 00126 *> \verbatim 00127 *> V is REAL array, dimension (LDV,N) 00128 *> If JOBV .EQ. 'V' then N rows of V are post-multipled by a 00129 *> sequence of Jacobi rotations. 00130 *> If JOBV .EQ. 'A' then MV rows of V are post-multipled by a 00131 *> sequence of Jacobi rotations. 00132 *> If JOBV = 'N', then V is not referenced. 00133 *> \endverbatim 00134 *> 00135 *> \param[in] LDV 00136 *> \verbatim 00137 *> LDV is INTEGER 00138 *> The leading dimension of the array V, LDV >= 1. 00139 *> If JOBV = 'V', LDV .GE. N. 00140 *> If JOBV = 'A', LDV .GE. MV. 00141 *> \endverbatim 00142 *> 00143 *> \param[in] EPS 00144 *> \verbatim 00145 *> EPS is INTEGER 00146 *> EPS = SLAMCH('Epsilon') 00147 *> \endverbatim 00148 *> 00149 *> \param[in] SFMIN 00150 *> \verbatim 00151 *> SFMIN is INTEGER 00152 *> SFMIN = SLAMCH('Safe Minimum') 00153 *> \endverbatim 00154 *> 00155 *> \param[in] TOL 00156 *> \verbatim 00157 *> TOL is REAL 00158 *> TOL is the threshold for Jacobi rotations. For a pair 00159 *> A(:,p), A(:,q) of pivot columns, the Jacobi rotation is 00160 *> applied only if ABS(COS(angle(A(:,p),A(:,q)))) .GT. TOL. 00161 *> \endverbatim 00162 *> 00163 *> \param[in] NSWEEP 00164 *> \verbatim 00165 *> NSWEEP is INTEGER 00166 *> NSWEEP is the number of sweeps of Jacobi rotations to be 00167 *> performed. 00168 *> \endverbatim 00169 *> 00170 *> \param[out] WORK 00171 *> \verbatim 00172 *> WORK is REAL array, dimension LWORK. 00173 *> \endverbatim 00174 *> 00175 *> \param[in] LWORK 00176 *> \verbatim 00177 *> LWORK is INTEGER 00178 *> LWORK is the dimension of WORK. LWORK .GE. M. 00179 *> \endverbatim 00180 *> 00181 *> \param[out] INFO 00182 *> \verbatim 00183 *> INFO is INTEGER 00184 *> = 0 : successful exit. 00185 *> < 0 : if INFO = -i, then the i-th argument had an illegal value 00186 *> \endverbatim 00187 * 00188 * Authors: 00189 * ======== 00190 * 00191 *> \author Univ. of Tennessee 00192 *> \author Univ. of California Berkeley 00193 *> \author Univ. of Colorado Denver 00194 *> \author NAG Ltd. 00195 * 00196 *> \date November 2011 00197 * 00198 *> \ingroup realOTHERcomputational 00199 * 00200 *> \par Further Details: 00201 * ===================== 00202 *> 00203 *> SGSVJ0 is used just to enable SGESVJ to call a simplified version of 00204 *> itself to work on a submatrix of the original matrix. 00205 *> 00206 *> \par Contributors: 00207 * ================== 00208 *> 00209 *> Zlatko Drmac (Zagreb, Croatia) and Kresimir Veselic (Hagen, Germany) 00210 *> 00211 *> \par Bugs, Examples and Comments: 00212 * ================================= 00213 *> 00214 *> Please report all bugs and send interesting test examples and comments to 00215 *> drmac@math.hr. Thank you. 00216 * 00217 * ===================================================================== 00218 SUBROUTINE SGSVJ0( JOBV, M, N, A, LDA, D, SVA, MV, V, LDV, EPS, 00219 $ SFMIN, TOL, NSWEEP, WORK, LWORK, INFO ) 00220 * 00221 * -- LAPACK computational routine (version 3.4.0) -- 00222 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00223 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00224 * November 2011 00225 * 00226 * .. Scalar Arguments .. 00227 INTEGER INFO, LDA, LDV, LWORK, M, MV, N, NSWEEP 00228 REAL EPS, SFMIN, TOL 00229 CHARACTER*1 JOBV 00230 * .. 00231 * .. Array Arguments .. 00232 REAL A( LDA, * ), SVA( N ), D( N ), V( LDV, * ), 00233 $ WORK( LWORK ) 00234 * .. 00235 * 00236 * ===================================================================== 00237 * 00238 * .. Local Parameters .. 00239 REAL ZERO, HALF, ONE, TWO 00240 PARAMETER ( ZERO = 0.0E0, HALF = 0.5E0, ONE = 1.0E0, 00241 $ TWO = 2.0E0 ) 00242 * .. 00243 * .. Local Scalars .. 00244 REAL AAPP, AAPP0, AAPQ, AAQQ, APOAQ, AQOAP, BIG, 00245 $ BIGTHETA, CS, MXAAPQ, MXSINJ, ROOTBIG, ROOTEPS, 00246 $ ROOTSFMIN, ROOTTOL, SMALL, SN, T, TEMP1, THETA, 00247 $ THSIGN 00248 INTEGER BLSKIP, EMPTSW, i, ibr, IERR, igl, IJBLSK, ir1, 00249 $ ISWROT, jbc, jgl, KBL, LKAHEAD, MVL, NBL, 00250 $ NOTROT, p, PSKIPPED, q, ROWSKIP, SWBAND 00251 LOGICAL APPLV, ROTOK, RSVEC 00252 * .. 00253 * .. Local Arrays .. 00254 REAL FASTR( 5 ) 00255 * .. 00256 * .. Intrinsic Functions .. 00257 INTRINSIC ABS, AMAX1, FLOAT, MIN0, SIGN, SQRT 00258 * .. 00259 * .. External Functions .. 00260 REAL SDOT, SNRM2 00261 INTEGER ISAMAX 00262 LOGICAL LSAME 00263 EXTERNAL ISAMAX, LSAME, SDOT, SNRM2 00264 * .. 00265 * .. External Subroutines .. 00266 EXTERNAL SAXPY, SCOPY, SLASCL, SLASSQ, SROTM, SSWAP 00267 * .. 00268 * .. Executable Statements .. 00269 * 00270 * Test the input parameters. 00271 * 00272 APPLV = LSAME( JOBV, 'A' ) 00273 RSVEC = LSAME( JOBV, 'V' ) 00274 IF( .NOT.( RSVEC .OR. APPLV .OR. LSAME( JOBV, 'N' ) ) ) THEN 00275 INFO = -1 00276 ELSE IF( M.LT.0 ) THEN 00277 INFO = -2 00278 ELSE IF( ( N.LT.0 ) .OR. ( N.GT.M ) ) THEN 00279 INFO = -3 00280 ELSE IF( LDA.LT.M ) THEN 00281 INFO = -5 00282 ELSE IF( ( RSVEC.OR.APPLV ) .AND. ( MV.LT.0 ) ) THEN 00283 INFO = -8 00284 ELSE IF( ( RSVEC.AND.( LDV.LT.N ) ).OR. 00285 $ ( APPLV.AND.( LDV.LT.MV ) ) ) THEN 00286 INFO = -10 00287 ELSE IF( TOL.LE.EPS ) THEN 00288 INFO = -13 00289 ELSE IF( NSWEEP.LT.0 ) THEN 00290 INFO = -14 00291 ELSE IF( LWORK.LT.M ) THEN 00292 INFO = -16 00293 ELSE 00294 INFO = 0 00295 END IF 00296 * 00297 * #:( 00298 IF( INFO.NE.0 ) THEN 00299 CALL XERBLA( 'SGSVJ0', -INFO ) 00300 RETURN 00301 END IF 00302 * 00303 IF( RSVEC ) THEN 00304 MVL = N 00305 ELSE IF( APPLV ) THEN 00306 MVL = MV 00307 END IF 00308 RSVEC = RSVEC .OR. APPLV 00309 00310 ROOTEPS = SQRT( EPS ) 00311 ROOTSFMIN = SQRT( SFMIN ) 00312 SMALL = SFMIN / EPS 00313 BIG = ONE / SFMIN 00314 ROOTBIG = ONE / ROOTSFMIN 00315 BIGTHETA = ONE / ROOTEPS 00316 ROOTTOL = SQRT( TOL ) 00317 * 00318 * .. Row-cyclic Jacobi SVD algorithm with column pivoting .. 00319 * 00320 EMPTSW = ( N*( N-1 ) ) / 2 00321 NOTROT = 0 00322 FASTR( 1 ) = ZERO 00323 * 00324 * .. Row-cyclic pivot strategy with de Rijk's pivoting .. 00325 * 00326 00327 SWBAND = 0 00328 *[TP] SWBAND is a tuning parameter. It is meaningful and effective 00329 * if SGESVJ is used as a computational routine in the preconditioned 00330 * Jacobi SVD algorithm SGESVJ. For sweeps i=1:SWBAND the procedure 00331 * ...... 00332 00333 KBL = MIN0( 8, N ) 00334 *[TP] KBL is a tuning parameter that defines the tile size in the 00335 * tiling of the p-q loops of pivot pairs. In general, an optimal 00336 * value of KBL depends on the matrix dimensions and on the 00337 * parameters of the computer's memory. 00338 * 00339 NBL = N / KBL 00340 IF( ( NBL*KBL ).NE.N )NBL = NBL + 1 00341 00342 BLSKIP = ( KBL**2 ) + 1 00343 *[TP] BLKSKIP is a tuning parameter that depends on SWBAND and KBL. 00344 00345 ROWSKIP = MIN0( 5, KBL ) 00346 *[TP] ROWSKIP is a tuning parameter. 00347 00348 LKAHEAD = 1 00349 *[TP] LKAHEAD is a tuning parameter. 00350 SWBAND = 0 00351 PSKIPPED = 0 00352 * 00353 DO 1993 i = 1, NSWEEP 00354 * .. go go go ... 00355 * 00356 MXAAPQ = ZERO 00357 MXSINJ = ZERO 00358 ISWROT = 0 00359 * 00360 NOTROT = 0 00361 PSKIPPED = 0 00362 * 00363 DO 2000 ibr = 1, NBL 00364 00365 igl = ( ibr-1 )*KBL + 1 00366 * 00367 DO 1002 ir1 = 0, MIN0( LKAHEAD, NBL-ibr ) 00368 * 00369 igl = igl + ir1*KBL 00370 * 00371 DO 2001 p = igl, MIN0( igl+KBL-1, N-1 ) 00372 00373 * .. de Rijk's pivoting 00374 q = ISAMAX( N-p+1, SVA( p ), 1 ) + p - 1 00375 IF( p.NE.q ) THEN 00376 CALL SSWAP( M, A( 1, p ), 1, A( 1, q ), 1 ) 00377 IF( RSVEC )CALL SSWAP( MVL, V( 1, p ), 1, 00378 $ V( 1, q ), 1 ) 00379 TEMP1 = SVA( p ) 00380 SVA( p ) = SVA( q ) 00381 SVA( q ) = TEMP1 00382 TEMP1 = D( p ) 00383 D( p ) = D( q ) 00384 D( q ) = TEMP1 00385 END IF 00386 * 00387 IF( ir1.EQ.0 ) THEN 00388 * 00389 * Column norms are periodically updated by explicit 00390 * norm computation. 00391 * Caveat: 00392 * Some BLAS implementations compute SNRM2(M,A(1,p),1) 00393 * as SQRT(SDOT(M,A(1,p),1,A(1,p),1)), which may result in 00394 * overflow for ||A(:,p)||_2 > SQRT(overflow_threshold), and 00395 * undeflow for ||A(:,p)||_2 < SQRT(underflow_threshold). 00396 * Hence, SNRM2 cannot be trusted, not even in the case when 00397 * the true norm is far from the under(over)flow boundaries. 00398 * If properly implemented SNRM2 is available, the IF-THEN-ELSE 00399 * below should read "AAPP = SNRM2( M, A(1,p), 1 ) * D(p)". 00400 * 00401 IF( ( SVA( p ).LT.ROOTBIG ) .AND. 00402 $ ( SVA( p ).GT.ROOTSFMIN ) ) THEN 00403 SVA( p ) = SNRM2( M, A( 1, p ), 1 )*D( p ) 00404 ELSE 00405 TEMP1 = ZERO 00406 AAPP = ONE 00407 CALL SLASSQ( M, A( 1, p ), 1, TEMP1, AAPP ) 00408 SVA( p ) = TEMP1*SQRT( AAPP )*D( p ) 00409 END IF 00410 AAPP = SVA( p ) 00411 ELSE 00412 AAPP = SVA( p ) 00413 END IF 00414 00415 * 00416 IF( AAPP.GT.ZERO ) THEN 00417 * 00418 PSKIPPED = 0 00419 * 00420 DO 2002 q = p + 1, MIN0( igl+KBL-1, N ) 00421 * 00422 AAQQ = SVA( q ) 00423 00424 IF( AAQQ.GT.ZERO ) THEN 00425 * 00426 AAPP0 = AAPP 00427 IF( AAQQ.GE.ONE ) THEN 00428 ROTOK = ( SMALL*AAPP ).LE.AAQQ 00429 IF( AAPP.LT.( BIG / AAQQ ) ) THEN 00430 AAPQ = ( SDOT( M, A( 1, p ), 1, A( 1, 00431 $ q ), 1 )*D( p )*D( q ) / AAQQ ) 00432 $ / AAPP 00433 ELSE 00434 CALL SCOPY( M, A( 1, p ), 1, WORK, 1 ) 00435 CALL SLASCL( 'G', 0, 0, AAPP, D( p ), 00436 $ M, 1, WORK, LDA, IERR ) 00437 AAPQ = SDOT( M, WORK, 1, A( 1, q ), 00438 $ 1 )*D( q ) / AAQQ 00439 END IF 00440 ELSE 00441 ROTOK = AAPP.LE.( AAQQ / SMALL ) 00442 IF( AAPP.GT.( SMALL / AAQQ ) ) THEN 00443 AAPQ = ( SDOT( M, A( 1, p ), 1, A( 1, 00444 $ q ), 1 )*D( p )*D( q ) / AAQQ ) 00445 $ / AAPP 00446 ELSE 00447 CALL SCOPY( M, A( 1, q ), 1, WORK, 1 ) 00448 CALL SLASCL( 'G', 0, 0, AAQQ, D( q ), 00449 $ M, 1, WORK, LDA, IERR ) 00450 AAPQ = SDOT( M, WORK, 1, A( 1, p ), 00451 $ 1 )*D( p ) / AAPP 00452 END IF 00453 END IF 00454 * 00455 MXAAPQ = AMAX1( MXAAPQ, ABS( AAPQ ) ) 00456 * 00457 * TO rotate or NOT to rotate, THAT is the question ... 00458 * 00459 IF( ABS( AAPQ ).GT.TOL ) THEN 00460 * 00461 * .. rotate 00462 * ROTATED = ROTATED + ONE 00463 * 00464 IF( ir1.EQ.0 ) THEN 00465 NOTROT = 0 00466 PSKIPPED = 0 00467 ISWROT = ISWROT + 1 00468 END IF 00469 * 00470 IF( ROTOK ) THEN 00471 * 00472 AQOAP = AAQQ / AAPP 00473 APOAQ = AAPP / AAQQ 00474 THETA = -HALF*ABS( AQOAP-APOAQ ) / AAPQ 00475 * 00476 IF( ABS( THETA ).GT.BIGTHETA ) THEN 00477 * 00478 T = HALF / THETA 00479 FASTR( 3 ) = T*D( p ) / D( q ) 00480 FASTR( 4 ) = -T*D( q ) / D( p ) 00481 CALL SROTM( M, A( 1, p ), 1, 00482 $ A( 1, q ), 1, FASTR ) 00483 IF( RSVEC )CALL SROTM( MVL, 00484 $ V( 1, p ), 1, 00485 $ V( 1, q ), 1, 00486 $ FASTR ) 00487 SVA( q ) = AAQQ*SQRT( AMAX1( ZERO, 00488 $ ONE+T*APOAQ*AAPQ ) ) 00489 AAPP = AAPP*SQRT( AMAX1( ZERO, 00490 $ ONE-T*AQOAP*AAPQ ) ) 00491 MXSINJ = AMAX1( MXSINJ, ABS( T ) ) 00492 * 00493 ELSE 00494 * 00495 * .. choose correct signum for THETA and rotate 00496 * 00497 THSIGN = -SIGN( ONE, AAPQ ) 00498 T = ONE / ( THETA+THSIGN* 00499 $ SQRT( ONE+THETA*THETA ) ) 00500 CS = SQRT( ONE / ( ONE+T*T ) ) 00501 SN = T*CS 00502 * 00503 MXSINJ = AMAX1( MXSINJ, ABS( SN ) ) 00504 SVA( q ) = AAQQ*SQRT( AMAX1( ZERO, 00505 $ ONE+T*APOAQ*AAPQ ) ) 00506 AAPP = AAPP*SQRT( AMAX1( ZERO, 00507 $ ONE-T*AQOAP*AAPQ ) ) 00508 * 00509 APOAQ = D( p ) / D( q ) 00510 AQOAP = D( q ) / D( p ) 00511 IF( D( p ).GE.ONE ) THEN 00512 IF( D( q ).GE.ONE ) THEN 00513 FASTR( 3 ) = T*APOAQ 00514 FASTR( 4 ) = -T*AQOAP 00515 D( p ) = D( p )*CS 00516 D( q ) = D( q )*CS 00517 CALL SROTM( M, A( 1, p ), 1, 00518 $ A( 1, q ), 1, 00519 $ FASTR ) 00520 IF( RSVEC )CALL SROTM( MVL, 00521 $ V( 1, p ), 1, V( 1, q ), 00522 $ 1, FASTR ) 00523 ELSE 00524 CALL SAXPY( M, -T*AQOAP, 00525 $ A( 1, q ), 1, 00526 $ A( 1, p ), 1 ) 00527 CALL SAXPY( M, CS*SN*APOAQ, 00528 $ A( 1, p ), 1, 00529 $ A( 1, q ), 1 ) 00530 D( p ) = D( p )*CS 00531 D( q ) = D( q ) / CS 00532 IF( RSVEC ) THEN 00533 CALL SAXPY( MVL, -T*AQOAP, 00534 $ V( 1, q ), 1, 00535 $ V( 1, p ), 1 ) 00536 CALL SAXPY( MVL, 00537 $ CS*SN*APOAQ, 00538 $ V( 1, p ), 1, 00539 $ V( 1, q ), 1 ) 00540 END IF 00541 END IF 00542 ELSE 00543 IF( D( q ).GE.ONE ) THEN 00544 CALL SAXPY( M, T*APOAQ, 00545 $ A( 1, p ), 1, 00546 $ A( 1, q ), 1 ) 00547 CALL SAXPY( M, -CS*SN*AQOAP, 00548 $ A( 1, q ), 1, 00549 $ A( 1, p ), 1 ) 00550 D( p ) = D( p ) / CS 00551 D( q ) = D( q )*CS 00552 IF( RSVEC ) THEN 00553 CALL SAXPY( MVL, T*APOAQ, 00554 $ V( 1, p ), 1, 00555 $ V( 1, q ), 1 ) 00556 CALL SAXPY( MVL, 00557 $ -CS*SN*AQOAP, 00558 $ V( 1, q ), 1, 00559 $ V( 1, p ), 1 ) 00560 END IF 00561 ELSE 00562 IF( D( p ).GE.D( q ) ) THEN 00563 CALL SAXPY( M, -T*AQOAP, 00564 $ A( 1, q ), 1, 00565 $ A( 1, p ), 1 ) 00566 CALL SAXPY( M, CS*SN*APOAQ, 00567 $ A( 1, p ), 1, 00568 $ A( 1, q ), 1 ) 00569 D( p ) = D( p )*CS 00570 D( q ) = D( q ) / CS 00571 IF( RSVEC ) THEN 00572 CALL SAXPY( MVL, 00573 $ -T*AQOAP, 00574 $ V( 1, q ), 1, 00575 $ V( 1, p ), 1 ) 00576 CALL SAXPY( MVL, 00577 $ CS*SN*APOAQ, 00578 $ V( 1, p ), 1, 00579 $ V( 1, q ), 1 ) 00580 END IF 00581 ELSE 00582 CALL SAXPY( M, T*APOAQ, 00583 $ A( 1, p ), 1, 00584 $ A( 1, q ), 1 ) 00585 CALL SAXPY( M, 00586 $ -CS*SN*AQOAP, 00587 $ A( 1, q ), 1, 00588 $ A( 1, p ), 1 ) 00589 D( p ) = D( p ) / CS 00590 D( q ) = D( q )*CS 00591 IF( RSVEC ) THEN 00592 CALL SAXPY( MVL, 00593 $ T*APOAQ, V( 1, p ), 00594 $ 1, V( 1, q ), 1 ) 00595 CALL SAXPY( MVL, 00596 $ -CS*SN*AQOAP, 00597 $ V( 1, q ), 1, 00598 $ V( 1, p ), 1 ) 00599 END IF 00600 END IF 00601 END IF 00602 END IF 00603 END IF 00604 * 00605 ELSE 00606 * .. have to use modified Gram-Schmidt like transformation 00607 CALL SCOPY( M, A( 1, p ), 1, WORK, 1 ) 00608 CALL SLASCL( 'G', 0, 0, AAPP, ONE, M, 00609 $ 1, WORK, LDA, IERR ) 00610 CALL SLASCL( 'G', 0, 0, AAQQ, ONE, M, 00611 $ 1, A( 1, q ), LDA, IERR ) 00612 TEMP1 = -AAPQ*D( p ) / D( q ) 00613 CALL SAXPY( M, TEMP1, WORK, 1, 00614 $ A( 1, q ), 1 ) 00615 CALL SLASCL( 'G', 0, 0, ONE, AAQQ, M, 00616 $ 1, A( 1, q ), LDA, IERR ) 00617 SVA( q ) = AAQQ*SQRT( AMAX1( ZERO, 00618 $ ONE-AAPQ*AAPQ ) ) 00619 MXSINJ = AMAX1( MXSINJ, SFMIN ) 00620 END IF 00621 * END IF ROTOK THEN ... ELSE 00622 * 00623 * In the case of cancellation in updating SVA(q), SVA(p) 00624 * recompute SVA(q), SVA(p). 00625 IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS ) 00626 $ THEN 00627 IF( ( AAQQ.LT.ROOTBIG ) .AND. 00628 $ ( AAQQ.GT.ROOTSFMIN ) ) THEN 00629 SVA( q ) = SNRM2( M, A( 1, q ), 1 )* 00630 $ D( q ) 00631 ELSE 00632 T = ZERO 00633 AAQQ = ONE 00634 CALL SLASSQ( M, A( 1, q ), 1, T, 00635 $ AAQQ ) 00636 SVA( q ) = T*SQRT( AAQQ )*D( q ) 00637 END IF 00638 END IF 00639 IF( ( AAPP / AAPP0 ).LE.ROOTEPS ) THEN 00640 IF( ( AAPP.LT.ROOTBIG ) .AND. 00641 $ ( AAPP.GT.ROOTSFMIN ) ) THEN 00642 AAPP = SNRM2( M, A( 1, p ), 1 )* 00643 $ D( p ) 00644 ELSE 00645 T = ZERO 00646 AAPP = ONE 00647 CALL SLASSQ( M, A( 1, p ), 1, T, 00648 $ AAPP ) 00649 AAPP = T*SQRT( AAPP )*D( p ) 00650 END IF 00651 SVA( p ) = AAPP 00652 END IF 00653 * 00654 ELSE 00655 * A(:,p) and A(:,q) already numerically orthogonal 00656 IF( ir1.EQ.0 )NOTROT = NOTROT + 1 00657 PSKIPPED = PSKIPPED + 1 00658 END IF 00659 ELSE 00660 * A(:,q) is zero column 00661 IF( ir1.EQ.0 )NOTROT = NOTROT + 1 00662 PSKIPPED = PSKIPPED + 1 00663 END IF 00664 * 00665 IF( ( i.LE.SWBAND ) .AND. 00666 $ ( PSKIPPED.GT.ROWSKIP ) ) THEN 00667 IF( ir1.EQ.0 )AAPP = -AAPP 00668 NOTROT = 0 00669 GO TO 2103 00670 END IF 00671 * 00672 2002 CONTINUE 00673 * END q-LOOP 00674 * 00675 2103 CONTINUE 00676 * bailed out of q-loop 00677 00678 SVA( p ) = AAPP 00679 00680 ELSE 00681 SVA( p ) = AAPP 00682 IF( ( ir1.EQ.0 ) .AND. ( AAPP.EQ.ZERO ) ) 00683 $ NOTROT = NOTROT + MIN0( igl+KBL-1, N ) - p 00684 END IF 00685 * 00686 2001 CONTINUE 00687 * end of the p-loop 00688 * end of doing the block ( ibr, ibr ) 00689 1002 CONTINUE 00690 * end of ir1-loop 00691 * 00692 *........................................................ 00693 * ... go to the off diagonal blocks 00694 * 00695 igl = ( ibr-1 )*KBL + 1 00696 * 00697 DO 2010 jbc = ibr + 1, NBL 00698 * 00699 jgl = ( jbc-1 )*KBL + 1 00700 * 00701 * doing the block at ( ibr, jbc ) 00702 * 00703 IJBLSK = 0 00704 DO 2100 p = igl, MIN0( igl+KBL-1, N ) 00705 * 00706 AAPP = SVA( p ) 00707 * 00708 IF( AAPP.GT.ZERO ) THEN 00709 * 00710 PSKIPPED = 0 00711 * 00712 DO 2200 q = jgl, MIN0( jgl+KBL-1, N ) 00713 * 00714 AAQQ = SVA( q ) 00715 * 00716 IF( AAQQ.GT.ZERO ) THEN 00717 AAPP0 = AAPP 00718 * 00719 * .. M x 2 Jacobi SVD .. 00720 * 00721 * .. Safe Gram matrix computation .. 00722 * 00723 IF( AAQQ.GE.ONE ) THEN 00724 IF( AAPP.GE.AAQQ ) THEN 00725 ROTOK = ( SMALL*AAPP ).LE.AAQQ 00726 ELSE 00727 ROTOK = ( SMALL*AAQQ ).LE.AAPP 00728 END IF 00729 IF( AAPP.LT.( BIG / AAQQ ) ) THEN 00730 AAPQ = ( SDOT( M, A( 1, p ), 1, A( 1, 00731 $ q ), 1 )*D( p )*D( q ) / AAQQ ) 00732 $ / AAPP 00733 ELSE 00734 CALL SCOPY( M, A( 1, p ), 1, WORK, 1 ) 00735 CALL SLASCL( 'G', 0, 0, AAPP, D( p ), 00736 $ M, 1, WORK, LDA, IERR ) 00737 AAPQ = SDOT( M, WORK, 1, A( 1, q ), 00738 $ 1 )*D( q ) / AAQQ 00739 END IF 00740 ELSE 00741 IF( AAPP.GE.AAQQ ) THEN 00742 ROTOK = AAPP.LE.( AAQQ / SMALL ) 00743 ELSE 00744 ROTOK = AAQQ.LE.( AAPP / SMALL ) 00745 END IF 00746 IF( AAPP.GT.( SMALL / AAQQ ) ) THEN 00747 AAPQ = ( SDOT( M, A( 1, p ), 1, A( 1, 00748 $ q ), 1 )*D( p )*D( q ) / AAQQ ) 00749 $ / AAPP 00750 ELSE 00751 CALL SCOPY( M, A( 1, q ), 1, WORK, 1 ) 00752 CALL SLASCL( 'G', 0, 0, AAQQ, D( q ), 00753 $ M, 1, WORK, LDA, IERR ) 00754 AAPQ = SDOT( M, WORK, 1, A( 1, p ), 00755 $ 1 )*D( p ) / AAPP 00756 END IF 00757 END IF 00758 * 00759 MXAAPQ = AMAX1( MXAAPQ, ABS( AAPQ ) ) 00760 * 00761 * TO rotate or NOT to rotate, THAT is the question ... 00762 * 00763 IF( ABS( AAPQ ).GT.TOL ) THEN 00764 NOTROT = 0 00765 * ROTATED = ROTATED + 1 00766 PSKIPPED = 0 00767 ISWROT = ISWROT + 1 00768 * 00769 IF( ROTOK ) THEN 00770 * 00771 AQOAP = AAQQ / AAPP 00772 APOAQ = AAPP / AAQQ 00773 THETA = -HALF*ABS( AQOAP-APOAQ ) / AAPQ 00774 IF( AAQQ.GT.AAPP0 )THETA = -THETA 00775 * 00776 IF( ABS( THETA ).GT.BIGTHETA ) THEN 00777 T = HALF / THETA 00778 FASTR( 3 ) = T*D( p ) / D( q ) 00779 FASTR( 4 ) = -T*D( q ) / D( p ) 00780 CALL SROTM( M, A( 1, p ), 1, 00781 $ A( 1, q ), 1, FASTR ) 00782 IF( RSVEC )CALL SROTM( MVL, 00783 $ V( 1, p ), 1, 00784 $ V( 1, q ), 1, 00785 $ FASTR ) 00786 SVA( q ) = AAQQ*SQRT( AMAX1( ZERO, 00787 $ ONE+T*APOAQ*AAPQ ) ) 00788 AAPP = AAPP*SQRT( AMAX1( ZERO, 00789 $ ONE-T*AQOAP*AAPQ ) ) 00790 MXSINJ = AMAX1( MXSINJ, ABS( T ) ) 00791 ELSE 00792 * 00793 * .. choose correct signum for THETA and rotate 00794 * 00795 THSIGN = -SIGN( ONE, AAPQ ) 00796 IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN 00797 T = ONE / ( THETA+THSIGN* 00798 $ SQRT( ONE+THETA*THETA ) ) 00799 CS = SQRT( ONE / ( ONE+T*T ) ) 00800 SN = T*CS 00801 MXSINJ = AMAX1( MXSINJ, ABS( SN ) ) 00802 SVA( q ) = AAQQ*SQRT( AMAX1( ZERO, 00803 $ ONE+T*APOAQ*AAPQ ) ) 00804 AAPP = AAPP*SQRT( AMAX1( ZERO, 00805 $ ONE-T*AQOAP*AAPQ ) ) 00806 * 00807 APOAQ = D( p ) / D( q ) 00808 AQOAP = D( q ) / D( p ) 00809 IF( D( p ).GE.ONE ) THEN 00810 * 00811 IF( D( q ).GE.ONE ) THEN 00812 FASTR( 3 ) = T*APOAQ 00813 FASTR( 4 ) = -T*AQOAP 00814 D( p ) = D( p )*CS 00815 D( q ) = D( q )*CS 00816 CALL SROTM( M, A( 1, p ), 1, 00817 $ A( 1, q ), 1, 00818 $ FASTR ) 00819 IF( RSVEC )CALL SROTM( MVL, 00820 $ V( 1, p ), 1, V( 1, q ), 00821 $ 1, FASTR ) 00822 ELSE 00823 CALL SAXPY( M, -T*AQOAP, 00824 $ A( 1, q ), 1, 00825 $ A( 1, p ), 1 ) 00826 CALL SAXPY( M, CS*SN*APOAQ, 00827 $ A( 1, p ), 1, 00828 $ A( 1, q ), 1 ) 00829 IF( RSVEC ) THEN 00830 CALL SAXPY( MVL, -T*AQOAP, 00831 $ V( 1, q ), 1, 00832 $ V( 1, p ), 1 ) 00833 CALL SAXPY( MVL, 00834 $ CS*SN*APOAQ, 00835 $ V( 1, p ), 1, 00836 $ V( 1, q ), 1 ) 00837 END IF 00838 D( p ) = D( p )*CS 00839 D( q ) = D( q ) / CS 00840 END IF 00841 ELSE 00842 IF( D( q ).GE.ONE ) THEN 00843 CALL SAXPY( M, T*APOAQ, 00844 $ A( 1, p ), 1, 00845 $ A( 1, q ), 1 ) 00846 CALL SAXPY( M, -CS*SN*AQOAP, 00847 $ A( 1, q ), 1, 00848 $ A( 1, p ), 1 ) 00849 IF( RSVEC ) THEN 00850 CALL SAXPY( MVL, T*APOAQ, 00851 $ V( 1, p ), 1, 00852 $ V( 1, q ), 1 ) 00853 CALL SAXPY( MVL, 00854 $ -CS*SN*AQOAP, 00855 $ V( 1, q ), 1, 00856 $ V( 1, p ), 1 ) 00857 END IF 00858 D( p ) = D( p ) / CS 00859 D( q ) = D( q )*CS 00860 ELSE 00861 IF( D( p ).GE.D( q ) ) THEN 00862 CALL SAXPY( M, -T*AQOAP, 00863 $ A( 1, q ), 1, 00864 $ A( 1, p ), 1 ) 00865 CALL SAXPY( M, CS*SN*APOAQ, 00866 $ A( 1, p ), 1, 00867 $ A( 1, q ), 1 ) 00868 D( p ) = D( p )*CS 00869 D( q ) = D( q ) / CS 00870 IF( RSVEC ) THEN 00871 CALL SAXPY( MVL, 00872 $ -T*AQOAP, 00873 $ V( 1, q ), 1, 00874 $ V( 1, p ), 1 ) 00875 CALL SAXPY( MVL, 00876 $ CS*SN*APOAQ, 00877 $ V( 1, p ), 1, 00878 $ V( 1, q ), 1 ) 00879 END IF 00880 ELSE 00881 CALL SAXPY( M, T*APOAQ, 00882 $ A( 1, p ), 1, 00883 $ A( 1, q ), 1 ) 00884 CALL SAXPY( M, 00885 $ -CS*SN*AQOAP, 00886 $ A( 1, q ), 1, 00887 $ A( 1, p ), 1 ) 00888 D( p ) = D( p ) / CS 00889 D( q ) = D( q )*CS 00890 IF( RSVEC ) THEN 00891 CALL SAXPY( MVL, 00892 $ T*APOAQ, V( 1, p ), 00893 $ 1, V( 1, q ), 1 ) 00894 CALL SAXPY( MVL, 00895 $ -CS*SN*AQOAP, 00896 $ V( 1, q ), 1, 00897 $ V( 1, p ), 1 ) 00898 END IF 00899 END IF 00900 END IF 00901 END IF 00902 END IF 00903 * 00904 ELSE 00905 IF( AAPP.GT.AAQQ ) THEN 00906 CALL SCOPY( M, A( 1, p ), 1, WORK, 00907 $ 1 ) 00908 CALL SLASCL( 'G', 0, 0, AAPP, ONE, 00909 $ M, 1, WORK, LDA, IERR ) 00910 CALL SLASCL( 'G', 0, 0, AAQQ, ONE, 00911 $ M, 1, A( 1, q ), LDA, 00912 $ IERR ) 00913 TEMP1 = -AAPQ*D( p ) / D( q ) 00914 CALL SAXPY( M, TEMP1, WORK, 1, 00915 $ A( 1, q ), 1 ) 00916 CALL SLASCL( 'G', 0, 0, ONE, AAQQ, 00917 $ M, 1, A( 1, q ), LDA, 00918 $ IERR ) 00919 SVA( q ) = AAQQ*SQRT( AMAX1( ZERO, 00920 $ ONE-AAPQ*AAPQ ) ) 00921 MXSINJ = AMAX1( MXSINJ, SFMIN ) 00922 ELSE 00923 CALL SCOPY( M, A( 1, q ), 1, WORK, 00924 $ 1 ) 00925 CALL SLASCL( 'G', 0, 0, AAQQ, ONE, 00926 $ M, 1, WORK, LDA, IERR ) 00927 CALL SLASCL( 'G', 0, 0, AAPP, ONE, 00928 $ M, 1, A( 1, p ), LDA, 00929 $ IERR ) 00930 TEMP1 = -AAPQ*D( q ) / D( p ) 00931 CALL SAXPY( M, TEMP1, WORK, 1, 00932 $ A( 1, p ), 1 ) 00933 CALL SLASCL( 'G', 0, 0, ONE, AAPP, 00934 $ M, 1, A( 1, p ), LDA, 00935 $ IERR ) 00936 SVA( p ) = AAPP*SQRT( AMAX1( ZERO, 00937 $ ONE-AAPQ*AAPQ ) ) 00938 MXSINJ = AMAX1( MXSINJ, SFMIN ) 00939 END IF 00940 END IF 00941 * END IF ROTOK THEN ... ELSE 00942 * 00943 * In the case of cancellation in updating SVA(q) 00944 * .. recompute SVA(q) 00945 IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS ) 00946 $ THEN 00947 IF( ( AAQQ.LT.ROOTBIG ) .AND. 00948 $ ( AAQQ.GT.ROOTSFMIN ) ) THEN 00949 SVA( q ) = SNRM2( M, A( 1, q ), 1 )* 00950 $ D( q ) 00951 ELSE 00952 T = ZERO 00953 AAQQ = ONE 00954 CALL SLASSQ( M, A( 1, q ), 1, T, 00955 $ AAQQ ) 00956 SVA( q ) = T*SQRT( AAQQ )*D( q ) 00957 END IF 00958 END IF 00959 IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN 00960 IF( ( AAPP.LT.ROOTBIG ) .AND. 00961 $ ( AAPP.GT.ROOTSFMIN ) ) THEN 00962 AAPP = SNRM2( M, A( 1, p ), 1 )* 00963 $ D( p ) 00964 ELSE 00965 T = ZERO 00966 AAPP = ONE 00967 CALL SLASSQ( M, A( 1, p ), 1, T, 00968 $ AAPP ) 00969 AAPP = T*SQRT( AAPP )*D( p ) 00970 END IF 00971 SVA( p ) = AAPP 00972 END IF 00973 * end of OK rotation 00974 ELSE 00975 NOTROT = NOTROT + 1 00976 PSKIPPED = PSKIPPED + 1 00977 IJBLSK = IJBLSK + 1 00978 END IF 00979 ELSE 00980 NOTROT = NOTROT + 1 00981 PSKIPPED = PSKIPPED + 1 00982 IJBLSK = IJBLSK + 1 00983 END IF 00984 * 00985 IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) ) 00986 $ THEN 00987 SVA( p ) = AAPP 00988 NOTROT = 0 00989 GO TO 2011 00990 END IF 00991 IF( ( i.LE.SWBAND ) .AND. 00992 $ ( PSKIPPED.GT.ROWSKIP ) ) THEN 00993 AAPP = -AAPP 00994 NOTROT = 0 00995 GO TO 2203 00996 END IF 00997 * 00998 2200 CONTINUE 00999 * end of the q-loop 01000 2203 CONTINUE 01001 * 01002 SVA( p ) = AAPP 01003 * 01004 ELSE 01005 IF( AAPP.EQ.ZERO )NOTROT = NOTROT + 01006 $ MIN0( jgl+KBL-1, N ) - jgl + 1 01007 IF( AAPP.LT.ZERO )NOTROT = 0 01008 END IF 01009 01010 2100 CONTINUE 01011 * end of the p-loop 01012 2010 CONTINUE 01013 * end of the jbc-loop 01014 2011 CONTINUE 01015 *2011 bailed out of the jbc-loop 01016 DO 2012 p = igl, MIN0( igl+KBL-1, N ) 01017 SVA( p ) = ABS( SVA( p ) ) 01018 2012 CONTINUE 01019 * 01020 2000 CONTINUE 01021 *2000 :: end of the ibr-loop 01022 * 01023 * .. update SVA(N) 01024 IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) ) 01025 $ THEN 01026 SVA( N ) = SNRM2( M, A( 1, N ), 1 )*D( N ) 01027 ELSE 01028 T = ZERO 01029 AAPP = ONE 01030 CALL SLASSQ( M, A( 1, N ), 1, T, AAPP ) 01031 SVA( N ) = T*SQRT( AAPP )*D( N ) 01032 END IF 01033 * 01034 * Additional steering devices 01035 * 01036 IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR. 01037 $ ( ISWROT.LE.N ) ) )SWBAND = i 01038 * 01039 IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.FLOAT( N )*TOL ) .AND. 01040 $ ( FLOAT( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN 01041 GO TO 1994 01042 END IF 01043 * 01044 IF( NOTROT.GE.EMPTSW )GO TO 1994 01045 01046 1993 CONTINUE 01047 * end i=1:NSWEEP loop 01048 * #:) Reaching this point means that the procedure has comleted the given 01049 * number of iterations. 01050 INFO = NSWEEP - 1 01051 GO TO 1995 01052 1994 CONTINUE 01053 * #:) Reaching this point means that during the i-th sweep all pivots were 01054 * below the given tolerance, causing early exit. 01055 * 01056 INFO = 0 01057 * #:) INFO = 0 confirms successful iterations. 01058 1995 CONTINUE 01059 * 01060 * Sort the vector D. 01061 DO 5991 p = 1, N - 1 01062 q = ISAMAX( N-p+1, SVA( p ), 1 ) + p - 1 01063 IF( p.NE.q ) THEN 01064 TEMP1 = SVA( p ) 01065 SVA( p ) = SVA( q ) 01066 SVA( q ) = TEMP1 01067 TEMP1 = D( p ) 01068 D( p ) = D( q ) 01069 D( q ) = TEMP1 01070 CALL SSWAP( M, A( 1, p ), 1, A( 1, q ), 1 ) 01071 IF( RSVEC )CALL SSWAP( MVL, V( 1, p ), 1, V( 1, q ), 1 ) 01072 END IF 01073 5991 CONTINUE 01074 * 01075 RETURN 01076 * .. 01077 * .. END OF SGSVJ0 01078 * .. 01079 END
1.7.6.1 
