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