src/gromacs/linearalgebra/gmx_arpack.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2004 David van der Spoel, Erik Lindahl, University of Groningen.
   5  * Copyright (c) 2012,2013,2014,2015,2016,2018,2019, by the GROMACS development team, led by
   6  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   7  * and including many others, as listed in the AUTHORS file in the
   8  * top-level source directory and at http://www.gromacs.org.
   9  *
  10  * GROMACS is free software; you can redistribute it and/or
  11  * modify it under the terms of the GNU Lesser General Public License
  12  * as published by the Free Software Foundation; either version 2.1
  13  * of the License, or (at your option) any later version.
  14  *
  15  * GROMACS is distributed in the hope that it will be useful,
  16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  18  * Lesser General Public License for more details.
  19  *
  20  * You should have received a copy of the GNU Lesser General Public
  21  * License along with GROMACS; if not, see
  22  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  23  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  24  *
  25  * If you want to redistribute modifications to GROMACS, please
  26  * consider that scientific software is very special. Version
  27  * control is crucial - bugs must be traceable. We will be happy to
  28  * consider code for inclusion in the official distribution, but
  29  * derived work must not be called official GROMACS. Details are found
  30  * in the README & COPYING files - if they are missing, get the
  31  * official version at http://www.gromacs.org.
  32  *
  33  * To help us fund GROMACS development, we humbly ask that you cite
  34  * the research papers on the package. Check out http://www.gromacs.org.
  35  */
  36 #include "gmxpre.h"
  37
  38 #include "gmx_arpack.h"
  39
  40 #include <cmath>
  41 #include <cstring>
  42
  43 #include "gromacs/utility/basedefinitions.h"
  44 #include "gromacs/utility/real.h"
  45
  46 #include "gmx_blas.h"
  47 #include "gmx_lapack.h"
  48
  49 static void
  50 F77_FUNC(dstqrb, DSTQRB) (int *      n,
  51                           double *   d__,
  52                           double *   e,
  53                           double *   z__,
  54                           double *   work,
  55                           int *      info)
  56 {
  57     int    i__1, i__2;
  58     double d__1, d__2;
  59     int    c__0  = 0;
  60     int    c__1  = 1;
  61     double c_b31 = 1.;
  62
  63     double b, c__, f, g;
  64     int    i__, j, k, l, m;
  65     double p, r__, s;
  66     int    l1, ii, mm, lm1, mm1, nm1;
  67     double rt1, rt2, eps;
  68     int    lsv;
  69     double tst, eps2;
  70     int    lend, jtot, lendm1, lendp1, iscale;
  71
  72     int    lendsv, nmaxit, icompz;
  73     double ssfmax, ssfmin, safmin, minval, safmax, anorm;
  74
  75
  76     --work;
  77     --z__;
  78     --e;
  79     --d__;
  80
  81     *info = 0;
  82
  83     icompz = 2;
  84
  85     if (*n == 0)
  86     {
  87         return;
  88     }
  89
  90     if (*n == 1)
  91     {
  92         z__[1] = 1.;
  93         return;
  94     }
  95
  96     eps = GMX_DOUBLE_EPS;
  97
  98     d__1   = eps;
  99     eps2   = d__1 * d__1;
 100     minval = GMX_DOUBLE_MIN;
 101     safmin = minval / GMX_DOUBLE_EPS;
 102     safmax = 1. / safmin;
 103     ssfmax = std::sqrt(safmax) / 3.;
 104     ssfmin = std::sqrt(safmin) / eps2;
 105
 106     i__1 = *n - 1;
 107     for (j = 1; j <= i__1; ++j)
 108     {
 109         z__[j] = 0.;
 110
 111     }
 112     z__[*n] = 1.;
 113
 114     nmaxit = *n * 30;
 115     jtot   = 0;
 116
 117     l1  = 1;
 118     nm1 = *n - 1;
 119
 120 L10:
 121     if (l1 > *n)
 122     {
 123         goto L160;
 124     }
 125     if (l1 > 1)
 126     {
 127         e[l1 - 1] = 0.;
 128     }
 129     if (l1 <= nm1)
 130     {
 131         i__1 = nm1;
 132         for (m = l1; m <= i__1; ++m)
 133         {
 134             tst = std::abs(e[m]);
 135             if (tst == 0.)
 136             {
 137                 goto L30;
 138             }
 139             if (tst <= std::sqrt(std::abs(d__[m])) * std::sqrt(std::abs(d__[m+1])) * eps)
 140             {
 141                 e[m] = 0.;
 142                 goto L30;
 143             }
 144         }
 145     }
 146     m = *n;
 147
 148 L30:
 149     l      = l1;
 150     lsv    = l;
 151     lend   = m;
 152     lendsv = lend;
 153     l1     = m + 1;
 154     if (lend == l)
 155     {
 156         goto L10;
 157     }
 158
 159     i__1   = lend - l + 1;
 160     anorm  = F77_FUNC(dlanst, DLANST) ("i", &i__1, &d__[l], &e[l]);
 161     iscale = 0;
 162     if (anorm == 0.)
 163     {
 164         goto L10;
 165     }
 166     if (anorm > ssfmax)
 167     {
 168         iscale = 1;
 169         i__1   = lend - l + 1;
 170         F77_FUNC(dlascl, DLASCL) ("g", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &d__[l], n,
 171                                   info);
 172         i__1 = lend - l;
 173         F77_FUNC(dlascl, DLASCL) ("g", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &e[l], n,
 174                                   info);
 175     }
 176     else if (anorm < ssfmin)
 177     {
 178         iscale = 2;
 179         i__1   = lend - l + 1;
 180         F77_FUNC(dlascl, DLASCL) ("g", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &d__[l], n,
 181                                   info);
 182         i__1 = lend - l;
 183         F77_FUNC(dlascl, DLASCL) ("g", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &e[l], n,
 184                                   info);
 185     }
 186
 187     if (std::abs(d__[lend]) < std::abs(d__[l]))
 188     {
 189         lend = lsv;
 190         l    = lendsv;
 191     }
 192
 193     if (lend > l)
 194     {
 195
 196 L40:
 197         if (l != lend)
 198         {
 199             lendm1 = lend - 1;
 200             i__1   = lendm1;
 201             for (m = l; m <= i__1; ++m)
 202             {
 203                 d__2 = std::abs(e[m]);
 204                 tst  = d__2 * d__2;
 205                 if (tst <= eps2 * std::abs(d__[m]) * std::abs(d__[m + 1]) + safmin)
 206                 {
 207                     goto L60;
 208                 }
 209             }
 210         }
 211
 212         m = lend;
 213
 214 L60:
 215         if (m < lend)
 216         {
 217             e[m] = 0.;
 218         }
 219         p = d__[l];
 220         if (m == l)
 221         {
 222             goto L80;
 223         }
 224
 225         if (m == l + 1)
 226         {
 227             if (icompz > 0)
 228             {
 229                 F77_FUNC(dlaev2, DLAEV2) (&d__[l], &e[l], &d__[l + 1], &rt1, &rt2, &c__, &s);
 230                 work[l]          = c__;
 231                 work[*n - 1 + l] = s;
 232
 233                 tst        = z__[l + 1];
 234                 z__[l + 1] = c__ * tst - s * z__[l];
 235                 z__[l]     = s * tst + c__ * z__[l];
 236             }
 237             else
 238             {
 239                 F77_FUNC(dlae2, DLAE2) (&d__[l], &e[l], &d__[l + 1], &rt1, &rt2);
 240             }
 241             d__[l]     = rt1;
 242             d__[l + 1] = rt2;
 243             e[l]       = 0.;
 244             l         += 2;
 245             if (l <= lend)
 246             {
 247                 goto L40;
 248             }
 249             goto L140;
 250         }
 251
 252         if (jtot == nmaxit)
 253         {
 254             goto L140;
 255         }
 256         ++jtot;
 257
 258         g   = (d__[l + 1] - p) / (e[l] * 2.);
 259         r__ = F77_FUNC(dlapy2, DLAPY2) (&g, &c_b31);
 260         g   = d__[m] - p + e[l] / (g + ((g > 0) ? r__ : -r__ ));
 261
 262         s   = 1.;
 263         c__ = 1.;
 264         p   = 0.;
 265
 266         mm1  = m - 1;
 267         i__1 = l;
 268         for (i__ = mm1; i__ >= i__1; --i__)
 269         {
 270             f = s * e[i__];
 271             b = c__ * e[i__];
 272             F77_FUNC(dlartg, DLARTG) (&g, &f, &c__, &s, &r__);
 273             if (i__ != m - 1)
 274             {
 275                 e[i__ + 1] = r__;
 276             }
 277             g            = d__[i__ + 1] - p;
 278             r__          = (d__[i__] - g) * s + c__ * 2. * b;
 279             p            = s * r__;
 280             d__[i__ + 1] = g + p;
 281             g            = c__ * r__ - b;
 282
 283             if (icompz > 0)
 284             {
 285                 work[i__]          = c__;
 286                 work[*n - 1 + i__] = -s;
 287             }
 288
 289         }
 290
 291         if (icompz > 0)
 292         {
 293             mm = m - l + 1;
 294
 295             F77_FUNC(dlasr, DLASR) ("r", "v", "b", &c__1, &mm, &work[l], &work[*n - 1 + l], &
 296                                     z__[l], &c__1);
 297         }
 298
 299         d__[l] -= p;
 300         e[l]    = g;
 301         goto L40;
 302
 303 L80:
 304         d__[l] = p;
 305
 306         ++l;
 307         if (l <= lend)
 308         {
 309             goto L40;
 310         }
 311         goto L140;
 312
 313     }
 314     else
 315     {
 316
 317 L90:
 318         if (l != lend)
 319         {
 320             lendp1 = lend + 1;
 321             i__1   = lendp1;
 322             for (m = l; m >= i__1; --m)
 323             {
 324                 d__2 = std::abs(e[m - 1]);
 325                 tst  = d__2 * d__2;
 326                 if (tst <= eps2 * std::abs(d__[m]) * std::abs(d__[m- 1]) + safmin)
 327                 {
 328                     goto L110;
 329                 }
 330             }
 331         }
 332
 333         m = lend;
 334
 335 L110:
 336         if (m > lend)
 337         {
 338             e[m - 1] = 0.;
 339         }
 340         p = d__[l];
 341         if (m == l)
 342         {
 343             goto L130;
 344         }
 345
 346         if (m == l - 1)
 347         {
 348             if (icompz > 0)
 349             {
 350                 F77_FUNC(dlaev2, DLAEV2) (&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2, &c__, &s)
 351                 ;
 352
 353                 tst        = z__[l];
 354                 z__[l]     = c__ * tst - s * z__[l - 1];
 355                 z__[l - 1] = s * tst + c__ * z__[l - 1];
 356             }
 357             else
 358             {
 359                 F77_FUNC(dlae2, DLAE2) (&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2);
 360             }
 361             d__[l - 1] = rt1;
 362             d__[l]     = rt2;
 363             e[l - 1]   = 0.;
 364             l         += -2;
 365             if (l >= lend)
 366             {
 367                 goto L90;
 368             }
 369             goto L140;
 370         }
 371
 372         if (jtot == nmaxit)
 373         {
 374             goto L140;
 375         }
 376         ++jtot;
 377
 378
 379         g   = (d__[l - 1] - p) / (e[l - 1] * 2.);
 380         r__ = F77_FUNC(dlapy2, DLAPY2) (&g, &c_b31);
 381         g   = d__[m] - p + e[l - 1] / (g + ((g > 0) ? r__ : -r__ ));
 382
 383         s   = 1.;
 384         c__ = 1.;
 385         p   = 0.;
 386
 387         lm1  = l - 1;
 388         i__1 = lm1;
 389         for (i__ = m; i__ <= i__1; ++i__)
 390         {
 391             f = s * e[i__];
 392             b = c__ * e[i__];
 393             F77_FUNC(dlartg, DLARTG) (&g, &f, &c__, &s, &r__);
 394             if (i__ != m)
 395             {
 396                 e[i__ - 1] = r__;
 397             }
 398             g        = d__[i__] - p;
 399             r__      = (d__[i__ + 1] - g) * s + c__ * 2. * b;
 400             p        = s * r__;
 401             d__[i__] = g + p;
 402             g        = c__ * r__ - b;
 403
 404             if (icompz > 0)
 405             {
 406                 work[i__]          = c__;
 407                 work[*n - 1 + i__] = s;
 408             }
 409
 410         }
 411
 412         if (icompz > 0)
 413         {
 414             mm = l - m + 1;
 415
 416             F77_FUNC(dlasr, DLASR) ("r", "v", "f", &c__1, &mm, &work[m], &work[*n - 1 + m], &
 417                                     z__[m], &c__1);
 418         }
 419
 420         d__[l] -= p;
 421         e[lm1]  = g;
 422         goto L90;
 423
 424 L130:
 425         d__[l] = p;
 426
 427         --l;
 428         if (l >= lend)
 429         {
 430             goto L90;
 431         }
 432         goto L140;
 433
 434     }
 435
 436 L140:
 437     if (iscale == 1)
 438     {
 439         i__1 = lendsv - lsv + 1;
 440         F77_FUNC(dlascl, DLASCL) ("g", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &d__[lsv],
 441                                   n, info);
 442         i__1 = lendsv - lsv;
 443         F77_FUNC(dlascl, DLASCL) ("g", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &e[lsv], n,
 444                                   info);
 445     }
 446     else if (iscale == 2)
 447     {
 448         i__1 = lendsv - lsv + 1;
 449         F77_FUNC(dlascl, DLASCL) ("g", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &d__[lsv],
 450                                   n, info);
 451         i__1 = lendsv - lsv;
 452         F77_FUNC(dlascl, DLASCL) ("g", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &e[lsv], n,
 453                                   info);
 454     }
 455
 456     if (jtot < nmaxit)
 457     {
 458         goto L10;
 459     }
 460     i__1 = *n - 1;
 461     for (i__ = 1; i__ <= i__1; ++i__)
 462     {
 463         if (e[i__] != 0.)
 464         {
 465             ++(*info);
 466         }
 467     }
 468     goto L190;
 469
 470 L160:
 471     if (icompz == 0)
 472     {
 473
 474         F77_FUNC(dlasrt, DLASRT) ("i", n, &d__[1], info);
 475
 476     }
 477     else
 478     {
 479
 480         i__1 = *n;
 481         for (ii = 2; ii <= i__1; ++ii)
 482         {
 483             i__  = ii - 1;
 484             k    = i__;
 485             p    = d__[i__];
 486             i__2 = *n;
 487             for (j = ii; j <= i__2; ++j)
 488             {
 489                 if (d__[j] < p)
 490                 {
 491                     k = j;
 492                     p = d__[j];
 493                 }
 494             }
 495             if (k != i__)
 496             {
 497                 d__[k]   = d__[i__];
 498                 d__[i__] = p;
 499
 500                 p        = z__[k];
 501                 z__[k]   = z__[i__];
 502                 z__[i__] = p;
 503             }
 504         }
 505     }
 506
 507 L190:
 508     return;
 509
 510 }
 511
 512 static void
 513 F77_FUNC(dgetv0, DGETV0) (int *     ido,
 514                           const char *    bmat,
 515                           int gmx_unused *     itry,
 516                           int *     initv,
 517                           int *     n,
 518                           int *     j,
 519                           double *  v,
 520                           int *     ldv,
 521                           double *  resid,
 522                           double *  rnorm,
 523                           int *     ipntr,
 524                           double *  workd,
 525                           int *     iwork,
 526                           int *     ierr)
 527 {
 528     int    c__1  = 1;
 529     double c_b22 = 1.;
 530     double c_b24 = 0.;
 531     double c_b27 = -1.;
 532     int    v_dim1, v_offset, i__1;
 533
 534     int    jj;
 535     int    idist;
 536
 537     --workd;
 538     --resid;
 539     v_dim1   = *ldv;
 540     v_offset = 1 + v_dim1;
 541     v       -= v_offset;
 542     --ipntr;
 543     --iwork;
 544
 545     if (*ido == 0)
 546     {
 547
 548         *ierr    = 0;
 549         iwork[7] = 0;
 550         iwork[5] = 0;
 551         iwork[6] = 0;
 552
 553         if (!(*initv))
 554         {
 555             idist = 2;
 556             F77_FUNC(dlarnv, DLARNV) (&idist, &iwork[1], n, &resid[1]);
 557         }
 558
 559         if (*bmat == 'G')
 560         {
 561             ipntr[1] = 1;
 562             ipntr[2] = *n + 1;
 563             F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[1], &c__1);
 564             *ido = -1;
 565             goto L9000;
 566         }
 567     }
 568
 569     if (iwork[5] == 1)
 570     {
 571         goto L20;
 572     }
 573
 574     if (iwork[6] == 1)
 575     {
 576         goto L40;
 577     }
 578
 579     iwork[5] = 1;
 580     if (*bmat == 'G')
 581     {
 582         F77_FUNC(dcopy, DCOPY) (n, &workd[*n + 1], &c__1, &resid[1], &c__1);
 583         ipntr[1] = *n + 1;
 584         ipntr[2] = 1;
 585         *ido     = 2;
 586         goto L9000;
 587     }
 588     else if (*bmat == 'I')
 589     {
 590         F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[1], &c__1);
 591     }
 592
 593 L20:
 594
 595
 596     iwork[5] = 0;
 597     if (*bmat == 'G')
 598     {
 599         workd[*n * 3 + 4] = F77_FUNC(ddot, DDOT) (n, &resid[1], &c__1, &workd[1], &c__1);
 600         workd[*n * 3 + 4] = std::sqrt(std::abs(workd[*n * 3 + 4]));
 601     }
 602     else if (*bmat == 'I')
 603     {
 604         workd[*n * 3 + 4] = F77_FUNC(dnrm2, DNRM2) (n, &resid[1], &c__1);
 605     }
 606     *rnorm = workd[*n * 3 + 4];
 607
 608     if (*j == 1)
 609     {
 610         goto L50;
 611     }
 612     iwork[6] = 1;
 613 L30:
 614
 615     i__1 = *j - 1;
 616     F77_FUNC(dgemv, DGEMV) ("T", n, &i__1, &c_b22, &v[v_offset], ldv, &workd[1], &c__1, &c_b24,
 617                             &workd[*n + 1], &c__1);
 618     i__1 = *j - 1;
 619     F77_FUNC(dgemv, DGEMV) ("N", n, &i__1, &c_b27, &v[v_offset], ldv, &workd[*n + 1], &c__1, &
 620                             c_b22, &resid[1], &c__1);
 621
 622     if (*bmat == 'G')
 623     {
 624         F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[*n + 1], &c__1);
 625         ipntr[1] = *n + 1;
 626         ipntr[2] = 1;
 627         *ido     = 2;
 628         goto L9000;
 629     }
 630     else if (*bmat == 'I')
 631     {
 632         F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[1], &c__1);
 633     }
 634
 635 L40:
 636
 637     if (*bmat == 'G')
 638     {
 639         *rnorm = F77_FUNC(ddot, DDOT) (n, &resid[1], &c__1, &workd[1], &c__1);
 640         *rnorm = std::sqrt(std::abs(*rnorm));
 641     }
 642     else if (*bmat == 'I')
 643     {
 644         *rnorm = F77_FUNC(dnrm2, DNRM2) (n, &resid[1], &c__1);
 645     }
 646
 647     if (*rnorm > workd[*n * 3 + 4] * .717F)
 648     {
 649         goto L50;
 650     }
 651
 652     ++iwork[7];
 653     if (iwork[7] <= 1)
 654     {
 655
 656         workd[*n * 3 + 4] = *rnorm;
 657         goto L30;
 658     }
 659     else
 660     {
 661
 662         i__1 = *n;
 663         for (jj = 1; jj <= i__1; ++jj)
 664         {
 665             resid[jj] = 0.;
 666         }
 667         *rnorm = 0.;
 668         *ierr  = -1;
 669     }
 670
 671 L50:
 672
 673     *ido = 99;
 674
 675 L9000:
 676     return;
 677 }
 678
 679
 680
 681
 682
 683 static void
 684 F77_FUNC(dsapps, DSAPPS) (int *     n,
 685                           int *     kev,
 686                           int *     np,
 687                           double *  shift,
 688                           double *  v,
 689                           int *     ldv,
 690                           double *  h__,
 691                           int *     ldh,
 692                           double *  resid,
 693                           double *  q,
 694                           int *     ldq,
 695                           double *  workd)
 696 {
 697     double c_b4  = 0.;
 698     double c_b5  = 1.;
 699     double c_b14 = -1.;
 700     int    c__1  = 1;
 701     int    h_dim1, h_offset, q_dim1, q_offset, v_dim1, v_offset, i__1, i__2,
 702            i__3, i__4;
 703     double c__, f, g;
 704     int    i__, j;
 705     double r__, s, a1, a2, a3, a4;
 706     int    jj;
 707     double big;
 708     int    iend, itop;
 709     double epsmch;
 710     int    istart, kplusp;
 711
 712     --workd;
 713     --resid;
 714     --shift;
 715     v_dim1   = *ldv;
 716     v_offset = 1 + v_dim1;
 717     v       -= v_offset;
 718     h_dim1   = *ldh;
 719     h_offset = 1 + h_dim1;
 720     h__     -= h_offset;
 721     q_dim1   = *ldq;
 722     q_offset = 1 + q_dim1;
 723     q       -= q_offset;
 724
 725     epsmch = GMX_DOUBLE_EPS;
 726     itop   = 1;
 727
 728
 729     kplusp = *kev + *np;
 730
 731     F77_FUNC(dlaset, DLASET) ("All", &kplusp, &kplusp, &c_b4, &c_b5, &q[q_offset], ldq);
 732
 733     if (*np == 0)
 734     {
 735         goto L9000;
 736     }
 737
 738     i__1 = *np;
 739     for (jj = 1; jj <= i__1; ++jj)
 740     {
 741
 742         istart = itop;
 743
 744 L20:
 745
 746         i__2 = kplusp - 1;
 747         for (i__ = istart; i__ <= i__2; ++i__)
 748         {
 749             big = std::abs(h__[i__ + (h_dim1*2)]) + std::abs(h__[i__ + 1 + (h_dim1*2)]);
 750             if (h__[i__ + 1 + h_dim1] <= epsmch * big)
 751             {
 752                 h__[i__ + 1 + h_dim1] = 0.;
 753                 iend                  = i__;
 754                 goto L40;
 755             }
 756         }
 757         iend = kplusp;
 758 L40:
 759
 760         if (istart < iend)
 761         {
 762
 763             f = h__[istart + (h_dim1 << 1)] - shift[jj];
 764             g = h__[istart + 1 + h_dim1];
 765             F77_FUNC(dlartg, DLARTG) (&f, &g, &c__, &s, &r__);
 766
 767             a1 = c__ * h__[istart + (h_dim1 << 1)] + s * h__[istart + 1 +
 768                                                              h_dim1];
 769             a2 = c__ * h__[istart + 1 + h_dim1] + s * h__[istart + 1 + (
 770                                                               h_dim1 << 1)];
 771             a4 = c__ * h__[istart + 1 + (h_dim1 << 1)] - s * h__[istart + 1 +
 772                                                                  h_dim1];
 773             a3 = c__ * h__[istart + 1 + h_dim1] - s * h__[istart + (h_dim1 <<
 774                                                                     1)];
 775             h__[istart + (h_dim1 << 1)]     = c__ * a1 + s * a2;
 776             h__[istart + 1 + (h_dim1 << 1)] = c__ * a4 - s * a3;
 777             h__[istart + 1 + h_dim1]        = c__ * a3 + s * a4;
 778
 779             i__3 = istart + jj;
 780             i__2 = (i__3 < kplusp) ? i__3 : kplusp;
 781             for (j = 1; j <= i__2; ++j)
 782             {
 783                 a1 = c__ * q[j + istart * q_dim1] + s * q[j + (istart + 1) *
 784                                                           q_dim1];
 785                 q[j + (istart + 1) * q_dim1] = -s * q[j + istart * q_dim1] +
 786                     c__ * q[j + (istart + 1) * q_dim1];
 787                 q[j + istart * q_dim1] = a1;
 788
 789             }
 790
 791             i__2 = iend - 1;
 792             for (i__ = istart + 1; i__ <= i__2; ++i__)
 793             {
 794
 795                 f = h__[i__ + h_dim1];
 796                 g = s * h__[i__ + 1 + h_dim1];
 797
 798                 h__[i__ + 1 + h_dim1] = c__ * h__[i__ + 1 + h_dim1];
 799                 F77_FUNC(dlartg, DLARTG) (&f, &g, &c__, &s, &r__);
 800
 801                 if (r__ < 0.)
 802                 {
 803                     r__ = -r__;
 804                     c__ = -c__;
 805                     s   = -s;
 806                 }
 807
 808                 h__[i__ + h_dim1] = r__;
 809
 810                 a1 = c__ * h__[i__ + (h_dim1 << 1)] + s * h__[i__ + 1 +
 811                                                               h_dim1];
 812                 a2 = c__ * h__[i__ + 1 + h_dim1] + s * h__[i__ + 1 + (h_dim1
 813                                                                       << 1)];
 814                 a3 = c__ * h__[i__ + 1 + h_dim1] - s * h__[i__ + (h_dim1 << 1)
 815                     ];
 816                 a4 = c__ * h__[i__ + 1 + (h_dim1 << 1)] - s * h__[i__ + 1 +
 817                                                                   h_dim1];
 818
 819                 h__[i__ + (h_dim1 << 1)]     = c__ * a1 + s * a2;
 820                 h__[i__ + 1 + (h_dim1 << 1)] = c__ * a4 - s * a3;
 821                 h__[i__ + 1 + h_dim1]        = c__ * a3 + s * a4;
 822
 823                 i__4 = j + jj;
 824                 i__3 = (i__4 < kplusp) ? i__4 : kplusp;
 825                 for (j = 1; j <= i__3; ++j)
 826                 {
 827                     a1 = c__ * q[j + i__ * q_dim1] + s * q[j + (i__ + 1) *
 828                                                            q_dim1];
 829                     q[j + (i__ + 1) * q_dim1] = -s * q[j + i__ * q_dim1] +
 830                         c__ * q[j + (i__ + 1) * q_dim1];
 831                     q[j + i__ * q_dim1] = a1;
 832                 }
 833
 834             }
 835
 836         }
 837
 838         istart = iend + 1;
 839
 840         if (h__[iend + h_dim1] < 0.)
 841         {
 842             h__[iend + h_dim1] = -h__[iend + h_dim1];
 843             F77_FUNC(dscal, DSCAL) (&kplusp, &c_b14, &q[iend * q_dim1 + 1], &c__1);
 844         }
 845
 846         if (iend < kplusp)
 847         {
 848             goto L20;
 849         }
 850
 851         i__2 = kplusp - 1;
 852         for (i__ = itop; i__ <= i__2; ++i__)
 853         {
 854             if (h__[i__ + 1 + h_dim1] > 0.)
 855             {
 856                 goto L90;
 857             }
 858             ++itop;
 859         }
 860
 861 L90:
 862         ;
 863     }
 864
 865     i__1 = kplusp - 1;
 866     for (i__ = itop; i__ <= i__1; ++i__)
 867     {
 868         big = std::abs(h__[i__ + (h_dim1*2)]) + std::abs(h__[i__+ 1 + (h_dim1*2)]);
 869         if (h__[i__ + 1 + h_dim1] <= epsmch * big)
 870         {
 871             h__[i__ + 1 + h_dim1] = 0.;
 872         }
 873
 874     }
 875
 876     if (h__[*kev + 1 + h_dim1] > 0.)
 877     {
 878         F77_FUNC(dgemv, DGEMV) ("N", n, &kplusp, &c_b5, &v[v_offset], ldv, &q[(*kev + 1) *
 879                                                                               q_dim1 + 1], &c__1, &c_b4, &workd[*n + 1], &c__1);
 880     }
 881
 882     i__1 = *kev;
 883     for (i__ = 1; i__ <= i__1; ++i__)
 884     {
 885         i__2 = kplusp - i__ + 1;
 886         F77_FUNC(dgemv, DGEMV) ("N", n, &i__2, &c_b5, &v[v_offset], ldv, &q[(*kev - i__ + 1) *
 887                                                                             q_dim1 + 1], &c__1, &c_b4, &workd[1], &c__1);
 888         F77_FUNC(dcopy, DCOPY) (n, &workd[1], &c__1, &v[(kplusp - i__ + 1) * v_dim1 + 1], &
 889                                 c__1);
 890
 891     }
 892
 893     F77_FUNC(dlacpy, DLACPY) ("All", n, kev, &v[(*np + 1) * v_dim1 + 1], ldv, &v[v_offset], ldv);
 894
 895     if (h__[*kev + 1 + h_dim1] > 0.)
 896     {
 897         F77_FUNC(dcopy, DCOPY) (n, &workd[*n + 1], &c__1, &v[(*kev + 1) * v_dim1 + 1], &c__1);
 898     }
 899
 900     F77_FUNC(dscal, DSCAL) (n, &q[kplusp + *kev * q_dim1], &resid[1], &c__1);
 901     if (h__[*kev + 1 + h_dim1] > 0.)
 902     {
 903         F77_FUNC(daxpy, DAXPY) (n, &h__[*kev + 1 + h_dim1], &v[(*kev + 1) * v_dim1 + 1], &c__1,
 904                                 &resid[1], &c__1);
 905     }
 906
 907
 908
 909 L9000:
 910     return;
 911
 912 }
 913
 914
 915
 916 static void
 917 F77_FUNC(dsortr, DSORTR) (const char *    which,
 918                           int *     apply,
 919                           int *     n,
 920                           double *  x1,
 921                           double *  x2)
 922 {
 923     int    i__1;
 924
 925     int    i__, j, igap;
 926     double temp;
 927
 928
 929
 930     igap = *n / 2;
 931
 932     if (!std::strncmp(which, "SA", 2))
 933     {
 934
 935 L10:
 936         if (igap == 0)
 937         {
 938             goto L9000;
 939         }
 940         i__1 = *n - 1;
 941         for (i__ = igap; i__ <= i__1; ++i__)
 942         {
 943             j = i__ - igap;
 944 L20:
 945
 946             if (j < 0)
 947             {
 948                 goto L30;
 949             }
 950
 951             if (x1[j] < x1[j + igap])
 952             {
 953                 temp         = x1[j];
 954                 x1[j]        = x1[j + igap];
 955                 x1[j + igap] = temp;
 956                 if (*apply)
 957                 {
 958                     temp         = x2[j];
 959                     x2[j]        = x2[j + igap];
 960                     x2[j + igap] = temp;
 961                 }
 962             }
 963             else
 964             {
 965                 goto L30;
 966             }
 967             j -= igap;
 968             goto L20;
 969 L30:
 970             ;
 971         }
 972         igap /= 2;
 973         goto L10;
 974
 975     }
 976     else if (!std::strncmp(which, "SM", 2))
 977     {
 978
 979 L40:
 980         if (igap == 0)
 981         {
 982             goto L9000;
 983         }
 984         i__1 = *n - 1;
 985         for (i__ = igap; i__ <= i__1; ++i__)
 986         {
 987             j = i__ - igap;
 988 L50:
 989
 990             if (j < 0)
 991             {
 992                 goto L60;
 993             }
 994
 995             if (std::abs(x1[j]) < std::abs(x1[j + igap]))
 996             {
 997                 temp         = x1[j];
 998                 x1[j]        = x1[j + igap];
 999                 x1[j + igap] = temp;
1000                 if (*apply)
1001                 {
1002                     temp         = x2[j];
1003                     x2[j]        = x2[j + igap];
1004                     x2[j + igap] = temp;
1005                 }
1006             }
1007             else
1008             {
1009                 goto L60;
1010             }
1011             j -= igap;
1012             goto L50;
1013 L60:
1014             ;
1015         }
1016         igap /= 2;
1017         goto L40;
1018
1019     }
1020     else if (!std::strncmp(which, "LA", 2))
1021     {
1022
1023 L70:
1024         if (igap == 0)
1025         {
1026             goto L9000;
1027         }
1028         i__1 = *n - 1;
1029         for (i__ = igap; i__ <= i__1; ++i__)
1030         {
1031             j = i__ - igap;
1032 L80:
1033
1034             if (j < 0)
1035             {
1036                 goto L90;
1037             }
1038
1039             if (x1[j] > x1[j + igap])
1040             {
1041                 temp         = x1[j];
1042                 x1[j]        = x1[j + igap];
1043                 x1[j + igap] = temp;
1044                 if (*apply)
1045                 {
1046                     temp         = x2[j];
1047                     x2[j]        = x2[j + igap];
1048                     x2[j + igap] = temp;
1049                 }
1050             }
1051             else
1052             {
1053                 goto L90;
1054             }
1055             j -= igap;
1056             goto L80;
1057 L90:
1058             ;
1059         }
1060         igap /= 2;
1061         goto L70;
1062
1063     }
1064     else if (!std::strncmp(which, "LM", 2))
1065     {
1066
1067
1068 L100:
1069         if (igap == 0)
1070         {
1071             goto L9000;
1072         }
1073         i__1 = *n - 1;
1074         for (i__ = igap; i__ <= i__1; ++i__)
1075         {
1076             j = i__ - igap;
1077 L110:
1078
1079             if (j < 0)
1080             {
1081                 goto L120;
1082             }
1083
1084             if (std::abs(x1[j]) > std::abs(x1[j + igap]))
1085             {
1086                 temp         = x1[j];
1087                 x1[j]        = x1[j + igap];
1088                 x1[j + igap] = temp;
1089                 if (*apply)
1090                 {
1091                     temp         = x2[j];
1092                     x2[j]        = x2[j + igap];
1093                     x2[j + igap] = temp;
1094                 }
1095             }
1096             else
1097             {
1098                 goto L120;
1099             }
1100             j -= igap;
1101             goto L110;
1102 L120:
1103             ;
1104         }
1105         igap /= 2;
1106         goto L100;
1107     }
1108
1109 L9000:
1110     return;
1111
1112 }
1113
1114
1115
1116
1117 static void
1118 F77_FUNC(dsesrt, DSESRT) (const char *    which,
1119                           int *     apply,
1120                           int *     n,
1121                           double *  x,
1122                           int *     na,
1123                           double *  a,
1124                           int *     lda)
1125 {
1126     int    a_dim1, a_offset, i__1;
1127     int    c__1 = 1;
1128
1129     int    i__, j, igap;
1130     double temp;
1131
1132     a_dim1   = *lda;
1133     a_offset = 1 + a_dim1 * 0;
1134     a       -= a_offset;
1135
1136     igap = *n / 2;
1137
1138     if (!std::strncmp(which, "SA", 2))
1139     {
1140
1141 L10:
1142         if (igap == 0)
1143         {
1144             goto L9000;
1145         }
1146         i__1 = *n - 1;
1147         for (i__ = igap; i__ <= i__1; ++i__)
1148         {
1149             j = i__ - igap;
1150 L20:
1151
1152             if (j < 0)
1153             {
1154                 goto L30;
1155             }
1156
1157             if (x[j] < x[j + igap])
1158             {
1159                 temp        = x[j];
1160                 x[j]        = x[j + igap];
1161                 x[j + igap] = temp;
1162                 if (*apply)
1163                 {
1164                     F77_FUNC(dswap, DSWAP) (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
1165                                                                               a_dim1 + 1], &c__1);
1166                 }
1167             }
1168             else
1169             {
1170                 goto L30;
1171             }
1172             j -= igap;
1173             goto L20;
1174 L30:
1175             ;
1176         }
1177         igap /= 2;
1178         goto L10;
1179
1180     }
1181     else if (!std::strncmp(which, "SM", 2))
1182     {
1183
1184 L40:
1185         if (igap == 0)
1186         {
1187             goto L9000;
1188         }
1189         i__1 = *n - 1;
1190         for (i__ = igap; i__ <= i__1; ++i__)
1191         {
1192             j = i__ - igap;
1193 L50:
1194
1195             if (j < 0)
1196             {
1197                 goto L60;
1198             }
1199
1200             if (std::abs(x[j]) < std::abs(x[j + igap]))
1201             {
1202                 temp        = x[j];
1203                 x[j]        = x[j + igap];
1204                 x[j + igap] = temp;
1205                 if (*apply)
1206                 {
1207                     F77_FUNC(dswap, DSWAP) (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
1208                                                                               a_dim1 + 1], &c__1);
1209                 }
1210             }
1211             else
1212             {
1213                 goto L60;
1214             }
1215             j -= igap;
1216             goto L50;
1217 L60:
1218             ;
1219         }
1220         igap /= 2;
1221         goto L40;
1222
1223     }
1224     else if (!std::strncmp(which, "LA", 2))
1225     {
1226
1227 L70:
1228         if (igap == 0)
1229         {
1230             goto L9000;
1231         }
1232         i__1 = *n - 1;
1233         for (i__ = igap; i__ <= i__1; ++i__)
1234         {
1235             j = i__ - igap;
1236 L80:
1237
1238             if (j < 0)
1239             {
1240                 goto L90;
1241             }
1242
1243             if (x[j] > x[j + igap])
1244             {
1245                 temp        = x[j];
1246                 x[j]        = x[j + igap];
1247                 x[j + igap] = temp;
1248                 if (*apply)
1249                 {
1250                     F77_FUNC(dswap, DSWAP) (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
1251                                                                               a_dim1 + 1], &c__1);
1252                 }
1253             }
1254             else
1255             {
1256                 goto L90;
1257             }
1258             j -= igap;
1259             goto L80;
1260 L90:
1261             ;
1262         }
1263         igap /= 2;
1264         goto L70;
1265
1266     }
1267     else if (!std::strncmp(which, "LM", 2))
1268     {
1269
1270 L100:
1271         if (igap == 0)
1272         {
1273             goto L9000;
1274         }
1275         i__1 = *n - 1;
1276         for (i__ = igap; i__ <= i__1; ++i__)
1277         {
1278             j = i__ - igap;
1279 L110:
1280
1281             if (j < 0)
1282             {
1283                 goto L120;
1284             }
1285
1286             if (std::abs(x[j]) > std::abs(x[j + igap]))
1287             {
1288                 temp        = x[j];
1289                 x[j]        = x[j + igap];
1290                 x[j + igap] = temp;
1291                 if (*apply)
1292                 {
1293                     F77_FUNC(dswap, DSWAP) (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
1294                                                                               a_dim1 + 1], &c__1);
1295                 }
1296             }
1297             else
1298             {
1299                 goto L120;
1300             }
1301             j -= igap;
1302             goto L110;
1303 L120:
1304             ;
1305         }
1306         igap /= 2;
1307         goto L100;
1308     }
1309
1310 L9000:
1311     return;
1312
1313 }
1314
1315
1316
1317
1318 static void
1319 F77_FUNC(dsgets, DSGETS) (int *     ishift,
1320                           const char *    which,
1321                           int *     kev,
1322                           int *     np,
1323                           double *  ritz,
1324                           double *  bounds,
1325                           double *  shifts)
1326 {
1327     int c__1 = 1;
1328     int i__1, i__2;
1329     int kevd2;
1330
1331     --shifts;
1332     --bounds;
1333     --ritz;
1334
1335     if (!std::strncmp(which, "BE", 2))
1336     {
1337         i__1 = *kev + *np;
1338         F77_FUNC(dsortr, DSORTR) ("LA", &c__1, &i__1, &ritz[1], &bounds[1]);
1339         kevd2 = *kev / 2;
1340         if (*kev > 1)
1341         {
1342             i__1 = (kevd2 < *np) ? kevd2 : *np;
1343             i__2 = (kevd2 > *np) ? kevd2 : *np;
1344             F77_FUNC(dswap, DSWAP) (&i__1, &ritz[1], &c__1,
1345                                     &ritz[i__2 + 1], &c__1);
1346             i__1 = (kevd2 < *np) ? kevd2 : *np;
1347             i__2 = (kevd2 > *np) ? kevd2 : *np;
1348             F77_FUNC(dswap, DSWAP) (&i__1, &bounds[1], &c__1,
1349                                     &bounds[i__2 + 1], &c__1);
1350         }
1351
1352     }
1353     else
1354     {
1355         i__1 = *kev + *np;
1356         F77_FUNC(dsortr, DSORTR) (which, &c__1, &i__1, &ritz[1], &bounds[1]);
1357     }
1358
1359     if (*ishift == 1 && *np > 0)
1360     {
1361
1362         F77_FUNC(dsortr, DSORTR) ("SM", &c__1, np, &bounds[1], &ritz[1]);
1363         F77_FUNC(dcopy, DCOPY) (np, &ritz[1], &c__1, &shifts[1], &c__1);
1364     }
1365
1366
1367     return;
1368 }
1369
1370
1371
1372 static void
1373 F77_FUNC(dsconv, DSCONV) (int *     n,
1374                           double *  ritz,
1375                           double *  bounds,
1376                           double *  tol,
1377                           int *     nconv)
1378 {
1379     double c_b3 = 2/3.;
1380     int    i__1;
1381     double d__2, d__3;
1382
1383     int    i__;
1384     double eps23, temp;
1385
1386     --bounds;
1387     --ritz;
1388
1389     eps23 = GMX_DOUBLE_EPS;
1390     eps23 = std::pow(eps23, c_b3);
1391
1392     *nconv = 0;
1393     i__1   = *n;
1394     for (i__ = 1; i__ <= i__1; ++i__)
1395     {
1396
1397         d__2 = eps23;
1398         d__3 = std::abs(ritz[i__]);
1399         temp = (d__2 > d__3) ? d__2 : d__3;
1400         if (bounds[i__] <= *tol * temp)
1401         {
1402             ++(*nconv);
1403         }
1404     }
1405
1406     return;
1407 }
1408
1409
1410 static void
1411 F77_FUNC(dseigt, DSEIGT) (double *  rnorm,
1412                           int *     n,
1413                           double *  h__,
1414                           int *     ldh,
1415                           double *  eig,
1416                           double *  bounds,
1417                           double *  workl,
1418                           int *     ierr)
1419 {
1420     int c__1 = 1;
1421     int h_dim1, h_offset, i__1;
1422
1423     int k;
1424
1425
1426     --workl;
1427     --bounds;
1428     --eig;
1429     h_dim1   = *ldh;
1430     h_offset = 1 + h_dim1;
1431     h__     -= h_offset;
1432
1433     F77_FUNC(dcopy, DCOPY) (n, &h__[(h_dim1 << 1) + 1], &c__1, &eig[1], &c__1);
1434     i__1 = *n - 1;
1435     F77_FUNC(dcopy, DCOPY) (&i__1, &h__[h_dim1 + 2], &c__1, &workl[1], &c__1);
1436     F77_FUNC(dstqrb, DSTQRB) (n, &eig[1], &workl[1], &bounds[1], &workl[*n + 1], ierr);
1437     if (*ierr != 0)
1438     {
1439         goto L9000;
1440     }
1441
1442     i__1 = *n;
1443     for (k = 1; k <= i__1; ++k)
1444     {
1445         bounds[k] = *rnorm * std::abs(bounds[k]);
1446
1447     }
1448
1449
1450 L9000:
1451     return;
1452 }
1453
1454
1455
1456
1457 static void
1458 F77_FUNC(dsaitr, DSAITR) (int *     ido,
1459                           const char *    bmat,
1460                           int *     n,
1461                           int *     k,
1462                           int *     np,
1463                           int *     mode,
1464                           double *  resid,
1465                           double *  rnorm,
1466                           double *  v,
1467                           int *     ldv,
1468                           double *  h__,
1469                           int *     ldh,
1470                           int *     ipntr,
1471                           double *  workd,
1472                           int *     iwork,
1473                           int *     info)
1474 {
1475
1476     int    c__0  = 0;
1477     int    c__1  = 1;
1478     double c_b18 = 1.;
1479     double c_b42 = 0.;
1480     double c_b50 = -1.;
1481
1482     int    h_dim1, h_offset, v_dim1, v_offset, i__1;
1483     int    i__, jj;
1484     double temp1;
1485     int    infol;
1486     double safmin, minval;
1487
1488
1489     --workd;
1490     --resid;
1491     v_dim1   = *ldv;
1492     v_offset = 1 + v_dim1;
1493     v       -= v_offset;
1494     h_dim1   = *ldh;
1495     h_offset = 1 + h_dim1;
1496     h__     -= h_offset;
1497     --ipntr;
1498     --iwork;
1499     minval = GMX_DOUBLE_MIN;
1500     safmin = minval / GMX_DOUBLE_EPS;
1501
1502     if (*ido == 0)
1503     {
1504         *info    = 0;
1505         iwork[5] = 0;
1506         iwork[6] = 0;
1507         iwork[4] = 0;
1508         iwork[2] = 0;
1509         iwork[3] = 0;
1510
1511         iwork[12] = *k + 1;
1512
1513         iwork[8]  = 1;
1514         iwork[9]  = iwork[8] + *n;
1515         iwork[10] = iwork[9] + *n;
1516     }
1517
1518     if (iwork[5] == 1)
1519     {
1520         goto L50;
1521     }
1522     if (iwork[6] == 1)
1523     {
1524         goto L60;
1525     }
1526     if (iwork[2] == 1)
1527     {
1528         goto L70;
1529     }
1530     if (iwork[3] == 1)
1531     {
1532         goto L90;
1533     }
1534     if (iwork[4] == 1)
1535     {
1536         goto L30;
1537     }
1538
1539 L1000:
1540
1541
1542     if (*rnorm > 0.)
1543     {
1544         goto L40;
1545     }
1546
1547     iwork[11] = 1;
1548 L20:
1549     iwork[4] = 1;
1550     *ido     = 0;
1551 L30:
1552
1553     F77_FUNC(dgetv0, DGETV0) (ido, bmat, &iwork[11], &c__0, n, &iwork[12], &v[v_offset], ldv,
1554                               &resid[1], rnorm, &ipntr[1], &workd[1], &iwork[21], &iwork[7]);
1555     if (*ido != 99)
1556     {
1557         goto L9000;
1558     }
1559     if (iwork[7] < 0)
1560     {
1561         ++iwork[11];
1562         if (iwork[11] <= 3)
1563         {
1564             goto L20;
1565         }
1566
1567         *info = iwork[12] - 1;
1568         *ido  = 99;
1569         goto L9000;
1570     }
1571
1572 L40:
1573
1574     F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &v[iwork[12] * v_dim1 + 1], &c__1);
1575     if (*rnorm >= safmin)
1576     {
1577         temp1 = 1. / *rnorm;
1578         F77_FUNC(dscal, DSCAL) (n, &temp1, &v[iwork[12] * v_dim1 + 1], &c__1);
1579         F77_FUNC(dscal, DSCAL) (n, &temp1, &workd[iwork[8]], &c__1);
1580     }
1581     else
1582     {
1583
1584         F77_FUNC(dlascl, DLASCL) ("General", &i__, &i__, rnorm, &c_b18, n, &c__1, &v[iwork[12] *
1585                                                                                      v_dim1 + 1], n, &infol);
1586         F77_FUNC(dlascl, DLASCL) ("General", &i__, &i__, rnorm, &c_b18, n, &c__1, &workd[iwork[
1587                                                                                              8]], n, &infol);
1588     }
1589
1590     iwork[5] = 1;
1591     F77_FUNC(dcopy, DCOPY) (n, &v[iwork[12] * v_dim1 + 1], &c__1, &workd[iwork[10]], &c__1);
1592     ipntr[1] = iwork[10];
1593     ipntr[2] = iwork[9];
1594     ipntr[3] = iwork[8];
1595     *ido     = 1;
1596
1597     goto L9000;
1598 L50:
1599
1600
1601     iwork[5] = 0;
1602
1603     F77_FUNC(dcopy, DCOPY) (n, &workd[iwork[9]], &c__1, &resid[1], &c__1);
1604
1605     if (*mode == 2)
1606     {
1607         goto L65;
1608     }
1609     if (*bmat == 'G')
1610     {
1611         iwork[6] = 1;
1612         ipntr[1] = iwork[9];
1613         ipntr[2] = iwork[8];
1614         *ido     = 2;
1615
1616         goto L9000;
1617     }
1618     else if (*bmat == 'I')
1619     {
1620         F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
1621     }
1622 L60:
1623
1624     iwork[6] = 0;
1625
1626 L65:
1627     if (*mode == 2)
1628     {
1629
1630         workd[*n * 3 + 3] = F77_FUNC(ddot, DDOT) (n, &resid[1], &c__1, &workd[iwork[10]], &
1631                                                   c__1);
1632         workd[*n * 3 + 3] = std::sqrt(std::abs(workd[*n * 3 + 3]));
1633     }
1634     else if (*bmat == 'G')
1635     {
1636         workd[*n * 3 + 3] = F77_FUNC(ddot, DDOT) (n, &resid[1], &c__1, &workd[iwork[8]], &
1637                                                   c__1);
1638         workd[*n * 3 + 3] = std::sqrt(std::abs(workd[*n * 3 + 3]));
1639     }
1640     else if (*bmat == 'I')
1641     {
1642         workd[*n * 3 + 3] = F77_FUNC(dnrm2, DNRM2) (n, &resid[1], &c__1);
1643     }
1644
1645     if (*mode != 2)
1646     {
1647         F77_FUNC(dgemv, DGEMV) ("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[8]],
1648                                 &c__1, &c_b42, &workd[iwork[9]], &c__1);
1649     }
1650     else
1651     {
1652         F77_FUNC(dgemv, DGEMV) ("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[10]
1653                                 ], &c__1, &c_b42, &workd[iwork[9]], &c__1);
1654     }
1655
1656     F77_FUNC(dgemv, DGEMV) ("N", n, &iwork[12], &c_b50, &v[v_offset], ldv, &workd[iwork[9]], &
1657                             c__1, &c_b18, &resid[1], &c__1);
1658
1659     h__[iwork[12] + (h_dim1 << 1)] = workd[iwork[9] + iwork[12] - 1];
1660     if (iwork[12] == 1 || iwork[4] == 1)
1661     {
1662         h__[iwork[12] + h_dim1] = 0.;
1663     }
1664     else
1665     {
1666         h__[iwork[12] + h_dim1] = *rnorm;
1667     }
1668
1669     iwork[2] = 1;
1670     iwork[1] = 0;
1671
1672     if (*bmat == 'G')
1673     {
1674         F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[iwork[9]], &c__1);
1675         ipntr[1] = iwork[9];
1676         ipntr[2] = iwork[8];
1677         *ido     = 2;
1678
1679         goto L9000;
1680     }
1681     else if (*bmat == 'I')
1682     {
1683         F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
1684     }
1685 L70:
1686
1687     iwork[2] = 0;
1688
1689     if (*bmat == 'G')
1690     {
1691         *rnorm = F77_FUNC(ddot, DDOT) (n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
1692         *rnorm = std::sqrt(std::abs(*rnorm));
1693     }
1694     else if (*bmat == 'I')
1695     {
1696         *rnorm = F77_FUNC(dnrm2, DNRM2) (n, &resid[1], &c__1);
1697     }
1698
1699     if (*rnorm > workd[*n * 3 + 3] * .717F)
1700     {
1701         goto L100;
1702     }
1703
1704 L80:
1705
1706     F77_FUNC(dgemv, DGEMV) ("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[8]], &
1707                             c__1, &c_b42, &workd[iwork[9]], &c__1);
1708
1709     F77_FUNC(dgemv, DGEMV) ("N", n, &iwork[12], &c_b50, &v[v_offset], ldv, &workd[iwork[9]], &
1710                             c__1, &c_b18, &resid[1], &c__1);
1711
1712     if (iwork[12] == 1 || iwork[4] == 1)
1713     {
1714         h__[iwork[12] + h_dim1] = 0.;
1715     }
1716     h__[iwork[12] + (h_dim1 << 1)] += workd[iwork[9] + iwork[12] - 1];
1717
1718     iwork[3] = 1;
1719     if (*bmat == 'G')
1720     {
1721         F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[iwork[9]], &c__1);
1722         ipntr[1] = iwork[9];
1723         ipntr[2] = iwork[8];
1724         *ido     = 2;
1725
1726         goto L9000;
1727     }
1728     else if (*bmat == 'I')
1729     {
1730         F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
1731     }
1732 L90:
1733
1734
1735     if (*bmat == 'G')
1736     {
1737         workd[*n * 3 + 2] = F77_FUNC(ddot, DDOT) (n, &resid[1], &c__1, &workd[iwork[8]], &
1738                                                   c__1);
1739         workd[*n * 3 + 2] = std::sqrt(std::abs(workd[*n * 3 + 2]));
1740     }
1741     else if (*bmat == 'I')
1742     {
1743         workd[*n * 3 + 2] = F77_FUNC(dnrm2, DNRM2) (n, &resid[1], &c__1);
1744     }
1745
1746
1747     if (workd[*n * 3 + 2] > *rnorm * .717F)
1748     {
1749
1750         *rnorm = workd[*n * 3 + 2];
1751
1752     }
1753     else
1754     {
1755
1756         *rnorm = workd[*n * 3 + 2];
1757         ++iwork[1];
1758         if (iwork[1] <= 1)
1759         {
1760             goto L80;
1761         }
1762
1763         i__1 = *n;
1764         for (jj = 1; jj <= i__1; ++jj)
1765         {
1766             resid[jj] = 0.;
1767         }
1768         *rnorm = 0.;
1769     }
1770
1771 L100:
1772
1773     iwork[4] = 0;
1774     iwork[3] = 0;
1775
1776     if (h__[iwork[12] + h_dim1] < 0.)
1777     {
1778         h__[iwork[12] + h_dim1] = -h__[iwork[12] + h_dim1];
1779         if (iwork[12] < *k + *np)
1780         {
1781             F77_FUNC(dscal, DSCAL) (n, &c_b50, &v[(iwork[12] + 1) * v_dim1 + 1], &c__1);
1782         }
1783         else
1784         {
1785             F77_FUNC(dscal, DSCAL) (n, &c_b50, &resid[1], &c__1);
1786         }
1787     }
1788
1789     ++iwork[12];
1790     if (iwork[12] > *k + *np)
1791     {
1792         *ido = 99;
1793
1794
1795         goto L9000;
1796     }
1797
1798     goto L1000;
1799
1800 L9000:
1801     return;
1802 }
1803
1804
1805
1806
1807
1808
1809 static void
1810 F77_FUNC(dsaup2, DSAUP2) (int *     ido,
1811                           const char *    bmat,
1812                           int *     n,
1813                           const char *    which,
1814                           int *     nev,
1815                           int *     np,
1816                           double *  tol,
1817                           double *  resid,
1818                           int *     mode,
1819                           int gmx_unused *     iupd,
1820                           int *     ishift,
1821                           int *     mxiter,
1822                           double *  v,
1823                           int *     ldv,
1824                           double *  h__,
1825                           int *     ldh,
1826                           double *  ritz,
1827                           double *  bounds,
1828                           double *  q,
1829                           int *     ldq,
1830                           double *  workl,
1831                           int *     ipntr,
1832                           double *  workd,
1833                           int *     iwork,
1834                           int *     info)
1835 {
1836     double c_b3 = 2/3.;
1837     int    c__1 = 1;
1838     int    c__0 = 0;
1839
1840     int    h_dim1, h_offset, q_dim1, q_offset, v_dim1, v_offset, i__1, i__2,
1841            i__3;
1842     double d__2, d__3;
1843     int    j;
1844     double eps23;
1845     int    ierr;
1846     double temp;
1847     int    nevd2;
1848     int    nevm2;
1849     int    nevbef;
1850     char   wprime[2];
1851     int    nptemp;
1852
1853     --workd;
1854     --resid;
1855     --workl;
1856     --bounds;
1857     --ritz;
1858     v_dim1   = *ldv;
1859     v_offset = 1 + v_dim1;
1860     v       -= v_offset;
1861     h_dim1   = *ldh;
1862     h_offset = 1 + h_dim1;
1863     h__     -= h_offset;
1864     q_dim1   = *ldq;
1865     q_offset = 1 + q_dim1;
1866     q       -= q_offset;
1867     --ipntr;
1868     --iwork;
1869     eps23 = GMX_DOUBLE_EPS;
1870     eps23 = std::pow(eps23, c_b3);
1871
1872     if (*ido == 0)
1873     {
1874
1875         iwork[41] = 1;
1876         iwork[42] = 3;
1877         iwork[43] = 5;
1878         iwork[44] = 7;
1879
1880         iwork[9]  = *nev;
1881         iwork[10] = *np;
1882
1883         iwork[7] = iwork[9] + iwork[10];
1884         iwork[8] = 0;
1885         iwork[6] = 0;
1886
1887         iwork[2] = 1;
1888         iwork[4] = 0;
1889         iwork[5] = 0;
1890         iwork[1] = 0;
1891
1892         if (*info != 0)
1893         {
1894
1895             iwork[3] = 1;
1896             *info    = 0;
1897         }
1898         else
1899         {
1900             iwork[3] = 0;
1901         }
1902     }
1903
1904     if (iwork[2] == 1)
1905     {
1906         F77_FUNC(dgetv0, DGETV0) (ido, bmat, &c__1, &iwork[3], n, &c__1, &v[v_offset], ldv, &
1907                                   resid[1], &workd[*n * 3 + 1], &ipntr[1], &workd[1], &iwork[41],
1908                                   info);
1909
1910         if (*ido != 99)
1911         {
1912             goto L9000;
1913         }
1914
1915         if (workd[*n * 3 + 1] == 0.)
1916         {
1917
1918             *info = -9;
1919             goto L1200;
1920         }
1921         iwork[2] = 0;
1922         *ido     = 0;
1923     }
1924
1925     if (iwork[4] == 1)
1926     {
1927         goto L20;
1928     }
1929
1930     if (iwork[5] == 1)
1931     {
1932         goto L50;
1933     }
1934
1935     if (iwork[1] == 1)
1936     {
1937         goto L100;
1938     }
1939
1940     F77_FUNC(dsaitr, DSAITR) (ido, bmat, n, &c__0, &iwork[9], mode, &resid[1], &workd[*n * 3 +
1941                                                                                       1], &v[v_offset], ldv, &h__[h_offset], ldh, &ipntr[1], &workd[1],
1942                               &iwork[21], info);
1943
1944     if (*ido != 99)
1945     {
1946         goto L9000;
1947     }
1948
1949     if (*info > 0)
1950     {
1951
1952         *np     = *info;
1953         *mxiter = iwork[6];
1954         *info   = -9999;
1955         goto L1200;
1956     }
1957
1958 L1000:
1959
1960     ++iwork[6];
1961
1962
1963     *ido = 0;
1964 L20:
1965     iwork[4] = 1;
1966
1967     F77_FUNC(dsaitr, DSAITR) (ido, bmat, n, nev, np, mode, &resid[1], &workd[*n * 3 + 1],
1968                               &v[v_offset], ldv, &h__[h_offset], ldh, &ipntr[1], &workd[1],
1969                               &iwork[21], info);
1970
1971     if (*ido != 99)
1972     {
1973         goto L9000;
1974     }
1975
1976     if (*info > 0)
1977     {
1978
1979         *np     = *info;
1980         *mxiter = iwork[6];
1981         *info   = -9999;
1982         goto L1200;
1983     }
1984     iwork[4] = 0;
1985
1986     F77_FUNC(dseigt, DSEIGT) (&workd[*n * 3 + 1], &iwork[7], &h__[h_offset], ldh, &ritz[1], &
1987                               bounds[1], &workl[1], &ierr);
1988
1989     if (ierr != 0)
1990     {
1991         *info = -8;
1992         goto L1200;
1993     }
1994
1995     F77_FUNC(dcopy, DCOPY) (&iwork[7], &ritz[1], &c__1, &workl[iwork[7] + 1], &c__1);
1996     F77_FUNC(dcopy, DCOPY) (&iwork[7], &bounds[1], &c__1, &workl[(iwork[7] << 1) + 1], &c__1);
1997
1998     *nev = iwork[9];
1999     *np  = iwork[10];
2000     F77_FUNC(dsgets, DSGETS) (ishift, which, nev, np, &ritz[1], &bounds[1], &workl[1]);
2001
2002     F77_FUNC(dcopy, DCOPY) (nev, &bounds[*np + 1], &c__1, &workl[*np + 1], &c__1);
2003     F77_FUNC(dsconv, DSCONV) (nev, &ritz[*np + 1], &workl[*np + 1], tol, &iwork[8]);
2004
2005     nptemp = *np;
2006     i__1   = nptemp;
2007     for (j = 1; j <= i__1; ++j)
2008     {
2009         if (bounds[j] == 0.)
2010         {
2011             --(*np);
2012             ++(*nev);
2013         }
2014     }
2015
2016     if (iwork[8] >= iwork[9] || iwork[6] > *mxiter || *np == 0)
2017     {
2018
2019         if (!std::strncmp(which, "BE", 2))
2020         {
2021
2022             std::strncpy(wprime, "SA", 2);
2023             F77_FUNC(dsortr, DSORTR) (wprime, &c__1, &iwork[7], &ritz[1], &bounds[1]);
2024             nevd2 = *nev / 2;
2025             nevm2 = *nev - nevd2;
2026             if (*nev > 1)
2027             {
2028                 i__1 = (nevd2 < *np) ? nevd2 : *np;
2029                 i__2 = iwork[7] - nevd2 + 1, i__3 = iwork[7] - *np + 1;
2030                 F77_FUNC(dswap, DSWAP) (&i__1, &ritz[nevm2 + 1], &c__1,
2031                                         &ritz[((i__2 > i__3) ? i__2 : i__3)],
2032                                         &c__1);
2033                 i__1 = (nevd2 < *np) ? nevd2 : *np;
2034                 i__2 = iwork[7] - nevd2 + 1, i__3 = iwork[7] - *np;
2035                 F77_FUNC(dswap, DSWAP) (&i__1, &bounds[nevm2 + 1], &c__1,
2036                                         &bounds[((i__2 > i__3) ? i__2 : i__3) + 1],
2037                                         &c__1);
2038             }
2039
2040         }
2041         else
2042         {
2043
2044             if (!std::strncmp(which, "LM", 2))
2045             {
2046                 std::strncpy(wprime, "SM", 2);
2047             }
2048             if (!std::strncmp(which, "SM", 2))
2049             {
2050                 std::strncpy(wprime, "LM", 2);
2051             }
2052             if (!std::strncmp(which, "LA", 2))
2053             {
2054                 std::strncpy(wprime, "SA", 2);
2055             }
2056             if (!std::strncmp(which, "SA", 2))
2057             {
2058                 std::strncpy(wprime, "LA", 2);
2059             }
2060
2061             F77_FUNC(dsortr, DSORTR) (wprime, &c__1, &iwork[7], &ritz[1], &bounds[1]);
2062
2063         }
2064
2065         i__1 = iwork[9];
2066         for (j = 1; j <= i__1; ++j)
2067         {
2068             d__2       = eps23;
2069             d__3       = std::abs(ritz[j]);
2070             temp       = (d__2 > d__3) ? d__2 : d__3;
2071             bounds[j] /= temp;
2072         }
2073
2074         std::strncpy(wprime, "LA", 2);
2075         F77_FUNC(dsortr, DSORTR) (wprime, &c__1, &iwork[9], &bounds[1], &ritz[1]);
2076
2077         i__1 = iwork[9];
2078         for (j = 1; j <= i__1; ++j)
2079         {
2080             d__2       = eps23;
2081             d__3       = std::abs(ritz[j]);
2082             temp       = (d__2 > d__3) ? d__2 : d__3;
2083             bounds[j] *= temp;
2084         }
2085
2086         if (!std::strncmp(which, "BE", 2))
2087         {
2088
2089             std::strncpy(wprime, "LA", 2);
2090             F77_FUNC(dsortr, DSORTR) (wprime, &c__1, &iwork[8], &ritz[1], &bounds[1]);
2091
2092         }
2093         else
2094         {
2095             F77_FUNC(dsortr, DSORTR) (which, &c__1, &iwork[8], &ritz[1], &bounds[1]);
2096
2097         }
2098
2099         h__[h_dim1 + 1] = workd[*n * 3 + 1];
2100
2101
2102         if (iwork[6] > *mxiter && iwork[8] < *nev)
2103         {
2104             *info = 1;
2105         }
2106         if (*np == 0 && iwork[8] < iwork[9])
2107         {
2108             *info = 2;
2109         }
2110
2111         *np = iwork[8];
2112         goto L1100;
2113
2114     }
2115     else if (iwork[8] < *nev && *ishift == 1)
2116     {
2117         nevbef = *nev;
2118         i__1   = iwork[8], i__2 = *np / 2;
2119         *nev  += (i__1 < i__2) ? i__1 : i__2;
2120         if (*nev == 1 && iwork[7] >= 6)
2121         {
2122             *nev = iwork[7] / 2;
2123         }
2124         else if (*nev == 1 && iwork[7] > 2)
2125         {
2126             *nev = 2;
2127         }
2128         *np = iwork[7] - *nev;
2129
2130
2131         if (nevbef < *nev)
2132         {
2133             F77_FUNC(dsgets, DSGETS) (ishift, which, nev, np, &ritz[1], &bounds[1], &workl[1]);
2134         }
2135
2136     }
2137
2138
2139     if (*ishift == 0)
2140     {
2141
2142         iwork[5] = 1;
2143         *ido     = 3;
2144         goto L9000;
2145     }
2146
2147 L50:
2148
2149     iwork[5] = 0;
2150
2151     if (*ishift == 0)
2152     {
2153         F77_FUNC(dcopy, DCOPY) (np, &workl[1], &c__1, &ritz[1], &c__1);
2154     }
2155
2156     F77_FUNC(dsapps, DSAPPS) (n, nev, np, &ritz[1], &v[v_offset], ldv, &h__[h_offset], ldh, &
2157                               resid[1], &q[q_offset], ldq, &workd[1]);
2158
2159     iwork[1] = 1;
2160     if (*bmat == 'G')
2161     {
2162         F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[*n + 1], &c__1);
2163         ipntr[1] = *n + 1;
2164         ipntr[2] = 1;
2165         *ido     = 2;
2166
2167         goto L9000;
2168     }
2169     else if (*bmat == 'I')
2170     {
2171         F77_FUNC(dcopy, DCOPY) (n, &resid[1], &c__1, &workd[1], &c__1);
2172     }
2173
2174 L100:
2175
2176     if (*bmat == 'G')
2177     {
2178         workd[*n * 3 + 1] = F77_FUNC(ddot, DDOT) (n, &resid[1], &c__1, &workd[1], &c__1);
2179         workd[*n * 3 + 1] = std::sqrt(std::abs(workd[*n * 3 + 1]));
2180     }
2181     else if (*bmat == 'I')
2182     {
2183         workd[*n * 3 + 1] = F77_FUNC(dnrm2, DNRM2) (n, &resid[1], &c__1);
2184     }
2185     iwork[1] = 0;
2186
2187     goto L1000;
2188
2189 L1100:
2190
2191     *mxiter = iwork[6];
2192     *nev    = iwork[8];
2193
2194 L1200:
2195     *ido = 99;
2196
2197 L9000:
2198     return;
2199
2200 }
2201
2202
2203
2204 void
2205 F77_FUNC(dsaupd, DSAUPD) (int *     ido,
2206                           const char *    bmat,
2207                           int *     n,
2208                           const char *      which,
2209                           int *     nev,
2210                           double *  tol,
2211                           double *  resid,
2212                           int *     ncv,
2213                           double *  v,
2214                           int *     ldv,
2215                           int *     iparam,
2216                           int *     ipntr,
2217                           double *  workd,
2218                           int *     iwork,
2219                           double *  workl,
2220                           int *     lworkl,
2221                           int *     info)
2222 {
2223     int v_dim1, v_offset, i__1, i__2;
2224     int j;
2225
2226     --workd;
2227     --resid;
2228     v_dim1   = *ldv;
2229     v_offset = 1 + v_dim1;
2230     v       -= v_offset;
2231     --iparam;
2232     --ipntr;
2233     --iwork;
2234     --workl;
2235
2236     if (*ido == 0)
2237     {
2238
2239
2240         iwork[2]  = 0;
2241         iwork[5]  = iparam[1];
2242         iwork[10] = iparam[3];
2243         iwork[12] = iparam[4];
2244
2245         iwork[6]  = 1;
2246         iwork[11] = iparam[7];
2247
2248
2249         if (*n <= 0)
2250         {
2251             iwork[2] = -1;
2252         }
2253         else if (*nev <= 0)
2254         {
2255             iwork[2] = -2;
2256         }
2257         else if (*ncv <= *nev || *ncv > *n)
2258         {
2259             iwork[2] = -3;
2260         }
2261
2262
2263         iwork[15] = *ncv - *nev;
2264
2265         if (iwork[10] <= 0)
2266         {
2267             iwork[2] = -4;
2268         }
2269         if (std::strncmp(which, "LM", 2) && std::strncmp(which, "SM", 2) &&
2270             std::strncmp(which, "LA", 2) && std::strncmp(which, "SA", 2) &&
2271             std::strncmp(which, "BE", 2))
2272         {
2273             iwork[2] = -5;
2274         }
2275         if (*bmat != 'I' && *bmat != 'G')
2276         {
2277             iwork[2] = -6;
2278         }
2279
2280         i__1 = *ncv;
2281         if (*lworkl < i__1 * i__1 + (*ncv << 3))
2282         {
2283             iwork[2] = -7;
2284         }
2285         if (iwork[11] < 1 || iwork[11] > 5)
2286         {
2287             iwork[2] = -10;
2288         }
2289         else if (iwork[11] == 1 && *bmat == 'G')
2290         {
2291             iwork[2] = -11;
2292         }
2293         else if (iwork[5] < 0 || iwork[5] > 1)
2294         {
2295             iwork[2] = -12;
2296         }
2297         else if (*nev == 1 && !std::strncmp(which, "BE", 2))
2298         {
2299             iwork[2] = -13;
2300         }
2301
2302         if (iwork[2] != 0)
2303         {
2304             *info = iwork[2];
2305             *ido  = 99;
2306             goto L9000;
2307         }
2308
2309         if (iwork[12] <= 0)
2310         {
2311             iwork[12] = 1;
2312         }
2313         if (*tol <= 0.)
2314         {
2315             *tol = GMX_DOUBLE_EPS;
2316         }
2317
2318         iwork[15] = *ncv - *nev;
2319         iwork[13] = *nev;
2320         i__2      = *ncv;
2321         i__1      = i__2 * i__2 + (*ncv << 3);
2322         for (j = 1; j <= i__1; ++j)
2323         {
2324             workl[j] = 0.;
2325         }
2326
2327         iwork[8]  = *ncv;
2328         iwork[9]  = *ncv;
2329         iwork[3]  = 1;
2330         iwork[16] = iwork[3] + (iwork[8] << 1);
2331         iwork[1]  = iwork[16] + *ncv;
2332         iwork[4]  = iwork[1] + *ncv;
2333         i__1      = *ncv;
2334         iwork[7]  = iwork[4] + i__1 * i__1;
2335         iwork[14] = iwork[7] + *ncv * 3;
2336
2337         ipntr[4]  = iwork[14];
2338         ipntr[5]  = iwork[3];
2339         ipntr[6]  = iwork[16];
2340         ipntr[7]  = iwork[1];
2341         ipntr[11] = iwork[7];
2342     }
2343
2344     F77_FUNC(dsaup2, DSAUP2) (ido, bmat, n, which, &iwork[13], &iwork[15], tol, &resid[1], &
2345                               iwork[11], &iwork[6], &iwork[5], &iwork[10], &v[v_offset], ldv, &
2346                               workl[iwork[3]], &iwork[8], &workl[iwork[16]], &workl[iwork[1]], &
2347                               workl[iwork[4]], &iwork[9], &workl[iwork[7]], &ipntr[1], &workd[1],
2348                               &iwork[21], info);
2349
2350     if (*ido == 3)
2351     {
2352         iparam[8] = iwork[15];
2353     }
2354     if (*ido != 99)
2355     {
2356         goto L9000;
2357     }
2358
2359     iparam[3] = iwork[10];
2360     iparam[5] = iwork[15];
2361
2362     if (*info < 0)
2363     {
2364         goto L9000;
2365     }
2366     if (*info == 2)
2367     {
2368         *info = 3;
2369     }
2370
2371 L9000:
2372
2373     return;
2374
2375 }
2376
2377
2378
2379 void
2380 F77_FUNC(dseupd, DSEUPD) (int *     rvec,
2381                           const char *    howmny,
2382                           int *     select,
2383                           double *  d__,
2384                           double *  z__,
2385                           int *     ldz,
2386                           double *  sigma,
2387                           const char *    bmat,
2388                           int *     n,
2389                           const char *    which,
2390                           int *     nev,
2391                           double *  tol,
2392                           double *  resid,
2393                           int *     ncv,
2394                           double *  v,
2395                           int *     ldv,
2396                           int *     iparam,
2397                           int *     ipntr,
2398                           double *  workd,
2399                           double *  workl,
2400                           int *     lworkl,
2401                           int *     info)
2402 {
2403     double c_b21  = 2/3.;
2404     int    c__1   = 1;
2405     double c_b102 = 1.;
2406     int    v_dim1, v_offset, z_dim1, z_offset, i__1;
2407     double d__1, d__2, d__3;
2408
2409     int    j, k, ih, iq, iw, ibd, ihb, ihd, ldh, ilg, ldq, ism, irz;
2410     int    mode;
2411     double eps23;
2412     int    ierr;
2413     double temp;
2414     int    next;
2415     char   type__[6];
2416     int    ritz;
2417     int    reord;
2418     int    nconv;
2419     double rnorm;
2420     double bnorm2;
2421     double thres1 = 0, thres2 = 0;
2422     int    bounds;
2423     int    ktrord;
2424     double tempbnd;
2425     int    leftptr, rghtptr;
2426
2427
2428     --workd;
2429     --resid;
2430     z_dim1   = *ldz;
2431     z_offset = 1 + z_dim1;
2432     z__     -= z_offset;
2433     --d__;
2434     --select;
2435     v_dim1   = *ldv;
2436     v_offset = 1 + v_dim1;
2437     v       -= v_offset;
2438     --iparam;
2439     --ipntr;
2440     --workl;
2441
2442     mode  = iparam[7];
2443     nconv = iparam[5];
2444     *info = 0;
2445
2446     if (nconv == 0)
2447     {
2448         goto L9000;
2449     }
2450     ierr = 0;
2451
2452     if (nconv <= 0)
2453     {
2454         ierr = -14;
2455     }
2456     if (*n <= 0)
2457     {
2458         ierr = -1;
2459     }
2460     if (*nev <= 0)
2461     {
2462         ierr = -2;
2463     }
2464     if (*ncv <= *nev || *ncv > *n)
2465     {
2466         ierr = -3;
2467     }
2468     if (std::strncmp(which, "LM", 2) && std::strncmp(which, "SM", 2) &&
2469         std::strncmp(which, "LA", 2) && std::strncmp(which, "SA", 2) &&
2470         std::strncmp(which, "BE", 2))
2471     {
2472         ierr = -5;
2473     }
2474     if (*bmat != 'I' && *bmat != 'G')
2475     {
2476         ierr = -6;
2477     }
2478     if (*howmny != 'A' && *howmny != 'P' &&
2479         *howmny != 'S' && *rvec)
2480     {
2481         ierr = -15;
2482     }
2483     if (*rvec && *howmny == 'S')
2484     {
2485         ierr = -16;
2486     }
2487     i__1 = *ncv;
2488     if (*rvec && *lworkl < i__1 * i__1 + (*ncv << 3))
2489     {
2490         ierr = -7;
2491     }
2492
2493     if (mode == 1 || mode == 2)
2494     {
2495         std::strncpy(type__, "REGULR", 6);
2496     }
2497     else if (mode == 3)
2498     {
2499         std::strncpy(type__, "SHIFTI", 6);
2500     }
2501     else if (mode == 4)
2502     {
2503         std::strncpy(type__, "BUCKLE", 6);
2504     }
2505     else if (mode == 5)
2506     {
2507         std::strncpy(type__, "CAYLEY", 6);
2508     }
2509     else
2510     {
2511         ierr = -10;
2512     }
2513     if (mode == 1 && *bmat == 'G')
2514     {
2515         ierr = -11;
2516     }
2517     if (*nev == 1 && !std::strncmp(which, "BE", 2))
2518     {
2519         ierr = -12;
2520     }
2521
2522     if (ierr != 0)
2523     {
2524         *info = ierr;
2525         goto L9000;
2526     }
2527
2528     ih        = ipntr[5];
2529     ritz      = ipntr[6];
2530     bounds    = ipntr[7];
2531     ldh       = *ncv;
2532     ldq       = *ncv;
2533     ihd       = bounds + ldh;
2534     ihb       = ihd + ldh;
2535     iq        = ihb + ldh;
2536     iw        = iq + ldh * *ncv;
2537     next      = iw + (*ncv << 1);
2538     ipntr[4]  = next;
2539     ipntr[8]  = ihd;
2540     ipntr[9]  = ihb;
2541     ipntr[10] = iq;
2542
2543     irz = ipntr[11] + *ncv;
2544     ibd = irz + *ncv;
2545
2546
2547     eps23 = GMX_DOUBLE_EPS;
2548     eps23 = std::pow(eps23, c_b21);
2549
2550     rnorm = workl[ih];
2551     if (*bmat == 'I')
2552     {
2553         bnorm2 = rnorm;
2554     }
2555     else if (*bmat == 'G')
2556     {
2557         bnorm2 = F77_FUNC(dnrm2, DNRM2) (n, &workd[1], &c__1);
2558     }
2559
2560     if (*rvec)
2561     {
2562
2563         if (!std::strncmp(which, "LM", 2) || !std::strncmp(which, "SM", 2) ||
2564             !std::strncmp(which, "LA", 2) || !std::strncmp(which, "SA", 2))
2565         {
2566
2567         }
2568         else if (!std::strncmp(which, "BE", 2))
2569         {
2570
2571
2572             ism    = (*nev > nconv) ? *nev : nconv;
2573             ism   /= 2;
2574             ilg    = ism + 1;
2575             thres1 = workl[ism];
2576             thres2 = workl[ilg];
2577
2578
2579         }
2580
2581         reord  = 0;
2582         ktrord = 0;
2583         i__1   = *ncv - 1;
2584         for (j = 0; j <= i__1; ++j)
2585         {
2586             select[j + 1] = 0;
2587             if (!std::strncmp(which, "LM", 2))
2588             {
2589                 if (std::abs(workl[irz + j]) >= std::abs(thres1))
2590                 {
2591                     d__2    = eps23;
2592                     d__3    = std::abs(workl[irz + j]);
2593                     tempbnd = (d__2 > d__3) ? d__2 : d__3;
2594                     if (workl[ibd + j] <= *tol * tempbnd)
2595                     {
2596                         select[j + 1] = 1;
2597                     }
2598                 }
2599             }
2600             else if (!std::strncmp(which, "SM", 2))
2601             {
2602                 if (std::abs(workl[irz + j]) <= std::abs(thres1))
2603                 {
2604                     d__2    = eps23;
2605                     d__3    = std::abs(workl[irz + j]);
2606                     tempbnd = (d__2 > d__3) ? d__2 : d__3;
2607                     if (workl[ibd + j] <= *tol * tempbnd)
2608                     {
2609                         select[j + 1] = 1;
2610                     }
2611                 }
2612             }
2613             else if (!std::strncmp(which, "LA", 2))
2614             {
2615                 if (workl[irz + j] >= thres1)
2616                 {
2617                     d__2    = eps23;
2618                     d__3    = std::abs(workl[irz + j]);
2619                     tempbnd = (d__2 > d__3) ? d__2 : d__3;
2620                     if (workl[ibd + j] <= *tol * tempbnd)
2621                     {
2622                         select[j + 1] = 1;
2623                     }
2624                 }
2625             }
2626             else if (!std::strncmp(which, "SA", 2))
2627             {
2628                 if (workl[irz + j] <= thres1)
2629                 {
2630                     d__2    = eps23;
2631                     d__3    = std::abs(workl[irz + j]);
2632                     tempbnd = (d__2 > d__3) ? d__2 : d__3;
2633                     if (workl[ibd + j] <= *tol * tempbnd)
2634                     {
2635                         select[j + 1] = 1;
2636                     }
2637                 }
2638             }
2639             else if (!std::strncmp(which, "BE", 2))
2640             {
2641                 if (workl[irz + j] <= thres1 || workl[irz + j] >= thres2)
2642                 {
2643                     d__2    = eps23;
2644                     d__3    = std::abs(workl[irz + j]);
2645                     tempbnd = (d__2 > d__3) ? d__2 : d__3;
2646                     if (workl[ibd + j] <= *tol * tempbnd)
2647                     {
2648                         select[j + 1] = 1;
2649                     }
2650                 }
2651             }
2652             if (j + 1 > nconv)
2653             {
2654                 reord = select[j + 1] || reord;
2655             }
2656             if (select[j + 1])
2657             {
2658                 ++ktrord;
2659             }
2660         }
2661
2662         i__1 = *ncv - 1;
2663         F77_FUNC(dcopy, DCOPY) (&i__1, &workl[ih + 1], &c__1, &workl[ihb], &c__1);
2664         F77_FUNC(dcopy, DCOPY) (ncv, &workl[ih + ldh], &c__1, &workl[ihd], &c__1);
2665
2666         F77_FUNC(dsteqr, DSTEQR) ("Identity", ncv, &workl[ihd], &workl[ihb], &workl[iq], &ldq, &
2667                                   workl[iw], &ierr);
2668
2669         if (ierr != 0)
2670         {
2671             *info = -8;
2672             goto L9000;
2673         }
2674
2675
2676         if (reord)
2677         {
2678
2679             leftptr = 1;
2680             rghtptr = *ncv;
2681
2682             if (*ncv == 1)
2683             {
2684                 goto L30;
2685             }
2686
2687 L20:
2688             if (select[leftptr])
2689             {
2690
2691                 ++leftptr;
2692
2693             }
2694             else if (!select[rghtptr])
2695             {
2696
2697                 --rghtptr;
2698
2699             }
2700             else
2701             {
2702
2703                 temp                     = workl[ihd + leftptr - 1];
2704                 workl[ihd + leftptr - 1] = workl[ihd + rghtptr - 1];
2705                 workl[ihd + rghtptr - 1] = temp;
2706                 F77_FUNC(dcopy, DCOPY) (ncv, &workl[iq + *ncv * (leftptr - 1)], &c__1, &workl[
2707                                             iw], &c__1);
2708                 F77_FUNC(dcopy, DCOPY) (ncv, &workl[iq + *ncv * (rghtptr - 1)], &c__1, &workl[
2709                                             iq + *ncv * (leftptr - 1)], &c__1);
2710                 F77_FUNC(dcopy, DCOPY) (ncv, &workl[iw], &c__1, &workl[iq + *ncv * (rghtptr -
2711                                                                                     1)], &c__1);
2712                 ++leftptr;
2713                 --rghtptr;
2714
2715             }
2716
2717             if (leftptr < rghtptr)
2718             {
2719                 goto L20;
2720             }
2721
2722 L30:
2723             ;
2724         }
2725
2726         F77_FUNC(dcopy, DCOPY) (&nconv, &workl[ihd], &c__1, &d__[1], &c__1);
2727
2728     }
2729     else
2730     {
2731
2732         F77_FUNC(dcopy, DCOPY) (&nconv, &workl[ritz], &c__1, &d__[1], &c__1);
2733         F77_FUNC(dcopy, DCOPY) (ncv, &workl[ritz], &c__1, &workl[ihd], &c__1);
2734
2735     }
2736     if (!std::strncmp(type__, "REGULR", 6))
2737     {
2738
2739         if (*rvec)
2740         {
2741             F77_FUNC(dsesrt, DSESRT) ("LA", rvec, &nconv, &d__[1], ncv, &workl[iq], &ldq);
2742         }
2743         else
2744         {
2745             F77_FUNC(dcopy, DCOPY) (ncv, &workl[bounds], &c__1, &workl[ihb], &c__1);
2746         }
2747
2748     }
2749     else
2750     {
2751
2752         F77_FUNC(dcopy, DCOPY) (ncv, &workl[ihd], &c__1, &workl[iw], &c__1);
2753         if (!std::strncmp(type__, "SHIFTI", 6))
2754         {
2755             i__1 = *ncv;
2756             for (k = 1; k <= i__1; ++k)
2757             {
2758                 workl[ihd + k - 1] = 1. / workl[ihd + k - 1] + *sigma;
2759             }
2760         }
2761         else if (!std::strncmp(type__, "BUCKLE", 6))
2762         {
2763             i__1 = *ncv;
2764             for (k = 1; k <= i__1; ++k)
2765             {
2766                 workl[ihd + k - 1] = *sigma * workl[ihd + k - 1] / (workl[ihd
2767                                                                           + k - 1] - 1.);
2768             }
2769         }
2770         else if (!std::strncmp(type__, "CAYLEY", 6))
2771         {
2772             i__1 = *ncv;
2773             for (k = 1; k <= i__1; ++k)
2774             {
2775                 workl[ihd + k - 1] = *sigma * (workl[ihd + k - 1] + 1.) / (
2776                         workl[ihd + k - 1] - 1.);
2777             }
2778         }
2779
2780         F77_FUNC(dcopy, DCOPY) (&nconv, &workl[ihd], &c__1, &d__[1], &c__1);
2781         F77_FUNC(dsortr, DSORTR) ("LA", &c__1, &nconv, &workl[ihd], &workl[iw]);
2782         if (*rvec)
2783         {
2784             F77_FUNC(dsesrt, DSESRT) ("LA", rvec, &nconv, &d__[1], ncv, &workl[iq], &ldq);
2785         }
2786         else
2787         {
2788             F77_FUNC(dcopy, DCOPY) (ncv, &workl[bounds], &c__1, &workl[ihb], &c__1);
2789             d__1 = bnorm2 / rnorm;
2790             F77_FUNC(dscal, DSCAL) (ncv, &d__1, &workl[ihb], &c__1);
2791             F77_FUNC(dsortr, DSORTR) ("LA", &c__1, &nconv, &d__[1], &workl[ihb]);
2792         }
2793
2794     }
2795
2796     if (*rvec && *howmny == 'A')
2797     {
2798
2799         F77_FUNC(dgeqr2, DGEQR2) (ncv, &nconv, &workl[iq], &ldq, &workl[iw + *ncv], &workl[ihb],
2800                                   &ierr);
2801
2802         F77_FUNC(dorm2r, DORM2R) ("Right", "Notranspose", n, ncv, &nconv, &workl[iq], &ldq, &
2803                                   workl[iw + *ncv], &v[v_offset], ldv, &workd[*n + 1], &ierr);
2804         F77_FUNC(dlacpy, DLACPY) ("All", n, &nconv, &v[v_offset], ldv, &z__[z_offset], ldz);
2805
2806         i__1 = *ncv - 1;
2807         for (j = 1; j <= i__1; ++j)
2808         {
2809             workl[ihb + j - 1] = 0.;
2810         }
2811         workl[ihb + *ncv - 1] = 1.;
2812         F77_FUNC(dorm2r, DORM2R) ("Left", "Transpose", ncv, &c__1, &nconv, &workl[iq], &ldq, &
2813                                   workl[iw + *ncv], &workl[ihb], ncv, &temp, &ierr);
2814
2815     }
2816     else if (*rvec && *howmny == 'S')
2817     {
2818
2819     }
2820
2821     if (!std::strncmp(type__, "REGULR", 6) && *rvec)
2822     {
2823
2824         i__1 = *ncv;
2825         for (j = 1; j <= i__1; ++j)
2826         {
2827             workl[ihb + j - 1] = rnorm * std::abs(workl[ihb + j - 1]);
2828         }
2829
2830     }
2831     else if (std::strncmp(type__, "REGULR", 6) && *rvec)
2832     {
2833
2834         F77_FUNC(dscal, DSCAL) (ncv, &bnorm2, &workl[ihb], &c__1);
2835         if (!std::strncmp(type__, "SHIFTI", 6))
2836         {
2837
2838             i__1 = *ncv;
2839             for (k = 1; k <= i__1; ++k)
2840             {
2841                 d__2               = workl[iw + k - 1];
2842                 workl[ihb + k - 1] = std::abs(workl[ihb + k - 1])/(d__2 * d__2);
2843             }
2844
2845         }
2846         else if (!std::strncmp(type__, "BUCKLE", 6))
2847         {
2848
2849             i__1 = *ncv;
2850             for (k = 1; k <= i__1; ++k)
2851             {
2852                 d__2               = workl[iw + k - 1] - 1.;
2853                 workl[ihb + k - 1] = *sigma * std::abs(workl[ihb + k - 1])/(d__2 * d__2);
2854             }
2855
2856         }
2857         else if (!std::strncmp(type__, "CAYLEY", 6))
2858         {
2859
2860             i__1 = *ncv;
2861             for (k = 1; k <= i__1; ++k)
2862             {
2863                 workl[ihb + k - 1] = std::abs(workl[ihb + k - 1] / workl[iw + k - 1] * (workl[iw + k - 1] - 1.));
2864
2865             }
2866
2867         }
2868
2869     }
2870
2871     if (*rvec && (!std::strncmp(type__, "SHIFTI", 6) || !std::strncmp(type__, "CAYLEY", 6)))
2872     {
2873
2874         i__1 = nconv - 1;
2875         for (k = 0; k <= i__1; ++k)
2876         {
2877             workl[iw + k] = workl[iq + k * ldq + *ncv - 1] / workl[iw + k];
2878         }
2879
2880     }
2881     else if (*rvec && !std::strncmp(type__, "BUCKLE", 6))
2882     {
2883
2884         i__1 = nconv - 1;
2885         for (k = 0; k <= i__1; ++k)
2886         {
2887             workl[iw + k] = workl[iq + k * ldq + *ncv - 1] / (workl[iw + k] -
2888                                                               1.);
2889         }
2890
2891     }
2892
2893     if (std::strncmp(type__, "REGULR", 6))
2894     {
2895         F77_FUNC(dger, DGER) (n, &nconv, &c_b102, &resid[1], &c__1, &workl[iw], &c__1, &z__[
2896                                   z_offset], ldz);
2897     }
2898
2899 L9000:
2900
2901     return;
2902
2903 }
2904
2905
2906
2907
2908
2909 /* Selected single precision arpack routines */
2910
2911
2912 static void
2913 F77_FUNC(sstqrb, SSTQRB) (int *      n,
2914                           float *   d__,
2915                           float *   e,
2916                           float *   z__,
2917                           float *   work,
2918                           int *      info)
2919 {
2920     int   i__1, i__2;
2921     float d__1, d__2;
2922     int   c__0  = 0;
2923     int   c__1  = 1;
2924     float c_b31 = 1.;
2925
2926     float b, c__, f, g;
2927     int   i__, j, k, l, m;
2928     float p, r__, s;
2929     int   l1, ii, mm, lm1, mm1, nm1;
2930     float rt1, rt2, eps;
2931     int   lsv;
2932     float tst, eps2;
2933     int   lend, jtot, lendm1, lendp1, iscale;
2934
2935     int   lendsv, nmaxit, icompz;
2936     float ssfmax, ssfmin, safmin, minval, safmax, anorm;
2937
2938
2939     --work;
2940     --z__;
2941     --e;
2942     --d__;
2943
2944     *info = 0;
2945
2946     icompz = 2;
2947
2948     if (*n == 0)
2949     {
2950         return;
2951     }
2952
2953     if (*n == 1)
2954     {
2955         z__[1] = 1.;
2956         return;
2957     }
2958
2959     eps = GMX_FLOAT_EPS;
2960
2961     d__1   = eps;
2962     eps2   = d__1 * d__1;
2963     minval = GMX_FLOAT_MIN;
2964     safmin = minval / GMX_FLOAT_EPS;
2965     safmax = 1. / safmin;
2966     ssfmax = std::sqrt(safmax) / 3.;
2967     ssfmin = std::sqrt(safmin) / eps2;
2968
2969     i__1 = *n - 1;
2970     for (j = 1; j <= i__1; ++j)
2971     {
2972         z__[j] = 0.;
2973
2974     }
2975     z__[*n] = 1.;
2976
2977     nmaxit = *n * 30;
2978     jtot   = 0;
2979
2980     l1  = 1;
2981     nm1 = *n - 1;
2982
2983 L10:
2984     if (l1 > *n)
2985     {
2986         goto L160;
2987     }
2988     if (l1 > 1)
2989     {
2990         e[l1 - 1] = 0.;
2991     }
2992     if (l1 <= nm1)
2993     {
2994         i__1 = nm1;
2995         for (m = l1; m <= i__1; ++m)
2996         {
2997             tst = std::abs(e[m]);
2998             if (tst == 0.)
2999             {
3000                 goto L30;
3001             }
3002             if (tst <= std::sqrt(std::abs(d__[m])) * std::sqrt(std::abs(d__[m+1])) * eps)
3003             {
3004                 e[m] = 0.;
3005                 goto L30;
3006             }
3007         }
3008     }
3009     m = *n;
3010
3011 L30:
3012     l      = l1;
3013     lsv    = l;
3014     lend   = m;
3015     lendsv = lend;
3016     l1     = m + 1;
3017     if (lend == l)
3018     {
3019         goto L10;
3020     }
3021
3022     i__1   = lend - l + 1;
3023     anorm  = F77_FUNC(slanst, SLANST) ("i", &i__1, &d__[l], &e[l]);
3024     iscale = 0;
3025     if (anorm == 0.)
3026     {
3027         goto L10;
3028     }
3029     if (anorm > ssfmax)
3030     {
3031         iscale = 1;
3032         i__1   = lend - l + 1;
3033         F77_FUNC(slascl, SLASCL) ("g", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &d__[l], n,
3034                                   info);
3035         i__1 = lend - l;
3036         F77_FUNC(slascl, SLASCL) ("g", &c__0, &c__0, &anorm, &ssfmax, &i__1, &c__1, &e[l], n,
3037                                   info);
3038     }
3039     else if (anorm < ssfmin)
3040     {
3041         iscale = 2;
3042         i__1   = lend - l + 1;
3043         F77_FUNC(slascl, SLASCL) ("g", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &d__[l], n,
3044                                   info);
3045         i__1 = lend - l;
3046         F77_FUNC(slascl, SLASCL) ("g", &c__0, &c__0, &anorm, &ssfmin, &i__1, &c__1, &e[l], n,
3047                                   info);
3048     }
3049
3050     if (std::abs(d__[lend]) < std::abs(d__[l]))
3051     {
3052         lend = lsv;
3053         l    = lendsv;
3054     }
3055
3056     if (lend > l)
3057     {
3058
3059 L40:
3060         if (l != lend)
3061         {
3062             lendm1 = lend - 1;
3063             i__1   = lendm1;
3064             for (m = l; m <= i__1; ++m)
3065             {
3066                 d__2 = std::abs(e[m]);
3067                 tst  = d__2 * d__2;
3068                 if (tst <= eps2 * std::abs(d__[m]) * std::abs(d__[m + 1]) + safmin)
3069                 {
3070                     goto L60;
3071                 }
3072             }
3073         }
3074
3075         m = lend;
3076
3077 L60:
3078         if (m < lend)
3079         {
3080             e[m] = 0.;
3081         }
3082         p = d__[l];
3083         if (m == l)
3084         {
3085             goto L80;
3086         }
3087
3088         if (m == l + 1)
3089         {
3090             if (icompz > 0)
3091             {
3092                 F77_FUNC(slaev2, SLAEV2) (&d__[l], &e[l], &d__[l + 1], &rt1, &rt2, &c__, &s);
3093                 work[l]          = c__;
3094                 work[*n - 1 + l] = s;
3095
3096                 tst        = z__[l + 1];
3097                 z__[l + 1] = c__ * tst - s * z__[l];
3098                 z__[l]     = s * tst + c__ * z__[l];
3099             }
3100             else
3101             {
3102                 F77_FUNC(slae2, SLAE2) (&d__[l], &e[l], &d__[l + 1], &rt1, &rt2);
3103             }
3104             d__[l]     = rt1;
3105             d__[l + 1] = rt2;
3106             e[l]       = 0.;
3107             l         += 2;
3108             if (l <= lend)
3109             {
3110                 goto L40;
3111             }
3112             goto L140;
3113         }
3114
3115         if (jtot == nmaxit)
3116         {
3117             goto L140;
3118         }
3119         ++jtot;
3120
3121         g   = (d__[l + 1] - p) / (e[l] * 2.);
3122         r__ = F77_FUNC(slapy2, SLAPY2) (&g, &c_b31);
3123         g   = d__[m] - p + e[l] / (g + ((g > 0) ? r__ : -r__ ));
3124
3125         s   = 1.;
3126         c__ = 1.;
3127         p   = 0.;
3128
3129         mm1  = m - 1;
3130         i__1 = l;
3131         for (i__ = mm1; i__ >= i__1; --i__)
3132         {
3133             f = s * e[i__];
3134             b = c__ * e[i__];
3135             F77_FUNC(slartg, SLARTG) (&g, &f, &c__, &s, &r__);
3136             if (i__ != m - 1)
3137             {
3138                 e[i__ + 1] = r__;
3139             }
3140             g            = d__[i__ + 1] - p;
3141             r__          = (d__[i__] - g) * s + c__ * 2. * b;
3142             p            = s * r__;
3143             d__[i__ + 1] = g + p;
3144             g            = c__ * r__ - b;
3145
3146             if (icompz > 0)
3147             {
3148                 work[i__]          = c__;
3149                 work[*n - 1 + i__] = -s;
3150             }
3151
3152         }
3153
3154         if (icompz > 0)
3155         {
3156             mm = m - l + 1;
3157
3158             F77_FUNC(slasr, SLASR) ("r", "v", "b", &c__1, &mm, &work[l], &work[*n - 1 + l], &
3159                                     z__[l], &c__1);
3160         }
3161
3162         d__[l] -= p;
3163         e[l]    = g;
3164         goto L40;
3165
3166 L80:
3167         d__[l] = p;
3168
3169         ++l;
3170         if (l <= lend)
3171         {
3172             goto L40;
3173         }
3174         goto L140;
3175
3176     }
3177     else
3178     {
3179
3180 L90:
3181         if (l != lend)
3182         {
3183             lendp1 = lend + 1;
3184             i__1   = lendp1;
3185             for (m = l; m >= i__1; --m)
3186             {
3187                 d__2 = std::abs(e[m - 1]);
3188                 tst  = d__2 * d__2;
3189                 if (tst <= eps2 * std::abs(d__[m]) * std::abs(d__[m- 1]) + safmin)
3190                 {
3191                     goto L110;
3192                 }
3193             }
3194         }
3195
3196         m = lend;
3197
3198 L110:
3199         if (m > lend)
3200         {
3201             e[m - 1] = 0.;
3202         }
3203         p = d__[l];
3204         if (m == l)
3205         {
3206             goto L130;
3207         }
3208
3209         if (m == l - 1)
3210         {
3211             if (icompz > 0)
3212             {
3213                 F77_FUNC(slaev2, SLAEV2) (&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2, &c__, &s)
3214                 ;
3215
3216                 tst        = z__[l];
3217                 z__[l]     = c__ * tst - s * z__[l - 1];
3218                 z__[l - 1] = s * tst + c__ * z__[l - 1];
3219             }
3220             else
3221             {
3222                 F77_FUNC(slae2, SLAE2) (&d__[l - 1], &e[l - 1], &d__[l], &rt1, &rt2);
3223             }
3224             d__[l - 1] = rt1;
3225             d__[l]     = rt2;
3226             e[l - 1]   = 0.;
3227             l         += -2;
3228             if (l >= lend)
3229             {
3230                 goto L90;
3231             }
3232             goto L140;
3233         }
3234
3235         if (jtot == nmaxit)
3236         {
3237             goto L140;
3238         }
3239         ++jtot;
3240
3241
3242         g   = (d__[l - 1] - p) / (e[l - 1] * 2.);
3243         r__ = F77_FUNC(slapy2, SLAPY2) (&g, &c_b31);
3244         g   = d__[m] - p + e[l - 1] / (g + ((g > 0) ? r__ : -r__ ));
3245
3246         s   = 1.;
3247         c__ = 1.;
3248         p   = 0.;
3249
3250         lm1  = l - 1;
3251         i__1 = lm1;
3252         for (i__ = m; i__ <= i__1; ++i__)
3253         {
3254             f = s * e[i__];
3255             b = c__ * e[i__];
3256             F77_FUNC(slartg, SLARTG) (&g, &f, &c__, &s, &r__);
3257             if (i__ != m)
3258             {
3259                 e[i__ - 1] = r__;
3260             }
3261             g        = d__[i__] - p;
3262             r__      = (d__[i__ + 1] - g) * s + c__ * 2. * b;
3263             p        = s * r__;
3264             d__[i__] = g + p;
3265             g        = c__ * r__ - b;
3266
3267             if (icompz > 0)
3268             {
3269                 work[i__]          = c__;
3270                 work[*n - 1 + i__] = s;
3271             }
3272
3273         }
3274
3275         if (icompz > 0)
3276         {
3277             mm = l - m + 1;
3278
3279             F77_FUNC(slasr, SLASR) ("r", "v", "f", &c__1, &mm, &work[m], &work[*n - 1 + m], &
3280                                     z__[m], &c__1);
3281         }
3282
3283         d__[l] -= p;
3284         e[lm1]  = g;
3285         goto L90;
3286
3287 L130:
3288         d__[l] = p;
3289
3290         --l;
3291         if (l >= lend)
3292         {
3293             goto L90;
3294         }
3295         goto L140;
3296
3297     }
3298
3299 L140:
3300     if (iscale == 1)
3301     {
3302         i__1 = lendsv - lsv + 1;
3303         F77_FUNC(slascl, SLASCL) ("g", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &d__[lsv],
3304                                   n, info);
3305         i__1 = lendsv - lsv;
3306         F77_FUNC(slascl, SLASCL) ("g", &c__0, &c__0, &ssfmax, &anorm, &i__1, &c__1, &e[lsv], n,
3307                                   info);
3308     }
3309     else if (iscale == 2)
3310     {
3311         i__1 = lendsv - lsv + 1;
3312         F77_FUNC(slascl, SLASCL) ("g", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &d__[lsv],
3313                                   n, info);
3314         i__1 = lendsv - lsv;
3315         F77_FUNC(slascl, SLASCL) ("g", &c__0, &c__0, &ssfmin, &anorm, &i__1, &c__1, &e[lsv], n,
3316                                   info);
3317     }
3318
3319     if (jtot < nmaxit)
3320     {
3321         goto L10;
3322     }
3323     i__1 = *n - 1;
3324     for (i__ = 1; i__ <= i__1; ++i__)
3325     {
3326         if (e[i__] != 0.)
3327         {
3328             ++(*info);
3329         }
3330     }
3331     goto L190;
3332
3333 L160:
3334     if (icompz == 0)
3335     {
3336
3337         F77_FUNC(slasrt, SLASRT) ("i", n, &d__[1], info);
3338
3339     }
3340     else
3341     {
3342
3343         i__1 = *n;
3344         for (ii = 2; ii <= i__1; ++ii)
3345         {
3346             i__  = ii - 1;
3347             k    = i__;
3348             p    = d__[i__];
3349             i__2 = *n;
3350             for (j = ii; j <= i__2; ++j)
3351             {
3352                 if (d__[j] < p)
3353                 {
3354                     k = j;
3355                     p = d__[j];
3356                 }
3357             }
3358             if (k != i__)
3359             {
3360                 d__[k]   = d__[i__];
3361                 d__[i__] = p;
3362
3363                 p        = z__[k];
3364                 z__[k]   = z__[i__];
3365                 z__[i__] = p;
3366             }
3367         }
3368     }
3369
3370 L190:
3371     return;
3372
3373 }
3374
3375 static void
3376 F77_FUNC(sgetv0, SGETV0) (int *     ido,
3377                           const char *    bmat,
3378                           int gmx_unused *     itry,
3379                           int *     initv,
3380                           int *     n,
3381                           int *     j,
3382                           float *  v,
3383                           int *     ldv,
3384                           float *  resid,
3385                           float *  rnorm,
3386                           int *     ipntr,
3387                           float *  workd,
3388                           int *     iwork,
3389                           int *     ierr)
3390 {
3391     int   c__1  = 1;
3392     float c_b22 = 1.;
3393     float c_b24 = 0.;
3394     float c_b27 = -1.;
3395     int   v_dim1, v_offset, i__1;
3396
3397     int   jj;
3398     int   idist;
3399
3400     --workd;
3401     --resid;
3402     v_dim1   = *ldv;
3403     v_offset = 1 + v_dim1;
3404     v       -= v_offset;
3405     --ipntr;
3406     --iwork;
3407
3408     if (*ido == 0)
3409     {
3410
3411         *ierr    = 0;
3412         iwork[7] = 0;
3413         iwork[5] = 0;
3414         iwork[6] = 0;
3415
3416         if (!(*initv))
3417         {
3418             idist = 2;
3419             F77_FUNC(slarnv, SLARNV) (&idist, &iwork[1], n, &resid[1]);
3420         }
3421
3422         if (*bmat == 'G')
3423         {
3424             ipntr[1] = 1;
3425             ipntr[2] = *n + 1;
3426             F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[1], &c__1);
3427             *ido = -1;
3428             goto L9000;
3429         }
3430     }
3431
3432     if (iwork[5] == 1)
3433     {
3434         goto L20;
3435     }
3436
3437     if (iwork[6] == 1)
3438     {
3439         goto L40;
3440     }
3441
3442     iwork[5] = 1;
3443     if (*bmat == 'G')
3444     {
3445         F77_FUNC(scopy, SCOPY) (n, &workd[*n + 1], &c__1, &resid[1], &c__1);
3446         ipntr[1] = *n + 1;
3447         ipntr[2] = 1;
3448         *ido     = 2;
3449         goto L9000;
3450     }
3451     else if (*bmat == 'I')
3452     {
3453         F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[1], &c__1);
3454     }
3455
3456 L20:
3457
3458
3459     iwork[5] = 0;
3460     if (*bmat == 'G')
3461     {
3462         workd[*n * 3 + 4] = F77_FUNC(sdot, SDOT) (n, &resid[1], &c__1, &workd[1], &c__1);
3463         workd[*n * 3 + 4] = std::sqrt(std::abs(workd[*n * 3 + 4]));
3464     }
3465     else if (*bmat == 'I')
3466     {
3467         workd[*n * 3 + 4] = F77_FUNC(snrm2, SNRM2) (n, &resid[1], &c__1);
3468     }
3469     *rnorm = workd[*n * 3 + 4];
3470
3471     if (*j == 1)
3472     {
3473         goto L50;
3474     }
3475     iwork[6] = 1;
3476 L30:
3477
3478     i__1 = *j - 1;
3479     F77_FUNC(sgemv, SGEMV) ("T", n, &i__1, &c_b22, &v[v_offset], ldv, &workd[1], &c__1, &c_b24,
3480                             &workd[*n + 1], &c__1);
3481     i__1 = *j - 1;
3482     F77_FUNC(sgemv, SGEMV) ("N", n, &i__1, &c_b27, &v[v_offset], ldv, &workd[*n + 1], &c__1, &
3483                             c_b22, &resid[1], &c__1);
3484
3485     if (*bmat == 'G')
3486     {
3487         F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[*n + 1], &c__1);
3488         ipntr[1] = *n + 1;
3489         ipntr[2] = 1;
3490         *ido     = 2;
3491         goto L9000;
3492     }
3493     else if (*bmat == 'I')
3494     {
3495         F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[1], &c__1);
3496     }
3497
3498 L40:
3499
3500     if (*bmat == 'G')
3501     {
3502         *rnorm = F77_FUNC(sdot, SDOT) (n, &resid[1], &c__1, &workd[1], &c__1);
3503         *rnorm = std::sqrt(std::abs(*rnorm));
3504     }
3505     else if (*bmat == 'I')
3506     {
3507         *rnorm = F77_FUNC(snrm2, SNRM2) (n, &resid[1], &c__1);
3508     }
3509
3510     if (*rnorm > workd[*n * 3 + 4] * .717F)
3511     {
3512         goto L50;
3513     }
3514
3515     ++iwork[7];
3516     if (iwork[7] <= 1)
3517     {
3518
3519         workd[*n * 3 + 4] = *rnorm;
3520         goto L30;
3521     }
3522     else
3523     {
3524
3525         i__1 = *n;
3526         for (jj = 1; jj <= i__1; ++jj)
3527         {
3528             resid[jj] = 0.;
3529         }
3530         *rnorm = 0.;
3531         *ierr  = -1;
3532     }
3533
3534 L50:
3535
3536     *ido = 99;
3537
3538 L9000:
3539     return;
3540 }
3541
3542
3543
3544
3545
3546 static void
3547 F77_FUNC(ssapps, SSAPPS) (int *     n,
3548                           int *     kev,
3549                           int *     np,
3550                           float *  shift,
3551                           float *  v,
3552                           int *     ldv,
3553                           float *  h__,
3554                           int *     ldh,
3555                           float *  resid,
3556                           float *  q,
3557                           int *     ldq,
3558                           float *  workd)
3559 {
3560     float c_b4  = 0.;
3561     float c_b5  = 1.;
3562     float c_b14 = -1.;
3563     int   c__1  = 1;
3564     int   h_dim1, h_offset, q_dim1, q_offset, v_dim1, v_offset, i__1, i__2,
3565           i__3, i__4;
3566     float c__, f, g;
3567     int   i__, j;
3568     float r__, s, a1, a2, a3, a4;
3569     int   jj;
3570     float big;
3571     int   iend, itop;
3572     float epsmch;
3573     int   istart, kplusp;
3574
3575     --workd;
3576     --resid;
3577     --shift;
3578     v_dim1   = *ldv;
3579     v_offset = 1 + v_dim1;
3580     v       -= v_offset;
3581     h_dim1   = *ldh;
3582     h_offset = 1 + h_dim1;
3583     h__     -= h_offset;
3584     q_dim1   = *ldq;
3585     q_offset = 1 + q_dim1;
3586     q       -= q_offset;
3587
3588     epsmch = GMX_FLOAT_EPS;
3589     itop   = 1;
3590
3591
3592     kplusp = *kev + *np;
3593
3594     F77_FUNC(slaset, SLASET) ("All", &kplusp, &kplusp, &c_b4, &c_b5, &q[q_offset], ldq);
3595
3596     if (*np == 0)
3597     {
3598         goto L9000;
3599     }
3600
3601     i__1 = *np;
3602     for (jj = 1; jj <= i__1; ++jj)
3603     {
3604
3605         istart = itop;
3606
3607 L20:
3608
3609         i__2 = kplusp - 1;
3610         for (i__ = istart; i__ <= i__2; ++i__)
3611         {
3612             big = std::abs(h__[i__ + (h_dim1*2)]) + std::abs(h__[i__ + 1 + (h_dim1*2)]);
3613             if (h__[i__ + 1 + h_dim1] <= epsmch * big)
3614             {
3615                 h__[i__ + 1 + h_dim1] = 0.;
3616                 iend                  = i__;
3617                 goto L40;
3618             }
3619         }
3620         iend = kplusp;
3621 L40:
3622
3623         if (istart < iend)
3624         {
3625
3626             f = h__[istart + (h_dim1 << 1)] - shift[jj];
3627             g = h__[istart + 1 + h_dim1];
3628             F77_FUNC(slartg, SLARTG) (&f, &g, &c__, &s, &r__);
3629
3630             a1 = c__ * h__[istart + (h_dim1 << 1)] + s * h__[istart + 1 +
3631                                                              h_dim1];
3632             a2 = c__ * h__[istart + 1 + h_dim1] + s * h__[istart + 1 + (
3633                                                               h_dim1 << 1)];
3634             a4 = c__ * h__[istart + 1 + (h_dim1 << 1)] - s * h__[istart + 1 +
3635                                                                  h_dim1];
3636             a3 = c__ * h__[istart + 1 + h_dim1] - s * h__[istart + (h_dim1 <<
3637                                                                     1)];
3638             h__[istart + (h_dim1 << 1)]     = c__ * a1 + s * a2;
3639             h__[istart + 1 + (h_dim1 << 1)] = c__ * a4 - s * a3;
3640             h__[istart + 1 + h_dim1]        = c__ * a3 + s * a4;
3641
3642             i__3 = istart + jj;
3643             i__2 = (i__3 < kplusp) ? i__3 : kplusp;
3644             for (j = 1; j <= i__2; ++j)
3645             {
3646                 a1 = c__ * q[j + istart * q_dim1] + s * q[j + (istart + 1) *
3647                                                           q_dim1];
3648                 q[j + (istart + 1) * q_dim1] = -s * q[j + istart * q_dim1] +
3649                     c__ * q[j + (istart + 1) * q_dim1];
3650                 q[j + istart * q_dim1] = a1;
3651
3652             }
3653
3654             i__2 = iend - 1;
3655             for (i__ = istart + 1; i__ <= i__2; ++i__)
3656             {
3657
3658                 f = h__[i__ + h_dim1];
3659                 g = s * h__[i__ + 1 + h_dim1];
3660
3661                 h__[i__ + 1 + h_dim1] = c__ * h__[i__ + 1 + h_dim1];
3662                 F77_FUNC(slartg, SLARTG) (&f, &g, &c__, &s, &r__);
3663
3664                 if (r__ < 0.)
3665                 {
3666                     r__ = -r__;
3667                     c__ = -c__;
3668                     s   = -s;
3669                 }
3670
3671                 h__[i__ + h_dim1] = r__;
3672
3673                 a1 = c__ * h__[i__ + (h_dim1 << 1)] + s * h__[i__ + 1 +
3674                                                               h_dim1];
3675                 a2 = c__ * h__[i__ + 1 + h_dim1] + s * h__[i__ + 1 + (h_dim1
3676                                                                       << 1)];
3677                 a3 = c__ * h__[i__ + 1 + h_dim1] - s * h__[i__ + (h_dim1 << 1)
3678                     ];
3679                 a4 = c__ * h__[i__ + 1 + (h_dim1 << 1)] - s * h__[i__ + 1 +
3680                                                                   h_dim1];
3681
3682                 h__[i__ + (h_dim1 << 1)]     = c__ * a1 + s * a2;
3683                 h__[i__ + 1 + (h_dim1 << 1)] = c__ * a4 - s * a3;
3684                 h__[i__ + 1 + h_dim1]        = c__ * a3 + s * a4;
3685
3686                 i__4 = j + jj;
3687                 i__3 = (i__4 < kplusp) ? i__4 : kplusp;
3688                 for (j = 1; j <= i__3; ++j)
3689                 {
3690                     a1 = c__ * q[j + i__ * q_dim1] + s * q[j + (i__ + 1) *
3691                                                            q_dim1];
3692                     q[j + (i__ + 1) * q_dim1] = -s * q[j + i__ * q_dim1] +
3693                         c__ * q[j + (i__ + 1) * q_dim1];
3694                     q[j + i__ * q_dim1] = a1;
3695                 }
3696
3697             }
3698
3699         }
3700
3701         istart = iend + 1;
3702
3703         if (h__[iend + h_dim1] < 0.)
3704         {
3705             h__[iend + h_dim1] = -h__[iend + h_dim1];
3706             F77_FUNC(sscal, SSCAL) (&kplusp, &c_b14, &q[iend * q_dim1 + 1], &c__1);
3707         }
3708
3709         if (iend < kplusp)
3710         {
3711             goto L20;
3712         }
3713
3714         i__2 = kplusp - 1;
3715         for (i__ = itop; i__ <= i__2; ++i__)
3716         {
3717             if (h__[i__ + 1 + h_dim1] > 0.)
3718             {
3719                 goto L90;
3720             }
3721             ++itop;
3722         }
3723
3724 L90:
3725         ;
3726     }
3727
3728     i__1 = kplusp - 1;
3729     for (i__ = itop; i__ <= i__1; ++i__)
3730     {
3731         big = std::abs(h__[i__ + (h_dim1*2)]) + std::abs(h__[i__+ 1 + (h_dim1*2)]);
3732         if (h__[i__ + 1 + h_dim1] <= epsmch * big)
3733         {
3734             h__[i__ + 1 + h_dim1] = 0.;
3735         }
3736
3737     }
3738
3739     if (h__[*kev + 1 + h_dim1] > 0.)
3740     {
3741         F77_FUNC(sgemv, SGEMV) ("N", n, &kplusp, &c_b5, &v[v_offset], ldv, &q[(*kev + 1) *
3742                                                                               q_dim1 + 1], &c__1, &c_b4, &workd[*n + 1], &c__1);
3743     }
3744
3745     i__1 = *kev;
3746     for (i__ = 1; i__ <= i__1; ++i__)
3747     {
3748         i__2 = kplusp - i__ + 1;
3749         F77_FUNC(sgemv, SGEMV) ("N", n, &i__2, &c_b5, &v[v_offset], ldv, &q[(*kev - i__ + 1) *
3750                                                                             q_dim1 + 1], &c__1, &c_b4, &workd[1], &c__1);
3751         F77_FUNC(scopy, SCOPY) (n, &workd[1], &c__1, &v[(kplusp - i__ + 1) * v_dim1 + 1], &
3752                                 c__1);
3753
3754     }
3755
3756     F77_FUNC(slacpy, SLACPY) ("All", n, kev, &v[(*np + 1) * v_dim1 + 1], ldv, &v[v_offset], ldv);
3757
3758     if (h__[*kev + 1 + h_dim1] > 0.)
3759     {
3760         F77_FUNC(scopy, SCOPY) (n, &workd[*n + 1], &c__1, &v[(*kev + 1) * v_dim1 + 1], &c__1);
3761     }
3762
3763     F77_FUNC(sscal, SSCAL) (n, &q[kplusp + *kev * q_dim1], &resid[1], &c__1);
3764     if (h__[*kev + 1 + h_dim1] > 0.)
3765     {
3766         F77_FUNC(saxpy, SAXPY) (n, &h__[*kev + 1 + h_dim1], &v[(*kev + 1) * v_dim1 + 1], &c__1,
3767                                 &resid[1], &c__1);
3768     }
3769
3770
3771
3772 L9000:
3773     return;
3774
3775 }
3776
3777
3778
3779 static void
3780 F77_FUNC(ssortr, SSORTR) (const char *    which,
3781                           int *     apply,
3782                           int *     n,
3783                           float *  x1,
3784                           float *  x2)
3785 {
3786     int   i__1;
3787
3788     int   i__, j, igap;
3789     float temp;
3790
3791
3792
3793     igap = *n / 2;
3794
3795     if (!std::strncmp(which, "SA", 2))
3796     {
3797
3798 L10:
3799         if (igap == 0)
3800         {
3801             goto L9000;
3802         }
3803         i__1 = *n - 1;
3804         for (i__ = igap; i__ <= i__1; ++i__)
3805         {
3806             j = i__ - igap;
3807 L20:
3808
3809             if (j < 0)
3810             {
3811                 goto L30;
3812             }
3813
3814             if (x1[j] < x1[j + igap])
3815             {
3816                 temp         = x1[j];
3817                 x1[j]        = x1[j + igap];
3818                 x1[j + igap] = temp;
3819                 if (*apply)
3820                 {
3821                     temp         = x2[j];
3822                     x2[j]        = x2[j + igap];
3823                     x2[j + igap] = temp;
3824                 }
3825             }
3826             else
3827             {
3828                 goto L30;
3829             }
3830             j -= igap;
3831             goto L20;
3832 L30:
3833             ;
3834         }
3835         igap /= 2;
3836         goto L10;
3837
3838     }
3839     else if (!std::strncmp(which, "SM", 2))
3840     {
3841
3842 L40:
3843         if (igap == 0)
3844         {
3845             goto L9000;
3846         }
3847         i__1 = *n - 1;
3848         for (i__ = igap; i__ <= i__1; ++i__)
3849         {
3850             j = i__ - igap;
3851 L50:
3852
3853             if (j < 0)
3854             {
3855                 goto L60;
3856             }
3857
3858             if (std::abs(x1[j]) < std::abs(x1[j + igap]))
3859             {
3860                 temp         = x1[j];
3861                 x1[j]        = x1[j + igap];
3862                 x1[j + igap] = temp;
3863                 if (*apply)
3864                 {
3865                     temp         = x2[j];
3866                     x2[j]        = x2[j + igap];
3867                     x2[j + igap] = temp;
3868                 }
3869             }
3870             else
3871             {
3872                 goto L60;
3873             }
3874             j -= igap;
3875             goto L50;
3876 L60:
3877             ;
3878         }
3879         igap /= 2;
3880         goto L40;
3881
3882     }
3883     else if (!std::strncmp(which, "LA", 2))
3884     {
3885
3886 L70:
3887         if (igap == 0)
3888         {
3889             goto L9000;
3890         }
3891         i__1 = *n - 1;
3892         for (i__ = igap; i__ <= i__1; ++i__)
3893         {
3894             j = i__ - igap;
3895 L80:
3896
3897             if (j < 0)
3898             {
3899                 goto L90;
3900             }
3901
3902             if (x1[j] > x1[j + igap])
3903             {
3904                 temp         = x1[j];
3905                 x1[j]        = x1[j + igap];
3906                 x1[j + igap] = temp;
3907                 if (*apply)
3908                 {
3909                     temp         = x2[j];
3910                     x2[j]        = x2[j + igap];
3911                     x2[j + igap] = temp;
3912                 }
3913             }
3914             else
3915             {
3916                 goto L90;
3917             }
3918             j -= igap;
3919             goto L80;
3920 L90:
3921             ;
3922         }
3923         igap /= 2;
3924         goto L70;
3925
3926     }
3927     else if (!std::strncmp(which, "LM", 2))
3928     {
3929
3930
3931 L100:
3932         if (igap == 0)
3933         {
3934             goto L9000;
3935         }
3936         i__1 = *n - 1;
3937         for (i__ = igap; i__ <= i__1; ++i__)
3938         {
3939             j = i__ - igap;
3940 L110:
3941
3942             if (j < 0)
3943             {
3944                 goto L120;
3945             }
3946
3947             if (std::abs(x1[j]) > std::abs(x1[j + igap]))
3948             {
3949                 temp         = x1[j];
3950                 x1[j]        = x1[j + igap];
3951                 x1[j + igap] = temp;
3952                 if (*apply)
3953                 {
3954                     temp         = x2[j];
3955                     x2[j]        = x2[j + igap];
3956                     x2[j + igap] = temp;
3957                 }
3958             }
3959             else
3960             {
3961                 goto L120;
3962             }
3963             j -= igap;
3964             goto L110;
3965 L120:
3966             ;
3967         }
3968         igap /= 2;
3969         goto L100;
3970     }
3971
3972 L9000:
3973     return;
3974
3975 }
3976
3977
3978
3979
3980 static void
3981 F77_FUNC(ssesrt, SSESRT) (const char *    which,
3982                           int *     apply,
3983                           int *     n,
3984                           float *  x,
3985                           int *     na,
3986                           float *  a,
3987                           int *     lda)
3988 {
3989     int   a_dim1, a_offset, i__1;
3990     int   c__1 = 1;
3991
3992     int   i__, j, igap;
3993     float temp;
3994
3995     a_dim1   = *lda;
3996     a_offset = 1 + a_dim1 * 0;
3997     a       -= a_offset;
3998
3999     igap = *n / 2;
4000
4001     if (!std::strncmp(which, "SA", 2))
4002     {
4003
4004 L10:
4005         if (igap == 0)
4006         {
4007             goto L9000;
4008         }
4009         i__1 = *n - 1;
4010         for (i__ = igap; i__ <= i__1; ++i__)
4011         {
4012             j = i__ - igap;
4013 L20:
4014
4015             if (j < 0)
4016             {
4017                 goto L30;
4018             }
4019
4020             if (x[j] < x[j + igap])
4021             {
4022                 temp        = x[j];
4023                 x[j]        = x[j + igap];
4024                 x[j + igap] = temp;
4025                 if (*apply)
4026                 {
4027                     F77_FUNC(sswap, SSWAP) (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
4028                                                                               a_dim1 + 1], &c__1);
4029                 }
4030             }
4031             else
4032             {
4033                 goto L30;
4034             }
4035             j -= igap;
4036             goto L20;
4037 L30:
4038             ;
4039         }
4040         igap /= 2;
4041         goto L10;
4042
4043     }
4044     else if (!std::strncmp(which, "SM", 2))
4045     {
4046
4047 L40:
4048         if (igap == 0)
4049         {
4050             goto L9000;
4051         }
4052         i__1 = *n - 1;
4053         for (i__ = igap; i__ <= i__1; ++i__)
4054         {
4055             j = i__ - igap;
4056 L50:
4057
4058             if (j < 0)
4059             {
4060                 goto L60;
4061             }
4062
4063             if (std::abs(x[j]) < std::abs(x[j + igap]))
4064             {
4065                 temp        = x[j];
4066                 x[j]        = x[j + igap];
4067                 x[j + igap] = temp;
4068                 if (*apply)
4069                 {
4070                     F77_FUNC(sswap, SSWAP) (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
4071                                                                               a_dim1 + 1], &c__1);
4072                 }
4073             }
4074             else
4075             {
4076                 goto L60;
4077             }
4078             j -= igap;
4079             goto L50;
4080 L60:
4081             ;
4082         }
4083         igap /= 2;
4084         goto L40;
4085
4086     }
4087     else if (!std::strncmp(which, "LA", 2))
4088     {
4089
4090 L70:
4091         if (igap == 0)
4092         {
4093             goto L9000;
4094         }
4095         i__1 = *n - 1;
4096         for (i__ = igap; i__ <= i__1; ++i__)
4097         {
4098             j = i__ - igap;
4099 L80:
4100
4101             if (j < 0)
4102             {
4103                 goto L90;
4104             }
4105
4106             if (x[j] > x[j + igap])
4107             {
4108                 temp        = x[j];
4109                 x[j]        = x[j + igap];
4110                 x[j + igap] = temp;
4111                 if (*apply)
4112                 {
4113                     F77_FUNC(sswap, SSWAP) (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
4114                                                                               a_dim1 + 1], &c__1);
4115                 }
4116             }
4117             else
4118             {
4119                 goto L90;
4120             }
4121             j -= igap;
4122             goto L80;
4123 L90:
4124             ;
4125         }
4126         igap /= 2;
4127         goto L70;
4128
4129     }
4130     else if (!std::strncmp(which, "LM", 2))
4131     {
4132
4133 L100:
4134         if (igap == 0)
4135         {
4136             goto L9000;
4137         }
4138         i__1 = *n - 1;
4139         for (i__ = igap; i__ <= i__1; ++i__)
4140         {
4141             j = i__ - igap;
4142 L110:
4143
4144             if (j < 0)
4145             {
4146                 goto L120;
4147             }
4148
4149             if (std::abs(x[j]) > std::abs(x[j + igap]))
4150             {
4151                 temp        = x[j];
4152                 x[j]        = x[j + igap];
4153                 x[j + igap] = temp;
4154                 if (*apply)
4155                 {
4156                     F77_FUNC(sswap, SSWAP) (na, &a[j * a_dim1 + 1], &c__1, &a[(j + igap) *
4157                                                                               a_dim1 + 1], &c__1);
4158                 }
4159             }
4160             else
4161             {
4162                 goto L120;
4163             }
4164             j -= igap;
4165             goto L110;
4166 L120:
4167             ;
4168         }
4169         igap /= 2;
4170         goto L100;
4171     }
4172
4173 L9000:
4174     return;
4175
4176 }
4177
4178
4179
4180
4181 static void
4182 F77_FUNC(ssgets, SSGETS) (int *     ishift,
4183                           const char *    which,
4184                           int *     kev,
4185                           int *     np,
4186                           float *  ritz,
4187                           float *  bounds,
4188                           float *  shifts)
4189 {
4190     int c__1 = 1;
4191     int i__1, i__2;
4192     int kevd2;
4193
4194     --shifts;
4195     --bounds;
4196     --ritz;
4197
4198     if (!std::strncmp(which, "BE", 2))
4199     {
4200         i__1 = *kev + *np;
4201         F77_FUNC(ssortr, SSORTR) ("LA", &c__1, &i__1, &ritz[1], &bounds[1]);
4202         kevd2 = *kev / 2;
4203         if (*kev > 1)
4204         {
4205             i__1 = (kevd2 < *np) ? kevd2 : *np;
4206             i__2 = (kevd2 > *np) ? kevd2 : *np;
4207             F77_FUNC(sswap, SSWAP) (&i__1, &ritz[1], &c__1,
4208                                     &ritz[i__2 + 1], &c__1);
4209             i__1 = (kevd2 < *np) ? kevd2 : *np;
4210             i__2 = (kevd2 > *np) ? kevd2 : *np;
4211             F77_FUNC(sswap, SSWAP) (&i__1, &bounds[1], &c__1,
4212                                     &bounds[i__2 + 1], &c__1);
4213         }
4214
4215     }
4216     else
4217     {
4218         i__1 = *kev + *np;
4219         F77_FUNC(ssortr, SSORTR) (which, &c__1, &i__1, &ritz[1], &bounds[1]);
4220     }
4221
4222     if (*ishift == 1 && *np > 0)
4223     {
4224
4225         F77_FUNC(ssortr, SSORTR) ("SM", &c__1, np, &bounds[1], &ritz[1]);
4226         F77_FUNC(scopy, SCOPY) (np, &ritz[1], &c__1, &shifts[1], &c__1);
4227     }
4228
4229
4230     return;
4231 }
4232
4233
4234
4235 static void
4236 F77_FUNC(ssconv, SSCONV) (int *     n,
4237                           float *  ritz,
4238                           float *  bounds,
4239                           float *  tol,
4240                           int *     nconv)
4241 {
4242     float c_b3 = 2/3.;
4243     int   i__1;
4244     float d__2, d__3;
4245
4246     int   i__;
4247     float eps23, temp;
4248
4249     --bounds;
4250     --ritz;
4251
4252     eps23 = GMX_FLOAT_EPS;
4253     eps23 = std::pow(eps23, c_b3);
4254
4255     *nconv = 0;
4256     i__1   = *n;
4257     for (i__ = 1; i__ <= i__1; ++i__)
4258     {
4259
4260         d__2 = eps23;
4261         d__3 = std::abs(ritz[i__]);
4262         temp = (d__2 > d__3) ? d__2 : d__3;
4263         if (bounds[i__] <= *tol * temp)
4264         {
4265             ++(*nconv);
4266         }
4267     }
4268
4269     return;
4270 }
4271
4272
4273 static void
4274 F77_FUNC(sseigt, SSEIGT) (float *  rnorm,
4275                           int *     n,
4276                           float *  h__,
4277                           int *     ldh,
4278                           float *  eig,
4279                           float *  bounds,
4280                           float *  workl,
4281                           int *     ierr)
4282 {
4283     int c__1 = 1;
4284     int h_dim1, h_offset, i__1;
4285
4286     int k;
4287
4288
4289     --workl;
4290     --bounds;
4291     --eig;
4292     h_dim1   = *ldh;
4293     h_offset = 1 + h_dim1;
4294     h__     -= h_offset;
4295
4296     F77_FUNC(scopy, SCOPY) (n, &h__[(h_dim1 << 1) + 1], &c__1, &eig[1], &c__1);
4297     i__1 = *n - 1;
4298     F77_FUNC(scopy, SCOPY) (&i__1, &h__[h_dim1 + 2], &c__1, &workl[1], &c__1);
4299     F77_FUNC(sstqrb, SSTQRB) (n, &eig[1], &workl[1], &bounds[1], &workl[*n + 1], ierr);
4300     if (*ierr != 0)
4301     {
4302         goto L9000;
4303     }
4304
4305     i__1 = *n;
4306     for (k = 1; k <= i__1; ++k)
4307     {
4308         bounds[k] = *rnorm * std::abs(bounds[k]);
4309
4310     }
4311
4312
4313 L9000:
4314     return;
4315 }
4316
4317
4318
4319
4320 static void
4321 F77_FUNC(ssaitr, SSAITR) (int *     ido,
4322                           const char *    bmat,
4323                           int *     n,
4324                           int *     k,
4325                           int *     np,
4326                           int *     mode,
4327                           float *  resid,
4328                           float *  rnorm,
4329                           float *  v,
4330                           int *     ldv,
4331                           float *  h__,
4332                           int *     ldh,
4333                           int *     ipntr,
4334                           float *  workd,
4335                           int *     iwork,
4336                           int *     info)
4337 {
4338
4339     int   c__0  = 0;
4340     int   c__1  = 1;
4341     float c_b18 = 1.;
4342     float c_b42 = 0.;
4343     float c_b50 = -1.;
4344
4345     int   h_dim1, h_offset, v_dim1, v_offset, i__1;
4346     int   i__, jj;
4347     float temp1;
4348     int   infol;
4349     float safmin, minval;
4350
4351
4352     --workd;
4353     --resid;
4354     v_dim1   = *ldv;
4355     v_offset = 1 + v_dim1;
4356     v       -= v_offset;
4357     h_dim1   = *ldh;
4358     h_offset = 1 + h_dim1;
4359     h__     -= h_offset;
4360     --ipntr;
4361     --iwork;
4362     minval = GMX_FLOAT_MIN;
4363     safmin = minval / GMX_FLOAT_EPS;
4364
4365     if (*ido == 0)
4366     {
4367         *info    = 0;
4368         iwork[5] = 0;
4369         iwork[6] = 0;
4370         iwork[4] = 0;
4371         iwork[2] = 0;
4372         iwork[3] = 0;
4373
4374         iwork[12] = *k + 1;
4375
4376         iwork[8]  = 1;
4377         iwork[9]  = iwork[8] + *n;
4378         iwork[10] = iwork[9] + *n;
4379     }
4380
4381     if (iwork[5] == 1)
4382     {
4383         goto L50;
4384     }
4385     if (iwork[6] == 1)
4386     {
4387         goto L60;
4388     }
4389     if (iwork[2] == 1)
4390     {
4391         goto L70;
4392     }
4393     if (iwork[3] == 1)
4394     {
4395         goto L90;
4396     }
4397     if (iwork[4] == 1)
4398     {
4399         goto L30;
4400     }
4401
4402 L1000:
4403
4404
4405     if (*rnorm > 0.)
4406     {
4407         goto L40;
4408     }
4409
4410     iwork[11] = 1;
4411 L20:
4412     iwork[4] = 1;
4413     *ido     = 0;
4414 L30:
4415
4416     F77_FUNC(sgetv0, sgetv0) (ido, bmat, &iwork[11], &c__0, n, &iwork[12], &v[v_offset], ldv,
4417                               &resid[1], rnorm, &ipntr[1], &workd[1], &iwork[21], &iwork[7]);
4418     if (*ido != 99)
4419     {
4420         goto L9000;
4421     }
4422     if (iwork[7] < 0)
4423     {
4424         ++iwork[11];
4425         if (iwork[11] <= 3)
4426         {
4427             goto L20;
4428         }
4429
4430         *info = iwork[12] - 1;
4431         *ido  = 99;
4432         goto L9000;
4433     }
4434
4435 L40:
4436
4437     F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &v[iwork[12] * v_dim1 + 1], &c__1);
4438     if (*rnorm >= safmin)
4439     {
4440         temp1 = 1. / *rnorm;
4441         F77_FUNC(sscal, SSCAL) (n, &temp1, &v[iwork[12] * v_dim1 + 1], &c__1);
4442         F77_FUNC(sscal, SSCAL) (n, &temp1, &workd[iwork[8]], &c__1);
4443     }
4444     else
4445     {
4446
4447         F77_FUNC(slascl, SLASCL) ("General", &i__, &i__, rnorm, &c_b18, n, &c__1, &v[iwork[12] *
4448                                                                                      v_dim1 + 1], n, &infol);
4449         F77_FUNC(slascl, SLASCL) ("General", &i__, &i__, rnorm, &c_b18, n, &c__1, &workd[iwork[
4450                                                                                              8]], n, &infol);
4451     }
4452
4453     iwork[5] = 1;
4454     F77_FUNC(scopy, SCOPY) (n, &v[iwork[12] * v_dim1 + 1], &c__1, &workd[iwork[10]], &c__1);
4455     ipntr[1] = iwork[10];
4456     ipntr[2] = iwork[9];
4457     ipntr[3] = iwork[8];
4458     *ido     = 1;
4459
4460     goto L9000;
4461 L50:
4462
4463
4464     iwork[5] = 0;
4465
4466     F77_FUNC(scopy, SCOPY) (n, &workd[iwork[9]], &c__1, &resid[1], &c__1);
4467
4468     if (*mode == 2)
4469     {
4470         goto L65;
4471     }
4472     if (*bmat == 'G')
4473     {
4474         iwork[6] = 1;
4475         ipntr[1] = iwork[9];
4476         ipntr[2] = iwork[8];
4477         *ido     = 2;
4478
4479         goto L9000;
4480     }
4481     else if (*bmat == 'I')
4482     {
4483         F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
4484     }
4485 L60:
4486
4487     iwork[6] = 0;
4488
4489 L65:
4490     if (*mode == 2)
4491     {
4492
4493         workd[*n * 3 + 3] = F77_FUNC(sdot, SDOT) (n, &resid[1], &c__1, &workd[iwork[10]], &
4494                                                   c__1);
4495         workd[*n * 3 + 3] = std::sqrt(std::abs(workd[*n * 3 + 3]));
4496     }
4497     else if (*bmat == 'G')
4498     {
4499         workd[*n * 3 + 3] = F77_FUNC(sdot, SDOT) (n, &resid[1], &c__1, &workd[iwork[8]], &
4500                                                   c__1);
4501         workd[*n * 3 + 3] = std::sqrt(std::abs(workd[*n * 3 + 3]));
4502     }
4503     else if (*bmat == 'I')
4504     {
4505         workd[*n * 3 + 3] = F77_FUNC(snrm2, SNRM2) (n, &resid[1], &c__1);
4506     }
4507
4508     if (*mode != 2)
4509     {
4510         F77_FUNC(sgemv, SGEMV) ("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[8]],
4511                                 &c__1, &c_b42, &workd[iwork[9]], &c__1);
4512     }
4513     else
4514     {
4515         F77_FUNC(sgemv, SGEMV) ("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[10]
4516                                 ], &c__1, &c_b42, &workd[iwork[9]], &c__1);
4517     }
4518
4519     F77_FUNC(sgemv, SGEMV) ("N", n, &iwork[12], &c_b50, &v[v_offset], ldv, &workd[iwork[9]], &
4520                             c__1, &c_b18, &resid[1], &c__1);
4521
4522     h__[iwork[12] + (h_dim1 << 1)] = workd[iwork[9] + iwork[12] - 1];
4523     if (iwork[12] == 1 || iwork[4] == 1)
4524     {
4525         h__[iwork[12] + h_dim1] = 0.;
4526     }
4527     else
4528     {
4529         h__[iwork[12] + h_dim1] = *rnorm;
4530     }
4531
4532     iwork[2] = 1;
4533     iwork[1] = 0;
4534
4535     if (*bmat == 'G')
4536     {
4537         F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[iwork[9]], &c__1);
4538         ipntr[1] = iwork[9];
4539         ipntr[2] = iwork[8];
4540         *ido     = 2;
4541
4542         goto L9000;
4543     }
4544     else if (*bmat == 'I')
4545     {
4546         F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
4547     }
4548 L70:
4549
4550     iwork[2] = 0;
4551
4552     if (*bmat == 'G')
4553     {
4554         *rnorm = F77_FUNC(sdot, SDOT) (n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
4555         *rnorm = std::sqrt(std::abs(*rnorm));
4556     }
4557     else if (*bmat == 'I')
4558     {
4559         *rnorm = F77_FUNC(snrm2, SNRM2) (n, &resid[1], &c__1);
4560     }
4561
4562     if (*rnorm > workd[*n * 3 + 3] * .717F)
4563     {
4564         goto L100;
4565     }
4566
4567 L80:
4568
4569     F77_FUNC(sgemv, SGEMV) ("T", n, &iwork[12], &c_b18, &v[v_offset], ldv, &workd[iwork[8]], &
4570                             c__1, &c_b42, &workd[iwork[9]], &c__1);
4571
4572     F77_FUNC(sgemv, SGEMV) ("N", n, &iwork[12], &c_b50, &v[v_offset], ldv, &workd[iwork[9]], &
4573                             c__1, &c_b18, &resid[1], &c__1);
4574
4575     if (iwork[12] == 1 || iwork[4] == 1)
4576     {
4577         h__[iwork[12] + h_dim1] = 0.;
4578     }
4579     h__[iwork[12] + (h_dim1 << 1)] += workd[iwork[9] + iwork[12] - 1];
4580
4581     iwork[3] = 1;
4582     if (*bmat == 'G')
4583     {
4584         F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[iwork[9]], &c__1);
4585         ipntr[1] = iwork[9];
4586         ipntr[2] = iwork[8];
4587         *ido     = 2;
4588
4589         goto L9000;
4590     }
4591     else if (*bmat == 'I')
4592     {
4593         F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[iwork[8]], &c__1);
4594     }
4595 L90:
4596
4597
4598     if (*bmat == 'G')
4599     {
4600         workd[*n * 3 + 2] = F77_FUNC(sdot, SDOT) (n, &resid[1], &c__1, &workd[iwork[8]], &
4601                                                   c__1);
4602         workd[*n * 3 + 2] = std::sqrt(std::abs(workd[*n * 3 + 2]));
4603     }
4604     else if (*bmat == 'I')
4605     {
4606         workd[*n * 3 + 2] = F77_FUNC(snrm2, SNRM2) (n, &resid[1], &c__1);
4607     }
4608
4609
4610     if (workd[*n * 3 + 2] > *rnorm * .717F)
4611     {
4612
4613         *rnorm = workd[*n * 3 + 2];
4614
4615     }
4616     else
4617     {
4618
4619         *rnorm = workd[*n * 3 + 2];
4620         ++iwork[1];
4621         if (iwork[1] <= 1)
4622         {
4623             goto L80;
4624         }
4625
4626         i__1 = *n;
4627         for (jj = 1; jj <= i__1; ++jj)
4628         {
4629             resid[jj] = 0.;
4630         }
4631         *rnorm = 0.;
4632     }
4633
4634 L100:
4635
4636     iwork[4] = 0;
4637     iwork[3] = 0;
4638
4639     if (h__[iwork[12] + h_dim1] < 0.)
4640     {
4641         h__[iwork[12] + h_dim1] = -h__[iwork[12] + h_dim1];
4642         if (iwork[12] < *k + *np)
4643         {
4644             F77_FUNC(sscal, SSCAL) (n, &c_b50, &v[(iwork[12] + 1) * v_dim1 + 1], &c__1);
4645         }
4646         else
4647         {
4648             F77_FUNC(sscal, SSCAL) (n, &c_b50, &resid[1], &c__1);
4649         }
4650     }
4651
4652     ++iwork[12];
4653     if (iwork[12] > *k + *np)
4654     {
4655         *ido = 99;
4656
4657
4658         goto L9000;
4659     }
4660
4661     goto L1000;
4662
4663 L9000:
4664     return;
4665 }
4666
4667
4668
4669
4670
4671
4672 static void
4673 F77_FUNC(ssaup2, SSAUP2) (int *     ido,
4674                           const char *    bmat,
4675                           int *     n,
4676                           const char *    which,
4677                           int *     nev,
4678                           int *     np,
4679                           float *  tol,
4680                           float *  resid,
4681                           int *     mode,
4682                           int gmx_unused *     iupd,
4683                           int *     ishift,
4684                           int *     mxiter,
4685                           float *  v,
4686                           int *     ldv,
4687                           float *  h__,
4688                           int *     ldh,
4689                           float *  ritz,
4690                           float *  bounds,
4691                           float *  q,
4692                           int *     ldq,
4693                           float *  workl,
4694                           int *     ipntr,
4695                           float *  workd,
4696                           int *     iwork,
4697                           int *     info)
4698 {
4699     float c_b3 = 2/3.;
4700     int   c__1 = 1;
4701     int   c__0 = 0;
4702
4703     int   h_dim1, h_offset, q_dim1, q_offset, v_dim1, v_offset, i__1, i__2,
4704           i__3;
4705     float d__2, d__3;
4706     int   j;
4707     float eps23;
4708     int   ierr;
4709     float temp;
4710     int   nevd2;
4711     int   nevm2;
4712     int   nevbef;
4713     char  wprime[2];
4714     int   nptemp;
4715
4716     --workd;
4717     --resid;
4718     --workl;
4719     --bounds;
4720     --ritz;
4721     v_dim1   = *ldv;
4722     v_offset = 1 + v_dim1;
4723     v       -= v_offset;
4724     h_dim1   = *ldh;
4725     h_offset = 1 + h_dim1;
4726     h__     -= h_offset;
4727     q_dim1   = *ldq;
4728     q_offset = 1 + q_dim1;
4729     q       -= q_offset;
4730     --ipntr;
4731     --iwork;
4732     eps23 = GMX_FLOAT_EPS;
4733     eps23 = std::pow(eps23, c_b3);
4734
4735     if (*ido == 0)
4736     {
4737
4738         iwork[41] = 1;
4739         iwork[42] = 3;
4740         iwork[43] = 5;
4741         iwork[44] = 7;
4742
4743         iwork[9]  = *nev;
4744         iwork[10] = *np;
4745
4746         iwork[7] = iwork[9] + iwork[10];
4747         iwork[8] = 0;
4748         iwork[6] = 0;
4749
4750         iwork[2] = 1;
4751         iwork[4] = 0;
4752         iwork[5] = 0;
4753         iwork[1] = 0;
4754
4755         if (*info != 0)
4756         {
4757
4758             iwork[3] = 1;
4759             *info    = 0;
4760         }
4761         else
4762         {
4763             iwork[3] = 0;
4764         }
4765     }
4766
4767     if (iwork[2] == 1)
4768     {
4769         F77_FUNC(sgetv0, SGETV0) (ido, bmat, &c__1, &iwork[3], n, &c__1, &v[v_offset], ldv, &
4770                                   resid[1], &workd[*n * 3 + 1], &ipntr[1], &workd[1], &iwork[41],
4771                                   info);
4772
4773         if (*ido != 99)
4774         {
4775             goto L9000;
4776         }
4777
4778         if (workd[*n * 3 + 1] == 0.)
4779         {
4780
4781             *info = -9;
4782             goto L1200;
4783         }
4784         iwork[2] = 0;
4785         *ido     = 0;
4786     }
4787
4788     if (iwork[4] == 1)
4789     {
4790         goto L20;
4791     }
4792
4793     if (iwork[5] == 1)
4794     {
4795         goto L50;
4796     }
4797
4798     if (iwork[1] == 1)
4799     {
4800         goto L100;
4801     }
4802
4803     F77_FUNC(ssaitr, SSAITR) (ido, bmat, n, &c__0, &iwork[9], mode, &resid[1], &workd[*n * 3 +
4804                                                                                       1], &v[v_offset], ldv, &h__[h_offset], ldh, &ipntr[1], &workd[1],
4805                               &iwork[21], info);
4806
4807     if (*ido != 99)
4808     {
4809         goto L9000;
4810     }
4811
4812     if (*info > 0)
4813     {
4814
4815         *np     = *info;
4816         *mxiter = iwork[6];
4817         *info   = -9999;
4818         goto L1200;
4819     }
4820
4821 L1000:
4822
4823     ++iwork[6];
4824
4825
4826     *ido = 0;
4827 L20:
4828     iwork[4] = 1;
4829
4830     F77_FUNC(ssaitr, SSAITR) (ido, bmat, n, nev, np, mode, &resid[1], &workd[*n * 3 + 1], &v[
4831                                   v_offset], ldv, &h__[h_offset], ldh, &ipntr[1], &workd[1], &iwork[
4832                                   21], info);
4833
4834     if (*ido != 99)
4835     {
4836         goto L9000;
4837     }
4838
4839     if (*info > 0)
4840     {
4841
4842         *np     = *info;
4843         *mxiter = iwork[6];
4844         *info   = -9999;
4845         goto L1200;
4846     }
4847     iwork[4] = 0;
4848
4849     F77_FUNC(sseigt, SSEIGT) (&workd[*n * 3 + 1], &iwork[7], &h__[h_offset], ldh, &ritz[1], &
4850                               bounds[1], &workl[1], &ierr);
4851
4852     if (ierr != 0)
4853     {
4854         *info = -8;
4855         goto L1200;
4856     }
4857
4858     F77_FUNC(scopy, SCOPY) (&iwork[7], &ritz[1], &c__1, &workl[iwork[7] + 1], &c__1);
4859     F77_FUNC(scopy, SCOPY) (&iwork[7], &bounds[1], &c__1, &workl[(iwork[7] << 1) + 1], &c__1);
4860
4861     *nev = iwork[9];
4862     *np  = iwork[10];
4863     F77_FUNC(ssgets, SSGETS) (ishift, which, nev, np, &ritz[1], &bounds[1], &workl[1]);
4864
4865     F77_FUNC(scopy, SCOPY) (nev, &bounds[*np + 1], &c__1, &workl[*np + 1], &c__1);
4866     F77_FUNC(ssconv, SSCONV) (nev, &ritz[*np + 1], &workl[*np + 1], tol, &iwork[8]);
4867
4868
4869     nptemp = *np;
4870     i__1   = nptemp;
4871     for (j = 1; j <= i__1; ++j)
4872     {
4873         if (bounds[j] == 0.)
4874         {
4875             --(*np);
4876             ++(*nev);
4877         }
4878     }
4879
4880     if (iwork[8] >= iwork[9] || iwork[6] > *mxiter || *np == 0)
4881     {
4882
4883         if (!std::strncmp(which, "BE", 2))
4884         {
4885
4886             std::strncpy(wprime, "SA", 2);
4887             F77_FUNC(ssortr, SSORTR) (wprime, &c__1, &iwork[7], &ritz[1], &bounds[1]);
4888             nevd2 = *nev / 2;
4889             nevm2 = *nev - nevd2;
4890             if (*nev > 1)
4891             {
4892                 i__1 = (nevd2 < *np) ? nevd2 : *np;
4893                 i__2 = iwork[7] - nevd2 + 1, i__3 = iwork[7] - *np + 1;
4894                 F77_FUNC(sswap, SSWAP) (&i__1, &ritz[nevm2 + 1], &c__1,
4895                                         &ritz[((i__2 > i__3) ? i__2 : i__3)],
4896                                         &c__1);
4897                 i__1 = (nevd2 < *np) ? nevd2 : *np;
4898                 i__2 = iwork[7] - nevd2 + 1, i__3 = iwork[7] - *np;
4899                 F77_FUNC(sswap, SSWAP) (&i__1, &bounds[nevm2 + 1], &c__1,
4900                                         &bounds[((i__2 > i__3) ? i__2 : i__3) + 1],
4901                                         &c__1);
4902             }
4903
4904         }
4905         else
4906         {
4907
4908             if (!std::strncmp(which, "LM", 2))
4909             {
4910                 std::strncpy(wprime, "SM", 2);
4911             }
4912             if (!std::strncmp(which, "SM", 2))
4913             {
4914                 std::strncpy(wprime, "LM", 2);
4915             }
4916             if (!std::strncmp(which, "LA", 2))
4917             {
4918                 std::strncpy(wprime, "SA", 2);
4919             }
4920             if (!std::strncmp(which, "SA", 2))
4921             {
4922                 std::strncpy(wprime, "LA", 2);
4923             }
4924
4925             F77_FUNC(ssortr, SSORTR) (wprime, &c__1, &iwork[7], &ritz[1], &bounds[1]);
4926
4927         }
4928
4929         i__1 = iwork[9];
4930         for (j = 1; j <= i__1; ++j)
4931         {
4932             d__2       = eps23;
4933             d__3       = std::abs(ritz[j]);
4934             temp       = (d__2 > d__3) ? d__2 : d__3;
4935             bounds[j] /= temp;
4936         }
4937
4938         std::strncpy(wprime, "LA", 2);
4939         F77_FUNC(ssortr, SSORTR) (wprime, &c__1, &iwork[9], &bounds[1], &ritz[1]);
4940
4941         i__1 = iwork[9];
4942         for (j = 1; j <= i__1; ++j)
4943         {
4944             d__2       = eps23;
4945             d__3       = std::abs(ritz[j]);
4946             temp       = (d__2 > d__3) ? d__2 : d__3;
4947             bounds[j] *= temp;
4948         }
4949
4950         if (!std::strncmp(which, "BE", 2))
4951         {
4952
4953             std::strncpy(wprime, "LA", 2);
4954             F77_FUNC(ssortr, SSORTR) (wprime, &c__1, &iwork[8], &ritz[1], &bounds[1]);
4955
4956         }
4957         else
4958         {
4959             F77_FUNC(ssortr, SSORTR) (which, &c__1, &iwork[8], &ritz[1], &bounds[1]);
4960
4961         }
4962
4963         h__[h_dim1 + 1] = workd[*n * 3 + 1];
4964
4965
4966         if (iwork[6] > *mxiter && iwork[8] < *nev)
4967         {
4968             *info = 1;
4969         }
4970         if (*np == 0 && iwork[8] < iwork[9])
4971         {
4972             *info = 2;
4973         }
4974
4975         *np = iwork[8];
4976         goto L1100;
4977
4978     }
4979     else if (iwork[8] < *nev && *ishift == 1)
4980     {
4981         nevbef = *nev;
4982         i__1   = iwork[8], i__2 = *np / 2;
4983         *nev  += (i__1 < i__2) ? i__1 : i__2;
4984         if (*nev == 1 && iwork[7] >= 6)
4985         {
4986             *nev = iwork[7] / 2;
4987         }
4988         else if (*nev == 1 && iwork[7] > 2)
4989         {
4990             *nev = 2;
4991         }
4992         *np = iwork[7] - *nev;
4993
4994
4995         if (nevbef < *nev)
4996         {
4997             F77_FUNC(ssgets, SSGETS) (ishift, which, nev, np, &ritz[1], &bounds[1], &workl[1]);
4998         }
4999
5000     }
5001
5002
5003     if (*ishift == 0)
5004     {
5005
5006         iwork[5] = 1;
5007         *ido     = 3;
5008         goto L9000;
5009     }
5010
5011 L50:
5012
5013     iwork[5] = 0;
5014
5015     if (*ishift == 0)
5016     {
5017         F77_FUNC(scopy, SCOPY) (np, &workl[1], &c__1, &ritz[1], &c__1);
5018     }
5019
5020     F77_FUNC(ssapps, SSAPPS) (n, nev, np, &ritz[1], &v[v_offset], ldv, &h__[h_offset], ldh, &
5021                               resid[1], &q[q_offset], ldq, &workd[1]);
5022
5023     iwork[1] = 1;
5024     if (*bmat == 'G')
5025     {
5026         F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[*n + 1], &c__1);
5027         ipntr[1] = *n + 1;
5028         ipntr[2] = 1;
5029         *ido     = 2;
5030
5031         goto L9000;
5032     }
5033     else if (*bmat == 'I')
5034     {
5035         F77_FUNC(scopy, SCOPY) (n, &resid[1], &c__1, &workd[1], &c__1);
5036     }
5037
5038 L100:
5039
5040     if (*bmat == 'G')
5041     {
5042         workd[*n * 3 + 1] = F77_FUNC(sdot, SDOT) (n, &resid[1], &c__1, &workd[1], &c__1);
5043         workd[*n * 3 + 1] = std::sqrt(std::abs(workd[*n * 3 + 1]));
5044     }
5045     else if (*bmat == 'I')
5046     {
5047         workd[*n * 3 + 1] = F77_FUNC(snrm2, SNRM2) (n, &resid[1], &c__1);
5048     }
5049     iwork[1] = 0;
5050
5051     goto L1000;
5052
5053 L1100:
5054
5055     *mxiter = iwork[6];
5056     *nev    = iwork[8];
5057
5058 L1200:
5059     *ido = 99;
5060
5061 L9000:
5062     return;
5063
5064 }
5065
5066
5067
5068 void
5069 F77_FUNC(ssaupd, SSAUPD) (int *     ido,
5070                           const char *    bmat,
5071                           int *     n,
5072                           const char *      which,
5073                           int *     nev,
5074                           float *  tol,
5075                           float *  resid,
5076                           int *     ncv,
5077                           float *  v,
5078                           int *     ldv,
5079                           int *     iparam,
5080                           int *     ipntr,
5081                           float *  workd,
5082                           int *     iwork,
5083                           float *  workl,
5084                           int *     lworkl,
5085                           int *     info)
5086 {
5087     int v_dim1, v_offset, i__1, i__2;
5088     int j;
5089
5090     --workd;
5091     --resid;
5092     v_dim1   = *ldv;
5093     v_offset = 1 + v_dim1;
5094     v       -= v_offset;
5095     --iparam;
5096     --ipntr;
5097     --iwork;
5098     --workl;
5099
5100     if (*ido == 0)
5101     {
5102
5103
5104         iwork[2]  = 0;
5105         iwork[5]  = iparam[1];
5106         iwork[10] = iparam[3];
5107         iwork[12] = iparam[4];
5108
5109         iwork[6]  = 1;
5110         iwork[11] = iparam[7];
5111
5112
5113         if (*n <= 0)
5114         {
5115             iwork[2] = -1;
5116         }
5117         else if (*nev <= 0)
5118         {
5119             iwork[2] = -2;
5120         }
5121         else if (*ncv <= *nev || *ncv > *n)
5122         {
5123             iwork[2] = -3;
5124         }
5125
5126
5127         iwork[15] = *ncv - *nev;
5128
5129         if (iwork[10] <= 0)
5130         {
5131             iwork[2] = -4;
5132         }
5133         if (std::strncmp(which, "LM", 2) && std::strncmp(which, "SM", 2) &&
5134             std::strncmp(which, "LA", 2) && std::strncmp(which, "SA", 2) &&
5135             std::strncmp(which, "BE", 2))
5136         {
5137             iwork[2] = -5;
5138         }
5139         if (*bmat != 'I' && *bmat != 'G')
5140         {
5141             iwork[2] = -6;
5142         }
5143
5144         i__1 = *ncv;
5145         if (*lworkl < i__1 * i__1 + (*ncv << 3))
5146         {
5147             iwork[2] = -7;
5148         }
5149         if (iwork[11] < 1 || iwork[11] > 5)
5150         {
5151             iwork[2] = -10;
5152         }
5153         else if (iwork[11] == 1 && *bmat == 'G')
5154         {
5155             iwork[2] = -11;
5156         }
5157         else if (iwork[5] < 0 || iwork[5] > 1)
5158         {
5159             iwork[2] = -12;
5160         }
5161         else if (*nev == 1 && !std::strncmp(which, "BE", 2))
5162         {
5163             iwork[2] = -13;
5164         }
5165
5166         if (iwork[2] != 0)
5167         {
5168             *info = iwork[2];
5169             *ido  = 99;
5170             goto L9000;
5171         }
5172
5173         if (iwork[12] <= 0)
5174         {
5175             iwork[12] = 1;
5176         }
5177         if (*tol <= 0.)
5178         {
5179             *tol = GMX_FLOAT_EPS;
5180         }
5181
5182         iwork[15] = *ncv - *nev;
5183         iwork[13] = *nev;
5184         i__2      = *ncv;
5185         i__1      = i__2 * i__2 + (*ncv << 3);
5186         for (j = 1; j <= i__1; ++j)
5187         {
5188             workl[j] = 0.;
5189         }
5190
5191         iwork[8]  = *ncv;
5192         iwork[9]  = *ncv;
5193         iwork[3]  = 1;
5194         iwork[16] = iwork[3] + (iwork[8] << 1);
5195         iwork[1]  = iwork[16] + *ncv;
5196         iwork[4]  = iwork[1] + *ncv;
5197         i__1      = *ncv;
5198         iwork[7]  = iwork[4] + i__1 * i__1;
5199         iwork[14] = iwork[7] + *ncv * 3;
5200
5201         ipntr[4]  = iwork[14];
5202         ipntr[5]  = iwork[3];
5203         ipntr[6]  = iwork[16];
5204         ipntr[7]  = iwork[1];
5205         ipntr[11] = iwork[7];
5206     }
5207
5208     F77_FUNC(ssaup2, SSAUP2) (ido, bmat, n, which, &iwork[13], &iwork[15], tol, &resid[1], &
5209                               iwork[11], &iwork[6], &iwork[5], &iwork[10], &v[v_offset], ldv, &
5210                               workl[iwork[3]], &iwork[8], &workl[iwork[16]], &workl[iwork[1]], &
5211                               workl[iwork[4]], &iwork[9], &workl[iwork[7]], &ipntr[1], &workd[1],
5212                               &iwork[21], info);
5213
5214     if (*ido == 3)
5215     {
5216         iparam[8] = iwork[15];
5217     }
5218     if (*ido != 99)
5219     {
5220         goto L9000;
5221     }
5222
5223     iparam[3] = iwork[10];
5224     iparam[5] = iwork[15];
5225
5226     if (*info < 0)
5227     {
5228         goto L9000;
5229     }
5230     if (*info == 2)
5231     {
5232         *info = 3;
5233     }
5234
5235 L9000:
5236
5237     return;
5238
5239 }
5240
5241
5242
5243 void
5244 F77_FUNC(sseupd, SSEUPD) (int *     rvec,
5245                           const char *    howmny,
5246                           int *     select,
5247                           float *  d__,
5248                           float *  z__,
5249                           int *     ldz,
5250                           float *  sigma,
5251                           const char *    bmat,
5252                           int *     n,
5253                           const char *    which,
5254                           int *     nev,
5255                           float *  tol,
5256                           float *  resid,
5257                           int *     ncv,
5258                           float *  v,
5259                           int *     ldv,
5260                           int *     iparam,
5261                           int *     ipntr,
5262                           float *  workd,
5263                           float *  workl,
5264                           int *     lworkl,
5265                           int *     info)
5266 {
5267     float c_b21  = 2/3.;
5268     int   c__1   = 1;
5269     float c_b102 = 1.;
5270     int   v_dim1, v_offset, z_dim1, z_offset, i__1;
5271     float d__1, d__2, d__3;
5272
5273     int   j, k, ih, iq, iw, ibd, ihb, ihd, ldh, ilg, ldq, ism, irz;
5274     int   mode;
5275     float eps23;
5276     int   ierr;
5277     float temp;
5278     int   next;
5279     char  type__[6];
5280     int   ritz;
5281     int   reord;
5282     int   nconv;
5283     float rnorm;
5284     float bnorm2;
5285     float thres1 = 0, thres2 = 0;
5286     int   bounds;
5287     int   ktrord;
5288     float tempbnd;
5289     int   leftptr, rghtptr;
5290
5291
5292     --workd;
5293     --resid;
5294     z_dim1   = *ldz;
5295     z_offset = 1 + z_dim1;
5296     z__     -= z_offset;
5297     --d__;
5298     --select;
5299     v_dim1   = *ldv;
5300     v_offset = 1 + v_dim1;
5301     v       -= v_offset;
5302     --iparam;
5303     --ipntr;
5304     --workl;
5305
5306     mode  = iparam[7];
5307     nconv = iparam[5];
5308     *info = 0;
5309
5310     if (nconv == 0)
5311     {
5312         goto L9000;
5313     }
5314     ierr = 0;
5315
5316     if (nconv <= 0)
5317     {
5318         ierr = -14;
5319     }
5320     if (*n <= 0)
5321     {
5322         ierr = -1;
5323     }
5324     if (*nev <= 0)
5325     {
5326         ierr = -2;
5327     }
5328     if (*ncv <= *nev || *ncv > *n)
5329     {
5330         ierr = -3;
5331     }
5332     if (std::strncmp(which, "LM", 2) && std::strncmp(which, "SM", 2) &&
5333         std::strncmp(which, "LA", 2) && std::strncmp(which, "SA", 2) &&
5334         std::strncmp(which, "BE", 2))
5335     {
5336         ierr = -5;
5337     }
5338     if (*bmat != 'I' && *bmat != 'G')
5339     {
5340         ierr = -6;
5341     }
5342     if (*howmny != 'A' && *howmny != 'P' &&
5343         *howmny != 'S' && *rvec)
5344     {
5345         ierr = -15;
5346     }
5347     if (*rvec && *howmny == 'S')
5348     {
5349         ierr = -16;
5350     }
5351     i__1 = *ncv;
5352     if (*rvec && *lworkl < i__1 * i__1 + (*ncv << 3))
5353     {
5354         ierr = -7;
5355     }
5356
5357     if (mode == 1 || mode == 2)
5358     {
5359         std::strncpy(type__, "REGULR", 6);
5360     }
5361     else if (mode == 3)
5362     {
5363         std::strncpy(type__, "SHIFTI", 6);
5364     }
5365     else if (mode == 4)
5366     {
5367         std::strncpy(type__, "BUCKLE", 6);
5368     }
5369     else if (mode == 5)
5370     {
5371         std::strncpy(type__, "CAYLEY", 6);
5372     }
5373     else
5374     {
5375         ierr = -10;
5376     }
5377     if (mode == 1 && *bmat == 'G')
5378     {
5379         ierr = -11;
5380     }
5381     if (*nev == 1 && !std::strncmp(which, "BE", 2))
5382     {
5383         ierr = -12;
5384     }
5385
5386     if (ierr != 0)
5387     {
5388         *info = ierr;
5389         goto L9000;
5390     }
5391
5392     ih        = ipntr[5];
5393     ritz      = ipntr[6];
5394     bounds    = ipntr[7];
5395     ldh       = *ncv;
5396     ldq       = *ncv;
5397     ihd       = bounds + ldh;
5398     ihb       = ihd + ldh;
5399     iq        = ihb + ldh;
5400     iw        = iq + ldh * *ncv;
5401     next      = iw + (*ncv << 1);
5402     ipntr[4]  = next;
5403     ipntr[8]  = ihd;
5404     ipntr[9]  = ihb;
5405     ipntr[10] = iq;
5406
5407     irz = ipntr[11] + *ncv;
5408     ibd = irz + *ncv;
5409
5410
5411     eps23 = GMX_FLOAT_EPS;
5412     eps23 = std::pow(eps23, c_b21);
5413
5414     rnorm = workl[ih];
5415     if (*bmat == 'I')
5416     {
5417         bnorm2 = rnorm;
5418     }
5419     else if (*bmat == 'G')
5420     {
5421         bnorm2 = F77_FUNC(snrm2, SNRM2) (n, &workd[1], &c__1);
5422     }
5423
5424     if (*rvec)
5425     {
5426
5427         if (!std::strncmp(which, "LM", 2) || !std::strncmp(which, "SM", 2) ||
5428             !std::strncmp(which, "LA", 2) || !std::strncmp(which, "SA", 2))
5429         {
5430
5431         }
5432         else if (!std::strncmp(which, "BE", 2))
5433         {
5434
5435
5436             ism    = (*nev > nconv) ? *nev : nconv;
5437             ism   /= 2;
5438             ilg    = ism + 1;
5439             thres1 = workl[ism];
5440             thres2 = workl[ilg];
5441
5442
5443         }
5444
5445         reord  = 0;
5446         ktrord = 0;
5447         i__1   = *ncv - 1;
5448         for (j = 0; j <= i__1; ++j)
5449         {
5450             select[j + 1] = 0;
5451             if (!std::strncmp(which, "LM", 2))
5452             {
5453                 if (std::abs(workl[irz + j]) >= std::abs(thres1))
5454                 {
5455                     d__2    = eps23;
5456                     d__3    = std::abs(workl[irz + j]);
5457                     tempbnd = (d__2 > d__3) ? d__2 : d__3;
5458                     if (workl[ibd + j] <= *tol * tempbnd)
5459                     {
5460                         select[j + 1] = 1;
5461                     }
5462                 }
5463             }
5464             else if (!std::strncmp(which, "SM", 2))
5465             {
5466                 if (std::abs(workl[irz + j]) <= std::abs(thres1))
5467                 {
5468                     d__2    = eps23;
5469                     d__3    = std::abs(workl[irz + j]);
5470                     tempbnd = (d__2 > d__3) ? d__2 : d__3;
5471                     if (workl[ibd + j] <= *tol * tempbnd)
5472                     {
5473                         select[j + 1] = 1;
5474                     }
5475                 }
5476             }
5477             else if (!std::strncmp(which, "LA", 2))
5478             {
5479                 if (workl[irz + j] >= thres1)
5480                 {
5481                     d__2    = eps23;
5482                     d__3    = std::abs(workl[irz + j]);
5483                     tempbnd = (d__2 > d__3) ? d__2 : d__3;
5484                     if (workl[ibd + j] <= *tol * tempbnd)
5485                     {
5486                         select[j + 1] = 1;
5487                     }
5488                 }
5489             }
5490             else if (!std::strncmp(which, "SA", 2))
5491             {
5492                 if (workl[irz + j] <= thres1)
5493                 {
5494                     d__2    = eps23;
5495                     d__3    = std::abs(workl[irz + j]);
5496                     tempbnd = (d__2 > d__3) ? d__2 : d__3;
5497                     if (workl[ibd + j] <= *tol * tempbnd)
5498                     {
5499                         select[j + 1] = 1;
5500                     }
5501                 }
5502             }
5503             else if (!std::strncmp(which, "BE", 2))
5504             {
5505                 if (workl[irz + j] <= thres1 || workl[irz + j] >= thres2)
5506                 {
5507                     d__2    = eps23;
5508                     d__3    = std::abs(workl[irz + j]);
5509                     tempbnd = (d__2 > d__3) ? d__2 : d__3;
5510                     if (workl[ibd + j] <= *tol * tempbnd)
5511                     {
5512                         select[j + 1] = 1;
5513                     }
5514                 }
5515             }
5516             if (j + 1 > nconv)
5517             {
5518                 reord = select[j + 1] || reord;
5519             }
5520             if (select[j + 1])
5521             {
5522                 ++ktrord;
5523             }
5524         }
5525
5526         i__1 = *ncv - 1;
5527         F77_FUNC(scopy, SCOPY) (&i__1, &workl[ih + 1], &c__1, &workl[ihb], &c__1);
5528         F77_FUNC(scopy, SCOPY) (ncv, &workl[ih + ldh], &c__1, &workl[ihd], &c__1);
5529
5530         F77_FUNC(ssteqr, SSTEQR) ("Identity", ncv, &workl[ihd], &workl[ihb], &workl[iq], &ldq, &
5531                                   workl[iw], &ierr);
5532
5533         if (ierr != 0)
5534         {
5535             *info = -8;
5536             goto L9000;
5537         }
5538
5539
5540         if (reord)
5541         {
5542
5543             leftptr = 1;
5544             rghtptr = *ncv;
5545
5546             if (*ncv == 1)
5547             {
5548                 goto L30;
5549             }
5550
5551 L20:
5552             if (select[leftptr])
5553             {
5554
5555                 ++leftptr;
5556
5557             }
5558             else if (!select[rghtptr])
5559             {
5560
5561                 --rghtptr;
5562
5563             }
5564             else
5565             {
5566
5567                 temp                     = workl[ihd + leftptr - 1];
5568                 workl[ihd + leftptr - 1] = workl[ihd + rghtptr - 1];
5569                 workl[ihd + rghtptr - 1] = temp;
5570                 F77_FUNC(scopy, SCOPY) (ncv, &workl[iq + *ncv * (leftptr - 1)], &c__1, &workl[
5571                                             iw], &c__1);
5572                 F77_FUNC(scopy, SCOPY) (ncv, &workl[iq + *ncv * (rghtptr - 1)], &c__1, &workl[
5573                                             iq + *ncv * (leftptr - 1)], &c__1);
5574                 F77_FUNC(scopy, SCOPY) (ncv, &workl[iw], &c__1, &workl[iq + *ncv * (rghtptr -
5575                                                                                     1)], &c__1);
5576                 ++leftptr;
5577                 --rghtptr;
5578
5579             }
5580
5581             if (leftptr < rghtptr)
5582             {
5583                 goto L20;
5584             }
5585
5586 L30:
5587             ;
5588         }
5589
5590         F77_FUNC(scopy, SCOPY) (&nconv, &workl[ihd], &c__1, &d__[1], &c__1);
5591
5592     }
5593     else
5594     {
5595
5596         F77_FUNC(scopy, SCOPY) (&nconv, &workl[ritz], &c__1, &d__[1], &c__1);
5597         F77_FUNC(scopy, SCOPY) (ncv, &workl[ritz], &c__1, &workl[ihd], &c__1);
5598
5599     }
5600     if (!std::strncmp(type__, "REGULR", 6))
5601     {
5602
5603         if (*rvec)
5604         {
5605             F77_FUNC(ssesrt, SSESRT) ("LA", rvec, &nconv, &d__[1], ncv, &workl[iq], &ldq);
5606         }
5607         else
5608         {
5609             F77_FUNC(scopy, SCOPY) (ncv, &workl[bounds], &c__1, &workl[ihb], &c__1);
5610         }
5611
5612     }
5613     else
5614     {
5615
5616         F77_FUNC(scopy, SCOPY) (ncv, &workl[ihd], &c__1, &workl[iw], &c__1);
5617         if (!std::strncmp(type__, "SHIFTI", 6))
5618         {
5619             i__1 = *ncv;
5620             for (k = 1; k <= i__1; ++k)
5621             {
5622                 workl[ihd + k - 1] = 1. / workl[ihd + k - 1] + *sigma;
5623             }
5624         }
5625         else if (!std::strncmp(type__, "BUCKLE", 6))
5626         {
5627             i__1 = *ncv;
5628             for (k = 1; k <= i__1; ++k)
5629             {
5630                 workl[ihd + k - 1] = *sigma * workl[ihd + k - 1] / (workl[ihd
5631                                                                           + k - 1] - 1.);
5632             }
5633         }
5634         else if (!std::strncmp(type__, "CAYLEY", 6))
5635         {
5636             i__1 = *ncv;
5637             for (k = 1; k <= i__1; ++k)
5638             {
5639                 workl[ihd + k - 1] = *sigma * (workl[ihd + k - 1] + 1.) / (
5640                         workl[ihd + k - 1] - 1.);
5641             }
5642         }
5643
5644         F77_FUNC(scopy, SCOPY) (&nconv, &workl[ihd], &c__1, &d__[1], &c__1);
5645         F77_FUNC(ssortr, SSORTR) ("LA", &c__1, &nconv, &workl[ihd], &workl[iw]);
5646         if (*rvec)
5647         {
5648             F77_FUNC(ssesrt, SSESRT) ("LA", rvec, &nconv, &d__[1], ncv, &workl[iq], &ldq);
5649         }
5650         else
5651         {
5652             F77_FUNC(scopy, SCOPY) (ncv, &workl[bounds], &c__1, &workl[ihb], &c__1);
5653             d__1 = bnorm2 / rnorm;
5654             F77_FUNC(sscal, SSCAL) (ncv, &d__1, &workl[ihb], &c__1);
5655             F77_FUNC(ssortr, SSORTR) ("LA", &c__1, &nconv, &d__[1], &workl[ihb]);
5656         }
5657
5658     }
5659
5660     if (*rvec && *howmny == 'A')
5661     {
5662
5663         F77_FUNC(sgeqr2, SGEQR2) (ncv, &nconv, &workl[iq], &ldq, &workl[iw + *ncv], &workl[ihb],
5664                                   &ierr);
5665
5666         F77_FUNC(sorm2r, SORM2R) ("Right", "Notranspose", n, ncv, &nconv, &workl[iq], &ldq, &
5667                                   workl[iw + *ncv], &v[v_offset], ldv, &workd[*n + 1], &ierr);
5668         F77_FUNC(slacpy, SLACPY) ("All", n, &nconv, &v[v_offset], ldv, &z__[z_offset], ldz);
5669
5670         i__1 = *ncv - 1;
5671         for (j = 1; j <= i__1; ++j)
5672         {
5673             workl[ihb + j - 1] = 0.;
5674         }
5675         workl[ihb + *ncv - 1] = 1.;
5676         F77_FUNC(sorm2r, SORM2R) ("Left", "Transpose", ncv, &c__1, &nconv, &workl[iq], &ldq, &
5677                                   workl[iw + *ncv], &workl[ihb], ncv, &temp, &ierr);
5678
5679     }
5680     else if (*rvec && *howmny == 'S')
5681     {
5682
5683     }
5684
5685     if (!std::strncmp(type__, "REGULR", 6) && *rvec)
5686     {
5687
5688         i__1 = *ncv;
5689         for (j = 1; j <= i__1; ++j)
5690         {
5691             workl[ihb + j - 1] = rnorm * std::abs(workl[ihb + j - 1]);
5692         }
5693
5694     }
5695     else if (std::strncmp(type__, "REGULR", 6) && *rvec)
5696     {
5697
5698         F77_FUNC(sscal, SSCAL) (ncv, &bnorm2, &workl[ihb], &c__1);
5699         if (!std::strncmp(type__, "SHIFTI", 6))
5700         {
5701
5702             i__1 = *ncv;
5703             for (k = 1; k <= i__1; ++k)
5704             {
5705                 d__2               = workl[iw + k - 1];
5706                 workl[ihb + k - 1] = std::abs(workl[ihb + k - 1])/(d__2 * d__2);
5707             }
5708
5709         }
5710         else if (!std::strncmp(type__, "BUCKLE", 6))
5711         {
5712
5713             i__1 = *ncv;
5714             for (k = 1; k <= i__1; ++k)
5715             {
5716                 d__2               = workl[iw + k - 1] - 1.;
5717                 workl[ihb + k - 1] = *sigma * std::abs(workl[ihb + k - 1])/(d__2 * d__2);
5718             }
5719
5720         }
5721         else if (!std::strncmp(type__, "CAYLEY", 6))
5722         {
5723
5724             i__1 = *ncv;
5725             for (k = 1; k <= i__1; ++k)
5726             {
5727                 workl[ihb + k - 1] = std::abs(workl[ihb + k - 1] / workl[iw + k - 1] * (workl[iw + k - 1] - 1.));
5728
5729             }
5730
5731         }
5732
5733     }
5734
5735     if (*rvec && (!std::strncmp(type__, "SHIFTI", 6) || !std::strncmp(type__, "CAYLEY", 6)))
5736     {
5737
5738         i__1 = nconv - 1;
5739         for (k = 0; k <= i__1; ++k)
5740         {
5741             workl[iw + k] = workl[iq + k * ldq + *ncv - 1] / workl[iw + k];
5742         }
5743
5744     }
5745     else if (*rvec && !std::strncmp(type__, "BUCKLE", 6))
5746     {
5747
5748         i__1 = nconv - 1;
5749         for (k = 0; k <= i__1; ++k)
5750         {
5751             workl[iw + k] = workl[iq + k * ldq + *ncv - 1] / (workl[iw + k] -
5752                                                               1.);
5753         }
5754
5755     }
5756
5757     if (std::strncmp(type__, "REGULR", 6))
5758     {
5759         F77_FUNC(sger, SGER) (n, &nconv, &c_b102, &resid[1], &c__1, &workl[iw], &c__1, &z__[
5760                                   z_offset], ldz);
5761     }
5762
5763 L9000:
5764
5765     return;
5766
5767 }