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