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