?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Merge branch 'master' of git://git.gromacs.org/gromacs
[gromacs/adressmacs.git] / src / gmxlib / gmx_arpack.c
blob3a4e4491195c0c61c3171fca8160e295f880acbd
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
got merge a long time ago