src/gromacs/linearalgebra/gmx_arpack.c

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