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