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