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