3rdparty/lapack/dlaebz.c

   1 #include "clapack.h"
   2
   3 /* Subroutine */ int dlaebz_(integer *ijob, integer *nitmax, integer *n,
   4         integer *mmax, integer *minp, integer *nbmin, doublereal *abstol,
   5         doublereal *reltol, doublereal *pivmin, doublereal *d__, doublereal *
   6         e, doublereal *e2, integer *nval, doublereal *ab, doublereal *c__,
   7         integer *mout, integer *nab, doublereal *work, integer *iwork,
   8         integer *info)
   9 {
  10     /* System generated locals */
  11     integer nab_dim1, nab_offset, ab_dim1, ab_offset, i__1, i__2, i__3, i__4,
  12             i__5, i__6;
  13     doublereal d__1, d__2, d__3, d__4;
  14
  15     /* Local variables */
  16     integer j, kf, ji, kl, jp, jit;
  17     doublereal tmp1, tmp2;
  18     integer itmp1, itmp2, kfnew, klnew;
  19
  20
  21 /*  -- LAPACK auxiliary routine (version 3.1) -- */
  22 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  23 /*     November 2006 */
  24
  25 /*     .. Scalar Arguments .. */
  26 /*     .. */
  27 /*     .. Array Arguments .. */
  28 /*     .. */
  29
  30 /*  Purpose */
  31 /*  ======= */
  32
  33 /*  DLAEBZ contains the iteration loops which compute and use the */
  34 /*  function N(w), which is the count of eigenvalues of a symmetric */
  35 /*  tridiagonal matrix T less than or equal to its argument  w.  It */
  36 /*  performs a choice of two types of loops: */
  37
  38 /*  IJOB=1, followed by */
  39 /*  IJOB=2: It takes as input a list of intervals and returns a list of */
  40 /*          sufficiently small intervals whose union contains the same */
  41 /*          eigenvalues as the union of the original intervals. */
  42 /*          The input intervals are (AB(j,1),AB(j,2)], j=1,...,MINP. */
  43 /*          The output interval (AB(j,1),AB(j,2)] will contain */
  44 /*          eigenvalues NAB(j,1)+1,...,NAB(j,2), where 1 <= j <= MOUT. */
  45
  46 /*  IJOB=3: It performs a binary search in each input interval */
  47 /*          (AB(j,1),AB(j,2)] for a point  w(j)  such that */
  48 /*          N(w(j))=NVAL(j), and uses  C(j)  as the starting point of */
  49 /*          the search.  If such a w(j) is found, then on output */
  50 /*          AB(j,1)=AB(j,2)=w.  If no such w(j) is found, then on output */
  51 /*          (AB(j,1),AB(j,2)] will be a small interval containing the */
  52 /*          point where N(w) jumps through NVAL(j), unless that point */
  53 /*          lies outside the initial interval. */
  54
  55 /*  Note that the intervals are in all cases half-open intervals, */
  56 /*  i.e., of the form  (a,b] , which includes  b  but not  a . */
  57
  58 /*  To avoid underflow, the matrix should be scaled so that its largest */
  59 /*  element is no greater than  overflow**(1/2) * underflow**(1/4) */
  60 /*  in absolute value.  To assure the most accurate computation */
  61 /*  of small eigenvalues, the matrix should be scaled to be */
  62 /*  not much smaller than that, either. */
  63
  64 /*  See W. Kahan "Accurate Eigenvalues of a Symmetric Tridiagonal */
  65 /*  Matrix", Report CS41, Computer Science Dept., Stanford */
  66 /*  University, July 21, 1966 */
  67
  68 /*  Note: the arguments are, in general, *not* checked for unreasonable */
  69 /*  values. */
  70
  71 /*  Arguments */
  72 /*  ========= */
  73
  74 /*  IJOB    (input) INTEGER */
  75 /*          Specifies what is to be done: */
  76 /*          = 1:  Compute NAB for the initial intervals. */
  77 /*          = 2:  Perform bisection iteration to find eigenvalues of T. */
  78 /*          = 3:  Perform bisection iteration to invert N(w), i.e., */
  79 /*                to find a point which has a specified number of */
  80 /*                eigenvalues of T to its left. */
  81 /*          Other values will cause DLAEBZ to return with INFO=-1. */
  82
  83 /*  NITMAX  (input) INTEGER */
  84 /*          The maximum number of "levels" of bisection to be */
  85 /*          performed, i.e., an interval of width W will not be made */
  86 /*          smaller than 2^(-NITMAX) * W.  If not all intervals */
  87 /*          have converged after NITMAX iterations, then INFO is set */
  88 /*          to the number of non-converged intervals. */
  89
  90 /*  N       (input) INTEGER */
  91 /*          The dimension n of the tridiagonal matrix T.  It must be at */
  92 /*          least 1. */
  93
  94 /*  MMAX    (input) INTEGER */
  95 /*          The maximum number of intervals.  If more than MMAX intervals */
  96 /*          are generated, then DLAEBZ will quit with INFO=MMAX+1. */
  97
  98 /*  MINP    (input) INTEGER */
  99 /*          The initial number of intervals.  It may not be greater than */
 100 /*          MMAX. */
 101
 102 /*  NBMIN   (input) INTEGER */
 103 /*          The smallest number of intervals that should be processed */
 104 /*          using a vector loop.  If zero, then only the scalar loop */
 105 /*          will be used. */
 106
 107 /*  ABSTOL  (input) DOUBLE PRECISION */
 108 /*          The minimum (absolute) width of an interval.  When an */
 109 /*          interval is narrower than ABSTOL, or than RELTOL times the */
 110 /*          larger (in magnitude) endpoint, then it is considered to be */
 111 /*          sufficiently small, i.e., converged.  This must be at least */
 112 /*          zero. */
 113
 114 /*  RELTOL  (input) DOUBLE PRECISION */
 115 /*          The minimum relative width of an interval.  When an interval */
 116 /*          is narrower than ABSTOL, or than RELTOL times the larger (in */
 117 /*          magnitude) endpoint, then it is considered to be */
 118 /*          sufficiently small, i.e., converged.  Note: this should */
 119 /*          always be at least radix*machine epsilon. */
 120
 121 /*  PIVMIN  (input) DOUBLE PRECISION */
 122 /*          The minimum absolute value of a "pivot" in the Sturm */
 123 /*          sequence loop.  This *must* be at least  max |e(j)**2| * */
 124 /*          safe_min  and at least safe_min, where safe_min is at least */
 125 /*          the smallest number that can divide one without overflow. */
 126
 127 /*  D       (input) DOUBLE PRECISION array, dimension (N) */
 128 /*          The diagonal elements of the tridiagonal matrix T. */
 129
 130 /*  E       (input) DOUBLE PRECISION array, dimension (N) */
 131 /*          The offdiagonal elements of the tridiagonal matrix T in */
 132 /*          positions 1 through N-1.  E(N) is arbitrary. */
 133
 134 /*  E2      (input) DOUBLE PRECISION array, dimension (N) */
 135 /*          The squares of the offdiagonal elements of the tridiagonal */
 136 /*          matrix T.  E2(N) is ignored. */
 137
 138 /*  NVAL    (input/output) INTEGER array, dimension (MINP) */
 139 /*          If IJOB=1 or 2, not referenced. */
 140 /*          If IJOB=3, the desired values of N(w).  The elements of NVAL */
 141 /*          will be reordered to correspond with the intervals in AB. */
 142 /*          Thus, NVAL(j) on output will not, in general be the same as */
 143 /*          NVAL(j) on input, but it will correspond with the interval */
 144 /*          (AB(j,1),AB(j,2)] on output. */
 145
 146 /*  AB      (input/output) DOUBLE PRECISION array, dimension (MMAX,2) */
 147 /*          The endpoints of the intervals.  AB(j,1) is  a(j), the left */
 148 /*          endpoint of the j-th interval, and AB(j,2) is b(j), the */
 149 /*          right endpoint of the j-th interval.  The input intervals */
 150 /*          will, in general, be modified, split, and reordered by the */
 151 /*          calculation. */
 152
 153 /*  C       (input/output) DOUBLE PRECISION array, dimension (MMAX) */
 154 /*          If IJOB=1, ignored. */
 155 /*          If IJOB=2, workspace. */
 156 /*          If IJOB=3, then on input C(j) should be initialized to the */
 157 /*          first search point in the binary search. */
 158
 159 /*  MOUT    (output) INTEGER */
 160 /*          If IJOB=1, the number of eigenvalues in the intervals. */
 161 /*          If IJOB=2 or 3, the number of intervals output. */
 162 /*          If IJOB=3, MOUT will equal MINP. */
 163
 164 /*  NAB     (input/output) INTEGER array, dimension (MMAX,2) */
 165 /*          If IJOB=1, then on output NAB(i,j) will be set to N(AB(i,j)). */
 166 /*          If IJOB=2, then on input, NAB(i,j) should be set.  It must */
 167 /*             satisfy the condition: */
 168 /*             N(AB(i,1)) <= NAB(i,1) <= NAB(i,2) <= N(AB(i,2)), */
 169 /*             which means that in interval i only eigenvalues */
 170 /*             NAB(i,1)+1,...,NAB(i,2) will be considered.  Usually, */
 171 /*             NAB(i,j)=N(AB(i,j)), from a previous call to DLAEBZ with */
 172 /*             IJOB=1. */
 173 /*             On output, NAB(i,j) will contain */
 174 /*             max(na(k),min(nb(k),N(AB(i,j)))), where k is the index of */
 175 /*             the input interval that the output interval */
 176 /*             (AB(j,1),AB(j,2)] came from, and na(k) and nb(k) are the */
 177 /*             the input values of NAB(k,1) and NAB(k,2). */
 178 /*          If IJOB=3, then on output, NAB(i,j) contains N(AB(i,j)), */
 179 /*             unless N(w) > NVAL(i) for all search points  w , in which */
 180 /*             case NAB(i,1) will not be modified, i.e., the output */
 181 /*             value will be the same as the input value (modulo */
 182 /*             reorderings -- see NVAL and AB), or unless N(w) < NVAL(i) */
 183 /*             for all search points  w , in which case NAB(i,2) will */
 184 /*             not be modified.  Normally, NAB should be set to some */
 185 /*             distinctive value(s) before DLAEBZ is called. */
 186
 187 /*  WORK    (workspace) DOUBLE PRECISION array, dimension (MMAX) */
 188 /*          Workspace. */
 189
 190 /*  IWORK   (workspace) INTEGER array, dimension (MMAX) */
 191 /*          Workspace. */
 192
 193 /*  INFO    (output) INTEGER */
 194 /*          = 0:       All intervals converged. */
 195 /*          = 1--MMAX: The last INFO intervals did not converge. */
 196 /*          = MMAX+1:  More than MMAX intervals were generated. */
 197
 198 /*  Further Details */
 199 /*  =============== */
 200
 201 /*      This routine is intended to be called only by other LAPACK */
 202 /*  routines, thus the interface is less user-friendly.  It is intended */
 203 /*  for two purposes: */
 204
 205 /*  (a) finding eigenvalues.  In this case, DLAEBZ should have one or */
 206 /*      more initial intervals set up in AB, and DLAEBZ should be called */
 207 /*      with IJOB=1.  This sets up NAB, and also counts the eigenvalues. */
 208 /*      Intervals with no eigenvalues would usually be thrown out at */
 209 /*      this point.  Also, if not all the eigenvalues in an interval i */
 210 /*      are desired, NAB(i,1) can be increased or NAB(i,2) decreased. */
 211 /*      For example, set NAB(i,1)=NAB(i,2)-1 to get the largest */
 212 /*      eigenvalue.  DLAEBZ is then called with IJOB=2 and MMAX */
 213 /*      no smaller than the value of MOUT returned by the call with */
 214 /*      IJOB=1.  After this (IJOB=2) call, eigenvalues NAB(i,1)+1 */
 215 /*      through NAB(i,2) are approximately AB(i,1) (or AB(i,2)) to the */
 216 /*      tolerance specified by ABSTOL and RELTOL. */
 217
 218 /*  (b) finding an interval (a',b'] containing eigenvalues w(f),...,w(l). */
 219 /*      In this case, start with a Gershgorin interval  (a,b).  Set up */
 220 /*      AB to contain 2 search intervals, both initially (a,b).  One */
 221 /*      NVAL element should contain  f-1  and the other should contain  l */
 222 /*      , while C should contain a and b, resp.  NAB(i,1) should be -1 */
 223 /*      and NAB(i,2) should be N+1, to flag an error if the desired */
 224 /*      interval does not lie in (a,b).  DLAEBZ is then called with */
 225 /*      IJOB=3.  On exit, if w(f-1) < w(f), then one of the intervals -- */
 226 /*      j -- will have AB(j,1)=AB(j,2) and NAB(j,1)=NAB(j,2)=f-1, while */
 227 /*      if, to the specified tolerance, w(f-k)=...=w(f+r), k > 0 and r */
 228 /*      >= 0, then the interval will have  N(AB(j,1))=NAB(j,1)=f-k and */
 229 /*      N(AB(j,2))=NAB(j,2)=f+r.  The cases w(l) < w(l+1) and */
 230 /*      w(l-r)=...=w(l+k) are handled similarly. */
 231
 232 /*  ===================================================================== */
 233
 234 /*     .. Parameters .. */
 235 /*     .. */
 236 /*     .. Local Scalars .. */
 237 /*     .. */
 238 /*     .. Intrinsic Functions .. */
 239 /*     .. */
 240 /*     .. Executable Statements .. */
 241
 242 /*     Check for Errors */
 243
 244     /* Parameter adjustments */
 245     nab_dim1 = *mmax;
 246     nab_offset = 1 + nab_dim1;
 247     nab -= nab_offset;
 248     ab_dim1 = *mmax;
 249     ab_offset = 1 + ab_dim1;
 250     ab -= ab_offset;
 251     --d__;
 252     --e;
 253     --e2;
 254     --nval;
 255     --c__;
 256     --work;
 257     --iwork;
 258
 259     /* Function Body */
 260     *info = 0;
 261     if (*ijob < 1 || *ijob > 3) {
 262         *info = -1;
 263         return 0;
 264     }
 265
 266 /*     Initialize NAB */
 267
 268     if (*ijob == 1) {
 269
 270 /*        Compute the number of eigenvalues in the initial intervals. */
 271
 272         *mout = 0;
 273 /* DIR$ NOVECTOR */
 274         i__1 = *minp;
 275         for (ji = 1; ji <= i__1; ++ji) {
 276             for (jp = 1; jp <= 2; ++jp) {
 277                 tmp1 = d__[1] - ab[ji + jp * ab_dim1];
 278                 if (abs(tmp1) < *pivmin) {
 279                     tmp1 = -(*pivmin);
 280                 }
 281                 nab[ji + jp * nab_dim1] = 0;
 282                 if (tmp1 <= 0.) {
 283                     nab[ji + jp * nab_dim1] = 1;
 284                 }
 285
 286                 i__2 = *n;
 287                 for (j = 2; j <= i__2; ++j) {
 288                     tmp1 = d__[j] - e2[j - 1] / tmp1 - ab[ji + jp * ab_dim1];
 289                     if (abs(tmp1) < *pivmin) {
 290                         tmp1 = -(*pivmin);
 291                     }
 292                     if (tmp1 <= 0.) {
 293                         ++nab[ji + jp * nab_dim1];
 294                     }
 295 /* L10: */
 296                 }
 297 /* L20: */
 298             }
 299             *mout = *mout + nab[ji + (nab_dim1 << 1)] - nab[ji + nab_dim1];
 300 /* L30: */
 301         }
 302         return 0;
 303     }
 304
 305 /*     Initialize for loop */
 306
 307 /*     KF and KL have the following meaning: */
 308 /*        Intervals 1,...,KF-1 have converged. */
 309 /*        Intervals KF,...,KL  still need to be refined. */
 310
 311     kf = 1;
 312     kl = *minp;
 313
 314 /*     If IJOB=2, initialize C. */
 315 /*     If IJOB=3, use the user-supplied starting point. */
 316
 317     if (*ijob == 2) {
 318         i__1 = *minp;
 319         for (ji = 1; ji <= i__1; ++ji) {
 320             c__[ji] = (ab[ji + ab_dim1] + ab[ji + (ab_dim1 << 1)]) * .5;
 321 /* L40: */
 322         }
 323     }
 324
 325 /*     Iteration loop */
 326
 327     i__1 = *nitmax;
 328     for (jit = 1; jit <= i__1; ++jit) {
 329
 330 /*        Loop over intervals */
 331
 332         if (kl - kf + 1 >= *nbmin && *nbmin > 0) {
 333
 334 /*           Begin of Parallel Version of the loop */
 335
 336             i__2 = kl;
 337             for (ji = kf; ji <= i__2; ++ji) {
 338
 339 /*              Compute N(c), the number of eigenvalues less than c */
 340
 341                 work[ji] = d__[1] - c__[ji];
 342                 iwork[ji] = 0;
 343                 if (work[ji] <= *pivmin) {
 344                     iwork[ji] = 1;
 345 /* Computing MIN */
 346                     d__1 = work[ji], d__2 = -(*pivmin);
 347                     work[ji] = min(d__1,d__2);
 348                 }
 349
 350                 i__3 = *n;
 351                 for (j = 2; j <= i__3; ++j) {
 352                     work[ji] = d__[j] - e2[j - 1] / work[ji] - c__[ji];
 353                     if (work[ji] <= *pivmin) {
 354                         ++iwork[ji];
 355 /* Computing MIN */
 356                         d__1 = work[ji], d__2 = -(*pivmin);
 357                         work[ji] = min(d__1,d__2);
 358                     }
 359 /* L50: */
 360                 }
 361 /* L60: */
 362             }
 363
 364             if (*ijob <= 2) {
 365
 366 /*              IJOB=2: Choose all intervals containing eigenvalues. */
 367
 368                 klnew = kl;
 369                 i__2 = kl;
 370                 for (ji = kf; ji <= i__2; ++ji) {
 371
 372 /*                 Insure that N(w) is monotone */
 373
 374 /* Computing MIN */
 375 /* Computing MAX */
 376                     i__5 = nab[ji + nab_dim1], i__6 = iwork[ji];
 377                     i__3 = nab[ji + (nab_dim1 << 1)], i__4 = max(i__5,i__6);
 378                     iwork[ji] = min(i__3,i__4);
 379
 380 /*                 Update the Queue -- add intervals if both halves */
 381 /*                 contain eigenvalues. */
 382
 383                     if (iwork[ji] == nab[ji + (nab_dim1 << 1)]) {
 384
 385 /*                    No eigenvalue in the upper interval: */
 386 /*                    just use the lower interval. */
 387
 388                         ab[ji + (ab_dim1 << 1)] = c__[ji];
 389
 390                     } else if (iwork[ji] == nab[ji + nab_dim1]) {
 391
 392 /*                    No eigenvalue in the lower interval: */
 393 /*                    just use the upper interval. */
 394
 395                         ab[ji + ab_dim1] = c__[ji];
 396                     } else {
 397                         ++klnew;
 398                         if (klnew <= *mmax) {
 399
 400 /*                       Eigenvalue in both intervals -- add upper to */
 401 /*                       queue. */
 402
 403                             ab[klnew + (ab_dim1 << 1)] = ab[ji + (ab_dim1 <<
 404                                     1)];
 405                             nab[klnew + (nab_dim1 << 1)] = nab[ji + (nab_dim1
 406                                     << 1)];
 407                             ab[klnew + ab_dim1] = c__[ji];
 408                             nab[klnew + nab_dim1] = iwork[ji];
 409                             ab[ji + (ab_dim1 << 1)] = c__[ji];
 410                             nab[ji + (nab_dim1 << 1)] = iwork[ji];
 411                         } else {
 412                             *info = *mmax + 1;
 413                         }
 414                     }
 415 /* L70: */
 416                 }
 417                 if (*info != 0) {
 418                     return 0;
 419                 }
 420                 kl = klnew;
 421             } else {
 422
 423 /*              IJOB=3: Binary search.  Keep only the interval containing */
 424 /*                      w   s.t. N(w) = NVAL */
 425
 426                 i__2 = kl;
 427                 for (ji = kf; ji <= i__2; ++ji) {
 428                     if (iwork[ji] <= nval[ji]) {
 429                         ab[ji + ab_dim1] = c__[ji];
 430                         nab[ji + nab_dim1] = iwork[ji];
 431                     }
 432                     if (iwork[ji] >= nval[ji]) {
 433                         ab[ji + (ab_dim1 << 1)] = c__[ji];
 434                         nab[ji + (nab_dim1 << 1)] = iwork[ji];
 435                     }
 436 /* L80: */
 437                 }
 438             }
 439
 440         } else {
 441
 442 /*           End of Parallel Version of the loop */
 443
 444 /*           Begin of Serial Version of the loop */
 445
 446             klnew = kl;
 447             i__2 = kl;
 448             for (ji = kf; ji <= i__2; ++ji) {
 449
 450 /*              Compute N(w), the number of eigenvalues less than w */
 451
 452                 tmp1 = c__[ji];
 453                 tmp2 = d__[1] - tmp1;
 454                 itmp1 = 0;
 455                 if (tmp2 <= *pivmin) {
 456                     itmp1 = 1;
 457 /* Computing MIN */
 458                     d__1 = tmp2, d__2 = -(*pivmin);
 459                     tmp2 = min(d__1,d__2);
 460                 }
 461
 462 /*              A series of compiler directives to defeat vectorization */
 463 /*              for the next loop */
 464
 465 /* $PL$ CMCHAR=' ' */
 466 /* DIR$          NEXTSCALAR */
 467 /* $DIR          SCALAR */
 468 /* DIR$          NEXT SCALAR */
 469 /* VD$L          NOVECTOR */
 470 /* DEC$          NOVECTOR */
 471 /* VD$           NOVECTOR */
 472 /* VDIR          NOVECTOR */
 473 /* VOCL          LOOP,SCALAR */
 474 /* IBM           PREFER SCALAR */
 475 /* $PL$ CMCHAR='*' */
 476
 477                 i__3 = *n;
 478                 for (j = 2; j <= i__3; ++j) {
 479                     tmp2 = d__[j] - e2[j - 1] / tmp2 - tmp1;
 480                     if (tmp2 <= *pivmin) {
 481                         ++itmp1;
 482 /* Computing MIN */
 483                         d__1 = tmp2, d__2 = -(*pivmin);
 484                         tmp2 = min(d__1,d__2);
 485                     }
 486 /* L90: */
 487                 }
 488
 489                 if (*ijob <= 2) {
 490
 491 /*                 IJOB=2: Choose all intervals containing eigenvalues. */
 492
 493 /*                 Insure that N(w) is monotone */
 494
 495 /* Computing MIN */
 496 /* Computing MAX */
 497                     i__5 = nab[ji + nab_dim1];
 498                     i__3 = nab[ji + (nab_dim1 << 1)], i__4 = max(i__5,itmp1);
 499                     itmp1 = min(i__3,i__4);
 500
 501 /*                 Update the Queue -- add intervals if both halves */
 502 /*                 contain eigenvalues. */
 503
 504                     if (itmp1 == nab[ji + (nab_dim1 << 1)]) {
 505
 506 /*                    No eigenvalue in the upper interval: */
 507 /*                    just use the lower interval. */
 508
 509                         ab[ji + (ab_dim1 << 1)] = tmp1;
 510
 511                     } else if (itmp1 == nab[ji + nab_dim1]) {
 512
 513 /*                    No eigenvalue in the lower interval: */
 514 /*                    just use the upper interval. */
 515
 516                         ab[ji + ab_dim1] = tmp1;
 517                     } else if (klnew < *mmax) {
 518
 519 /*                    Eigenvalue in both intervals -- add upper to queue. */
 520
 521                         ++klnew;
 522                         ab[klnew + (ab_dim1 << 1)] = ab[ji + (ab_dim1 << 1)];
 523                         nab[klnew + (nab_dim1 << 1)] = nab[ji + (nab_dim1 <<
 524                                 1)];
 525                         ab[klnew + ab_dim1] = tmp1;
 526                         nab[klnew + nab_dim1] = itmp1;
 527                         ab[ji + (ab_dim1 << 1)] = tmp1;
 528                         nab[ji + (nab_dim1 << 1)] = itmp1;
 529                     } else {
 530                         *info = *mmax + 1;
 531                         return 0;
 532                     }
 533                 } else {
 534
 535 /*                 IJOB=3: Binary search.  Keep only the interval */
 536 /*                         containing  w  s.t. N(w) = NVAL */
 537
 538                     if (itmp1 <= nval[ji]) {
 539                         ab[ji + ab_dim1] = tmp1;
 540                         nab[ji + nab_dim1] = itmp1;
 541                     }
 542                     if (itmp1 >= nval[ji]) {
 543                         ab[ji + (ab_dim1 << 1)] = tmp1;
 544                         nab[ji + (nab_dim1 << 1)] = itmp1;
 545                     }
 546                 }
 547 /* L100: */
 548             }
 549             kl = klnew;
 550
 551 /*           End of Serial Version of the loop */
 552
 553         }
 554
 555 /*        Check for convergence */
 556
 557         kfnew = kf;
 558         i__2 = kl;
 559         for (ji = kf; ji <= i__2; ++ji) {
 560             tmp1 = (d__1 = ab[ji + (ab_dim1 << 1)] - ab[ji + ab_dim1], abs(
 561                     d__1));
 562 /* Computing MAX */
 563             d__3 = (d__1 = ab[ji + (ab_dim1 << 1)], abs(d__1)), d__4 = (d__2 =
 564                      ab[ji + ab_dim1], abs(d__2));
 565             tmp2 = max(d__3,d__4);
 566 /* Computing MAX */
 567             d__1 = max(*abstol,*pivmin), d__2 = *reltol * tmp2;
 568             if (tmp1 < max(d__1,d__2) || nab[ji + nab_dim1] >= nab[ji + (
 569                     nab_dim1 << 1)]) {
 570
 571 /*              Converged -- Swap with position KFNEW, */
 572 /*                           then increment KFNEW */
 573
 574                 if (ji > kfnew) {
 575                     tmp1 = ab[ji + ab_dim1];
 576                     tmp2 = ab[ji + (ab_dim1 << 1)];
 577                     itmp1 = nab[ji + nab_dim1];
 578                     itmp2 = nab[ji + (nab_dim1 << 1)];
 579                     ab[ji + ab_dim1] = ab[kfnew + ab_dim1];
 580                     ab[ji + (ab_dim1 << 1)] = ab[kfnew + (ab_dim1 << 1)];
 581                     nab[ji + nab_dim1] = nab[kfnew + nab_dim1];
 582                     nab[ji + (nab_dim1 << 1)] = nab[kfnew + (nab_dim1 << 1)];
 583                     ab[kfnew + ab_dim1] = tmp1;
 584                     ab[kfnew + (ab_dim1 << 1)] = tmp2;
 585                     nab[kfnew + nab_dim1] = itmp1;
 586                     nab[kfnew + (nab_dim1 << 1)] = itmp2;
 587                     if (*ijob == 3) {
 588                         itmp1 = nval[ji];
 589                         nval[ji] = nval[kfnew];
 590                         nval[kfnew] = itmp1;
 591                     }
 592                 }
 593                 ++kfnew;
 594             }
 595 /* L110: */
 596         }
 597         kf = kfnew;
 598
 599 /*        Choose Midpoints */
 600
 601         i__2 = kl;
 602         for (ji = kf; ji <= i__2; ++ji) {
 603             c__[ji] = (ab[ji + ab_dim1] + ab[ji + (ab_dim1 << 1)]) * .5;
 604 /* L120: */
 605         }
 606
 607 /*        If no more intervals to refine, quit. */
 608
 609         if (kf > kl) {
 610             goto L140;
 611         }
 612 /* L130: */
 613     }
 614
 615 /*     Converged */
 616
 617 L140:
 618 /* Computing MAX */
 619     i__1 = kl + 1 - kf;
 620     *info = max(i__1,0);
 621     *mout = kl;
 622
 623     return 0;
 624
 625 /*     End of DLAEBZ */
 626
 627 } /* dlaebz_ */