src/gmxlib/gmx_arpack.c

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