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