src/gromacs/linearalgebra/gmx_arpack.c

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