src/gromacs/linearalgebra/gmx_arpack.c

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