src/gmxlib/gmx_arpack.c

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