1 /* -- translated by f2c (version 20191129). 2 You must link the resulting object file with libf2c: 3 on Microsoft Windows system, link with libf2c.lib; 4 on Linux or Unix systems, link with .../path/to/libf2c.a -lm 5 or, if you install libf2c.a in a standard place, with -lf2c -lm 6 -- in that order, at the end of the command line, as in 7 cc *.o -lf2c -lm 8 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., 9 10 http://www.netlib.org/f2c/libf2c.zip 11 */ 12 13 #include "f2c.h" 14 15 /* Table of constant values */ 16 17 static doublereal c_b5 = 0.; 18 static integer c__1 = 1; 19 static integer c__2 = 2; 20 21 /* > \brief \b DLARRV computes the eigenvectors of the tridiagonal matrix T = L D LT given L, D and the eigenv 22 alues of L D LT. 23 24 =========== DOCUMENTATION =========== 25 26 Online html documentation available at 27 http://www.netlib.org/lapack/explore-html/ 28 29 > \htmlonly 30 > Download DLARRV + dependencies 31 > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarrv. 32 f"> 33 > [TGZ]</a> 34 > <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarrv. 35 f"> 36 > [ZIP]</a> 37 > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarrv. 38 f"> 39 > [TXT]</a> 40 > \endhtmlonly 41 42 Definition: 43 =========== 44 45 SUBROUTINE DLARRV( N, VL, VU, D, L, PIVMIN, 46 ISPLIT, M, DOL, DOU, MINRGP, 47 RTOL1, RTOL2, W, WERR, WGAP, 48 IBLOCK, INDEXW, GERS, Z, LDZ, ISUPPZ, 49 WORK, IWORK, INFO ) 50 51 INTEGER DOL, DOU, INFO, LDZ, M, N 52 DOUBLE PRECISION MINRGP, PIVMIN, RTOL1, RTOL2, VL, VU 53 INTEGER IBLOCK( * ), INDEXW( * ), ISPLIT( * ), 54 $ ISUPPZ( * ), IWORK( * ) 55 DOUBLE PRECISION D( * ), GERS( * ), L( * ), W( * ), WERR( * ), 56 $ WGAP( * ), WORK( * ) 57 DOUBLE PRECISION Z( LDZ, * ) 58 59 60 > \par Purpose: 61 ============= 62 > 63 > \verbatim 64 > 65 > DLARRV computes the eigenvectors of the tridiagonal matrix 66 > T = L D L**T given L, D and APPROXIMATIONS to the eigenvalues of L D L**T. 67 > The input eigenvalues should have been computed by DLARRE. 68 > \endverbatim 69 70 Arguments: 71 ========== 72 73 > \param[in] N 74 > \verbatim 75 > N is INTEGER 76 > The order of the matrix. N >= 0. 77 > \endverbatim 78 > 79 > \param[in] VL 80 > \verbatim 81 > VL is DOUBLE PRECISION 82 > \endverbatim 83 > 84 > \param[in] VU 85 > \verbatim 86 > VU is DOUBLE PRECISION 87 > Lower and upper bounds of the interval that contains the desired 88 > eigenvalues. VL < VU. Needed to compute gaps on the left or right 89 > end of the extremal eigenvalues in the desired RANGE. 90 > \endverbatim 91 > 92 > \param[in,out] D 93 > \verbatim 94 > D is DOUBLE PRECISION array, dimension (N) 95 > On entry, the N diagonal elements of the diagonal matrix D. 96 > On exit, D may be overwritten. 97 > \endverbatim 98 > 99 > \param[in,out] L 100 > \verbatim 101 > L is DOUBLE PRECISION array, dimension (N) 102 > On entry, the (N-1) subdiagonal elements of the unit 103 > bidiagonal matrix L are in elements 1 to N-1 of L 104 > (if the matrix is not splitted.) At the end of each block 105 > is stored the corresponding shift as given by DLARRE. 106 > On exit, L is overwritten. 107 > \endverbatim 108 > 109 > \param[in] PIVMIN 110 > \verbatim 111 > PIVMIN is DOUBLE PRECISION 112 > The minimum pivot allowed in the Sturm sequence. 113 > \endverbatim 114 > 115 > \param[in] ISPLIT 116 > \verbatim 117 > ISPLIT is INTEGER array, dimension (N) 118 > The splitting points, at which T breaks up into blocks. 119 > The first block consists of rows/columns 1 to 120 > ISPLIT( 1 ), the second of rows/columns ISPLIT( 1 )+1 121 > through ISPLIT( 2 ), etc. 122 > \endverbatim 123 > 124 > \param[in] M 125 > \verbatim 126 > M is INTEGER 127 > The total number of input eigenvalues. 0 <= M <= N. 128 > \endverbatim 129 > 130 > \param[in] DOL 131 > \verbatim 132 > DOL is INTEGER 133 > \endverbatim 134 > 135 > \param[in] DOU 136 > \verbatim 137 > DOU is INTEGER 138 > If the user wants to compute only selected eigenvectors from all 139 > the eigenvalues supplied, he can specify an index range DOL:DOU. 140 > Or else the setting DOL=1, DOU=M should be applied. 141 > Note that DOL and DOU refer to the order in which the eigenvalues 142 > are stored in W. 143 > If the user wants to compute only selected eigenpairs, then 144 > the columns DOL-1 to DOU+1 of the eigenvector space Z contain the 145 > computed eigenvectors. All other columns of Z are set to zero. 146 > \endverbatim 147 > 148 > \param[in] MINRGP 149 > \verbatim 150 > MINRGP is DOUBLE PRECISION 151 > \endverbatim 152 > 153 > \param[in] RTOL1 154 > \verbatim 155 > RTOL1 is DOUBLE PRECISION 156 > \endverbatim 157 > 158 > \param[in] RTOL2 159 > \verbatim 160 > RTOL2 is DOUBLE PRECISION 161 > Parameters for bisection. 162 > An interval [LEFT,RIGHT] has converged if 163 > RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) ) 164 > \endverbatim 165 > 166 > \param[in,out] W 167 > \verbatim 168 > W is DOUBLE PRECISION array, dimension (N) 169 > The first M elements of W contain the APPROXIMATE eigenvalues for 170 > which eigenvectors are to be computed. The eigenvalues 171 > should be grouped by split-off block and ordered from 172 > smallest to largest within the block ( The output array 173 > W from DLARRE is expected here ). Furthermore, they are with 174 > respect to the shift of the corresponding root representation 175 > for their block. On exit, W holds the eigenvalues of the 176 > UNshifted matrix. 177 > \endverbatim 178 > 179 > \param[in,out] WERR 180 > \verbatim 181 > WERR is DOUBLE PRECISION array, dimension (N) 182 > The first M elements contain the semiwidth of the uncertainty 183 > interval of the corresponding eigenvalue in W 184 > \endverbatim 185 > 186 > \param[in,out] WGAP 187 > \verbatim 188 > WGAP is DOUBLE PRECISION array, dimension (N) 189 > The separation from the right neighbor eigenvalue in W. 190 > \endverbatim 191 > 192 > \param[in] IBLOCK 193 > \verbatim 194 > IBLOCK is INTEGER array, dimension (N) 195 > The indices of the blocks (submatrices) associated with the 196 > corresponding eigenvalues in W; IBLOCK(i)=1 if eigenvalue 197 > W(i) belongs to the first block from the top, =2 if W(i) 198 > belongs to the second block, etc. 199 > \endverbatim 200 > 201 > \param[in] INDEXW 202 > \verbatim 203 > INDEXW is INTEGER array, dimension (N) 204 > The indices of the eigenvalues within each block (submatrix); 205 > for example, INDEXW(i)= 10 and IBLOCK(i)=2 imply that the 206 > i-th eigenvalue W(i) is the 10-th eigenvalue in the second block. 207 > \endverbatim 208 > 209 > \param[in] GERS 210 > \verbatim 211 > GERS is DOUBLE PRECISION array, dimension (2*N) 212 > The N Gerschgorin intervals (the i-th Gerschgorin interval 213 > is (GERS(2*i-1), GERS(2*i)). The Gerschgorin intervals should 214 > be computed from the original UNshifted matrix. 215 > \endverbatim 216 > 217 > \param[out] Z 218 > \verbatim 219 > Z is DOUBLE PRECISION array, dimension (LDZ, max(1,M) ) 220 > If INFO = 0, the first M columns of Z contain the 221 > orthonormal eigenvectors of the matrix T 222 > corresponding to the input eigenvalues, with the i-th 223 > column of Z holding the eigenvector associated with W(i). 224 > Note: the user must ensure that at least max(1,M) columns are 225 > supplied in the array Z. 226 > \endverbatim 227 > 228 > \param[in] LDZ 229 > \verbatim 230 > LDZ is INTEGER 231 > The leading dimension of the array Z. LDZ >= 1, and if 232 > JOBZ = 'V', LDZ >= max(1,N). 233 > \endverbatim 234 > 235 > \param[out] ISUPPZ 236 > \verbatim 237 > ISUPPZ is INTEGER array, dimension ( 2*max(1,M) ) 238 > The support of the eigenvectors in Z, i.e., the indices 239 > indicating the nonzero elements in Z. The I-th eigenvector 240 > is nonzero only in elements ISUPPZ( 2*I-1 ) through 241 > ISUPPZ( 2*I ). 242 > \endverbatim 243 > 244 > \param[out] WORK 245 > \verbatim 246 > WORK is DOUBLE PRECISION array, dimension (12*N) 247 > \endverbatim 248 > 249 > \param[out] IWORK 250 > \verbatim 251 > IWORK is INTEGER array, dimension (7*N) 252 > \endverbatim 253 > 254 > \param[out] INFO 255 > \verbatim 256 > INFO is INTEGER 257 > = 0: successful exit 258 > 259 > > 0: A problem occured in DLARRV. 260 > < 0: One of the called subroutines signaled an internal problem. 261 > Needs inspection of the corresponding parameter IINFO 262 > for further information. 263 > 264 > =-1: Problem in DLARRB when refining a child's eigenvalues. 265 > =-2: Problem in DLARRF when computing the RRR of a child. 266 > When a child is inside a tight cluster, it can be difficult 267 > to find an RRR. A partial remedy from the user's point of 268 > view is to make the parameter MINRGP smaller and recompile. 269 > However, as the orthogonality of the computed vectors is 270 > proportional to 1/MINRGP, the user should be aware that 271 > he might be trading in precision when he decreases MINRGP. 272 > =-3: Problem in DLARRB when refining a single eigenvalue 273 > after the Rayleigh correction was rejected. 274 > = 5: The Rayleigh Quotient Iteration failed to converge to 275 > full accuracy in MAXITR steps. 276 > \endverbatim 277 278 Authors: 279 ======== 280 281 > \author Univ. of Tennessee 282 > \author Univ. of California Berkeley 283 > \author Univ. of Colorado Denver 284 > \author NAG Ltd. 285 286 > \date September 2012 287 288 > \ingroup doubleOTHERauxiliary 289 290 > \par Contributors: 291 ================== 292 > 293 > Beresford Parlett, University of California, Berkeley, USA \n 294 > Jim Demmel, University of California, Berkeley, USA \n 295 > Inderjit Dhillon, University of Texas, Austin, USA \n 296 > Osni Marques, LBNL/NERSC, USA \n 297 > Christof Voemel, University of California, Berkeley, USA 298 299 ===================================================================== igraphdlarrv_(integer * n,doublereal * vl,doublereal * vu,doublereal * d__,doublereal * l,doublereal * pivmin,integer * isplit,integer * m,integer * dol,integer * dou,doublereal * minrgp,doublereal * rtol1,doublereal * rtol2,doublereal * w,doublereal * werr,doublereal * wgap,integer * iblock,integer * indexw,doublereal * gers,doublereal * z__,integer * ldz,integer * isuppz,doublereal * work,integer * iwork,integer * info)300 Subroutine */ int igraphdlarrv_(integer *n, doublereal *vl, doublereal *vu, 301 doublereal *d__, doublereal *l, doublereal *pivmin, integer *isplit, 302 integer *m, integer *dol, integer *dou, doublereal *minrgp, 303 doublereal *rtol1, doublereal *rtol2, doublereal *w, doublereal *werr, 304 doublereal *wgap, integer *iblock, integer *indexw, doublereal *gers, 305 doublereal *z__, integer *ldz, integer *isuppz, doublereal *work, 306 integer *iwork, integer *info) 307 { 308 /* System generated locals */ 309 integer z_dim1, z_offset, i__1, i__2, i__3, i__4, i__5; 310 doublereal d__1, d__2; 311 logical L__1; 312 313 /* Builtin functions */ 314 double log(doublereal); 315 316 /* Local variables */ 317 integer minwsize, i__, j, k, p, q, miniwsize, ii; 318 doublereal gl; 319 integer im, in; 320 doublereal gu, gap, eps, tau, tol, tmp; 321 integer zto; 322 doublereal ztz; 323 integer iend, jblk; 324 doublereal lgap; 325 integer done; 326 doublereal rgap, left; 327 integer wend, iter; 328 doublereal bstw; 329 integer itmp1; 330 extern /* Subroutine */ int igraphdscal_(integer *, doublereal *, doublereal *, 331 integer *); 332 integer indld; 333 doublereal fudge; 334 integer idone; 335 doublereal sigma; 336 integer iinfo, iindr; 337 doublereal resid; 338 logical eskip; 339 doublereal right; 340 extern /* Subroutine */ int igraphdcopy_(integer *, doublereal *, integer *, 341 doublereal *, integer *); 342 integer nclus, zfrom; 343 doublereal rqtol; 344 integer iindc1, iindc2; 345 extern /* Subroutine */ int igraphdlar1v_(integer *, integer *, integer *, 346 doublereal *, doublereal *, doublereal *, doublereal *, 347 doublereal *, doublereal *, doublereal *, doublereal *, logical *, 348 integer *, doublereal *, doublereal *, integer *, integer *, 349 doublereal *, doublereal *, doublereal *, doublereal *); 350 logical stp2ii; 351 doublereal lambda; 352 extern doublereal igraphdlamch_(char *); 353 integer ibegin, indeig; 354 logical needbs; 355 integer indlld; 356 doublereal sgndef, mingma; 357 extern /* Subroutine */ int igraphdlarrb_(integer *, doublereal *, doublereal *, 358 integer *, integer *, doublereal *, doublereal *, integer *, 359 doublereal *, doublereal *, doublereal *, doublereal *, integer *, 360 doublereal *, doublereal *, integer *, integer *); 361 integer oldien, oldncl, wbegin; 362 doublereal spdiam; 363 integer negcnt; 364 extern /* Subroutine */ int igraphdlarrf_(integer *, doublereal *, doublereal *, 365 doublereal *, integer *, integer *, doublereal *, doublereal *, 366 doublereal *, doublereal *, doublereal *, doublereal *, 367 doublereal *, doublereal *, doublereal *, doublereal *, 368 doublereal *, integer *); 369 integer oldcls; 370 doublereal savgap; 371 integer ndepth; 372 doublereal ssigma; 373 extern /* Subroutine */ int igraphdlaset_(char *, integer *, integer *, 374 doublereal *, doublereal *, doublereal *, integer *); 375 logical usedbs; 376 integer iindwk, offset; 377 doublereal gaptol; 378 integer newcls, oldfst, indwrk, windex, oldlst; 379 logical usedrq; 380 integer newfst, newftt, parity, windmn, windpl, isupmn, newlst, zusedl; 381 doublereal bstres; 382 integer newsiz, zusedu, zusedw; 383 doublereal nrminv, rqcorr; 384 logical tryrqc; 385 integer isupmx; 386 387 388 /* -- LAPACK auxiliary routine (version 3.4.2) -- 389 -- LAPACK is a software package provided by Univ. of Tennessee, -- 390 -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 391 September 2012 392 393 394 ===================================================================== 395 396 The first N entries of WORK are reserved for the eigenvalues 397 Parameter adjustments */ 398 --d__; 399 --l; 400 --isplit; 401 --w; 402 --werr; 403 --wgap; 404 --iblock; 405 --indexw; 406 --gers; 407 z_dim1 = *ldz; 408 z_offset = 1 + z_dim1; 409 z__ -= z_offset; 410 --isuppz; 411 --work; 412 --iwork; 413 414 /* Function Body */ 415 indld = *n + 1; 416 indlld = (*n << 1) + 1; 417 indwrk = *n * 3 + 1; 418 minwsize = *n * 12; 419 i__1 = minwsize; 420 for (i__ = 1; i__ <= i__1; ++i__) { 421 work[i__] = 0.; 422 /* L5: */ 423 } 424 /* IWORK(IINDR+1:IINDR+N) hold the twist indices R for the 425 factorization used to compute the FP vector */ 426 iindr = 0; 427 /* IWORK(IINDC1+1:IINC2+N) are used to store the clusters of the current 428 layer and the one above. */ 429 iindc1 = *n; 430 iindc2 = *n << 1; 431 iindwk = *n * 3 + 1; 432 miniwsize = *n * 7; 433 i__1 = miniwsize; 434 for (i__ = 1; i__ <= i__1; ++i__) { 435 iwork[i__] = 0; 436 /* L10: */ 437 } 438 zusedl = 1; 439 if (*dol > 1) { 440 /* Set lower bound for use of Z */ 441 zusedl = *dol - 1; 442 } 443 zusedu = *m; 444 if (*dou < *m) { 445 /* Set lower bound for use of Z */ 446 zusedu = *dou + 1; 447 } 448 /* The width of the part of Z that is used */ 449 zusedw = zusedu - zusedl + 1; 450 igraphdlaset_("Full", n, &zusedw, &c_b5, &c_b5, &z__[zusedl * z_dim1 + 1], ldz); 451 eps = igraphdlamch_("Precision"); 452 rqtol = eps * 2.; 453 454 /* Set expert flags for standard code. */ 455 tryrqc = TRUE_; 456 if (*dol == 1 && *dou == *m) { 457 } else { 458 /* Only selected eigenpairs are computed. Since the other evalues 459 are not refined by RQ iteration, bisection has to compute to full 460 accuracy. */ 461 *rtol1 = eps * 4.; 462 *rtol2 = eps * 4.; 463 } 464 /* The entries WBEGIN:WEND in W, WERR, WGAP correspond to the 465 desired eigenvalues. The support of the nonzero eigenvector 466 entries is contained in the interval IBEGIN:IEND. 467 Remark that if k eigenpairs are desired, then the eigenvectors 468 are stored in k contiguous columns of Z. 469 DONE is the number of eigenvectors already computed */ 470 done = 0; 471 ibegin = 1; 472 wbegin = 1; 473 i__1 = iblock[*m]; 474 for (jblk = 1; jblk <= i__1; ++jblk) { 475 iend = isplit[jblk]; 476 sigma = l[iend]; 477 /* Find the eigenvectors of the submatrix indexed IBEGIN 478 through IEND. */ 479 wend = wbegin - 1; 480 L15: 481 if (wend < *m) { 482 if (iblock[wend + 1] == jblk) { 483 ++wend; 484 goto L15; 485 } 486 } 487 if (wend < wbegin) { 488 ibegin = iend + 1; 489 goto L170; 490 } else if (wend < *dol || wbegin > *dou) { 491 ibegin = iend + 1; 492 wbegin = wend + 1; 493 goto L170; 494 } 495 /* Find local spectral diameter of the block */ 496 gl = gers[(ibegin << 1) - 1]; 497 gu = gers[ibegin * 2]; 498 i__2 = iend; 499 for (i__ = ibegin + 1; i__ <= i__2; ++i__) { 500 /* Computing MIN */ 501 d__1 = gers[(i__ << 1) - 1]; 502 gl = min(d__1,gl); 503 /* Computing MAX */ 504 d__1 = gers[i__ * 2]; 505 gu = max(d__1,gu); 506 /* L20: */ 507 } 508 spdiam = gu - gl; 509 /* OLDIEN is the last index of the previous block */ 510 oldien = ibegin - 1; 511 /* Calculate the size of the current block */ 512 in = iend - ibegin + 1; 513 /* The number of eigenvalues in the current block */ 514 im = wend - wbegin + 1; 515 /* This is for a 1x1 block */ 516 if (ibegin == iend) { 517 ++done; 518 z__[ibegin + wbegin * z_dim1] = 1.; 519 isuppz[(wbegin << 1) - 1] = ibegin; 520 isuppz[wbegin * 2] = ibegin; 521 w[wbegin] += sigma; 522 work[wbegin] = w[wbegin]; 523 ibegin = iend + 1; 524 ++wbegin; 525 goto L170; 526 } 527 /* The desired (shifted) eigenvalues are stored in W(WBEGIN:WEND) 528 Note that these can be approximations, in this case, the corresp. 529 entries of WERR give the size of the uncertainty interval. 530 The eigenvalue approximations will be refined when necessary as 531 high relative accuracy is required for the computation of the 532 corresponding eigenvectors. */ 533 igraphdcopy_(&im, &w[wbegin], &c__1, &work[wbegin], &c__1); 534 /* We store in W the eigenvalue approximations w.r.t. the original 535 matrix T. */ 536 i__2 = im; 537 for (i__ = 1; i__ <= i__2; ++i__) { 538 w[wbegin + i__ - 1] += sigma; 539 /* L30: */ 540 } 541 /* NDEPTH is the current depth of the representation tree */ 542 ndepth = 0; 543 /* PARITY is either 1 or 0 */ 544 parity = 1; 545 /* NCLUS is the number of clusters for the next level of the 546 representation tree, we start with NCLUS = 1 for the root */ 547 nclus = 1; 548 iwork[iindc1 + 1] = 1; 549 iwork[iindc1 + 2] = im; 550 /* IDONE is the number of eigenvectors already computed in the current 551 block */ 552 idone = 0; 553 /* loop while( IDONE.LT.IM ) 554 generate the representation tree for the current block and 555 compute the eigenvectors */ 556 L40: 557 if (idone < im) { 558 /* This is a crude protection against infinitely deep trees */ 559 if (ndepth > *m) { 560 *info = -2; 561 return 0; 562 } 563 /* breadth first processing of the current level of the representation 564 tree: OLDNCL = number of clusters on current level */ 565 oldncl = nclus; 566 /* reset NCLUS to count the number of child clusters */ 567 nclus = 0; 568 569 parity = 1 - parity; 570 if (parity == 0) { 571 oldcls = iindc1; 572 newcls = iindc2; 573 } else { 574 oldcls = iindc2; 575 newcls = iindc1; 576 } 577 /* Process the clusters on the current level */ 578 i__2 = oldncl; 579 for (i__ = 1; i__ <= i__2; ++i__) { 580 j = oldcls + (i__ << 1); 581 /* OLDFST, OLDLST = first, last index of current cluster. 582 cluster indices start with 1 and are relative 583 to WBEGIN when accessing W, WGAP, WERR, Z */ 584 oldfst = iwork[j - 1]; 585 oldlst = iwork[j]; 586 if (ndepth > 0) { 587 /* Retrieve relatively robust representation (RRR) of cluster 588 that has been computed at the previous level 589 The RRR is stored in Z and overwritten once the eigenvectors 590 have been computed or when the cluster is refined */ 591 if (*dol == 1 && *dou == *m) { 592 /* Get representation from location of the leftmost evalue 593 of the cluster */ 594 j = wbegin + oldfst - 1; 595 } else { 596 if (wbegin + oldfst - 1 < *dol) { 597 /* Get representation from the left end of Z array */ 598 j = *dol - 1; 599 } else if (wbegin + oldfst - 1 > *dou) { 600 /* Get representation from the right end of Z array */ 601 j = *dou; 602 } else { 603 j = wbegin + oldfst - 1; 604 } 605 } 606 igraphdcopy_(&in, &z__[ibegin + j * z_dim1], &c__1, &d__[ibegin] 607 , &c__1); 608 i__3 = in - 1; 609 igraphdcopy_(&i__3, &z__[ibegin + (j + 1) * z_dim1], &c__1, &l[ 610 ibegin], &c__1); 611 sigma = z__[iend + (j + 1) * z_dim1]; 612 /* Set the corresponding entries in Z to zero */ 613 igraphdlaset_("Full", &in, &c__2, &c_b5, &c_b5, &z__[ibegin + j 614 * z_dim1], ldz); 615 } 616 /* Compute DL and DLL of current RRR */ 617 i__3 = iend - 1; 618 for (j = ibegin; j <= i__3; ++j) { 619 tmp = d__[j] * l[j]; 620 work[indld - 1 + j] = tmp; 621 work[indlld - 1 + j] = tmp * l[j]; 622 /* L50: */ 623 } 624 if (ndepth > 0) { 625 /* P and Q are index of the first and last eigenvalue to compute 626 within the current block */ 627 p = indexw[wbegin - 1 + oldfst]; 628 q = indexw[wbegin - 1 + oldlst]; 629 /* Offset for the arrays WORK, WGAP and WERR, i.e., the P-OFFSET 630 through the Q-OFFSET elements of these arrays are to be used. 631 OFFSET = P-OLDFST */ 632 offset = indexw[wbegin] - 1; 633 /* perform limited bisection (if necessary) to get approximate 634 eigenvalues to the precision needed. */ 635 igraphdlarrb_(&in, &d__[ibegin], &work[indlld + ibegin - 1], &p, 636 &q, rtol1, rtol2, &offset, &work[wbegin], &wgap[ 637 wbegin], &werr[wbegin], &work[indwrk], &iwork[ 638 iindwk], pivmin, &spdiam, &in, &iinfo); 639 if (iinfo != 0) { 640 *info = -1; 641 return 0; 642 } 643 /* We also recompute the extremal gaps. W holds all eigenvalues 644 of the unshifted matrix and must be used for computation 645 of WGAP, the entries of WORK might stem from RRRs with 646 different shifts. The gaps from WBEGIN-1+OLDFST to 647 WBEGIN-1+OLDLST are correctly computed in DLARRB. 648 However, we only allow the gaps to become greater since 649 this is what should happen when we decrease WERR */ 650 if (oldfst > 1) { 651 /* Computing MAX */ 652 d__1 = wgap[wbegin + oldfst - 2], d__2 = w[wbegin + 653 oldfst - 1] - werr[wbegin + oldfst - 1] - w[ 654 wbegin + oldfst - 2] - werr[wbegin + oldfst - 655 2]; 656 wgap[wbegin + oldfst - 2] = max(d__1,d__2); 657 } 658 if (wbegin + oldlst - 1 < wend) { 659 /* Computing MAX */ 660 d__1 = wgap[wbegin + oldlst - 1], d__2 = w[wbegin + 661 oldlst] - werr[wbegin + oldlst] - w[wbegin + 662 oldlst - 1] - werr[wbegin + oldlst - 1]; 663 wgap[wbegin + oldlst - 1] = max(d__1,d__2); 664 } 665 /* Each time the eigenvalues in WORK get refined, we store 666 the newly found approximation with all shifts applied in W */ 667 i__3 = oldlst; 668 for (j = oldfst; j <= i__3; ++j) { 669 w[wbegin + j - 1] = work[wbegin + j - 1] + sigma; 670 /* L53: */ 671 } 672 } 673 /* Process the current node. */ 674 newfst = oldfst; 675 i__3 = oldlst; 676 for (j = oldfst; j <= i__3; ++j) { 677 if (j == oldlst) { 678 /* we are at the right end of the cluster, this is also the 679 boundary of the child cluster */ 680 newlst = j; 681 } else if (wgap[wbegin + j - 1] >= *minrgp * (d__1 = work[ 682 wbegin + j - 1], abs(d__1))) { 683 /* the right relative gap is big enough, the child cluster 684 (NEWFST,..,NEWLST) is well separated from the following */ 685 newlst = j; 686 } else { 687 /* inside a child cluster, the relative gap is not 688 big enough. */ 689 goto L140; 690 } 691 /* Compute size of child cluster found */ 692 newsiz = newlst - newfst + 1; 693 /* NEWFTT is the place in Z where the new RRR or the computed 694 eigenvector is to be stored */ 695 if (*dol == 1 && *dou == *m) { 696 /* Store representation at location of the leftmost evalue 697 of the cluster */ 698 newftt = wbegin + newfst - 1; 699 } else { 700 if (wbegin + newfst - 1 < *dol) { 701 /* Store representation at the left end of Z array */ 702 newftt = *dol - 1; 703 } else if (wbegin + newfst - 1 > *dou) { 704 /* Store representation at the right end of Z array */ 705 newftt = *dou; 706 } else { 707 newftt = wbegin + newfst - 1; 708 } 709 } 710 if (newsiz > 1) { 711 712 /* Current child is not a singleton but a cluster. 713 Compute and store new representation of child. 714 715 716 Compute left and right cluster gap. 717 718 LGAP and RGAP are not computed from WORK because 719 the eigenvalue approximations may stem from RRRs 720 different shifts. However, W hold all eigenvalues 721 of the unshifted matrix. Still, the entries in WGAP 722 have to be computed from WORK since the entries 723 in W might be of the same order so that gaps are not 724 exhibited correctly for very close eigenvalues. */ 725 if (newfst == 1) { 726 /* Computing MAX */ 727 d__1 = 0., d__2 = w[wbegin] - werr[wbegin] - *vl; 728 lgap = max(d__1,d__2); 729 } else { 730 lgap = wgap[wbegin + newfst - 2]; 731 } 732 rgap = wgap[wbegin + newlst - 1]; 733 734 /* Compute left- and rightmost eigenvalue of child 735 to high precision in order to shift as close 736 as possible and obtain as large relative gaps 737 as possible */ 738 739 for (k = 1; k <= 2; ++k) { 740 if (k == 1) { 741 p = indexw[wbegin - 1 + newfst]; 742 } else { 743 p = indexw[wbegin - 1 + newlst]; 744 } 745 offset = indexw[wbegin] - 1; 746 igraphdlarrb_(&in, &d__[ibegin], &work[indlld + ibegin 747 - 1], &p, &p, &rqtol, &rqtol, &offset, & 748 work[wbegin], &wgap[wbegin], &werr[wbegin] 749 , &work[indwrk], &iwork[iindwk], pivmin, & 750 spdiam, &in, &iinfo); 751 /* L55: */ 752 } 753 754 if (wbegin + newlst - 1 < *dol || wbegin + newfst - 1 755 > *dou) { 756 /* if the cluster contains no desired eigenvalues 757 skip the computation of that branch of the rep. tree 758 759 We could skip before the refinement of the extremal 760 eigenvalues of the child, but then the representation 761 tree could be different from the one when nothing is 762 skipped. For this reason we skip at this place. */ 763 idone = idone + newlst - newfst + 1; 764 goto L139; 765 } 766 767 /* Compute RRR of child cluster. 768 Note that the new RRR is stored in Z 769 770 DLARRF needs LWORK = 2*N */ 771 igraphdlarrf_(&in, &d__[ibegin], &l[ibegin], &work[indld + 772 ibegin - 1], &newfst, &newlst, &work[wbegin], 773 &wgap[wbegin], &werr[wbegin], &spdiam, &lgap, 774 &rgap, pivmin, &tau, &z__[ibegin + newftt * 775 z_dim1], &z__[ibegin + (newftt + 1) * z_dim1], 776 &work[indwrk], &iinfo); 777 if (iinfo == 0) { 778 /* a new RRR for the cluster was found by DLARRF 779 update shift and store it */ 780 ssigma = sigma + tau; 781 z__[iend + (newftt + 1) * z_dim1] = ssigma; 782 /* WORK() are the midpoints and WERR() the semi-width 783 Note that the entries in W are unchanged. */ 784 i__4 = newlst; 785 for (k = newfst; k <= i__4; ++k) { 786 fudge = eps * 3. * (d__1 = work[wbegin + k - 787 1], abs(d__1)); 788 work[wbegin + k - 1] -= tau; 789 fudge += eps * 4. * (d__1 = work[wbegin + k - 790 1], abs(d__1)); 791 /* Fudge errors */ 792 werr[wbegin + k - 1] += fudge; 793 /* Gaps are not fudged. Provided that WERR is small 794 when eigenvalues are close, a zero gap indicates 795 that a new representation is needed for resolving 796 the cluster. A fudge could lead to a wrong decision 797 of judging eigenvalues 'separated' which in 798 reality are not. This could have a negative impact 799 on the orthogonality of the computed eigenvectors. 800 L116: */ 801 } 802 ++nclus; 803 k = newcls + (nclus << 1); 804 iwork[k - 1] = newfst; 805 iwork[k] = newlst; 806 } else { 807 *info = -2; 808 return 0; 809 } 810 } else { 811 812 /* Compute eigenvector of singleton */ 813 814 iter = 0; 815 816 tol = log((doublereal) in) * 4. * eps; 817 818 k = newfst; 819 windex = wbegin + k - 1; 820 /* Computing MAX */ 821 i__4 = windex - 1; 822 windmn = max(i__4,1); 823 /* Computing MIN */ 824 i__4 = windex + 1; 825 windpl = min(i__4,*m); 826 lambda = work[windex]; 827 ++done; 828 /* Check if eigenvector computation is to be skipped */ 829 if (windex < *dol || windex > *dou) { 830 eskip = TRUE_; 831 goto L125; 832 } else { 833 eskip = FALSE_; 834 } 835 left = work[windex] - werr[windex]; 836 right = work[windex] + werr[windex]; 837 indeig = indexw[windex]; 838 /* Note that since we compute the eigenpairs for a child, 839 all eigenvalue approximations are w.r.t the same shift. 840 In this case, the entries in WORK should be used for 841 computing the gaps since they exhibit even very small 842 differences in the eigenvalues, as opposed to the 843 entries in W which might "look" the same. */ 844 if (k == 1) { 845 /* In the case RANGE='I' and with not much initial 846 accuracy in LAMBDA and VL, the formula 847 LGAP = MAX( ZERO, (SIGMA - VL) + LAMBDA ) 848 can lead to an overestimation of the left gap and 849 thus to inadequately early RQI 'convergence'. 850 Prevent this by forcing a small left gap. 851 Computing MAX */ 852 d__1 = abs(left), d__2 = abs(right); 853 lgap = eps * max(d__1,d__2); 854 } else { 855 lgap = wgap[windmn]; 856 } 857 if (k == im) { 858 /* In the case RANGE='I' and with not much initial 859 accuracy in LAMBDA and VU, the formula 860 can lead to an overestimation of the right gap and 861 thus to inadequately early RQI 'convergence'. 862 Prevent this by forcing a small right gap. 863 Computing MAX */ 864 d__1 = abs(left), d__2 = abs(right); 865 rgap = eps * max(d__1,d__2); 866 } else { 867 rgap = wgap[windex]; 868 } 869 gap = min(lgap,rgap); 870 if (k == 1 || k == im) { 871 /* The eigenvector support can become wrong 872 because significant entries could be cut off due to a 873 large GAPTOL parameter in LAR1V. Prevent this. */ 874 gaptol = 0.; 875 } else { 876 gaptol = gap * eps; 877 } 878 isupmn = in; 879 isupmx = 1; 880 /* Update WGAP so that it holds the minimum gap 881 to the left or the right. This is crucial in the 882 case where bisection is used to ensure that the 883 eigenvalue is refined up to the required precision. 884 The correct value is restored afterwards. */ 885 savgap = wgap[windex]; 886 wgap[windex] = gap; 887 /* We want to use the Rayleigh Quotient Correction 888 as often as possible since it converges quadratically 889 when we are close enough to the desired eigenvalue. 890 However, the Rayleigh Quotient can have the wrong sign 891 and lead us away from the desired eigenvalue. In this 892 case, the best we can do is to use bisection. */ 893 usedbs = FALSE_; 894 usedrq = FALSE_; 895 /* Bisection is initially turned off unless it is forced */ 896 needbs = ! tryrqc; 897 L120: 898 /* Check if bisection should be used to refine eigenvalue */ 899 if (needbs) { 900 /* Take the bisection as new iterate */ 901 usedbs = TRUE_; 902 itmp1 = iwork[iindr + windex]; 903 offset = indexw[wbegin] - 1; 904 d__1 = eps * 2.; 905 igraphdlarrb_(&in, &d__[ibegin], &work[indlld + ibegin 906 - 1], &indeig, &indeig, &c_b5, &d__1, & 907 offset, &work[wbegin], &wgap[wbegin], & 908 werr[wbegin], &work[indwrk], &iwork[ 909 iindwk], pivmin, &spdiam, &itmp1, &iinfo); 910 if (iinfo != 0) { 911 *info = -3; 912 return 0; 913 } 914 lambda = work[windex]; 915 /* Reset twist index from inaccurate LAMBDA to 916 force computation of true MINGMA */ 917 iwork[iindr + windex] = 0; 918 } 919 /* Given LAMBDA, compute the eigenvector. */ 920 L__1 = ! usedbs; 921 igraphdlar1v_(&in, &c__1, &in, &lambda, &d__[ibegin], &l[ 922 ibegin], &work[indld + ibegin - 1], &work[ 923 indlld + ibegin - 1], pivmin, &gaptol, &z__[ 924 ibegin + windex * z_dim1], &L__1, &negcnt, & 925 ztz, &mingma, &iwork[iindr + windex], &isuppz[ 926 (windex << 1) - 1], &nrminv, &resid, &rqcorr, 927 &work[indwrk]); 928 if (iter == 0) { 929 bstres = resid; 930 bstw = lambda; 931 } else if (resid < bstres) { 932 bstres = resid; 933 bstw = lambda; 934 } 935 /* Computing MIN */ 936 i__4 = isupmn, i__5 = isuppz[(windex << 1) - 1]; 937 isupmn = min(i__4,i__5); 938 /* Computing MAX */ 939 i__4 = isupmx, i__5 = isuppz[windex * 2]; 940 isupmx = max(i__4,i__5); 941 ++iter; 942 /* sin alpha <= |resid|/gap 943 Note that both the residual and the gap are 944 proportional to the matrix, so ||T|| doesn't play 945 a role in the quotient 946 947 Convergence test for Rayleigh-Quotient iteration 948 (omitted when Bisection has been used) */ 949 950 if (resid > tol * gap && abs(rqcorr) > rqtol * abs( 951 lambda) && ! usedbs) { 952 /* We need to check that the RQCORR update doesn't 953 move the eigenvalue away from the desired one and 954 towards a neighbor. -> protection with bisection */ 955 if (indeig <= negcnt) { 956 /* The wanted eigenvalue lies to the left */ 957 sgndef = -1.; 958 } else { 959 /* The wanted eigenvalue lies to the right */ 960 sgndef = 1.; 961 } 962 /* We only use the RQCORR if it improves the 963 the iterate reasonably. */ 964 if (rqcorr * sgndef >= 0. && lambda + rqcorr <= 965 right && lambda + rqcorr >= left) { 966 usedrq = TRUE_; 967 /* Store new midpoint of bisection interval in WORK */ 968 if (sgndef == 1.) { 969 /* The current LAMBDA is on the left of the true 970 eigenvalue */ 971 left = lambda; 972 /* We prefer to assume that the error estimate 973 is correct. We could make the interval not 974 as a bracket but to be modified if the RQCORR 975 chooses to. In this case, the RIGHT side should 976 be modified as follows: 977 RIGHT = MAX(RIGHT, LAMBDA + RQCORR) */ 978 } else { 979 /* The current LAMBDA is on the right of the true 980 eigenvalue */ 981 right = lambda; 982 /* See comment about assuming the error estimate is 983 correct above. 984 LEFT = MIN(LEFT, LAMBDA + RQCORR) */ 985 } 986 work[windex] = (right + left) * .5; 987 /* Take RQCORR since it has the correct sign and 988 improves the iterate reasonably */ 989 lambda += rqcorr; 990 /* Update width of error interval */ 991 werr[windex] = (right - left) * .5; 992 } else { 993 needbs = TRUE_; 994 } 995 if (right - left < rqtol * abs(lambda)) { 996 /* The eigenvalue is computed to bisection accuracy 997 compute eigenvector and stop */ 998 usedbs = TRUE_; 999 goto L120; 1000 } else if (iter < 10) { 1001 goto L120; 1002 } else if (iter == 10) { 1003 needbs = TRUE_; 1004 goto L120; 1005 } else { 1006 *info = 5; 1007 return 0; 1008 } 1009 } else { 1010 stp2ii = FALSE_; 1011 if (usedrq && usedbs && bstres <= resid) { 1012 lambda = bstw; 1013 stp2ii = TRUE_; 1014 } 1015 if (stp2ii) { 1016 /* improve error angle by second step */ 1017 L__1 = ! usedbs; 1018 igraphdlar1v_(&in, &c__1, &in, &lambda, &d__[ibegin] 1019 , &l[ibegin], &work[indld + ibegin - 1020 1], &work[indlld + ibegin - 1], 1021 pivmin, &gaptol, &z__[ibegin + windex 1022 * z_dim1], &L__1, &negcnt, &ztz, & 1023 mingma, &iwork[iindr + windex], & 1024 isuppz[(windex << 1) - 1], &nrminv, & 1025 resid, &rqcorr, &work[indwrk]); 1026 } 1027 work[windex] = lambda; 1028 } 1029 1030 /* Compute FP-vector support w.r.t. whole matrix */ 1031 1032 isuppz[(windex << 1) - 1] += oldien; 1033 isuppz[windex * 2] += oldien; 1034 zfrom = isuppz[(windex << 1) - 1]; 1035 zto = isuppz[windex * 2]; 1036 isupmn += oldien; 1037 isupmx += oldien; 1038 /* Ensure vector is ok if support in the RQI has changed */ 1039 if (isupmn < zfrom) { 1040 i__4 = zfrom - 1; 1041 for (ii = isupmn; ii <= i__4; ++ii) { 1042 z__[ii + windex * z_dim1] = 0.; 1043 /* L122: */ 1044 } 1045 } 1046 if (isupmx > zto) { 1047 i__4 = isupmx; 1048 for (ii = zto + 1; ii <= i__4; ++ii) { 1049 z__[ii + windex * z_dim1] = 0.; 1050 /* L123: */ 1051 } 1052 } 1053 i__4 = zto - zfrom + 1; 1054 igraphdscal_(&i__4, &nrminv, &z__[zfrom + windex * z_dim1], 1055 &c__1); 1056 L125: 1057 /* Update W */ 1058 w[windex] = lambda + sigma; 1059 /* Recompute the gaps on the left and right 1060 But only allow them to become larger and not 1061 smaller (which can only happen through "bad" 1062 cancellation and doesn't reflect the theory 1063 where the initial gaps are underestimated due 1064 to WERR being too crude.) */ 1065 if (! eskip) { 1066 if (k > 1) { 1067 /* Computing MAX */ 1068 d__1 = wgap[windmn], d__2 = w[windex] - werr[ 1069 windex] - w[windmn] - werr[windmn]; 1070 wgap[windmn] = max(d__1,d__2); 1071 } 1072 if (windex < wend) { 1073 /* Computing MAX */ 1074 d__1 = savgap, d__2 = w[windpl] - werr[windpl] 1075 - w[windex] - werr[windex]; 1076 wgap[windex] = max(d__1,d__2); 1077 } 1078 } 1079 ++idone; 1080 } 1081 /* here ends the code for the current child */ 1082 1083 L139: 1084 /* Proceed to any remaining child nodes */ 1085 newfst = j + 1; 1086 L140: 1087 ; 1088 } 1089 /* L150: */ 1090 } 1091 ++ndepth; 1092 goto L40; 1093 } 1094 ibegin = iend + 1; 1095 wbegin = wend + 1; 1096 L170: 1097 ; 1098 } 1099 1100 return 0; 1101 1102 /* End of DLARRV */ 1103 1104 } /* igraphdlarrv_ */ 1105 1106