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