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