|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
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.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 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 ) / 00379 + AAPQ 00380 IF( AAQQ.GT.AAPP0 )THETA = -THETA 00381 00382 IF( DABS( THETA ).GT.BIGTHETA ) THEN 00383 T = HALF / THETA 00384 FASTR( 3 ) = T*D( p ) / D( q ) 00385 FASTR( 4 ) = -T*D( q ) / D( p ) 00386 CALL DROTM( M, A( 1, p ), 1, 00387 + A( 1, q ), 1, FASTR ) 00388 IF( RSVEC )CALL DROTM( MVL, 00389 + V( 1, p ), 1, 00390 + V( 1, q ), 1, 00391 + FASTR ) 00392 SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, 00393 + ONE+T*APOAQ*AAPQ ) ) 00394 AAPP = AAPP*DSQRT( DMAX1( ZERO, 00395 + ONE-T*AQOAP*AAPQ ) ) 00396 MXSINJ = DMAX1( MXSINJ, DABS( T ) ) 00397 ELSE 00398 * 00399 * .. choose correct signum for THETA and rotate 00400 * 00401 THSIGN = -DSIGN( ONE, AAPQ ) 00402 IF( AAQQ.GT.AAPP0 )THSIGN = -THSIGN 00403 T = ONE / ( THETA+THSIGN* 00404 + DSQRT( ONE+THETA*THETA ) ) 00405 CS = DSQRT( ONE / ( ONE+T*T ) ) 00406 SN = T*CS 00407 MXSINJ = DMAX1( MXSINJ, DABS( SN ) ) 00408 SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, 00409 + ONE+T*APOAQ*AAPQ ) ) 00410 AAPP = AAPP*DSQRT( DMAX1( ZERO, 00411 + ONE-T*AQOAP*AAPQ ) ) 00412 00413 APOAQ = D( p ) / D( q ) 00414 AQOAP = D( q ) / D( p ) 00415 IF( D( p ).GE.ONE ) THEN 00416 * 00417 IF( D( q ).GE.ONE ) THEN 00418 FASTR( 3 ) = T*APOAQ 00419 FASTR( 4 ) = -T*AQOAP 00420 D( p ) = D( p )*CS 00421 D( q ) = D( q )*CS 00422 CALL DROTM( M, A( 1, p ), 1, 00423 + A( 1, q ), 1, 00424 + FASTR ) 00425 IF( RSVEC )CALL DROTM( MVL, 00426 + V( 1, p ), 1, V( 1, q ), 00427 + 1, FASTR ) 00428 ELSE 00429 CALL DAXPY( M, -T*AQOAP, 00430 + A( 1, q ), 1, 00431 + A( 1, p ), 1 ) 00432 CALL DAXPY( M, CS*SN*APOAQ, 00433 + A( 1, p ), 1, 00434 + A( 1, q ), 1 ) 00435 IF( RSVEC ) THEN 00436 CALL DAXPY( MVL, -T*AQOAP, 00437 + V( 1, q ), 1, 00438 + V( 1, p ), 1 ) 00439 CALL DAXPY( MVL, 00440 + CS*SN*APOAQ, 00441 + V( 1, p ), 1, 00442 + V( 1, q ), 1 ) 00443 END IF 00444 D( p ) = D( p )*CS 00445 D( q ) = D( q ) / CS 00446 END IF 00447 ELSE 00448 IF( D( q ).GE.ONE ) THEN 00449 CALL DAXPY( M, T*APOAQ, 00450 + A( 1, p ), 1, 00451 + A( 1, q ), 1 ) 00452 CALL DAXPY( M, -CS*SN*AQOAP, 00453 + A( 1, q ), 1, 00454 + A( 1, p ), 1 ) 00455 IF( RSVEC ) THEN 00456 CALL DAXPY( MVL, T*APOAQ, 00457 + V( 1, p ), 1, 00458 + V( 1, q ), 1 ) 00459 CALL DAXPY( MVL, 00460 + -CS*SN*AQOAP, 00461 + V( 1, q ), 1, 00462 + V( 1, p ), 1 ) 00463 END IF 00464 D( p ) = D( p ) / CS 00465 D( q ) = D( q )*CS 00466 ELSE 00467 IF( D( p ).GE.D( q ) ) THEN 00468 CALL DAXPY( M, -T*AQOAP, 00469 + A( 1, q ), 1, 00470 + A( 1, p ), 1 ) 00471 CALL DAXPY( M, CS*SN*APOAQ, 00472 + A( 1, p ), 1, 00473 + A( 1, q ), 1 ) 00474 D( p ) = D( p )*CS 00475 D( q ) = D( q ) / CS 00476 IF( RSVEC ) THEN 00477 CALL DAXPY( MVL, 00478 + -T*AQOAP, 00479 + V( 1, q ), 1, 00480 + V( 1, p ), 1 ) 00481 CALL DAXPY( MVL, 00482 + CS*SN*APOAQ, 00483 + V( 1, p ), 1, 00484 + V( 1, q ), 1 ) 00485 END IF 00486 ELSE 00487 CALL DAXPY( M, T*APOAQ, 00488 + A( 1, p ), 1, 00489 + A( 1, q ), 1 ) 00490 CALL DAXPY( M, 00491 + -CS*SN*AQOAP, 00492 + A( 1, q ), 1, 00493 + A( 1, p ), 1 ) 00494 D( p ) = D( p ) / CS 00495 D( q ) = D( q )*CS 00496 IF( RSVEC ) THEN 00497 CALL DAXPY( MVL, 00498 + T*APOAQ, V( 1, p ), 00499 + 1, V( 1, q ), 1 ) 00500 CALL DAXPY( MVL, 00501 + -CS*SN*AQOAP, 00502 + V( 1, q ), 1, 00503 + V( 1, p ), 1 ) 00504 END IF 00505 END IF 00506 END IF 00507 END IF 00508 END IF 00509 00510 ELSE 00511 IF( AAPP.GT.AAQQ ) THEN 00512 CALL DCOPY( M, A( 1, p ), 1, WORK, 00513 + 1 ) 00514 CALL DLASCL( 'G', 0, 0, AAPP, ONE, 00515 + M, 1, WORK, LDA, IERR ) 00516 CALL DLASCL( 'G', 0, 0, AAQQ, ONE, 00517 + M, 1, A( 1, q ), LDA, 00518 + IERR ) 00519 TEMP1 = -AAPQ*D( p ) / D( q ) 00520 CALL DAXPY( M, TEMP1, WORK, 1, 00521 + A( 1, q ), 1 ) 00522 CALL DLASCL( 'G', 0, 0, ONE, AAQQ, 00523 + M, 1, A( 1, q ), LDA, 00524 + IERR ) 00525 SVA( q ) = AAQQ*DSQRT( DMAX1( ZERO, 00526 + ONE-AAPQ*AAPQ ) ) 00527 MXSINJ = DMAX1( MXSINJ, SFMIN ) 00528 ELSE 00529 CALL DCOPY( M, A( 1, q ), 1, WORK, 00530 + 1 ) 00531 CALL DLASCL( 'G', 0, 0, AAQQ, ONE, 00532 + M, 1, WORK, LDA, IERR ) 00533 CALL DLASCL( 'G', 0, 0, AAPP, ONE, 00534 + M, 1, A( 1, p ), LDA, 00535 + IERR ) 00536 TEMP1 = -AAPQ*D( q ) / D( p ) 00537 CALL DAXPY( M, TEMP1, WORK, 1, 00538 + A( 1, p ), 1 ) 00539 CALL DLASCL( 'G', 0, 0, ONE, AAPP, 00540 + M, 1, A( 1, p ), LDA, 00541 + IERR ) 00542 SVA( p ) = AAPP*DSQRT( DMAX1( ZERO, 00543 + ONE-AAPQ*AAPQ ) ) 00544 MXSINJ = DMAX1( MXSINJ, SFMIN ) 00545 END IF 00546 END IF 00547 * END IF ROTOK THEN ... ELSE 00548 * 00549 * In the case of cancellation in updating SVA(q) 00550 * .. recompute SVA(q) 00551 IF( ( SVA( q ) / AAQQ )**2.LE.ROOTEPS ) 00552 + THEN 00553 IF( ( AAQQ.LT.ROOTBIG ) .AND. 00554 + ( AAQQ.GT.ROOTSFMIN ) ) THEN 00555 SVA( q ) = DNRM2( M, A( 1, q ), 1 )* 00556 + D( q ) 00557 ELSE 00558 T = ZERO 00559 AAQQ = ONE 00560 CALL DLASSQ( M, A( 1, q ), 1, T, 00561 + AAQQ ) 00562 SVA( q ) = T*DSQRT( AAQQ )*D( q ) 00563 END IF 00564 END IF 00565 IF( ( AAPP / AAPP0 )**2.LE.ROOTEPS ) THEN 00566 IF( ( AAPP.LT.ROOTBIG ) .AND. 00567 + ( AAPP.GT.ROOTSFMIN ) ) THEN 00568 AAPP = DNRM2( M, A( 1, p ), 1 )* 00569 + D( p ) 00570 ELSE 00571 T = ZERO 00572 AAPP = ONE 00573 CALL DLASSQ( M, A( 1, p ), 1, T, 00574 + AAPP ) 00575 AAPP = T*DSQRT( AAPP )*D( p ) 00576 END IF 00577 SVA( p ) = AAPP 00578 END IF 00579 * end of OK rotation 00580 ELSE 00581 NOTROT = NOTROT + 1 00582 * SKIPPED = SKIPPED + 1 00583 PSKIPPED = PSKIPPED + 1 00584 IJBLSK = IJBLSK + 1 00585 END IF 00586 ELSE 00587 NOTROT = NOTROT + 1 00588 PSKIPPED = PSKIPPED + 1 00589 IJBLSK = IJBLSK + 1 00590 END IF 00591 00592 * IF ( NOTROT .GE. EMPTSW ) GO TO 2011 00593 IF( ( i.LE.SWBAND ) .AND. ( IJBLSK.GE.BLSKIP ) ) 00594 + THEN 00595 SVA( p ) = AAPP 00596 NOTROT = 0 00597 GO TO 2011 00598 END IF 00599 IF( ( i.LE.SWBAND ) .AND. 00600 + ( PSKIPPED.GT.ROWSKIP ) ) THEN 00601 AAPP = -AAPP 00602 NOTROT = 0 00603 GO TO 2203 00604 END IF 00605 00606 * 00607 2200 CONTINUE 00608 * end of the q-loop 00609 2203 CONTINUE 00610 00611 SVA( p ) = AAPP 00612 * 00613 ELSE 00614 IF( AAPP.EQ.ZERO )NOTROT = NOTROT + 00615 + MIN0( jgl+KBL-1, N ) - jgl + 1 00616 IF( AAPP.LT.ZERO )NOTROT = 0 00617 *** IF ( NOTROT .GE. EMPTSW ) GO TO 2011 00618 END IF 00619 00620 2100 CONTINUE 00621 * end of the p-loop 00622 2010 CONTINUE 00623 * end of the jbc-loop 00624 2011 CONTINUE 00625 *2011 bailed out of the jbc-loop 00626 DO 2012 p = igl, MIN0( igl+KBL-1, N ) 00627 SVA( p ) = DABS( SVA( p ) ) 00628 2012 CONTINUE 00629 *** IF ( NOTROT .GE. EMPTSW ) GO TO 1994 00630 2000 CONTINUE 00631 *2000 :: end of the ibr-loop 00632 * 00633 * .. update SVA(N) 00634 IF( ( SVA( N ).LT.ROOTBIG ) .AND. ( SVA( N ).GT.ROOTSFMIN ) ) 00635 + THEN 00636 SVA( N ) = DNRM2( M, A( 1, N ), 1 )*D( N ) 00637 ELSE 00638 T = ZERO 00639 AAPP = ONE 00640 CALL DLASSQ( M, A( 1, N ), 1, T, AAPP ) 00641 SVA( N ) = T*DSQRT( AAPP )*D( N ) 00642 END IF 00643 * 00644 * Additional steering devices 00645 * 00646 IF( ( i.LT.SWBAND ) .AND. ( ( MXAAPQ.LE.ROOTTOL ) .OR. 00647 + ( ISWROT.LE.N ) ) )SWBAND = i 00648 00649 IF( ( i.GT.SWBAND+1 ) .AND. ( MXAAPQ.LT.DBLE( N )*TOL ) .AND. 00650 + ( DBLE( N )*MXAAPQ*MXSINJ.LT.TOL ) ) THEN 00651 GO TO 1994 00652 END IF 00653 00654 * 00655 IF( NOTROT.GE.EMPTSW )GO TO 1994 00656 00657 1993 CONTINUE 00658 * end i=1:NSWEEP loop 00659 * #:) Reaching this point means that the procedure has completed the given 00660 * number of sweeps. 00661 INFO = NSWEEP - 1 00662 GO TO 1995 00663 1994 CONTINUE 00664 * #:) Reaching this point means that during the i-th sweep all pivots were 00665 * below the given threshold, causing early exit. 00666 00667 INFO = 0 00668 * #:) INFO = 0 confirms successful iterations. 00669 1995 CONTINUE 00670 * 00671 * Sort the vector D 00672 * 00673 DO 5991 p = 1, N - 1 00674 q = IDAMAX( N-p+1, SVA( p ), 1 ) + p - 1 00675 IF( p.NE.q ) THEN 00676 TEMP1 = SVA( p ) 00677 SVA( p ) = SVA( q ) 00678 SVA( q ) = TEMP1 00679 TEMP1 = D( p ) 00680 D( p ) = D( q ) 00681 D( q ) = TEMP1 00682 CALL DSWAP( M, A( 1, p ), 1, A( 1, q ), 1 ) 00683 IF( RSVEC )CALL DSWAP( MVL, V( 1, p ), 1, V( 1, q ), 1 ) 00684 END IF 00685 5991 CONTINUE 00686 * 00687 RETURN 00688 * .. 00689 * .. END OF DGSVJ1 00690 * .. 00691 END
1.7.3 
