src/gromacs/linearalgebra/gmx_arpack.cpp

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