1 /* Homotopies.c */
2
3 /*
4 This is the implementation code from the various files that used
5 to reside in the Cont subdirectory of Pelican0.80/source.
6 */
7
8 #include "pelhomot.h"
9
10 /**********************************************************************/
11 /*************** implementation code from Hom_params.c ****************/
12 /**********************************************************************/
13
14 double Hom_tol=.000001; /* tolerence for path tracking */
15
print_Hom_params(FILE * outfile)16 void print_Hom_params(FILE *outfile){
17 #ifdef HOM_PRINT
18 fprintf(outfile," Hominuation Parameters:\n");
19 fprintf(outfile," Hom_tol=%g\n",Hom_tol)
20 #endif
21 ;
22 }
23
24 /* end Hom_params.c */
25
26 /************************************************************************/
27 /***************************** code from fixpnf.c ***********************/
28 /************************************************************************/
29
30 #define TRUE_ (1)
31 #define FALSE_ (0)
32 double get_abs_homog();
33 double coord_r, coord_i;
34
35 /* Table of constant values */
36
37 static integer c__4 = 4;
38 static integer c__1 = 1;
39
fixpnf_0_(int n__,int * n,double * y,int * iflag,double * arcre,double * arcae,double * ansre,double * ansae,int * trace,double * a,int * nfe,double * arclen,double * yp,double * yold,double * ypold,double * qr,double * alpha,double * tz,int * pivot,double * w,double * wp,double * z0,double * z1,double * sspar,double * par,int * ipar,int tweak)40 /* Subroutine */ int fixpnf_0_(int n__,
41 int *n,
42 double *y,
43 int *iflag,
44 double *arcre,
45 double *arcae,
46 double *ansre,
47 double *ansae,
48 int *trace,
49 double *a,
50 int *nfe,
51 double *arclen,
52 double *yp,
53 double *yold,
54 double *ypold,
55 double *qr,
56 double *alpha,
57 double *tz,
58 int *pivot,
59 double *w,
60 double *wp,
61 double *z0,
62 double *z1,
63 double *sspar,
64 double *par,
65 int *ipar,
66
67 int tweak) /* added to adjust step size */
68 {
69 /* System generated locals */
70 int qr_dim1, qr_offset, i__1, i__2;
71 double d__1 /* ,abx UNUSED */;
72
73 /* Builtin functions */
74 /* double sqrt(double); CANT DECLARE BUILTINS UNDER C++ */
75 // integer s_wsfe(), do_fio(), e_wsfe(); /* in Hom_params.c these are int's */
76
77 /* Local variables */
78 static int nfec;
79 static double hold;
80 static int iter;
81 extern double dnrm2_(integer *n,
82 doublereal *dx,
83 integer *incx);
84 static double h, s;
85 static long int crash;
86 static int limit;
87 extern double d1mach_(integer *);
88 static long int start;
89 static int nc, iflagc, jw;
90 static double abserr, relerr;
91 extern /* Subroutine */ int stepnf_(integer *n,
92 integer *nfe,
93 integer *iflag,
94 logical *start,
95 logical *crash,
96 doublereal *hold,
97 doublereal *h,
98 doublereal *relerr,
99 doublereal *abserr,
100 doublereal *s,
101 doublereal *y,
102 doublereal *yp,
103 doublereal *yold,
104 doublereal *ypold,
105 doublereal *a,
106 doublereal *qr,
107 doublereal *alpha,
108 doublereal *tz,
109 integer *pivot,
110 doublereal *w,
111 doublereal *wp,
112 doublereal *z0,
113 doublereal *z1,
114 doublereal *sspar,
115 doublereal *par,
116 integer *ipar);
117
118 static double curtol;
119
120 extern /* Subroutine */ int rootnf_(int *n,
121 int *nfe,
122 int *iflag,
123 double *relerr,
124 double *abserr,
125 double *y,
126 double *yp,
127 double *yold,
128 double *ypold,
129 double *a,
130 double *qr,
131 double *alpha,
132 double *tz,
133 int *pivot,
134 double *w,
135 double *wp,
136 double *par,
137 int *ipar);
138 static int np1;
139 static long int polsys;
140
141
142
143 /* SUBROUTINE FIXPNF FINDS A FIXED POINT OR ZERO OF THE */
144 /* N-DIMENSIONAL VECTOR FUNCTION F(X), OR TRACKS A ZERO CURVE */
145 /* OF A GENERAL HOMOTOPY MAP RHO(A,LAMBDA,X). FOR THE FIXED */
146 /* POINT PROBLEM F(X) IS ASSUMED TO BE A C2 MAP OF SOME BALL */
147 /* INTO ITSELF. THE EQUATION X = F(X) IS SOLVED BY */
148 /* FOLLOWING THE ZERO CURVE OF THE HOMOTOPY MAP */
149
150 /* LAMBDA*(X - F(X)) + (1 - LAMBDA)*(X - A) , */
151
152 /* STARTING FROM LAMBDA = 0, X = A. THE CURVE IS PARAMETERIZED */
153 /* BY ARC LENGTH S, AND IS FOLLOWED BY SOLVING THE ORDINARY */
154 /* DIFFERENTIAL EQUATION D(HOMOTOPY MAP)/DS = 0 FOR */
155 /* Y(S) = (LAMBDA(S), X(S)) USING A HERMITE CUBIC PREDICTOR AND A */
156 /* CORRECTOR WHICH RETURNS TO THE ZERO CURVE ALONG THE FLOW NORMAL */
157 /* TO THE DAVIDENKO FLOW (WHICH CONSISTS OF THE INTEGRAL CURVES OF */
158 /* D(HOMOTOPY MAP)/DS ). */
159
160 /* FOR THE ZERO FINDING PROBLEM F(X) IS ASSUMED TO BE A C2 MAP */
161 /* SUCH THAT FOR SOME R > 0, X*F(X) >= 0 WHENEVER NORM(X) = R. */
162 /* THE EQUATION F(X) = 0 IS SOLVED BY FOLLOWING THE ZERO CURVE */
163 /* OF THE HOMOTOPY MAP */
164
165 /* LAMBDA*F(X) + (1 - LAMBDA)*(X - A) */
166
167 /* EMANATING FROM LAMBDA = 0, X = A. */
168
169 /* A MUST BE AN INTERIOR POINT OF THE ABOVE MENTIONED BALLS. */
170
171 /* FOR THE CURVE TRACKING PROBLEM RHO(A,LAMBDA,X) IS ASSUMED TO */
172 /* BE A C2 MAP FROM E**M X [0,1) X E**N INTO E**N, WHICH FOR */
173 /* ALMOST ALL PARAMETER VECTORS A IN SOME NONEMPTY OPEN SUBSET */
174 /* OF E**M SATISFIES */
175
176 /* RANK [D RHO(A,LAMBDA,X)/D LAMBDA , D RHO(A,LAMBDA,X)/DX] = N */
177
178 /* FOR ALL POINTS (LAMBDA,X) SUCH THAT RHO(A,LAMBDA,X)=0. IT IS */
179 /* FURTHER ASSUMED THAT */
180
181 /* RANK [ D RHO(A,0,X0)/DX ] = N . */
182
183 /* WITH A FIXED, THE ZERO CURVE OF RHO(A,LAMBDA,X) EMANATING */
184 /* FROM LAMBDA = 0, X = X0 IS TRACKED UNTIL LAMBDA = 1 BY */
185 /* SOLVING THE ORDINARY DIFFERENTIAL EQUATION */
186 /* D RHO(A,LAMBDA(S),X(S))/DS = 0 FOR Y(S) = (LAMBDA(S), X(S)), */
187 /* WHERE S IS ARC LENGTH ALONG THE ZERO CURVE. ALSO THE HOMOTOPY */
188 /* MAP RHO(A,LAMBDA,X) IS ASSUMED TO BE CONSTRUCTED SUCH THAT */
189
190 /* D LAMBDA(0)/DS > 0 . */
191
192
193 /* FOR THE FIXED POINT AND ZERO FINDING PROBLEMS, THE USER MUST SUPPLY */
194 /* A SUBROUTINE F(X,V) WHICH EVALUATES F(X) AT X AND RETURNS THE */
195 /* VECTOR F(X) IN V, AND A SUBROUTINE FJAC(X,V,K) WHICH RETURNS IN V */
196 /* THE KTH COLUMN OF THE JACOBIAN MATRIX OF F(X) EVALUATED AT X. FOR */
197 /* THE CURVE TRACKING PROBLEM, THE USER MUST SUPPLY A SUBROUTINE */
198 /* RHO(A,LAMBDA,X,V,PAR,IPAR) WHICH EVALUATES THE HOMOTOPY MAP RHO AT */
199 /* (A,LAMBDA,X) AND RETURNS THE VECTOR RHO(A,LAMBDA,X) IN V, AND A */
200 /* SUBROUTINE RHOJAC(A,LAMBDA,X,V,K,PAR,IPAR) WHICH RETURNS IN V THE KTH
201 */
202 /* COLUMN OF THE N X (N+1) JACOBIAN MATRIX [D RHO/D LAMBDA, D RHO/DX] */
203 /* EVALUATED AT (A,LAMBDA,X). FIXPNF DIRECTLY OR INDIRECTLY USES */
204 /* THE SUBROUTINES STEPNF , TANGNF , ROOTNF , ROOT , F (OR RHO ), */
205 /* FJAC (OR RHOJAC ), D1MACH , AND THE BLAS FUNCTIONS DDOT AND */
206 /* DNRM2 . ONLY D1MACH CONTAINS MACHINE DEPENDENT CONSTANTS. */
207 /* NO OTHER MODIFICATIONS BY THE USER ARE REQUIRED. */
208
209
210 /* ON INPUT: */
211
212 /* N IS THE DIMENSION OF X, F(X), AND RHO(A,LAMBDA,X). */
213
214 /* Y IS AN ARRRAY OF LENGTH N + 1. (Y(2),...,Y(N+1)) = A IS THE */
215 /* STARTING POINT FOR THE ZERO CURVE FOR THE FIXED POINT AND */
216 /* ZERO FINDING PROBLEMS. (Y(2),...,Y(N+1)) = X0 FOR THE CURVE */
217 /* TRACKING PROBLEM. */
218
219 /* IFLAG CAN BE -2, -1, 0, 2, OR 3. IFLAG SHOULD BE 0 ON THE */
220 /* FIRST CALL TO FIXPNF FOR THE PROBLEM X=F(X), -1 FOR THE */
221 /* PROBLEM F(X)=0, AND -2 FOR THE PROBLEM RHO(A,LAMBDA,X)=0. */
222 /* IN CERTAIN SITUATIONS IFLAG IS SET TO 2 OR 3 BY FIXPNF, */
223 /* AND FIXPNF CAN BE CALLED AGAIN WITHOUT CHANGING IFLAG. */
224
225 /* ARCRE , ARCAE ARE THE RELATIVE AND ABSOLUTE ERRORS, RESPECTIVELY, */
226 /* ALLOWED THE NORMAL FLOW ITERATION ALONG THE ZERO CURVE. IF */
227 /* ARC?E .LE. 0.0 ON INPUT IT IS RESET TO .5*SQRT(ANS?E) . */
228 /* NORMALLY ARC?E SHOULD BE CONSIDERABLY LARGER THAN ANS?E . */
229
230 /* ANSRE , ANSAE ARE THE RELATIVE AND ABSOLUTE ERROR VALUES USED FOR */
231 /* THE ANSWER AT LAMBDA = 1. THE ACCEPTED ANSWER Y = (LAMBDA, X) */
232 /* SATISFIES */
233
234 /* |Y(1) - 1| .LE. ANSRE + ANSAE .AND. */
235
236 /* ||Z|| .LE. ANSRE*||X|| + ANSAE WHERE */
237
238 /* (.,Z) IS THE NEWTON STEP TO Y. */
239
240 /* TRACE IS AN INTEGER SPECIFYING THE LOGICAL I/O UNIT FOR */
241 /* INTERMEDIATE OUTPUT. IF TRACE .GT. 0 THE POINTS COMPUTED ON */
242 /* THE ZERO CURVE ARE WRITTEN TO I/O UNIT TRACE . */
243
244 /* A(1:*) CONTAINS THE PARAMETER VECTOR A . FOR THE FIXED POINT */
245 /* AND ZERO FINDING PROBLEMS, A NEED NOT BE INITIALIZED BY THE */
246 /* USER, AND IS ASSUMED TO HAVE LENGTH N. FOR THE CURVE */
247 /* TRACKING PROBLEM, A MUST BE INITIALIZED BY THE USER. */
248
249 /* YP(1:N+1) IS A WORK ARRAY CONTAINING THE TANGENT VECTOR TO */
250 /* THE ZERO CURVE AT THE CURRENT POINT Y . */
251
252 /* YOLD(1:N+1) IS A WORK ARRAY CONTAINING THE PREVIOUS POINT FOUND */
253 /* ON THE ZERO CURVE. */
254
255 /* YPOLD(1:N+1) IS A WORK ARRAY CONTAINING THE TANGENT VECTOR TO */
256 /* THE ZERO CURVE AT YOLD . */
257
258 /* QR(1:N,1:N+2), ALPHA(1:N), TZ(1:N+1), PIVOT(1:N+1) , W(1:N+1) , */
259 /* WP(1:N+1) , Z0(1:N+1) , Z1(1:N+1) ARE ALL WORK ARRAYS USED BY */
260 /* STEPNF TO CALCULATE THE TANGENT VECTORS AND NEWTON STEPS. */
261
262 /* SSPAR(1:8) = (LIDEAL, RIDEAL, DIDEAL, HMIN, HMAX, BMIN, BMAX, P) IS */
263 /* A VECTOR OF PARAMETERS USED FOR THE OPTIMAL STEP SIZE ESTIMATION. */
264 /* IF SSPAR(J) .LE. 0.0 ON INPUT, IT IS RESET TO A DEFAULT VALUE */
265 /* BY FIXPNF . OTHERWISE THE INPUT VALUE OF SSPAR(J) IS USED. */
266 /* SEE THE COMMENTS BELOW AND IN STEPNF FOR MORE INFORMATION ABOUT */
267 /* THESE CONSTANTS. */
268
269 /* PAR(1:*) AND IPAR(1:*) ARE ARRAYS FOR (OPTIONAL) USER PARAMETERS, */
270 /* WHICH ARE SIMPLY PASSED THROUGH TO THE USER WRITTEN SUBROUTINES */
271 /* RHO, RHOJAC. */
272
273
274 /* ON OUTPUT: */
275
276 /* N , TRACE , A ARE UNCHANGED. */
277
278 /* Y(1) = LAMBDA, (Y(2),...,Y(N+1)) = X, AND Y IS AN APPROXIMATE */
279 /* ZERO OF THE HOMOTOPY MAP. NORMALLY LAMBDA = 1 AND X IS A */
280 /* FIXED POINT(ZERO) OF F(X). IN ABNORMAL SITUATIONS LAMBDA */
281 /* MAY ONLY BE NEAR 1 AND X IS NEAR A FIXED POINT(ZERO). */
282
283 /* IFLAG = */
284 /* -2 CAUSES FIXPNF TO INITIALIZE EVERYTHING FOR THE PROBLEM */
285 /* RHO(A,LAMBDA,X) = 0 (USE ON FIRST CALL). */
286
287 /* -1 CAUSES FIXPNF TO INITIALIZE EVERYTHING FOR THE PROBLEM */
288 /* F(X) = 0 (USE ON FIRST CALL). */
289
290 /* 0 CAUSES FIXPNF TO INITIALIZE EVERYTHING FOR THE PROBLEM */
291 /* X = F(X) (USE ON FIRST CALL). */
292
293 /* 1 NORMAL RETURN. */
294
295 /* 2 SPECIFIED ERROR TOLERANCE CANNOT BE MET. SOME OR ALL OF */
296 /* ARCRE , ARCAE , ANSRE , ANSAE HAVE BEEN INCREASED TO */
297 /* SUITABLE VALUES. TO CONTINUE, JUST CALL FIXPNF AGAIN */
298 /* WITHOUT CHANGING ANY PARAMETERS. */
299
300 /* 3 STEPNF HAS BEEN CALLED 1000 TIMES. TO CONTINUE, CALL */
301 /* FIXPNF AGAIN WITHOUT CHANGING ANY PARAMETERS. */
302
303 /* 4 JACOBIAN MATRIX DOES NOT HAVE FULL RANK. THE ALGORITHM */
304 /* HAS FAILED (THE ZERO CURVE OF THE HOMOTOPY MAP CANNOT BE */
305 /* FOLLOWED ANY FURTHER). */
306
307 /* 5 THE TRACKING ALGORITHM HAS LOST THE ZERO CURVE OF THE */
308 /* HOMOTOPY MAP AND IS NOT MAKING PROGRESS. THE ERROR TOLERANCES */
309 /* ARC?E AND ANS?E WERE TOO LENIENT. THE PROBLEM SHOULD BE */
310 /* RESTARTED BY CALLING FIXPNF WITH SMALLER ERROR TOLERANCES */
311 /* AND IFLAG = 0 (-1, -2). */
312
313 /* 6 THE NORMAL FLOW NEWTON ITERATION IN STEPNF OR ROOTNF */
314 /* FAILED TO CONVERGE. THE ERROR TOLERANCES ANS?E MAY BE TOO */
315 /* STRINGENT. */
316
317 /* 7 ILLEGAL INPUT PARAMETERS, A FATAL ERROR. */
318
319 /* ARCRE , ARCAE , ANSRE , ANSAE ARE UNCHANGED AFTER A NORMAL RETURN */
320 /* (IFLAG = 1). THEY ARE INCREASED TO APPROPRIATE VALUES ON THE */
321 /* RETURN IFLAG = 2 . */
322
323 /* NFE IS THE NUMBER OF FUNCTION EVALUATIONS (= NUMBER OF */
324 /* JACOBIAN EVALUATIONS). */
325
326 /* ARCLEN IS THE LENGTH OF THE PATH FOLLOWED. */
327
328
329
330
331 /* ***** ARRAY DECLARATIONS. ***** */
332
333
334 /* ***** END OF DIMENSIONAL INFORMATION. ***** */
335
336
337 /* LIMITD IS AN UPPER BOUND ON THE NUMBER OF STEPS. IT MAY BE */
338 /* CHANGED BY CHANGING THE FOLLOWING PARAMETER STATEMENT: */
339
340 /* SWITCH FROM THE TOLERANCE ARC?E TO THE (FINER) TOLERANCE ANS?E IF
341 */
342 /* THE CURVATURE OF ANY COMPONENT OF Y EXCEEDS CURSW. */
343
344
345
346 /* : : : : : : : : : : : : : : : : : : : : : : : :
347 */
348 /* SET LOGICAL SWITCH TO REFLECT ENTRY POINT. */
349 /* Parameter adjustments */
350 --y;
351 --a;
352 --yp;
353 --yold;
354 --ypold;
355 qr_dim1 = *n;
356 qr_offset = qr_dim1 + 1;
357 qr -= qr_offset;
358 --alpha;
359 --tz;
360 --pivot;
361 --w;
362 --wp;
363 --z0;
364 --z1;
365 --sspar;
366 --par;
367 --ipar;
368
369 /* Function Body */
370 switch(n__) {
371 case 1: goto L_polynf;
372 }
373
374 polsys = FALSE_;
375 goto L11;
376
377 L_polynf:
378 polsys = TRUE_;
379 L11:
380
381 if (*n <= 0 || *ansre <= (float)0. || *ansae < (float)0.) {
382 *iflag = 7;
383 }
384 if (*iflag >= -2 && *iflag <= 0) {
385 goto L20;
386 }
387 if (*iflag == 2) {
388 goto L120;
389 }
390 if (*iflag == 3) {
391 goto L90;
392 }
393 /* ONLY VALID INPUT FOR IFLAG IS -2, -1, 0, 2, 3. */
394 *iflag = 7;
395 return 0;
396
397 /* ***** INITIALIZATION BLOCK. ***** */
398
399 L20:
400 *arclen = (float)0.;
401 if (*arcre <= (float)0.) {
402 *arcre = sqrt(*ansre) * (float).5;
403 }
404 if (*arcae <= (float)0.) {
405 *arcae = sqrt(*ansae) * (float).5;
406 }
407 nc = *n;
408 nfec = 0;
409 iflagc = *iflag;
410 np1 = *n + 1;
411 /* SET INITIAL CONDITIONS FOR FIRST CALL TO STEPNF . */
412 start = TRUE_;
413 crash = FALSE_;
414 hold = (float)1.;
415 h = (float).1;
416 s = (float)0.;
417 ypold[1] = (float)1.;
418 yp[1] = (float)1.;
419 y[1] = (float)0.;
420 i__1 = np1;
421 for (jw = 2; jw <= i__1; ++jw) {
422 ypold[jw] = (float)0.;
423 yp[jw] = (float)0.;
424 /* L40: */
425 }
426 /* SET OPTIMAL STEP SIZE ESTIMATION PARAMETERS. */
427 /* LET Z[K] DENOTE THE NEWTON ITERATES ALONG THE FLOW NORMAL TO THE */
428 /* DAVIDENKO FLOW AND Y THEIR LIMIT. */
429 /* IDEAL CONTRACTION FACTOR: ||Z[2] - Z[1]|| / ||Z[1] - Z[0]|| */
430 if (sspar[1] <= (float)0.) {
431 sspar[1] = (float).5;
432 }
433 /* IDEAL RESIDUAL FACTOR: ||RHO(A, Z[1])|| / ||RHO(A, Z[0])|| */
434 if (sspar[2] <= (float)0.) {
435 sspar[2] = (float).01;
436 }
437 /* IDEAL DISTANCE FACTOR: ||Z[1] - Y|| / ||Z[0] - Y|| */
438 if (sspar[3] <= (float)0.) {
439 sspar[3] = (float).5;
440 }
441 /* MINIMUM STEP SIZE HMIN . */
442 if (sspar[4] <= (float)0.) {
443 sspar[4] = (sqrt(*n + (float)1.) + (float)4.) * d1mach_(&c__4);
444 }
445 /* MAXIMUM STEP SIZE HMAX . */
446 if (sspar[5] <= (float)0.) {
447 sspar[5] = (float)1.;
448 }
449 /* MINIMUM STEP SIZE REDUCTION FACTOR BMIN . */
450 if (sspar[6] <= (float)0.) {
451 sspar[6] = (float).1;
452 }
453 /* MAXIMUM STEP SIZE EXPANSION FACTOR BMAX . */
454 if (sspar[7] <= (float)0.) {
455 sspar[7] = (float)3.;
456 }
457 /* ASSUMED OPERATING ORDER P . */
458 if (sspar[8] <= (float)0.) {
459 sspar[8] = (float)2.;
460 }
461
462 /* Adjustments of step size */
463 if (tweak == 1)
464 sspar[5] = (float).1;
465 else if (tweak == 2)
466 sspar[5] = (float).01;
467 else if (tweak == 3)
468 sspar[5] = (float).001;
469
470
471 /* LOAD A FOR THE FIXED POINT AND ZERO FINDING PROBLEMS. */
472 if (iflagc >= -1) {
473 i__1 = np1;
474 for (jw = 2; jw <= i__1; ++jw) {
475 a[jw - 1] = y[jw];
476 /* L60: */
477 }
478 }
479 L90:
480 limit = 1000;
481
482 /* Adjustment of maximum number of steps for smaller step size */
483 if (tweak == 1)
484 limit = 10000;
485 else if (tweak == 2)
486 limit = 100000;
487 else if (tweak == 3)
488 limit = 1000000;
489
490 /* ***** END OF INITIALIZATION BLOCK. ***** */
491
492
493 /* ***** MAIN LOOP. ***** */
494
495 L120:
496 i__1 = limit;
497 for (iter = 1; iter <= i__1; ++iter) {
498 if (y[1] < (float)0.) {
499 *arclen = s;
500 *iflag = 5;
501 return 0;
502 }
503
504 /* SET DIFFERENT ERROR TOLERANCE IF THE TRAJECTORY Y(S) HAS ANY HIGH
505 */
506 /* CURVATURE COMPONENTS. */
507 /* L140: */
508 curtol = hold * 10.;
509 relerr = *arcre;
510 abserr = *arcae;
511 i__2 = np1;
512 for (jw = 1; jw <= i__2; ++jw) {
513 if ((d__1 = yp[jw] - ypold[jw], abs((int)d__1)) > curtol) {
514 relerr = *ansre;
515 abserr = *ansae;
516 goto L200;
517 }
518 /* L160: */
519 }
520
521 /* TAKE A STEP ALONG THE CURVE. */
522 L200:
523 stepnf_((integer *)&nc, (integer *)&nfec, (integer *)&iflagc, &start,
524 &crash, &hold, &h, &relerr, &abserr,
525 &s, &y[1], &yp[1], &yold[1], &ypold[1], &a[1], &qr[
526 qr_offset], &alpha[1], &tz[1], (integer *)&pivot[1],
527 &w[1], &wp[1], &z0[
528 1], &z1[1], &sspar[1], &par[1], (integer *)&ipar[1]);
529 /* PRINT LATEST POINT ON CURVE IF REQUESTED. */
530 if (*trace > 0) {
531 print_homog(y+2,&coord_r, &coord_i);
532 #ifdef HOM_PRINT
533 fprintf(Hom_LogFile,"C %g %g",coord_r,coord_i)
534 #endif
535 ;
536 for(jw=2; jw<=np1; jw++){
537 #ifdef HOM_PRINT
538 fprintf(Hom_LogFile," %d = %g",jw,y[jw])
539 #endif
540 ;
541 }
542 #ifdef HOM_PRINT
543 fprintf(Hom_LogFile," %g",y[1])
544 #endif
545 ;
546 }
547 #ifdef HOM_PRINT
548 fprintf(Hom_LogFile," 4 0 %d %d %f\n",iter,nfec,s)
549 #endif
550 ;
551 *nfe = nfec;
552 /* CHECK IF THE STEP WAS SUCCESSFUL. */
553 if (iflagc > 0) {
554 *arclen = s;
555 *iflag = iflagc;
556 return 0;
557 }
558 if (crash) {
559 /* RETURN CODE FOR ERROR TOLERANCE TOO SMALL. */
560 *iflag = 2;
561 /* CHANGE ERROR TOLERANCES. */
562 if (*arcre < relerr) {
563 *arcre = relerr;
564 }
565 if (*ansre < relerr) {
566 *ansre = relerr;
567 }
568 if (*arcae < abserr) {
569 *arcae = abserr;
570 }
571 if (*ansae < abserr) {
572 *ansae = abserr;
573 }
574 /* CHANGE LIMIT ON NUMBER OF ITERATIONS. */
575 limit -= iter;
576 return 0;
577 }
578
579 if (y[1] >= (float)1.) {
580
581 /* USE HERMITE CUBIC INTERPOLATION AND NEWTON ITERATION TO GET THE
582 */
583 /* ANSWER AT LAMBDA = 1.0 . */
584
585 /* SAVE YOLD FOR ARC LENGTH CALCULATION LATER. */
586 i__2 = np1;
587 for (jw = 1; jw <= i__2; ++jw) {
588 z0[jw] = yold[jw];
589 /* L260: */
590 }
591 rootnf_(&nc, &nfec, &iflagc, ansre, ansae, &y[1], &yp[1], &yold[1],
592 &ypold[1], &a[1], &qr[qr_offset], &alpha[1], &tz[1],
593 &pivot[1], &w[1], &wp[1], &par[1], &ipar[1]);
594
595 *nfe = nfec;
596 *iflag = 1;
597 /* SET ERROR FLAG IF ROOTNF COULD NOT GET THE POINT ON THE ZERO
598 */
599 /* CURVE AT LAMBDA = 1.0 . */
600 if (iflagc > 0) {
601 *iflag = iflagc;
602 }
603 /* CALCULATE FINAL ARC LENGTH. */
604 i__2 = np1;
605 for (jw = 1; jw <= i__2; ++jw) {
606 w[jw] = y[jw] - z0[jw];
607 /* L290: */
608 }
609 *arclen = s - hold + dnrm2_((integer *)&np1, &w[1], &c__1);
610 return 0;
611 }
612
613 /* FOR POLYNOMIAL SYSTEMS AND THE POLSYS HOMOTOPY MAP, */
614 /* D LAMBDA/DS .GE. 0 NECESSARILY. THIS CONDITION IS FORCED HERE IF
615 */
616 /* THE ENTRY POINT WAS POLYNF . */
617
618 /* if (polsys) { */
619 if (yp[1] < (float)0.) {
620 /* REVERSE TANGENT DIRECTION SO D LAMBDA/DS = YP(1) > 0 . */
621 i__2 = np1;
622 for (jw = 1; jw <= i__2; ++jw) {
623 yp[jw] = -yp[jw];
624 ypold[jw] = yp[jw];
625 /* L310: */
626 /* } */
627 /* FORCE STEPNF TO USE THE LINEAR PREDICTOR FOR THE NEXT STE
628 P ONLY. */
629 start = TRUE_;
630 }
631 }
632
633 /* L400: */
634 }
635
636 /* ***** END OF MAIN LOOP. ***** */
637
638 /* LAMBDA HAS NOT REACHED 1 IN 1000 STEPS. */
639 *iflag = 3;
640 *arclen = s;
641 return 0;
642
643 } /* fixpnf_ */
644
fixpnf_(int * n,double * y,int * iflag,double * arcre,double * arcae,double * ansre,double * ansae,int * trace,double * a,int * nfe,double * arclen,double * yp,double * yold,double * ypold,double * qr,double * alpha,double * tz,int * pivot,double * w,double * wp,double * z0,double * z1,double * sspar,double * par,int * ipar,int tweak)645 /* Subroutine */ int fixpnf_(int *n,
646 double *y,
647 int *iflag,
648 double *arcre,
649 double *arcae,
650 double *ansre,
651 double *ansae,
652 int *trace,
653 double *a,
654 int *nfe,
655 double *arclen,
656 double *yp,
657 double *yold,
658 double *ypold,
659 double *qr,
660 double *alpha,
661 double *tz,
662 int *pivot,
663 double *w,
664 double *wp,
665 double *z0,
666 double *z1,
667 double *sspar,
668 double *par,
669 int *ipar,
670
671 int tweak)
672 {
673 return fixpnf_0_(0, n, y, iflag, arcre, arcae, ansre, ansae, trace, a,
674 nfe, arclen, yp, yold, ypold, qr, alpha, tz, pivot, w, wp, z0, z1,
675 sspar, par, ipar, tweak);
676 }
677
polynf_(int * n,double * y,int * iflag,double * arcre,double * arcae,double * ansre,double * ansae,int * trace,double * a,int * nfe,double * arclen,double * yp,double * yold,double * ypold,double * qr,double * alpha,double * tz,int * pivot,double * w,double * wp,double * z0,double * z1,double * sspar,double * par,int * ipar,int tweak)678 /* Subroutine */ int polynf_(int *n,
679 double *y,
680 int *iflag,
681 double *arcre,
682 double *arcae,
683 double *ansre,
684 double *ansae,
685 int *trace,
686 double *a,
687 int *nfe,
688 double *arclen,
689 double *yp,
690 double *yold,
691 double *ypold,
692 double *qr,
693 double *alpha,
694 double *tz,
695 int *pivot,
696 double *w,
697 double *wp,
698 double *z0,
699 double *z1,
700 double *sspar,
701 double *par,
702 int *ipar,
703
704 int tweak)
705 {
706 return fixpnf_0_(1, n, y, iflag, arcre, arcae, ansre, ansae, trace, a,
707 nfe, arclen, yp, yold, ypold, qr, alpha, tz, pivot, w, wp, z0, z1,
708 sspar, par, ipar, tweak);
709 }
710
711
712
713 /* end fixpnf.c */
714
715 /**********************************************************************/
716 /***************** implementation code from Hom_Mem.c *****************/
717 /**********************************************************************/
718
719 /*--------------------------------------------------------------------
720 Workspace and Workspace management macros
721 -------------------------------------------------------------------*/
722 #define TOPN 20 /*Assume no more than 20 variables*/
723 #define TOPM 500 /* Assume no more than 500 monomials total*/
724 static int Istore[TOPM*(TOPN+1)+3*TOPN+2+TOPM];
725 static double Dstore[TOPM*2+12*(TOPN+1)+8+TOPN*(TOPN+2)+1];
726 static int didx=0,iidx=0;
727
728 /* initialization for storage currently static should be
729 made dynamic*/
730
731 /* access funtions to double storage */
Dres(int sz)732 double *Dres(int sz){ int v=didx; didx+=sz; return Dstore+v;}
Dtop()733 int Dtop(){return didx;}
Dfree(int ntop)734 int Dfree(int ntop){ didx=ntop; return 0;}
735
736 /* access funtions to int storage */
Ires(int sz)737 int *Ires(int sz){ int v=iidx; iidx+=sz; return Istore+v;}
Itop()738 int Itop(){return iidx;}
Ifree(int ntop)739 int Ifree(int ntop){ iidx=ntop; return 0;}
740
741 /* end Hom_Mem.c */
742
743 /**********************************************************************/
744 /***************** implementation code from d1mach.c ******************/
745 /**********************************************************************/
746
747 /* Table of constant values */
748
d1mach_(integer * i)749 doublereal d1mach_(integer *i)
750 {
751 /* Initialized data */
752
753 static integer sc = 987;
754 static struct {
755 integer e_1[10];
756 doublereal e_2;
757 } equiv_4 = { { 0, 1048576, -1, 2146435071, 0, 1017118720, 0,
758 1018167296, 1352628735, 1070810131 }, 0.0 };
759
760
761
762 /* System generated locals */
763 doublereal ret_val = 0.0;
764
765 /* Builtin functions */
766 /* Subroutine */ // int s_stop();
767 // integer s_wsfe(), do_fio(), e_wsfe();
768
769 /* Local variables */
770 #define log10 ((integer *)&equiv_4 + 8)
771 #define dmach ((doublereal *)&equiv_4)
772 #define large ((integer *)&equiv_4 + 2)
773 #define small ((integer *)&equiv_4)
774 #define diver ((integer *)&equiv_4 + 6)
775 #define right ((integer *)&equiv_4 + 4)
776
777
778
779
780 /* DOUBLE-PRECISION MACHINE CONSTANTS */
781
782 /* D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE. */
783
784 /* D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE. */
785
786 /* D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING. */
787
788 /* D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING. */
789
790 /* D1MACH( 5) = LOG10(B) */
791
792 /* TO ALTER THIS FUNCTION FOR A PARTICULAR ENVIRONMENT, */
793 /* THE DESIRED SET OF DATA STATEMENTS SHOULD BE ACTIVATED BY */
794 /* REMOVING THE C FROM COLUMN 1. */
795 /* ON RARE MACHINES A STATIC STATEMENT MAY NEED TO BE ADDED. */
796 /* (BUT PROBABLY MORE SYSTEMS PROHIBIT IT THAN REQUIRE IT.) */
797
798 /* FOR IEEE-ARITHMETIC MACHINES (BINARY STANDARD), ONE OF THE FIRST */
799 /* TWO SETS OF CONSTANTS BELOW SHOULD BE APPROPRIATE. IF YOU DO NOT */
800 /* KNOW WHICH SET TO USE, TRY BOTH AND SEE WHICH GIVES PLAUSIBLE */
801 /* VALUES. */
802
803 /* WHERE POSSIBLE, DECIMAL, OCTAL OR HEXADECIMAL CONSTANTS ARE USED */
804 /* TO SPECIFY THE CONSTANTS EXACTLY. SOMETIMES THIS REQUIRES USING */
805 /* EQUIVALENT INTEGER ARRAYS. IF YOUR COMPILER USES HALF-WORD */
806 /* INTEGERS BY DEFAULT (SOMETIMES CALLED INTEGER*2), YOU MAY NEED TO */
807 /* CHANGE INTEGER TO INTEGER*4 OR OTHERWISE INSTRUCT YOUR COMPILER */
808 /* TO USE FULL-WORD INTEGERS IN THE NEXT 5 DECLARATIONS. */
809
810 /* COMMENTS JUST BEFORE THE END STATEMENT (LINES STARTING WITH *) */
811 /* GIVE C SOURCE FOR D1MACH. */
812
813
814
815
816 /* MACHINE CONSTANTS FOR BIG-ENDIAN IEEE ARITHMETIC (BINARY FORMAT) */
817 /* MACHINES IN WHICH THE MOST SIGNIFICANT BYTE IS STORED FIRST, */
818 /* SUCH AS THE AT&T 3B SERIES, MOTOROLA 68000 BASED MACHINES (E.G. */
819 /* SUN 3), AND MACHINES THAT USE SPARC, HP, OR IBM RISC CHIPS. */
820
821 /* DATA SMALL(1),SMALL(2) / 1048576, 0 / */
822 /* DATA LARGE(1),LARGE(2) / 2146435071, -1 / */
823 /* DATA RIGHT(1),RIGHT(2) / 1017118720, 0 / */
824 /* DATA DIVER(1),DIVER(2) / 1018167296, 0 / */
825 /* DATA LOG10(1),LOG10(2) / 1070810131, 1352628735 /, SC/987/ */
826
827 /* MACHINE CONSTANTS FOR LITTLE-ENDIAN (BINARY) IEEE ARITHMETIC */
828 /* MACHINES IN WHICH THE LEAST SIGNIFICANT BYTE IS STORED FIRST, */
829 /* E.G. IBM PCS AND OTHER MACHINES THAT USE INTEL 80X87 OR DEC */
830 /* ALPHA CHIPS. */
831
832
833 /* MACHINE CONSTANTS FOR AMDAHL MACHINES. */
834
835 /* DATA SMALL(1),SMALL(2) / 1048576, 0 / */
836 /* DATA LARGE(1),LARGE(2) / 2147483647, -1 / */
837 /* DATA RIGHT(1),RIGHT(2) / 856686592, 0 / */
838 /* DATA DIVER(1),DIVER(2) / 873463808, 0 / */
839 /* DATA LOG10(1),LOG10(2) / 1091781651, 1352628735 /, SC/987/ */
840
841 /* MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. */
842
843 /* DATA SMALL(1) / ZC00800000 / */
844 /* DATA SMALL(2) / Z000000000 / */
845
846 /* DATA LARGE(1) / ZDFFFFFFFF / */
847 /* DATA LARGE(2) / ZFFFFFFFFF / */
848
849 /* DATA RIGHT(1) / ZCC5800000 / */
850 /* DATA RIGHT(2) / Z000000000 / */
851
852 /* DATA DIVER(1) / ZCC6800000 / */
853 /* DATA DIVER(2) / Z000000000 / */
854
855 /* DATA LOG10(1) / ZD00E730E7 / */
856 /* DATA LOG10(2) / ZC77800DC0 /, SC/987/ */
857
858 /* MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM. */
859
860 /* DATA SMALL(1) / O1771000000000000 / */
861 /* DATA SMALL(2) / O0000000000000000 / */
862
863 /* DATA LARGE(1) / O0777777777777777 / */
864 /* DATA LARGE(2) / O0007777777777777 / */
865
866 /* DATA RIGHT(1) / O1461000000000000 / */
867 /* DATA RIGHT(2) / O0000000000000000 / */
868
869 /* DATA DIVER(1) / O1451000000000000 / */
870 /* DATA DIVER(2) / O0000000000000000 / */
871
872 /* DATA LOG10(1) / O1157163034761674 / */
873 /* DATA LOG10(2) / O0006677466732724 /, SC/987/ */
874
875 /* MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS. */
876
877 /* DATA SMALL(1) / O1771000000000000 / */
878 /* DATA SMALL(2) / O7770000000000000 / */
879
880 /* DATA LARGE(1) / O0777777777777777 / */
881 /* DATA LARGE(2) / O7777777777777777 / */
882
883 /* DATA RIGHT(1) / O1461000000000000 / */
884 /* DATA RIGHT(2) / O0000000000000000 / */
885
886 /* DATA DIVER(1) / O1451000000000000 / */
887 /* DATA DIVER(2) / O0000000000000000 / */
888
889 /* DATA LOG10(1) / O1157163034761674 / */
890 /* DATA LOG10(2) / O0006677466732724 /, SC/987/ */
891
892 /* MACHINE CONSTANTS FOR FTN4 ON THE CDC 6000/7000 SERIES. */
893
894 /* DATA SMALL(1) / 00564000000000000000B / */
895 /* DATA SMALL(2) / 00000000000000000000B / */
896
897 /* DATA LARGE(1) / 37757777777777777777B / */
898 /* DATA LARGE(2) / 37157777777777777774B / */
899
900 /* DATA RIGHT(1) / 15624000000000000000B / */
901 /* DATA RIGHT(2) / 00000000000000000000B / */
902
903 /* DATA DIVER(1) / 15634000000000000000B / */
904 /* DATA DIVER(2) / 00000000000000000000B / */
905
906 /* DATA LOG10(1) / 17164642023241175717B / */
907 /* DATA LOG10(2) / 16367571421742254654B /, SC/987/ */
908
909 /* MACHINE CONSTANTS FOR FTN5 ON THE CDC 6000/7000 SERIES. */
910
911 /* DATA SMALL(1) / O"00564000000000000000" / */
912 /* DATA SMALL(2) / O"00000000000000000000" / */
913
914 /* DATA LARGE(1) / O"37757777777777777777" / */
915 /* DATA LARGE(2) / O"37157777777777777774" / */
916
917 /* DATA RIGHT(1) / O"15624000000000000000" / */
918 /* DATA RIGHT(2) / O"00000000000000000000" / */
919
920 /* DATA DIVER(1) / O"15634000000000000000" / */
921 /* DATA DIVER(2) / O"00000000000000000000" / */
922
923 /* DATA LOG10(1) / O"17164642023241175717" / */
924 /* DATA LOG10(2) / O"16367571421742254654" /, SC/987/ */
925
926 /* MACHINE CONSTANTS FOR CONVEX C-1 */
927
928 /* DATA SMALL(1),SMALL(2) / '00100000'X, '00000000'X / */
929 /* DATA LARGE(1),LARGE(2) / '7FFFFFFF'X, 'FFFFFFFF'X / */
930 /* DATA RIGHT(1),RIGHT(2) / '3CC00000'X, '00000000'X / */
931 /* DATA DIVER(1),DIVER(2) / '3CD00000'X, '00000000'X / */
932 /* DATA LOG10(1),LOG10(2) / '3FF34413'X, '509F79FF'X /, SC/987/ */
933
934 /* MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3. */
935
936 /* DATA SMALL(1) / 201354000000000000000B / */
937 /* DATA SMALL(2) / 000000000000000000000B / */
938
939 /* DATA LARGE(1) / 577767777777777777777B / */
940 /* DATA LARGE(2) / 000007777777777777776B / */
941
942 /* DATA RIGHT(1) / 376434000000000000000B / */
943 /* DATA RIGHT(2) / 000000000000000000000B / */
944
945 /* DATA DIVER(1) / 376444000000000000000B / */
946 /* DATA DIVER(2) / 000000000000000000000B / */
947
948 /* DATA LOG10(1) / 377774642023241175717B / */
949 /* DATA LOG10(2) / 000007571421742254654B /, SC/987/ */
950
951 /* MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200 */
952
953 /* SMALL, LARGE, RIGHT, DIVER, LOG10 SHOULD BE DECLARED */
954 /* INTEGER SMALL(4), LARGE(4), RIGHT(4), DIVER(4), LOG10(4) */
955
956 /* NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING LINE - */
957 /* STATIC DMACH(5) */
958
959 /* DATA SMALL/20K,3*0/,LARGE/77777K,3*177777K/ */
960 /* DATA RIGHT/31420K,3*0/,DIVER/32020K,3*0/ */
961 /* DATA LOG10/40423K,42023K,50237K,74776K/, SC/987/ */
962
963 /* MACHINE CONSTANTS FOR THE HARRIS SLASH 6 AND SLASH 7 */
964
965 /* DATA SMALL(1),SMALL(2) / '20000000, '00000201 / */
966 /* DATA LARGE(1),LARGE(2) / '37777777, '37777577 / */
967 /* DATA RIGHT(1),RIGHT(2) / '20000000, '00000333 / */
968 /* DATA DIVER(1),DIVER(2) / '20000000, '00000334 / */
969 /* DATA LOG10(1),LOG10(2) / '23210115, '10237777 /, SC/987/ */
970
971 /* MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES. */
972
973 /* DATA SMALL(1),SMALL(2) / O402400000000, O000000000000 / */
974 /* DATA LARGE(1),LARGE(2) / O376777777777, O777777777777 / */
975 /* DATA RIGHT(1),RIGHT(2) / O604400000000, O000000000000 / */
976 /* DATA DIVER(1),DIVER(2) / O606400000000, O000000000000 / */
977 /* DATA LOG10(1),LOG10(2) / O776464202324, O117571775714 /, SC/987/
978 */
979
980 /* MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, */
981 /* THE XEROX SIGMA 5/7/9 AND THE SEL SYSTEMS 85/86. */
982
983 /* DATA SMALL(1),SMALL(2) / Z00100000, Z00000000 / */
984 /* DATA LARGE(1),LARGE(2) / Z7FFFFFFF, ZFFFFFFFF / */
985 /* DATA RIGHT(1),RIGHT(2) / Z33100000, Z00000000 / */
986 /* DATA DIVER(1),DIVER(2) / Z34100000, Z00000000 / */
987 /* DATA LOG10(1),LOG10(2) / Z41134413, Z509F79FF /, SC/987/ */
988
989 /* MACHINE CONSTANTS FOR THE INTERDATA 8/32 */
990 /* WITH THE UNIX SYSTEM FORTRAN 77 COMPILER. */
991
992 /* FOR THE INTERDATA FORTRAN VII COMPILER REPLACE */
993 /* THE Z'S SPECIFYING HEX CONSTANTS WITH Y'S. */
994
995 /* DATA SMALL(1),SMALL(2) / Z'00100000', Z'00000000' / */
996 /* DATA LARGE(1),LARGE(2) / Z'7EFFFFFF', Z'FFFFFFFF' / */
997 /* DATA RIGHT(1),RIGHT(2) / Z'33100000', Z'00000000' / */
998 /* DATA DIVER(1),DIVER(2) / Z'34100000', Z'00000000' / */
999 /* DATA LOG10(1),LOG10(2) / Z'41134413', Z'509F79FF' /, SC/987/ */
1000
1001 /* MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR). */
1002
1003 /* DATA SMALL(1),SMALL(2) / "033400000000, "000000000000 / */
1004 /* DATA LARGE(1),LARGE(2) / "377777777777, "344777777777 / */
1005 /* DATA RIGHT(1),RIGHT(2) / "113400000000, "000000000000 / */
1006 /* DATA DIVER(1),DIVER(2) / "114400000000, "000000000000 / */
1007 /* DATA LOG10(1),LOG10(2) / "177464202324, "144117571776 /, SC/987/
1008 */
1009
1010 /* MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR). */
1011
1012 /* DATA SMALL(1),SMALL(2) / "000400000000, "000000000000 / */
1013 /* DATA LARGE(1),LARGE(2) / "377777777777, "377777777777 / */
1014 /* DATA RIGHT(1),RIGHT(2) / "103400000000, "000000000000 / */
1015 /* DATA DIVER(1),DIVER(2) / "104400000000, "000000000000 / */
1016 /* DATA LOG10(1),LOG10(2) / "177464202324, "047674776746 /, SC/987/
1017 */
1018
1019 /* MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING */
1020 /* 32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). */
1021
1022 /* DATA SMALL(1),SMALL(2) / 8388608, 0 / */
1023 /* DATA LARGE(1),LARGE(2) / 2147483647, -1 / */
1024 /* DATA RIGHT(1),RIGHT(2) / 612368384, 0 / */
1025 /* DATA DIVER(1),DIVER(2) / 620756992, 0 / */
1026 /* DATA LOG10(1),LOG10(2) / 1067065498, -2063872008 /, SC/987/ */
1027
1028 /* DATA SMALL(1),SMALL(2) / O00040000000, O00000000000 / */
1029 /* DATA LARGE(1),LARGE(2) / O17777777777, O37777777777 / */
1030 /* DATA RIGHT(1),RIGHT(2) / O04440000000, O00000000000 / */
1031 /* DATA DIVER(1),DIVER(2) / O04500000000, O00000000000 / */
1032 /* DATA LOG10(1),LOG10(2) / O07746420232, O20476747770 /, SC/987/ */
1033
1034 /* MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING */
1035 /* 16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). */
1036
1037 /* SMALL, LARGE, RIGHT, DIVER, LOG10 SHOULD BE DECLARED */
1038 /* INTEGER SMALL(4), LARGE(4), RIGHT(4), DIVER(4), LOG10(4) */
1039
1040 /* DATA SMALL(1),SMALL(2) / 128, 0 / */
1041 /* DATA SMALL(3),SMALL(4) / 0, 0 / */
1042
1043 /* DATA LARGE(1),LARGE(2) / 32767, -1 / */
1044 /* DATA LARGE(3),LARGE(4) / -1, -1 / */
1045
1046 /* DATA RIGHT(1),RIGHT(2) / 9344, 0 / */
1047 /* DATA RIGHT(3),RIGHT(4) / 0, 0 / */
1048
1049 /* DATA DIVER(1),DIVER(2) / 9472, 0 / */
1050 /* DATA DIVER(3),DIVER(4) / 0, 0 / */
1051
1052 /* DATA LOG10(1),LOG10(2) / 16282, 8346 / */
1053 /* DATA LOG10(3),LOG10(4) / -31493, -12296 /, SC/987/ */
1054
1055 /* DATA SMALL(1),SMALL(2) / O000200, O000000 / */
1056 /* DATA SMALL(3),SMALL(4) / O000000, O000000 / */
1057
1058 /* DATA LARGE(1),LARGE(2) / O077777, O177777 / */
1059 /* DATA LARGE(3),LARGE(4) / O177777, O177777 / */
1060
1061 /* DATA RIGHT(1),RIGHT(2) / O022200, O000000 / */
1062 /* DATA RIGHT(3),RIGHT(4) / O000000, O000000 / */
1063
1064 /* DATA DIVER(1),DIVER(2) / O022400, O000000 / */
1065 /* DATA DIVER(3),DIVER(4) / O000000, O000000 / */
1066
1067 /* DATA LOG10(1),LOG10(2) / O037632, O020232 / */
1068 /* DATA LOG10(3),LOG10(4) / O102373, O147770 /, SC/987/ */
1069
1070 /* MACHINE CONSTANTS FOR THE PRIME 50 SERIES SYSTEMS */
1071 /* WITH 32-BIT INTEGERS AND 64V MODE INSTRUCTIONS, */
1072 /* SUPPLIED BY IGOR BRAY. */
1073
1074 /* DATA SMALL(1),SMALL(2) / :10000000000, :00000100001 / */
1075 /* DATA LARGE(1),LARGE(2) / :17777777777, :37777677775 / */
1076 /* DATA RIGHT(1),RIGHT(2) / :10000000000, :00000000122 / */
1077 /* DATA DIVER(1),DIVER(2) / :10000000000, :00000000123 / */
1078 /* DATA LOG10(1),LOG10(2) / :11504046501, :07674600177 /, SC/987/ */
1079
1080 /* MACHINE CONSTANTS FOR THE SEQUENT BALANCE 8000 */
1081
1082 /* DATA SMALL(1),SMALL(2) / $00000000, $00100000 / */
1083 /* DATA LARGE(1),LARGE(2) / $FFFFFFFF, $7FEFFFFF / */
1084 /* DATA RIGHT(1),RIGHT(2) / $00000000, $3CA00000 / */
1085 /* DATA DIVER(1),DIVER(2) / $00000000, $3CB00000 / */
1086 /* DATA LOG10(1),LOG10(2) / $509F79FF, $3FD34413 /, SC/987/ */
1087
1088 /* MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. */
1089
1090 /* DATA SMALL(1),SMALL(2) / O000040000000, O000000000000 / */
1091 /* DATA LARGE(1),LARGE(2) / O377777777777, O777777777777 / */
1092 /* DATA RIGHT(1),RIGHT(2) / O170540000000, O000000000000 / */
1093 /* DATA DIVER(1),DIVER(2) / O170640000000, O000000000000 / */
1094 /* DATA LOG10(1),LOG10(2) / O177746420232, O411757177572 /, SC/987/
1095 */
1096
1097 /* MACHINE CONSTANTS FOR THE VAX UNIX F77 COMPILER */
1098
1099 /* DATA SMALL(1),SMALL(2) / 128, 0 / */
1100 /* DATA LARGE(1),LARGE(2) / -32769, -1 / */
1101 /* DATA RIGHT(1),RIGHT(2) / 9344, 0 / */
1102 /* DATA DIVER(1),DIVER(2) / 9472, 0 / */
1103 /* DATA LOG10(1),LOG10(2) / 546979738, -805796613 /, SC/987/ */
1104
1105 /* MACHINE CONSTANTS FOR THE VAX-11 WITH */
1106 /* FORTRAN IV-PLUS COMPILER */
1107
1108 /* DATA SMALL(1),SMALL(2) / Z00000080, Z00000000 / */
1109 /* DATA LARGE(1),LARGE(2) / ZFFFF7FFF, ZFFFFFFFF / */
1110 /* DATA RIGHT(1),RIGHT(2) / Z00002480, Z00000000 / */
1111 /* DATA DIVER(1),DIVER(2) / Z00002500, Z00000000 / */
1112 /* DATA LOG10(1),LOG10(2) / Z209A3F9A, ZCFF884FB /, SC/987/ */
1113
1114 /* MACHINE CONSTANTS FOR VAX/VMS VERSION 2.2 */
1115
1116 /* DATA SMALL(1),SMALL(2) / '80'X, '0'X / */
1117 /* DATA LARGE(1),LARGE(2) / 'FFFF7FFF'X, 'FFFFFFFF'X / */
1118 /* DATA RIGHT(1),RIGHT(2) / '2480'X, '0'X / */
1119 /* DATA DIVER(1),DIVER(2) / '2500'X, '0'X / */
1120 /* DATA LOG10(1),LOG10(2) / '209A3F9A'X, 'CFF884FB'X /, SC/987/ */
1121
1122 /* *** ISSUE STOP 779 IF ALL DATA STATEMENTS ARE COMMENTED... */
1123 if (sc != 987) {
1124 /* s_stop("779", 3L); */exit(0);
1125 }
1126 /* *** ISSUE STOP 778 IF ALL DATA STATEMENTS ARE OBVIOUSLY WRONG... */
1127 if (dmach[3] >= 1.) {
1128 /* s_stop("778", 3L); */ exit(0);
1129 }
1130 if (*i < 1 || *i > 5) {
1131 goto L999;
1132 }
1133 ret_val = dmach[*i - 1];
1134 return ret_val;
1135 L999:
1136 /* s_wsfe(&io___8);
1137 do_fio(&c__1, (char *)&(*i), (ftnlen)sizeof(integer));
1138 e_wsfe();
1139 s_stop("", 0L); */
1140 return ret_val;
1141 } /* d1mach_ */
1142
1143 #undef right
1144 #undef diver
1145 #undef small
1146 #undef large
1147 #undef dmach
1148 #undef log10
1149
1150 /* end d1mach.c */
1151
1152 /**********************************************************************/
1153 /***************** implementation code from daxpy.c *******************/
1154 /**********************************************************************/
1155
daxpy_(integer * n,doublereal * da,doublereal * dx,integer * incx,doublereal * dy,integer * incy)1156 /* Subroutine */ int daxpy_(integer *n,
1157 doublereal *da,
1158 doublereal *dx,
1159 integer *incx,
1160 doublereal *dy,
1161 integer *incy)
1162 {
1163 /* System generated locals */
1164 integer i__1;
1165
1166 /* Local variables */
1167 static integer i, m, ix, iy, mp1;
1168
1169
1170 /* constant times a vector plus a vector. */
1171 /* uses unrolled loops for increments equal to one. */
1172 /* jack dongarra, linpack, 3/11/78. */
1173
1174
1175 /* Parameter adjustments */
1176 --dy;
1177 --dx;
1178
1179 /* Function Body */
1180 if (*n <= 0) {
1181 return 0;
1182 }
1183 if (*da == 0.) {
1184 return 0;
1185 }
1186 if (*incx == 1 && *incy == 1) {
1187 goto L20;
1188 }
1189
1190 /* code for unequal increments or equal increments */
1191 /* not equal to 1 */
1192
1193 ix = 1;
1194 iy = 1;
1195 if (*incx < 0) {
1196 ix = (-(*n) + 1) * *incx + 1;
1197 }
1198 if (*incy < 0) {
1199 iy = (-(*n) + 1) * *incy + 1;
1200 }
1201 i__1 = *n;
1202 for (i = 1; i <= i__1; ++i) {
1203 dy[iy] += *da * dx[ix];
1204 ix += *incx;
1205 iy += *incy;
1206 /* L10: */
1207 }
1208 return 0;
1209
1210 /* code for both increments equal to 1 */
1211
1212
1213 /* clean-up loop */
1214
1215 L20:
1216 m = *n % 4;
1217 if (m == 0) {
1218 goto L40;
1219 }
1220 i__1 = m;
1221 for (i = 1; i <= i__1; ++i) {
1222 dy[i] += *da * dx[i];
1223 /* L30: */
1224 }
1225 if (*n < 4) {
1226 return 0;
1227 }
1228 L40:
1229 mp1 = m + 1;
1230 i__1 = *n;
1231 for (i = mp1; i <= i__1; i += 4) {
1232 dy[i] += *da * dx[i];
1233 dy[i + 1] += *da * dx[i + 1];
1234 dy[i + 2] += *da * dx[i + 2];
1235 dy[i + 3] += *da * dx[i + 3];
1236 /* L50: */
1237 }
1238 return 0;
1239 } /* daxpy_ */
1240
1241 /* end daxpy.c */
1242
1243 /**********************************************************************/
1244 /***************** implementation code from dcopy.c *******************/
1245 /**********************************************************************/
1246
dcopy_(integer * n,doublereal * dx,integer * incx,doublereal * dy,integer * incy)1247 /* Subroutine */ int dcopy_(integer *n,
1248 doublereal *dx,
1249 integer *incx,
1250 doublereal *dy,
1251 integer *incy)
1252 {
1253 /* System generated locals */
1254 integer i__1;
1255
1256 /* Local variables */
1257 static integer i, m, ix, iy, mp1;
1258
1259
1260 /* copies a vector, x, to a vector, y. */
1261 /* uses unrolled loops for increments equal to one. */
1262 /* jack dongarra, linpack, 3/11/78. */
1263
1264
1265 /* Parameter adjustments */
1266 --dy;
1267 --dx;
1268
1269 /* Function Body */
1270 if (*n <= 0) {
1271 return 0;
1272 }
1273 if (*incx == 1 && *incy == 1) {
1274 goto L20;
1275 }
1276
1277 /* code for unequal increments or equal increments */
1278 /* not equal to 1 */
1279
1280 ix = 1;
1281 iy = 1;
1282 if (*incx < 0) {
1283 ix = (-(*n) + 1) * *incx + 1;
1284 }
1285 if (*incy < 0) {
1286 iy = (-(*n) + 1) * *incy + 1;
1287 }
1288 i__1 = *n;
1289 for (i = 1; i <= i__1; ++i) {
1290 dy[iy] = dx[ix];
1291 ix += *incx;
1292 iy += *incy;
1293 /* L10: */
1294 }
1295 return 0;
1296
1297 /* code for both increments equal to 1 */
1298
1299
1300 /* clean-up loop */
1301
1302 L20:
1303 m = *n % 7;
1304 if (m == 0) {
1305 goto L40;
1306 }
1307 i__1 = m;
1308 for (i = 1; i <= i__1; ++i) {
1309 dy[i] = dx[i];
1310 /* L30: */
1311 }
1312 if (*n < 7) {
1313 return 0;
1314 }
1315 L40:
1316 mp1 = m + 1;
1317 i__1 = *n;
1318 for (i = mp1; i <= i__1; i += 7) {
1319 dy[i] = dx[i];
1320 dy[i + 1] = dx[i + 1];
1321 dy[i + 2] = dx[i + 2];
1322 dy[i + 3] = dx[i + 3];
1323 dy[i + 4] = dx[i + 4];
1324 dy[i + 5] = dx[i + 5];
1325 dy[i + 6] = dx[i + 6];
1326 /* L50: */
1327 }
1328 return 0;
1329 } /* dcopy_ */
1330
1331 /* end dcopy.c */
1332
1333 /**********************************************************************/
1334 /***************** implementation code from dcpose.c ******************/
1335 /**********************************************************************/
1336
1337 /* Table of constant values */
1338
1339 /* static integer c__1 = 1; NOW REDUNDANT */
1340
dcpose_(integer * ndim,integer * n,doublereal * qr,doublereal * alpha,integer * pivot,integer * ierr,doublereal * y,doublereal * sum)1341 /* Subroutine */ int dcpose_(integer *ndim,
1342 integer *n,
1343 doublereal *qr,
1344 doublereal *alpha,
1345 integer *pivot,
1346 integer *ierr,
1347 doublereal *y,
1348 doublereal *sum)
1349 {
1350 /* System generated locals */
1351 integer qr_dim1, qr_offset, i__1, i__2, i__3;
1352 doublereal d__1;
1353
1354 /* Builtin functions */
1355 /* double sqrt(double); CANT DECLARE BUILTINS UNDER C++ */
1356
1357 /* Local variables */
1358 static doublereal beta;
1359 static integer jbar;
1360 extern doublereal ddot_(integer *n,
1361 doublereal *dx,
1362 integer *incx,
1363 doublereal *dy,
1364 integer *incy);
1365 static doublereal qrkk;
1366 static integer i, j, k;
1367 static doublereal sigma, alphak;
1368 static integer kp1, np1;
1369
1370
1371 /* SUBROUTINE DCPOSE IS A MODIFICATION OF THE ALGOL PROCEDURE */
1372 /* DECOMPOSE IN P. BUSINGER AND G. H. GOLUB, LINEAR LEAST */
1373 /* SQUARES SOLUTIONS BY HOUSEHOLDER TRANSFORMATIONS, */
1374 /* NUMER. MATH. 7 (1965) 269-276. */
1375
1376 /* Parameter adjustments */
1377 qr_dim1 = *ndim;
1378 qr_offset = qr_dim1 + 1;
1379 qr -= qr_offset;
1380 --alpha;
1381 --pivot;
1382 --y;
1383 --sum;
1384
1385 /* Function Body */
1386 *ierr = 0;
1387 np1 = *n + 1;
1388 i__1 = np1;
1389 for (j = 1; j <= i__1; ++j) {
1390 sum[j] = ddot_(n, &qr[j * qr_dim1 + 1], &c__1, &qr[j * qr_dim1 + 1], &
1391 c__1);
1392 /* L20: */
1393 pivot[j] = j;
1394 }
1395 i__1 = *n;
1396 for (k = 1; k <= i__1; ++k) {
1397 sigma = sum[k];
1398 jbar = k;
1399 kp1 = k + 1;
1400 i__2 = np1;
1401 for (j = kp1; j <= i__2; ++j) {
1402 if (sigma >= sum[j]) {
1403 goto L40;
1404 }
1405 sigma = sum[j];
1406 jbar = j;
1407 L40:
1408 ;
1409 }
1410 if (jbar == k) {
1411 goto L70;
1412 }
1413 i = pivot[k];
1414 pivot[k] = pivot[jbar];
1415 pivot[jbar] = i;
1416 sum[jbar] = sum[k];
1417 sum[k] = sigma;
1418 i__2 = *n;
1419 for (i = 1; i <= i__2; ++i) {
1420 sigma = qr[i + k * qr_dim1];
1421 qr[i + k * qr_dim1] = qr[i + jbar * qr_dim1];
1422 qr[i + jbar * qr_dim1] = sigma;
1423 /* L50: */
1424 }
1425 /* END OF COLUMN INTERCHANGE. */
1426 L70:
1427 i__2 = *n - k + 1;
1428 sigma = ddot_(&i__2, &qr[k + k * qr_dim1], &c__1, &qr[k + k * qr_dim1]
1429 , &c__1);
1430 if (sigma != (float)0.) {
1431 goto L60;
1432 }
1433 *ierr = 1;
1434 return 0;
1435 L60:
1436 if (k == *n) {
1437 goto L500;
1438 }
1439 qrkk = qr[k + k * qr_dim1];
1440 alphak = -sqrt(sigma);
1441 if (qrkk < (float)0.) {
1442 alphak = -alphak;
1443 }
1444 alpha[k] = alphak;
1445 beta = (float)1. / (sigma - qrkk * alphak);
1446 qr[k + k * qr_dim1] = qrkk - alphak;
1447 i__2 = np1;
1448 for (j = kp1; j <= i__2; ++j) {
1449 /* L80: */
1450 i__3 = *n - k + 1;
1451 y[j] = beta * ddot_(&i__3, &qr[k + k * qr_dim1], &c__1, &qr[k + j
1452 * qr_dim1], &c__1);
1453 }
1454 i__3 = np1;
1455 for (j = kp1; j <= i__3; ++j) {
1456 i__2 = *n;
1457 for (i = k; i <= i__2; ++i) {
1458 qr[i + j * qr_dim1] -= qr[i + k * qr_dim1] * y[j];
1459 /* L90: */
1460 }
1461 /* Computing 2nd power */
1462 d__1 = qr[k + j * qr_dim1];
1463 sum[j] -= d__1 * d__1;
1464 /* L100: */
1465 }
1466 L500:
1467 ;
1468 }
1469 alpha[*n] = qr[*n + *n * qr_dim1];
1470 return 0;
1471 } /* dcpose_ */
1472
1473 /* end dcpose.c */
1474
1475 /**********************************************************************/
1476 /****************** implementation code from ddot.c *******************/
1477 /**********************************************************************/
1478
ddot_(integer * n,doublereal * dx,integer * incx,doublereal * dy,integer * incy)1479 doublereal ddot_(integer *n,
1480 doublereal *dx,
1481 integer *incx,
1482 doublereal *dy,
1483 integer *incy)
1484 {
1485 /* System generated locals */
1486 integer i__1;
1487 doublereal ret_val;
1488
1489 /* Local variables */
1490 static integer i, m;
1491 static doublereal dtemp;
1492 static integer ix, iy, mp1;
1493
1494
1495 /* forms the dot product of two vectors. */
1496 /* uses unrolled loops for increments equal to one. */
1497 /* jack dongarra, linpack, 3/11/78. */
1498
1499
1500 /* Parameter adjustments */
1501 --dy;
1502 --dx;
1503
1504 /* Function Body */
1505 ret_val = 0.;
1506 dtemp = 0.;
1507 if (*n <= 0) {
1508 return ret_val;
1509 }
1510 if (*incx == 1 && *incy == 1) {
1511 goto L20;
1512 }
1513
1514 /* code for unequal increments or equal increments */
1515 /* not equal to 1 */
1516
1517 ix = 1;
1518 iy = 1;
1519 if (*incx < 0) {
1520 ix = (-(*n) + 1) * *incx + 1;
1521 }
1522 if (*incy < 0) {
1523 iy = (-(*n) + 1) * *incy + 1;
1524 }
1525 i__1 = *n;
1526 for (i = 1; i <= i__1; ++i) {
1527 dtemp += dx[ix] * dy[iy];
1528 ix += *incx;
1529 iy += *incy;
1530 /* L10: */
1531 }
1532 ret_val = dtemp;
1533 return ret_val;
1534
1535 /* code for both increments equal to 1 */
1536
1537
1538 /* clean-up loop */
1539
1540 L20:
1541 m = *n % 5;
1542 if (m == 0) {
1543 goto L40;
1544 }
1545 i__1 = m;
1546 for (i = 1; i <= i__1; ++i) {
1547 dtemp += dx[i] * dy[i];
1548 /* L30: */
1549 }
1550 if (*n < 5) {
1551 goto L60;
1552 }
1553 L40:
1554 mp1 = m + 1;
1555 i__1 = *n;
1556 for (i = mp1; i <= i__1; i += 5) {
1557 dtemp = dtemp + dx[i] * dy[i] + dx[i + 1] * dy[i + 1] + dx[i + 2] *
1558 dy[i + 2] + dx[i + 3] * dy[i + 3] + dx[i + 4] * dy[i + 4];
1559 /* L50: */
1560 }
1561 L60:
1562 ret_val = dtemp;
1563 return ret_val;
1564 } /* ddot_ */
1565
1566 /* end ddot.c */
1567
1568 /**********************************************************************/
1569 /******************* implementation code from divp.c *********************/
1570 /**********************************************************************/
1571
1572 /* Table of constant values */
1573
1574 static integer c__2 = 2;
1575
divp_(doublereal * xxxx,doublereal * yyyy,doublereal * zzzz,integer * ierr)1576 /* Subroutine */ int divp_(doublereal *xxxx,
1577 doublereal *yyyy,
1578 doublereal *zzzz,
1579 integer *ierr)
1580 {
1581 static doublereal xnum, denom;
1582 extern doublereal d1mach_(integer *i);
1583
1584
1585 /* THIS SUBROUTINE PERFORMS DIVISION OF COMPLEX NUMBERS: */
1586 /* ZZZZ = XXXX/YYYY */
1587
1588 /* ON INPUT: */
1589
1590 /* XXXX IS AN ARRAY OF LENGTH TWO REPRESENTING THE FIRST COMPLEX */
1591 /* NUMBER, WHERE XXXX(1) = REAL PART OF XXXX AND XXXX(2) = */
1592 /* IMAGINARY PART OF XXXX. */
1593
1594 /* YYYY IS AN ARRAY OF LENGTH TWO REPRESENTING THE SECOND COMPLEX */
1595 /* NUMBER, WHERE YYYY(1) = REAL PART OF YYYY AND YYYY(2) = */
1596 /* IMAGINARY PART OF YYYY. */
1597
1598 /* ON OUTPUT: */
1599
1600 /* ZZZZ IS AN ARRAY OF LENGTH TWO REPRESENTING THE RESULT OF */
1601 /* THE DIVISION, ZZZZ = XXXX/YYYY, WHERE ZZZZ(1) = */
1602 /* REAL PART OF ZZZZ AND ZZZZ(2) = IMAGINARY PART OF ZZZZ. */
1603
1604 /* IERR = */
1605 /* 1 IF DIVISION WOULD HAVE CAUSED OVERFLOW. IN THIS CASE, THE */
1606 /* APPROPRIATE PARTS OF ZZZZ ARE SET EQUAL TO THE LARGEST */
1607 /* FLOATING POINT NUMBER, AS GIVEN BY FUNCTION D1MACH . */
1608
1609 /* 0 IF DIVISION DOES NOT CAUSE OVERFLOW. */
1610
1611 /* DECLARATION OF INPUT */
1612
1613 /* DECLARATION OF OUTPUT */
1614
1615 /* DECLARATION OF VARIABLES */
1616
1617 /* Parameter adjustments */
1618 --zzzz;
1619 --yyyy;
1620 --xxxx;
1621
1622 /* Function Body */
1623 *ierr = 0;
1624 denom = yyyy[1] * yyyy[1] + yyyy[2] * yyyy[2];
1625 xnum = xxxx[1] * yyyy[1] + xxxx[2] * yyyy[2];
1626 if (dblabs(denom) >= (float)1. || (dblabs(denom) < (float)1. && dblabs(xnum) /
1627 d1mach_(&c__2) < dblabs(denom)) ) {
1628 zzzz[1] = xnum / denom;
1629 } else {
1630 zzzz[1] = d1mach_(&c__2);
1631 *ierr = 1;
1632 }
1633 xnum = xxxx[2] * yyyy[1] - xxxx[1] * yyyy[2];
1634 if (dblabs(denom) >= (float)1. || (dblabs(denom) < (float)1. && dblabs(xnum) /
1635 d1mach_(&c__2) < dblabs(denom))) {
1636 zzzz[2] = xnum / denom;
1637 } else {
1638 zzzz[2] = d1mach_(&c__2);
1639 *ierr = 1;
1640 }
1641 return 0;
1642 } /* divp_ */
1643
1644 /* end divp.c */
1645
1646 /**********************************************************************/
1647 /****************** implementation code from dnrm2.c ******************/
1648 /**********************************************************************/
1649
dnrm2_(integer * n,doublereal * dx,integer * incx)1650 doublereal dnrm2_(integer *n,
1651 doublereal *dx,
1652 integer *incx)
1653 {
1654 /* Initialized data */
1655
1656 static doublereal zero = 0.;
1657 static doublereal one = 1.;
1658 static doublereal cutlo = 8.232e-11;
1659 static doublereal cuthi = 1.304e19;
1660
1661 /* Format strings */
1662 static char fmt_30[] = "";
1663 static char fmt_50[] = "";
1664 static char fmt_70[] = "";
1665 static char fmt_110[] = "";
1666
1667 /* System generated locals */
1668 integer i__1;
1669 doublereal ret_val, d__1;
1670
1671 /* Builtin functions */
1672 /* double sqrt(double); CANT DECLARE BUILTINS UNDER C++ */
1673
1674 /* Local variables */
1675 static doublereal xmax;
1676 static integer next, i, j, ix;
1677 static doublereal hitest, sum;
1678
1679 /* Assigned format variables */
1680 char *next_fmt;
1681
1682 /* Parameter adjustments */
1683 --dx;
1684
1685 /* Function Body */
1686
1687 /* euclidean norm of the n-vector stored in dx() with storage */
1688 /* increment incx . */
1689 /* if n .le. 0 return with result = 0. */
1690 /* if n .ge. 1 then incx must be .ge. 1 */
1691
1692 /* c.l.lawson, 1978 jan 08 */
1693 /* modified to correct failure to update ix, 1/25/92. */
1694 /* modified 3/93 to return if incx .le. 0. */
1695
1696 /* four phase method using two built-in constants that are */
1697 /* hopefully applicable to all machines. */
1698 /* cutlo = maximum of dsqrt(u/eps) over all known machines. */
1699 /* cuthi = minimum of dsqrt(v) over all known machines. */
1700 /* where */
1701 /* eps = smallest no. such that eps + 1. .gt. 1. */
1702 /* u = smallest positive no. (underflow limit) */
1703 /* v = largest no. (overflow limit) */
1704
1705 /* brief outline of algorithm.. */
1706
1707 /* phase 1 scans zero components. */
1708 /* move to phase 2 when a component is nonzero and .le. cutlo */
1709 /* move to phase 3 when a component is .gt. cutlo */
1710 /* move to phase 4 when a component is .ge. cuthi/m */
1711 /* where m = n for x() real and m = 2*n for complex. */
1712
1713 /* values for cutlo and cuthi.. */
1714 /* from the environmental parameters listed in the imsl converter */
1715 /* document the limiting values are as follows.. */
1716 /* cutlo, s.p. u/eps = 2**(-102) for honeywell. close seconds are
1717 */
1718 /* univac and dec at 2**(-103) */
1719 /* thus cutlo = 2**(-51) = 4.44089e-16 */
1720 /* cuthi, s.p. v = 2**127 for univac, honeywell, and dec. */
1721 /* thus cuthi = 2**(63.5) = 1.30438e19 */
1722 /* cutlo, d.p. u/eps = 2**(-67) for honeywell and dec. */
1723 /* thus cutlo = 2**(-33.5) = 8.23181d-11 */
1724 /* cuthi, d.p. same as s.p. cuthi = 1.30438d19 */
1725 /* data cutlo, cuthi / 8.232d-11, 1.304d19 / */
1726 /* data cutlo, cuthi / 4.441e-16, 1.304e19 / */
1727
1728 if (*n > 0 && *incx > 0) {
1729 goto L10;
1730 }
1731 ret_val = zero;
1732 goto L300;
1733
1734 L10:
1735 next = 0;
1736 next_fmt = fmt_30;
1737 sum = zero;
1738 i = 1;
1739 ix = 1;
1740 /* begin main loop */
1741 L20:
1742 switch ((int)next) {
1743 case 0: goto L30;
1744 case 1: goto L50;
1745 case 2: goto L70;
1746 case 3: goto L110;
1747 }
1748 L30:
1749 if ((d__1 = dx[i], dblabs(d__1)) > cutlo) {
1750 goto L85;
1751 }
1752 next = 1;
1753 next_fmt = fmt_50;
1754 xmax = zero;
1755
1756 /* phase 1. sum is zero */
1757
1758 L50:
1759 if (dx[i] == zero) {
1760 goto L200;
1761 }
1762 if ((d__1 = dx[i], dblabs(d__1)) > cutlo) {
1763 goto L85;
1764 }
1765
1766 /* prepare for phase 2. */
1767 next = 2;
1768 next_fmt = fmt_70;
1769 goto L105;
1770
1771 /* prepare for phase 4. */
1772
1773 L100:
1774 ix = j;
1775 next = 3;
1776 next_fmt = fmt_110;
1777 sum = sum / dx[i] / dx[i];
1778 L105:
1779 xmax = (d__1 = dx[i], dblabs(d__1));
1780 goto L115;
1781
1782 /* phase 2. sum is small. */
1783 /* scale to avoid destructive underflow. */
1784
1785 L70:
1786 if ((d__1 = dx[i], dblabs(d__1)) > cutlo) {
1787 goto L75;
1788 }
1789
1790 /* common code for phases 2 and 4. */
1791 /* in phase 4 sum is large. scale to avoid overflow.
1792 */
1793
1794 L110:
1795 if ((d__1 = dx[i], dblabs(d__1)) <= xmax) {
1796 goto L115;
1797 }
1798 /* Computing 2nd power */
1799 d__1 = xmax / dx[i];
1800 sum = one + sum * (d__1 * d__1);
1801 xmax = (d__1 = dx[i], dblabs(d__1));
1802 goto L200;
1803
1804 L115:
1805 /* Computing 2nd power */
1806 d__1 = dx[i] / xmax;
1807 sum += d__1 * d__1;
1808 goto L200;
1809
1810
1811 /* prepare for phase 3. */
1812
1813 L75:
1814 sum = sum * xmax * xmax;
1815
1816
1817 /* for real or d.p. set hitest = cuthi/n */
1818 /* for complex set hitest = cuthi/(2*n) */
1819
1820 L85:
1821 hitest = cuthi / (real) (*n);
1822
1823 /* phase 3. sum is mid-range. no scaling. */
1824
1825 i__1 = *n;
1826 for (j = ix; j <= i__1; ++j) {
1827 if ((d__1 = dx[i], dblabs(d__1)) >= hitest) {
1828 goto L100;
1829 }
1830 /* Computing 2nd power */
1831 d__1 = dx[i];
1832 sum += d__1 * d__1;
1833 i += *incx;
1834 /* L95: */
1835 }
1836 ret_val = sqrt(sum);
1837 goto L300;
1838
1839 L200:
1840 ++ix;
1841 i += *incx;
1842 if (ix <= *n) {
1843 goto L20;
1844 }
1845
1846 /* end of main loop. */
1847
1848 /* compute square root and adjust for scaling. */
1849
1850 ret_val = xmax * sqrt(sum);
1851 L300:
1852 return ret_val;
1853 } /* dnrm2_ */
1854
1855 /* end dnrm2.c */
1856
1857 /**********************************************************************/
1858 /***************** implementation code from dscal.c *******************/
1859 /**********************************************************************/
1860
dscal_(integer * n,doublereal * da,doublereal * dx,integer * incx)1861 /* Subroutine */ int dscal_(integer *n,
1862 doublereal *da,
1863 doublereal *dx,
1864 integer *incx)
1865 {
1866 /* System generated locals */
1867 integer i__1, i__2;
1868
1869 /* Local variables */
1870 static integer i, m, nincx, mp1;
1871
1872
1873 /* scales a vector by a constant. */
1874 /* uses unrolled loops for increment equal to one. */
1875 /* jack dongarra, linpack, 3/11/78. */
1876 /* modified 3/93 to return if incx .le. 0. */
1877
1878
1879 /* Parameter adjustments */
1880 --dx;
1881
1882 /* Function Body */
1883 if (*n <= 0 || *incx <= 0) {
1884 return 0;
1885 }
1886 if (*incx == 1) {
1887 goto L20;
1888 }
1889
1890 /* code for increment not equal to 1 */
1891
1892 nincx = *n * *incx;
1893 i__1 = nincx;
1894 i__2 = *incx;
1895 for (i = 1; i__2 < 0 ? i >= i__1 : i <= i__1; i += i__2) {
1896 dx[i] = *da * dx[i];
1897 /* L10: */
1898 }
1899 return 0;
1900
1901 /* code for increment equal to 1 */
1902
1903
1904 /* clean-up loop */
1905
1906 L20:
1907 m = *n % 5;
1908 if (m == 0) {
1909 goto L40;
1910 }
1911 i__2 = m;
1912 for (i = 1; i <= i__2; ++i) {
1913 dx[i] = *da * dx[i];
1914 /* L30: */
1915 }
1916 if (*n < 5) {
1917 return 0;
1918 }
1919 L40:
1920 mp1 = m + 1;
1921 i__2 = *n;
1922 for (i = mp1; i <= i__2; i += 5) {
1923 dx[i] = *da * dx[i];
1924 dx[i + 1] = *da * dx[i + 1];
1925 dx[i + 2] = *da * dx[i + 2];
1926 dx[i + 3] = *da * dx[i + 3];
1927 dx[i + 4] = *da * dx[i + 4];
1928 /* L50: */
1929 }
1930 return 0;
1931 } /* dscal_ */
1932
1933 /* end dscal.c */
1934
1935 /**********************************************************************/
1936 /******************* implementation code from f.c *********************/
1937 /**********************************************************************/
1938
f_(doublereal * x,doublereal * v)1939 /* Subroutine */ int f_(doublereal *x,
1940 doublereal *v)
1941 {
1942
1943 /* EVALUATE F(X) AND RETURN IN THE VECTOR V . */
1944
1945 /* Parameter adjustments */
1946 --v;
1947 --x;
1948
1949 /* Function Body */
1950 return 0;
1951 } /* f_ */
1952
1953 /* end f.c */
1954
1955 /**********************************************************************/
1956 /***************** implementation code from fjac.c ********************/
1957 /**********************************************************************/
1958
fjac_(doublereal * x,doublereal * v,integer * k)1959 /* Subroutine */ int fjac_(doublereal *x,
1960 doublereal *v,
1961 integer *k)
1962 {
1963
1964 /* RETURN IN V THE KTH COLUMN OF THE JACOBIAN MATRIX OF */
1965 /* F(X) EVALUATED AT X . */
1966
1967 /* Parameter adjustments */
1968 --v;
1969 --x;
1970
1971 /* Function Body */
1972 return 0;
1973 } /* fjac_ */
1974
1975 /* end fjac.c */
1976
1977 /**********************************************************************/
1978 /***************** implementation code from idamax.c ******************/
1979 /**********************************************************************/
1980
idamax_(integer * n,doublereal * dx,integer * incx)1981 integer idamax_(integer *n,
1982 doublereal *dx,
1983 integer *incx)
1984 {
1985 /* System generated locals */
1986 integer ret_val, i__1;
1987 doublereal d__1;
1988
1989 /* Local variables */
1990 static doublereal dmax_;
1991 static integer i, ix;
1992
1993
1994 /* finds the index of element having max. absolute value. */
1995 /* jack dongarra, linpack, 3/11/78. */
1996 /* modified 3/93 to return if incx .le. 0. */
1997
1998
1999 /* Parameter adjustments */
2000 --dx;
2001
2002 /* Function Body */
2003 ret_val = 0;
2004 if (*n < 1 || *incx <= 0) {
2005 return ret_val;
2006 }
2007 ret_val = 1;
2008 if (*n == 1) {
2009 return ret_val;
2010 }
2011 if (*incx == 1) {
2012 goto L20;
2013 }
2014
2015 /* code for increment not equal to 1 */
2016
2017 ix = 1;
2018 dmax_ = dblabs(dx[1]);
2019 ix += *incx;
2020 i__1 = *n;
2021 for (i = 2; i <= i__1; ++i) {
2022 if ((d__1 = dx[ix], dblabs(d__1)) <= dmax_) {
2023 goto L5;
2024 }
2025 ret_val = i;
2026 dmax_ = (d__1 = dx[ix], dblabs(d__1));
2027 L5:
2028 ix += *incx;
2029 /* L10: */
2030 }
2031 return ret_val;
2032
2033 /* code for increment equal to 1 */
2034
2035 L20:
2036 dmax_ = dblabs(dx[1]);
2037 i__1 = *n;
2038 for (i = 2; i <= i__1; ++i) {
2039 if ((d__1 = dx[i], dblabs(d__1)) <= dmax_) {
2040 goto L30;
2041 }
2042 ret_val = i;
2043 dmax_ = (d__1 = dx[i], dblabs(d__1));
2044 L30:
2045 ;
2046 }
2047 return ret_val;
2048 } /* idamax_ */
2049
2050 /* end idamax.c */
2051
2052 /**********************************************************************/
2053 /****************** implementation code from root.c *******************/
2054 /**********************************************************************/
2055
d_sign(double * arg1,double * arg2)2056 double d_sign(double *arg1, double *arg2){
2057 if (*arg2>=0.0) return dblabs(*arg1);
2058 else return -1.0*dblabs(*arg1);
2059 }
2060
2061 /* Table of constant values */
2062
2063 /* static integer c__4 = 4; NOW REDUNDANT */
2064 static doublereal c_b17 = 1.;
2065
root_(doublereal * t,doublereal * ft,doublereal * b,doublereal * c,doublereal * relerr,doublereal * abserr,integer * iflag)2066 /* Subroutine */ int root_(doublereal *t,
2067 doublereal *ft,
2068 doublereal *b,
2069 doublereal *c,
2070 doublereal *relerr,
2071 doublereal *abserr,
2072 integer *iflag)
2073 {
2074 /* System generated locals */
2075 doublereal d__1, d__2;
2076
2077 /* Builtin functions */
2078 double d_sign(double *arg1, double *arg2);
2079
2080 /* Local variables */
2081 static doublereal acmb, acbs, a, p, q, u;
2082 extern doublereal d1mach_(integer *i);
2083 static integer kount;
2084 static doublereal ae, fa, fb, fc;
2085 static integer ic;
2086 static doublereal re, fx, cmb, tol;
2087
2088
2089 /* ROOT COMPUTES A ROOT OF THE NONLINEAR EQUATION F(X)=0 */
2090 /* WHERE F(X) IS A CONTINOUS REAL FUNCTION OF A SINGLE REAL */
2091 /* VARIABLE X. THE METHOD USED IS A COMBINATION OF BISECTION */
2092 /* AND THE SECANT RULE. */
2093
2094 /* NORMAL INPUT CONSISTS OF A CONTINUOS FUNCTION F AND AN */
2095 /* INTERVAL (B,C) SUCH THAT F(B)*F(C).LE.0.0. EACH ITERATION */
2096 /* FINDS NEW VALUES OF B AND C SUCH THAT THE INTERVAL(B,C) IS */
2097 /* SHRUNK AND F(B)*F(C).LE.0.0. THE STOPPING CRITERION IS */
2098
2099 /* DABS(B-C).LE.2.0*(RELERR*DABS(B)+ABSERR) */
2100
2101 /* WHERE RELERR=RELATIVE ERROR AND ABSERR=ABSOLUTE ERROR ARE */
2102 /* INPUT QUANTITIES. SET THE FLAG, IFLAG, POSITIVE TO INITIALIZE */
2103 /* THE COMPUTATION. AS B,C AND IFLAG ARE USED FOR BOTH INPUT AND */
2104 /* OUTPUT, THEY MUST BE VARIABLES IN THE CALLING PROGRAM. */
2105
2106 /* IF 0 IS A POSSIBLE ROOT, ONE SHOULD NOT CHOOSE ABSERR=0.0. */
2107
2108 /* THE OUTPUT VALUE OF B IS THE BETTER APPROXIMATION TO A ROOT */
2109 /* AS B AND C ARE ALWAYS REDEFINED SO THAT DABS(F(B)).LE.DABS(F(C)). */
2110
2111 /* TO SOLVE THE EQUATION, ROOT MUST EVALUATE F(X) REPEATEDLY. THIS */
2112 /* IS DONE IN THE CALLING PROGRAM. WHEN AN EVALUATION OF F IS */
2113 /* NEEDED AT T, ROOT RETURNS WITH IFLAG NEGATIVE. EVALUATE FT=F(T) */
2114 /* AND CALL ROOT AGAIN. DO NOT ALTER IFLAG. */
2115
2116 /* WHEN THE COMPUTATION IS COMPLETE, ROOT RETURNS TO THE CALLING */
2117 /* PROGRAM WITH IFLAG POSITIVE= */
2118
2119 /* IFLAG=1 IF F(B)*F(C).LT.0 AND THE STOPPING CRITERION IS MET. */
2120
2121 /* =2 IF A VALUE B IS FOUND SUCH THAT THE COMPUTED VALUE */
2122 /* F(B) IS EXACTLY ZERO. THE INTERVAL (B,C) MAY NOT */
2123 /* SATISFY THE STOPPING CRITERION. */
2124
2125 /* =3 IF DABS(F(B)) EXCEEDS THE INPUT VALUES DABS(F(B)), */
2126 /* DABS(F(C)). IN THIS CASE IT IS LIKELY THAT B IS CLOSE */
2127 /* TO A POLE OF F. */
2128
2129 /* =4 IF NO ODD ORDER ROOT WAS FOUND IN THE INTERVAL. A */
2130 /* LOCAL MINIMUM MAY HAVE BEEN OBTAINED. */
2131
2132 /* =5 IF TOO MANY FUNCTION EVALUATIONS WERE MADE. */
2133 /* (AS PROGRAMMED, 500 ARE ALLOWED.) */
2134
2135 /* THIS CODE IS A MODIFICATION OF THE CODE ZEROIN WHICH IS COMPLETELY */
2136 /* EXPLAINED AND DOCUMENTED IN THE TEXT NUMERICAL COMPUTING: AN */
2137 /* INTRODUCTION, BY L. F. SHAMPINE AND R. C. ALLEN. */
2138
2139 /* CALLS D1MACH . */
2140
2141
2142 if (*iflag >= 0) {
2143 goto L100;
2144 }
2145 *iflag = dblabs(*iflag);
2146 switch ((int)*iflag) {
2147 case 1: goto L200;
2148 case 2: goto L300;
2149 case 3: goto L400;
2150 }
2151 L100:
2152 u = d1mach_(&c__4);
2153 re = max(*relerr,u);
2154 ae = max(*abserr,0.);
2155 ic = 0;
2156 acbs = (d__1 = *b - *c, dblabs(d__1));
2157 a = *c;
2158 *t = a;
2159 *iflag = -1;
2160 return 0;
2161 L200:
2162 fa = *ft;
2163 *t = *b;
2164 *iflag = -2;
2165 return 0;
2166 L300:
2167 fb = *ft;
2168 fc = fa;
2169 kount = 2;
2170 /* Computing MAX */
2171 d__1 = dblabs(fb), d__2 = dblabs(fc);
2172 fx = max(d__1,d__2);
2173 L1:
2174 if (dblabs(fc) >= dblabs(fb)) {
2175 goto L2;
2176 }
2177
2178 /* INTERCHANGE B AND C SO THAT ABS(F(B)).LE.ABS(F(C)). */
2179
2180 a = *b;
2181 fa = fb;
2182 *b = *c;
2183 fb = fc;
2184 *c = a;
2185 fc = fa;
2186 L2:
2187 cmb = (*c - *b) * (float).5;
2188 acmb = dblabs(cmb);
2189 tol = re * dblabs(*b) + ae;
2190
2191 /* TEST STOPPING CRITERION AND FUNCTION COUNT. */
2192
2193 if (acmb <= tol) {
2194 goto L8;
2195 }
2196 if (kount >= 500) {
2197 goto L12;
2198 }
2199
2200 /* CALCULATE NEW ITERATE EXPLICITLY AS B+P/Q */
2201 /* WHERE WE ARRANGE P.GE.0. THE IMPLICIT */
2202 /* FORM IS USED TO PREVENT OVERFLOW. */
2203
2204 p = (*b - a) * fb;
2205 q = fa - fb;
2206 if (p >= (float)0.) {
2207 goto L3;
2208 }
2209 p = -p;
2210 q = -q;
2211
2212 /* UPDATE A, CHECK IF REDUCTION IN THE SIZE OF BRACKETING */
2213 /* INTERVAL IS SATISFACTORY. IF NOT BISECT UNTIL IT IS. */
2214
2215 L3:
2216 a = *b;
2217 fa = fb;
2218 ++ic;
2219 if (ic < 4) {
2220 goto L4;
2221 }
2222 if (acmb * (float)8. >= acbs) {
2223 goto L6;
2224 }
2225 ic = 0;
2226 acbs = acmb;
2227
2228 /* TEST FOR TOO SMALL A CHANGE. */
2229
2230 L4:
2231 if (p > dblabs(q) * tol) {
2232 goto L5;
2233 }
2234
2235 /* INCREMENT BY TOLERANCE */
2236
2237 /* sign gives sign(arg2)*|arg1| */
2238 *b += d_sign(&tol, &cmb);
2239 goto L7;
2240
2241 /* ROOT OUGHT TO BE BETWEEN B AND (C+B)/2 */
2242
2243 L5:
2244 if (p >= cmb * q) {
2245 goto L6;
2246 }
2247
2248 /* USE SECANT RULE. */
2249
2250 *b += p / q;
2251 goto L7;
2252
2253 /* USE BISECTION. */
2254
2255 L6:
2256 *b = (*c + *b) * (float).5;
2257
2258 /* HAVE COMPLETED COMPUTATION FOR NEW ITERATE B. */
2259
2260 L7:
2261 *t = *b;
2262 *iflag = -3;
2263 return 0;
2264 L400:
2265 fb = *ft;
2266 if (fb == (float)0.) {
2267 goto L9;
2268 }
2269 ++kount;
2270 if (d_sign(&c_b17, &fb) != d_sign(&c_b17, &fc)) {
2271 goto L1;
2272 }
2273 *c = a;
2274 fc = fa;
2275 goto L1;
2276
2277 /* FINISHED. SET IFLAG. */
2278
2279 L8:
2280 if (d_sign(&c_b17, &fb) == d_sign(&c_b17, &fc)) {
2281 goto L11;
2282 }
2283 if (dblabs(fb) > fx) {
2284 goto L10;
2285 }
2286 *iflag = 1;
2287 return 0;
2288 L9:
2289 *iflag = 2;
2290 return 0;
2291 L10:
2292 *iflag = 3;
2293 return 0;
2294 L11:
2295 *iflag = 4;
2296 return 0;
2297 L12:
2298 *iflag = 5;
2299 return 0;
2300 } /* root_ */
2301
2302 /* end root.c */
2303
2304 /**********************************************************************/
2305 /***************** implementation code from rootnf.c ******************/
2306 /**********************************************************************/
2307
2308 /* Table of constant values */
2309
2310 /*
2311 static int c__4 = 4;
2312 static int c__1 = 1;
2313 */
2314
rootnf_(int * n,int * nfe,int * iflag,double * relerr,double * abserr,double * y,double * yp,double * yold,double * ypold,double * a,double * qr,double * alpha,double * tz,int * pivot,double * w,double * wp,double * par,int * ipar)2315 /* Subroutine */ int rootnf_(int *n,
2316 int *nfe,
2317 int *iflag,
2318 double *relerr,
2319 double *abserr,
2320 double *y,
2321 double *yp,
2322 double *yold,
2323 double *ypold,
2324 double *a,
2325 double *qr,
2326 double *alpha,
2327 double *tz,
2328 int *pivot,
2329 double *w,
2330 double *wp,
2331 double *par,
2332 int *ipar)
2333 {
2334 /* System generated locals */
2335 int qr_dim1, qr_offset, i__1;
2336 double d__1;
2337
2338 /* Local variables */
2339 static double dels, aerr, rerr;
2340 static int judy;
2341 extern /* Subroutine */ int root_(doublereal *t,
2342 doublereal *ft,
2343 doublereal *b,
2344 doublereal *c,
2345 doublereal *relerr,
2346 doublereal *abserr,
2347 integer *iflag);
2348 static double sout;
2349 extern double dnrm2_(integer *n,
2350 doublereal *dx,
2351 integer *incx);
2352 static double u;
2353 static int lcode;
2354 extern double d1mach_(integer *);
2355 static double qsout, sa, sb;
2356 static int jw;
2357 /* extern */ /* Subroutine */ /* int tangnf_(); IN pelutils.h */
2358 static int np1;
2359
2360
2361 /* ROOTNF FINDS THE POINT YBAR = (1, XBAR) ON THE ZERO CURVE OF THE */
2362 /* HOMOTOPY MAP. IT STARTS WITH TWO POINTS YOLD=(LAMBDAOLD,XOLD) AND */
2363 /* Y=(LAMBDA,X) SUCH THAT LAMBDAOLD < 1 <= LAMBDA , AND ALTERNATES */
2364 /* BETWEEN HERMITE CUBIC INTERPOLATION AND NEWTON ITERATION UNTIL */
2365 /* CONVERGENCE. */
2366
2367 /* ON INPUT: */
2368
2369 /* N = DIMENSION OF X AND THE HOMOTOPY MAP. */
2370
2371 /* NFE = NUMBER OF JACOBIAN MATRIX EVALUATIONS. */
2372
2373 /* IFLAG = -2, -1, OR 0, INDICATING THE PROBLEM TYPE. */
2374
2375 /* RELERR, ABSERR = RELATIVE AND ABSOLUTE ERROR VALUES. THE ITERATION IS
2376 */
2377 /* CONSIDERED TO HAVE CONVERGED WHEN A POINT Y=(LAMBDA,X) IS FOUND */
2378 /* SUCH THAT */
2379
2380 /* |Y(1) - 1| <= RELERR + ABSERR AND */
2381
2382 /* ||Z|| <= RELERR*||X|| + ABSERR , WHERE */
2383
2384 /* (?,Z) IS THE NEWTON STEP TO Y=(LAMBDA,X). */
2385
2386 /* Y(1:N+1) = POINT (LAMBDA(S), X(S)) ON ZERO CURVE OF HOMOTOPY MAP. */
2387
2388 /* YP(1:N+1) = UNIT TANGENT VECTOR TO THE ZERO CURVE OF THE HOMOTOPY MAP
2389 */
2390 /* AT Y . */
2391
2392 /* YOLD(1:N+1) = A POINT DIFFERENT FROM Y ON THE ZERO CURVE. */
2393
2394 /* YPOLD(1:N+1) = UNIT TANGENT VECTOR TO THE ZERO CURVE OF THE HOMOTOPY */
2395 /* MAP AT YOLD . */
2396
2397 /* A(1:*) = PARAMETER VECTOR IN THE HOMOTOPY MAP. */
2398
2399 /* QR(1:N,1:N+2), ALPHA(1:N), TZ(1:N+1), PIVOT(1:N+1), W(1:N+1), */
2400 /* WP(1:N+1) ARE WORK ARRAYS USED FOR THE QR FACTORIZATION (IN THE */
2401 /* NEWTON STEP CALCULATION) AND THE INTERPOLATION. */
2402
2403 /* PAR(1:*) AND IPAR(1:*) ARE ARRAYS FOR (OPTIONAL) USER PARAMETERS, */
2404 /* WHICH ARE SIMPLY PASSED THROUGH TO THE USER WRITTEN SUBROUTINES */
2405 /* RHO, RHOJAC. */
2406
2407 /* ON OUTPUT: */
2408
2409 /* N , RELERR , ABSERR , A ARE UNCHANGED. */
2410
2411 /* NFE HAS BEEN UPDATED. */
2412
2413 /* IFLAG */
2414 /* = -2, -1, OR 0 (UNCHANGED) ON A NORMAL RETURN. */
2415
2416 /* = 4 IF A JACOBIAN MATRIX WITH RANK < N HAS OCCURRED. THE */
2417 /* ITERATION WAS NOT COMPLETED. */
2418
2419 /* = 6 IF THE ITERATION FAILED TO CONVERGE. Y AND YOLD CONTAIN */
2420 /* THE LAST TWO POINTS FOUND ON THE ZERO CURVE. */
2421
2422 /* Y IS THE POINT ON THE ZERO CURVE OF THE HOMOTOPY MAP AT LAMBDA = 1 .
2423 */
2424
2425
2426 /* CALLS D1MACH , DNRM2 , ROOT , TANGNF . */
2427
2428
2429 /* ***** ARRAY DECLARATIONS. ***** */
2430
2431
2432 /* ***** END OF DIMENSIONAL INFORMATION. ***** */
2433
2434 /* THE LIMIT ON THE NUMBER OF ITERATIONS ALLOWED MAY BE CHANGED BY */
2435 /* CHANGING THE FOLLOWING PARAMETER STATEMENT: */
2436
2437 /* DEFINITION OF HERMITE CUBIC INTERPOLANT VIA DIVIDED DIFFERENCES. */
2438
2439
2440
2441 /* Parameter adjustments */
2442 --wp;
2443 --w;
2444 --pivot;
2445 --tz;
2446 --alpha;
2447 qr_dim1 = *n;
2448 qr_offset = qr_dim1 + 1;
2449 qr -= qr_offset;
2450 --a;
2451 --ypold;
2452 --yold;
2453 --yp;
2454 --y;
2455 --par;
2456 --ipar;
2457
2458 /* Function Body */
2459 u = d1mach_(&c__4);
2460 rerr = max(*relerr,u);
2461 aerr = max(*abserr,0.);
2462 np1 = *n + 1;
2463
2464 /* ***** MAIN LOOP. ***** */
2465
2466 /* L100: */
2467 for (judy = 1; judy <= 20; ++judy) {
2468 i__1 = np1;
2469 for (jw = 1; jw <= i__1; ++jw) {
2470 tz[jw] = y[jw] - yold[jw];
2471 /* L110: */
2472 }
2473 dels = dnrm2_((integer *)&np1, &tz[1], &c__1);
2474
2475 /* USING TWO POINTS AND TANGENTS ON THE HOMOTOPY ZERO CURVE, CONSTRUCT
2476 */
2477 /* THE HERMITE CUBIC INTERPOLANT Q(S). THEN USE ROOT TO FIND THE S
2478 */
2479 /* CORRESPONDING TO LAMBDA = 1 . THE TWO POINTS ON THE ZERO CURVE AR
2480 E */
2481 /* ALWAYS CHOSEN TO BRACKET LAMBDA=1, WITH THE BRACKETING INTERVAL */
2482 /* ALWAYS BEING [0, DELS]. */
2483
2484 sa = (float)0.;
2485 sb = dels;
2486 lcode = 1;
2487 L130:
2488 root_(&sout, &qsout, &sa, &sb, &rerr, &aerr, (integer *)&lcode);
2489 if (lcode > 0) {
2490 goto L140;
2491 }
2492 qsout = ((((yp[1] - (y[1] - yold[1]) / dels) / dels - ((y[1] - yold[1]
2493 ) / dels - ypold[1]) / dels) / dels * (sout - dels) + ((y[1]
2494 - yold[1]) / dels - ypold[1]) / dels) * sout + ypold[1]) *
2495 sout + yold[1] - (float)1.;
2496 goto L130;
2497 /* IF LAMBDA = 1 WERE BRACKETED, ROOT CANNOT FAIL. */
2498 L140:
2499 if (lcode > 2) {
2500 *iflag = 6;
2501 return 0;
2502 }
2503
2504 /* CALCULATE Q(SA) AS THE INITIAL POINT FOR A NEWTON ITERATION. */
2505 i__1 = np1;
2506 for (jw = 1; jw <= i__1; ++jw) {
2507 w[jw] = ((((yp[jw] - (y[jw] - yold[jw]) / dels) / dels - ((y[jw]
2508 - yold[jw]) / dels - ypold[jw]) / dels) / dels * (sa -
2509 dels) + ((y[jw] - yold[jw]) / dels - ypold[jw]) / dels) *
2510 sa + ypold[jw]) * sa + yold[jw];
2511 /* L150: */
2512 }
2513 /* CALCULATE NEWTON STEP AT Q(SA). */
2514 tangnf_(&sa, &w[1], &wp[1], &ypold[1], &a[1], &qr[qr_offset], &alpha[
2515 1], &tz[1], &pivot[1], nfe, n, iflag, &par[1], &ipar[1]);
2516 if (*iflag > 0) {
2517 return 0;
2518 }
2519 /* NEXT POINT = CURRENT POINT + NEWTON STEP. */
2520 i__1 = np1;
2521 for (jw = 1; jw <= i__1; ++jw) {
2522 w[jw] += tz[jw];
2523 /* L160: */
2524 }
2525 /* GET THE TANGENT WP AT W AND THE NEXT NEWTON STEP IN TZ . */
2526 tangnf_(&sa, &w[1], &wp[1], &ypold[1], &a[1], &qr[qr_offset], &alpha[
2527 1], &tz[1], &pivot[1], nfe, n, iflag, &par[1], &ipar[1]);
2528 if (*iflag > 0) {
2529 return 0;
2530 }
2531 /* TAKE NEWTON STEP AND CHECK CONVERGENCE. */
2532 i__1 = np1;
2533 for (jw = 1; jw <= i__1; ++jw) {
2534 w[jw] += tz[jw];
2535 /* L170: */
2536 }
2537 if ((d__1 = w[1] - (float)1., dblabs(d__1)) <= rerr + aerr && dnrm2_((integer *)n, &
2538 tz[2], &c__1) <= rerr * dnrm2_((integer*)n, &w[2], &c__1) + aerr) {
2539 i__1 = np1;
2540 for (jw = 1; jw <= i__1; ++jw) {
2541 y[jw] = w[jw];
2542 /* L180: */
2543 }
2544 return 0;
2545 }
2546 /* IF THE ITERATION HAS NOT CONVERGED, DISCARD ONE OF THE OLD POINTS
2547 */
2548 /* SUCH THAT LAMBDA = 1 IS STILL BRACKETED. */
2549 if ((yold[1] - (float)1.) * (w[1] - (float)1.) > (float)0.) {
2550 i__1 = np1;
2551 for (jw = 1; jw <= i__1; ++jw) {
2552 yold[jw] = w[jw];
2553 ypold[jw] = wp[jw];
2554 /* L200: */
2555 }
2556 } else {
2557 i__1 = np1;
2558 for (jw = 1; jw <= i__1; ++jw) {
2559 y[jw] = w[jw];
2560 yp[jw] = wp[jw];
2561 /* L210: */
2562 }
2563 }
2564 /* L300: */
2565 }
2566
2567 /* ***** END OF MAIN LOOP. ***** */
2568
2569 /* THE ALTERNATING OSCULATORY CUBIC INTERPOLATION AND NEWTON ITERATION */
2570 /* HAS NOT CONVERGED IN LIMIT STEPS. ERROR RETURN. */
2571 *iflag = 6;
2572 return 0;
2573 } /* rootnf_ */
2574
2575 /* end rootnf.c */
2576
2577 /**********************************************************************/
2578 /******************* implementation code from stepnf.c *********************/
2579 /**********************************************************************/
2580
2581 /* Table of constant values */
2582
2583 /*
2584 static integer c__4 = 4;
2585 static integer c__1 = 1;
2586 */
2587
stepnf_(integer * n,integer * nfe,integer * iflag,logical * start,logical * crash,doublereal * hold,doublereal * h,doublereal * relerr,doublereal * abserr,doublereal * s,doublereal * y,doublereal * yp,doublereal * yold,doublereal * ypold,doublereal * a,doublereal * qr,doublereal * alpha,doublereal * tz,integer * pivot,doublereal * w,doublereal * wp,doublereal * z0,doublereal * z1,doublereal * sspar,doublereal * par,integer * ipar)2588 /* Subroutine */ int stepnf_(integer *n,
2589 integer *nfe,
2590 integer *iflag,
2591 logical *start,
2592 logical *crash,
2593 doublereal *hold,
2594 doublereal *h,
2595 doublereal *relerr,
2596 doublereal *abserr,
2597 doublereal *s,
2598 doublereal *y,
2599 doublereal *yp,
2600 doublereal *yold,
2601 doublereal *ypold,
2602 doublereal *a,
2603 doublereal *qr,
2604 doublereal *alpha,
2605 doublereal *tz,
2606 integer *pivot,
2607 doublereal *w,
2608 doublereal *wp,
2609 doublereal *z0,
2610 doublereal *z1,
2611 doublereal *sspar,
2612 doublereal *par,
2613 integer *ipar)
2614 {
2615 /* System generated locals */
2616 integer qr_dim1, qr_offset, i__1;
2617 doublereal d__1, d__2, d__3, d__4;
2618
2619 /* Builtin functions */
2620 /* double sqrt(double), pow(double,double); CANT DECLARE BUILTINS UNDER C++ */
2621
2622 /* Local variables */
2623 static logical fail;
2624 static doublereal temp;
2625 static integer judy;
2626 static doublereal twou;
2627 extern doublereal dnrm2_(integer *n,
2628 doublereal *dx,
2629 integer *incx);
2630 static doublereal dcalc;
2631 static integer j;
2632 static doublereal lcalc, hfail, rcalc;
2633 static integer itnum;
2634 extern doublereal d1mach_(integer *i);
2635 static doublereal fouru, ht;
2636 /*extern */ /* Subroutine */ /* int tangnf_(); IN pelutils.h */
2637 static doublereal rholen;
2638 static integer np1;
2639
2640
2641 /* STEPNF TAKES ONE STEP ALONG THE ZERO CURVE OF THE HOMOTOPY MAP */
2642 /* USING A PREDICTOR-CORRECTOR ALGORITHM. THE PREDICTOR USES A HERMITE */
2643 /* CUBIC INTERPOLANT, AND THE CORRECTOR RETURNS TO THE ZERO CURVE ALONG */
2644 /* THE FLOW NORMAL TO THE DAVIDENKO FLOW. STEPNF ALSO ESTIMATES A */
2645 /* STEP SIZE H FOR THE NEXT STEP ALONG THE ZERO CURVE. NORMALLY */
2646 /* STEPNF IS USED INDIRECTLY THROUGH FIXPNF , AND SHOULD BE CALLED */
2647 /* DIRECTLY ONLY IF IT IS NECESSARY TO MODIFY THE STEPPING ALGORITHM'S */
2648 /* PARAMETERS. */
2649
2650 /* ON INPUT: */
2651
2652 /* N = DIMENSION OF X AND THE HOMOTOPY MAP. */
2653
2654 /* NFE = NUMBER OF JACOBIAN MATRIX EVALUATIONS. */
2655
2656 /* IFLAG = -2, -1, OR 0, INDICATING THE PROBLEM TYPE. */
2657
2658 /* START = .TRUE. ON FIRST CALL TO STEPNF , .FALSE. OTHERWISE. */
2659
2660 /* HOLD = ||Y - YOLD||; SHOULD NOT BE MODIFIED BY THE USER. */
2661
2662 /* H = UPPER LIMIT ON LENGTH OF STEP THAT WILL BE ATTEMPTED. H MUST BE
2663 */
2664 /* SET TO A POSITIVE NUMBER ON THE FIRST CALL TO STEPNF . */
2665 /* THEREAFTER STEPNF CALCULATES AN OPTIMAL VALUE FOR H , AND H */
2666 /* SHOULD NOT BE MODIFIED BY THE USER. */
2667
2668 /* RELERR, ABSERR = RELATIVE AND ABSOLUTE ERROR VALUES. THE ITERATION IS
2669 */
2670 /* CONSIDERED TO HAVE CONVERGED WHEN A POINT W=(LAMBDA,X) IS FOUND */
2671 /* SUCH THAT */
2672
2673 /* ||Z|| <= RELERR*||W|| + ABSERR , WHERE */
2674
2675 /* Z IS THE NEWTON STEP TO W=(LAMBDA,X). */
2676
2677 /* S = (APPROXIMATE) ARC LENGTH ALONG THE HOMOTOPY ZERO CURVE UP TO */
2678 /* Y(S) = (LAMBDA(S), X(S)). */
2679
2680 /* Y(1:N+1) = PREVIOUS POINT (LAMBDA(S), X(S)) FOUND ON THE ZERO CURVE OF
2681 */
2682 /* THE HOMOTOPY MAP. */
2683
2684 /* YP(1:N+1) = UNIT TANGENT VECTOR TO THE ZERO CURVE OF THE HOMOTOPY MAP
2685 */
2686 /* AT Y . */
2687
2688 /* YOLD(1:N+1) = A POINT BEFORE Y ON THE ZERO CURVE OF THE HOMOTOPY MAP.
2689 */
2690
2691 /* YPOLD(1:N+1) = UNIT TANGENT VECTOR TO THE ZERO CURVE OF THE HOMOTOPY */
2692 /* MAP AT YOLD . */
2693
2694 /* A(1:*) = PARAMETER VECTOR IN THE HOMOTOPY MAP. */
2695
2696 /* QR(1:N,1:N+2), ALPHA(1:N), TZ(1:N+1), PIVOT(1:N+1), W(1:N+1), */
2697 /* WP(1:N+1) ARE WORK ARRAYS USED FOR THE QR FACTORIZATION (IN THE */
2698 /* NEWTON STEP CALCULATION) AND THE INTERPOLATION. */
2699
2700 /* Z0(1:N+1), Z1(1:N+1) ARE WORK ARRAYS USED FOR THE ESTIMATION OF THE */
2701 /* NEXT STEP SIZE H . */
2702
2703 /* SSPAR(1:8) = (LIDEAL, RIDEAL, DIDEAL, HMIN, HMAX, BMIN, BMAX, P) IS */
2704 /* A VECTOR OF PARAMETERS USED FOR THE OPTIMAL STEP SIZE ESTIMATION. */
2705
2706 /* PAR(1:*) AND IPAR(1:*) ARE ARRAYS FOR (OPTIONAL) USER PARAMETERS, */
2707 /* WHICH ARE SIMPLY PASSED THROUGH TO THE USER WRITTEN SUBROUTINES */
2708 /* RHO, RHOJAC. */
2709
2710 /* ON OUTPUT: */
2711
2712 /* N , A , SSPAR ARE UNCHANGED. */
2713
2714 /* NFE HAS BEEN UPDATED. */
2715
2716 /* IFLAG */
2717 /* = -2, -1, OR 0 (UNCHANGED) ON A NORMAL RETURN. */
2718
2719 /* = 4 IF A JACOBIAN MATRIX WITH RANK < N HAS OCCURRED. THE */
2720 /* ITERATION WAS NOT COMPLETED. */
2721
2722 /* = 6 IF THE ITERATION FAILED TO CONVERGE. W CONTAINS THE LAST */
2723 /* NEWTON ITERATE. */
2724
2725 /* START = .FALSE. ON A NORMAL RETURN. */
2726
2727 /* CRASH */
2728 /* = .FALSE. ON A NORMAL RETURN. */
2729
2730 /* = .TRUE. IF THE STEP SIZE H WAS TOO SMALL. H HAS BEEN */
2731 /* INCREASED TO AN ACCEPTABLE VALUE, WITH WHICH STEPNF MAY BE */
2732 /* CALLED AGAIN. */
2733
2734 /* = .TRUE. IF RELERR AND/OR ABSERR WERE TOO SMALL. THEY HAVE */
2735 /* BEEN INCREASED TO ACCEPTABLE VALUES, WITH WHICH STEPNF MAY */
2736 /* BE CALLED AGAIN. */
2737
2738 /* HOLD = ||Y - YOLD||. */
2739
2740 /* H = OPTIMAL VALUE FOR NEXT STEP TO BE ATTEMPTED. NORMALLY H SHOULD
2741 */
2742 /* NOT BE MODIFIED BY THE USER. */
2743
2744 /* RELERR, ABSERR ARE UNCHANGED ON A NORMAL RETURN. */
2745
2746 /* S = (APPROXIMATE) ARC LENGTH ALONG THE ZERO CURVE OF THE HOMOTOPY MAP
2747 */
2748 /* UP TO THE LATEST POINT FOUND, WHICH IS RETURNED IN Y . */
2749
2750 /* Y, YP, YOLD, YPOLD CONTAIN THE TWO MOST RECENT POINTS AND TANGENT */
2751 /* VECTORS FOUND ON THE ZERO CURVE OF THE HOMOTOPY MAP. */
2752
2753
2754 /* CALLS D1MACH , DNRM2 , TANGNF . */
2755
2756
2757 /* ***** ARRAY DECLARATIONS. ***** */
2758
2759
2760 /* ***** END OF DIMENSIONAL INFORMATION. ***** */
2761
2762 /* THE LIMIT ON THE NUMBER OF NEWTON ITERATIONS ALLOWED BEFORE REDUCING */
2763 /* THE STEP SIZE H MAY BE CHANGED BY CHANGING THE FOLLOWING PARAMETER */
2764 /* STATEMENT: */
2765
2766 /* DEFINITION OF HERMITE CUBIC INTERPOLANT VIA DIVIDED DIFFERENCES. */
2767
2768
2769
2770 /* Parameter adjustments */
2771 --z1;
2772 --z0;
2773 --wp;
2774 --w;
2775 --pivot;
2776 --tz;
2777 --alpha;
2778 qr_dim1 = *n;
2779 qr_offset = qr_dim1 + 1;
2780 qr -= qr_offset;
2781 --a;
2782 --ypold;
2783 --yold;
2784 --yp;
2785 --y;
2786 --sspar;
2787 --par;
2788 --ipar;
2789
2790 /* Function Body */
2791 twou = d1mach_(&c__4) * (float)2.;
2792 fouru = twou + twou;
2793 np1 = *n + 1;
2794 *crash = TRUE_;
2795 /* THE ARCLENGTH S MUST BE NONNEGATIVE. */
2796 if (*s < (float)0.) {
2797 return 0;
2798 }
2799 /* IF STEP SIZE IS TOO SMALL, DETERMINE AN ACCEPTABLE ONE. */
2800 if (*h < fouru * (*s + (float)1.)) {
2801 *h = fouru * (*s + (float)1.);
2802 return 0;
2803 }
2804 /* IF ERROR TOLERANCES ARE TOO SMALL, INCREASE THEM TO ACCEPTABLE VALUES.
2805 */
2806 temp = dnrm2_(&np1, &y[1], &c__1);
2807 if ((*relerr * temp + *abserr) * (float).5 >= twou * temp) {
2808 goto L40;
2809 }
2810 if (*relerr != (float)0.) {
2811 *relerr = fouru * (fouru + (float)1.);
2812 *abserr = max(*abserr,0.);
2813 } else {
2814 *abserr = fouru * temp;
2815 }
2816 return 0;
2817 L40:
2818 *crash = FALSE_;
2819 if (! (*start)) {
2820 goto L300;
2821 }
2822
2823 /* ***** STARTUP SECTION(FIRST STEP ALONG ZERO CURVE. ***** */
2824
2825 fail = FALSE_;
2826 *start = FALSE_;
2827 /* DETERMINE SUITABLE INITIAL STEP SIZE. */
2828 /* Computing MIN */
2829 d__1 = min(*h,.1), d__2 = sqrt(sqrt(*relerr * temp + *abserr));
2830 *h = min(d__1,d__2);
2831 /* USE LINEAR PREDICTOR ALONG TANGENT DIRECTION TO START NEWTON ITERATION.
2832 */
2833 ypold[1] = (float)1.;
2834 i__1 = np1;
2835 for (j = 2; j <= i__1; ++j) {
2836 ypold[j] = (float)0.;
2837 /* L50: */
2838 }
2839 tangnf_(s, &y[1], &yp[1], &ypold[1], &a[1], &qr[qr_offset], &alpha[1],
2840 &tz[1], (int *)&pivot[1], (int *)nfe,
2841 (int *)n, (int *)iflag, &par[1], (int *)&ipar[1]);
2842 if (*iflag > 0) {
2843 return 0;
2844 }
2845 L70:
2846 i__1 = np1;
2847 for (j = 1; j <= i__1; ++j) {
2848 temp = y[j] + *h * yp[j];
2849 w[j] = temp;
2850 z0[j] = temp;
2851 /* L80: */
2852 }
2853 for (judy = 1; judy <= 4; ++judy) {
2854 rholen = (float)-1.;
2855 /* CALCULATE THE NEWTON STEP TZ AT THE CURRENT POINT W . */
2856 tangnf_(&rholen, &w[1], &wp[1], &ypold[1], &a[1], &qr[qr_offset],
2857 &alpha[1], &tz[1], (int *)&pivot[1], (int *)nfe,
2858 (int *)n, (int *)iflag, &par[1], (int *)&ipar[1])
2859 ;
2860 if (*iflag > 0) {
2861 return 0;
2862 }
2863
2864 /* TAKE NEWTON STEP AND CHECK CONVERGENCE. */
2865 i__1 = np1;
2866 for (j = 1; j <= i__1; ++j) {
2867 w[j] += tz[j];
2868 /* L90: */
2869 }
2870 itnum = judy;
2871 /* COMPUTE QUANTITIES USED FOR OPTIMAL STEP SIZE ESTIMATION. */
2872 if (judy == 1) {
2873 lcalc = dnrm2_(&np1, &tz[1], &c__1);
2874 rcalc = rholen;
2875 i__1 = np1;
2876 for (j = 1; j <= i__1; ++j) {
2877 z1[j] = w[j];
2878 /* L110: */
2879 }
2880 } else if (judy == 2) {
2881 lcalc = dnrm2_(&np1, &tz[1], &c__1) / lcalc;
2882 rcalc = rholen / rcalc;
2883 }
2884 /* GO TO MOP-UP SECTION AFTER CONVERGENCE. */
2885 if (dnrm2_(&np1, &tz[1], &c__1) <= *relerr * dnrm2_(&np1, &w[1], &
2886 c__1) + *abserr) {
2887 goto L600;
2888 }
2889
2890 /* L200: */
2891 }
2892
2893 /* NO CONVERGENCE IN LITFH ITERATIONS. REDUCE H AND TRY AGAIN. */
2894 if (*h <= fouru * (*s + (float)1.)) {
2895 *iflag = 6;
2896 return 0;
2897 }
2898 *h *= (float).5;
2899 goto L70;
2900
2901 /* ***** END OF STARTUP SECTION. ***** */
2902
2903 /* ***** PREDICTOR SECTION. ***** */
2904
2905 L300:
2906 fail = FALSE_;
2907 /* COMPUTE POINT PREDICTED BY HERMITE INTERPOLANT. USE STEP SIZE H */
2908 /* COMPUTED ON LAST CALL TO STEPNF . */
2909 L320:
2910 i__1 = np1;
2911 for (j = 1; j <= i__1; ++j) {
2912 d__1 = *hold + *h;
2913 temp = ((((yp[j] - (y[j] - yold[j]) / *hold) / *hold - ((y[j] - yold[
2914 j]) / *hold - ypold[j]) / *hold) / *hold * (d__1 - *hold) + ((
2915 y[j] - yold[j]) / *hold - ypold[j]) / *hold) * d__1 + ypold[j]
2916 ) * d__1 + yold[j];
2917 w[j] = temp;
2918 z0[j] = temp;
2919 /* L330: */
2920 }
2921
2922 /* ***** END OF PREDICTOR SECTION. ***** */
2923
2924 /* ***** CORRECTOR SECTION. ***** */
2925
2926 for (judy = 1; judy <= 4; ++judy) {
2927 rholen = (float)-1.;
2928 /* CALCULATE THE NEWTON STEP TZ AT THE CURRENT POINT W . */
2929 tangnf_(&rholen, &w[1], &wp[1], &yp[1], &a[1], &qr[qr_offset],
2930 &alpha[1], &tz[1], (int *)&pivot[1], (int *)nfe,
2931 (int *)n, (int *)iflag, &par[1], (int *)&ipar[1]);
2932 if (*iflag > 0) {
2933 return 0;
2934 }
2935
2936 /* TAKE NEWTON STEP AND CHECK CONVERGENCE. */
2937 i__1 = np1;
2938 for (j = 1; j <= i__1; ++j) {
2939 w[j] += tz[j];
2940 /* L420: */
2941 }
2942 itnum = judy;
2943 /* COMPUTE QUANTITIES USED FOR OPTIMAL STEP SIZE ESTIMATION. */
2944 if (judy == 1) {
2945 lcalc = dnrm2_(&np1, &tz[1], &c__1);
2946 rcalc = rholen;
2947 i__1 = np1;
2948 for (j = 1; j <= i__1; ++j) {
2949 z1[j] = w[j];
2950 /* L440: */
2951 }
2952 } else if (judy == 2) {
2953 lcalc = dnrm2_(&np1, &tz[1], &c__1) / lcalc;
2954 rcalc = rholen / rcalc;
2955 }
2956 /* GO TO MOP-UP SECTION AFTER CONVERGENCE. */
2957 if (dnrm2_(&np1, &tz[1], &c__1) <= *relerr * dnrm2_(&np1, &w[1], &
2958 c__1) + *abserr) {
2959 goto L600;
2960 }
2961
2962 /* L500: */
2963 }
2964
2965 /*NO CONVERGENCE IN LITFH ITERATIONS. RECORD FAILURE AT CALCULATED H ,
2966 */
2967 /* SAVE THIS STEP SIZE, REDUCE H AND TRY AGAIN. */
2968 fail = TRUE_;
2969 hfail = *h;
2970 if (*h <= fouru * (*s + (float)1.)) {
2971 *iflag = 6;
2972 return 0;
2973 }
2974 *h *= (float).5;
2975 goto L320;
2976
2977 /* ***** END OF CORRECTOR SECTION. ***** */
2978
2979 /* ***** MOP-UP SECTION. ***** */
2980
2981 /* YOLD AND Y ALWAYS CONTAIN THE LAST TWO POINTS FOUND ON THE ZERO */
2982 /* CURVE OF THE HOMOTOPY MAP. YPOLD AND YP CONTAIN THE TANGENT */
2983 /* VECTORS TO THE ZERO CURVE AT YOLD AND Y , RESPECTIVELY. */
2984
2985 L600:
2986 i__1 = np1;
2987 for (j = 1; j <= i__1; ++j) {
2988 yold[j] = y[j];
2989 ypold[j] = yp[j];
2990 y[j] = w[j];
2991 yp[j] = wp[j];
2992 w[j] = y[j] - yold[j];
2993 /* L620: */
2994 }
2995 /* UPDATE ARC LENGTH. */
2996 *hold = dnrm2_(&np1, &w[1], &c__1);
2997 *s += *hold;
2998
2999 /* ***** END OF MOP-UP SECTION. ***** */
3000
3001 /* ***** OPTIMAL STEP SIZE ESTIMATION SECTION. ***** */
3002
3003 /* CALCULATE THE DISTANCE FACTOR DCALC . */
3004 /* L700: */
3005 i__1 = np1;
3006 for (j = 1; j <= i__1; ++j) {
3007 tz[j] = z0[j] - y[j];
3008 w[j] = z1[j] - y[j];
3009 /* L710: */
3010 }
3011 dcalc = dnrm2_(&np1, &tz[1], &c__1);
3012 if (dcalc != (float)0.) {
3013 dcalc = dnrm2_(&np1, &w[1], &c__1) / dcalc;
3014 }
3015
3016 /* THE OPTIMAL STEP SIZE HBAR IS DEFINED BY */
3017
3018 /* HT=HOLD * [MIN(LIDEAL/LCALC, RIDEAL/RCALC, DIDEAL/DCALC)]**(1/P) */
3019
3020 /* HBAR = MIN [ MAX(HT, BMIN*HOLD, HMIN), BMAX*HOLD, HMAX ] */
3021
3022 /* IF CONVERGENCE HAD OCCURRED AFTER 1 ITERATION, SET THE CONTRACTION */
3023 /* FACTOR LCALC TO ZERO. */
3024 if (itnum == 1) {
3025 lcalc = (float)0.;
3026 }
3027 /* FORMULA FOR OPTIMAL STEP SIZE. */
3028 if (lcalc + rcalc + dcalc == (float)0.) {
3029 ht = sspar[7] * *hold;
3030 } else {
3031 /* Computing MAX */
3032 d__2 = lcalc / sspar[1], d__3 = rcalc / sspar[2], d__2 = max(d__2,
3033 d__3), d__3 = dcalc / sspar[3];
3034 d__1 = (float)1. / max(d__2,d__3);
3035 d__4 = (float)1. / sspar[8];
3036 /* ht = pow_dd(&d__1, &d__4) * *hold; */ ht=pow(d__1,d__4) * *hold;
3037 }
3038 /* HT CONTAINS THE ESTIMATED OPTIMAL STEP SIZE. NOW PUT IT WITHIN */
3039 /* REASONABLE BOUNDS. */
3040 /* Computing MIN */
3041 /* Computing MAX */
3042 d__3 = ht, d__4 = sspar[6] * *hold, d__3 = max(d__3,d__4);
3043 d__1 = max(d__3,sspar[4]), d__2 = sspar[7] * *hold, d__1 = min(d__1,d__2);
3044 *h = min(d__1,sspar[5]);
3045 if (itnum == 1) {
3046 /* IF CONVERGENCE HAD OCCURRED AFTER 1 ITERATION, DON'T DECREASE H .
3047 */
3048 *h = max(*h,*hold);
3049 } else if (itnum == 4) {
3050 /* IF CONVERGENCE REQUIRED THE MAXIMUM LITFH ITERATIONS, DON'T */
3051 /* INCREASE H . */
3052 *h = min(*h,*hold);
3053 }
3054 /* IF CONVERGENCE DID NOT OCCUR IN LITFH ITERATIONS FOR A PARTICULAR */
3055 /* H = HFAIL , DON'T CHOOSE THE NEW STEP SIZE LARGER THAN HFAIL . */
3056 if (fail) {
3057 *h = min(*h,hfail);
3058 }
3059
3060
3061 return 0;
3062 } /* stepnf_ */
3063
3064 /* end stepnf.c */
3065
3066 /**********************************************************************/
3067 /***************** implementation code from Pmatrix.c *****************/
3068 /**********************************************************************/
3069
3070 #define min(i,j) ((i) < (j) ? (i): (j))
3071 #define FALSE 0
3072 #define TRUE 1
3073 /*---------------------------------------------------------------
3074 vector/matrix type a linear array of int, whith auxilary info.
3075 *) the number of elements that can be stored is in elt[0]
3076 *) the current number of rows is in elt[1]
3077 *) the current number of collumbs is in elt[2]
3078 The actual data are then stored in row major order from elt[3] on
3079 ---------------------------------------------------------------*/
3080 struct Pmatrix_t {
3081 int store; /* maximum number of helts reserved */
3082 int topc; /* effective number of columbs */
3083 int topr; /* effective number of rows */
3084 int ncols; /* number of elts between first elts of rows */
3085 polynomial1 *coords;
3086 };
3087
3088 /*
3089 ** a couple of functions which are defined in files that
3090 ** I don't want to include right now, when things stableize
3091 ** bad_error -- used mainly to stop things when a function
3092 ** is passed bad input, as things stablizes
3093 ** less catastrophic reactions may be introduced.
3094 ** mem_alloc,
3095 ** mem_free -- call malloc and free with some extra book-keeping.
3096 */
3097 void bad_error(const char *);
3098
3099
3100 /* WEIRD - The first pair of declarations below are the global one,
3101 which are now (6/13/99) included from pelutils.h. They replace the
3102 second pair, which had been declared at this point, specifically
3103 for this file, EVEN THOUGH THE ARGUMENT LIST IS DIFFERENT.
3104 I do not understand. - AMM */
3105 /*
3106 void *mem_malloc(int);
3107 void mem_free(void *);
3108 */
3109 /*
3110 char *mem_malloc(size_t);
3111 void mem_free(char *);
3112 */
3113
3114
3115
3116 /*-------------------------------------------------------------
3117 vector access macroes (which ignore any rows except for first)
3118 -------------------------------------------------------------*/
3119 #define Vstore(V) ((V->store)) /* maximum #elts available */
3120 #define Vlength(V) ((V->topr)) /* actual #elts stored */
3121 #define Vref1(V,i) (&(((V->coords)[i-1]))) /*acces ith elt (starting at 1) */
3122 #define Vref0(V,i) (&(((V->coords)[i]))) /*acces ith elt (starting at 0) */
3123 #define Vref(V,i) Vref1(V,i)
3124
3125 /*------------------------------------------------------------
3126 matrix access macroes
3127 -------------------------------------------------------------*/
3128 #define Mstore(V) ((V->store)) /* maximum #elts available */
3129 #define MMrows(V) ((V->store/V->ncols)) /* maximum #rows */
3130 #define Mrows(V) ((V->topr)) /* number rows stored */
3131 #define Mcols(V) ((V->topc)) /* number cols stored */
3132 #define MNcols(V) ((V->ncols)) /* number cols stored */
3133
3134 #define Mref1(V,i,j) (&(((V->coords)[(i-1)*(V->ncols)+j-1])))
3135 #define Mref0(V,i,j) (&(((V->coords)[i*(V->ncols)+j])))
3136 #define Mref(V,i,j) Mref1((V),i,j)
3137
3138
3139
PMstore(Pmatrix M)3140 int PMstore(Pmatrix M)
3141 {
3142 return Mstore(M);
3143 }
PMMrows(Pmatrix M)3144 int PMMrows(Pmatrix M)
3145 {
3146 return MMrows(M);
3147 }
PMrows(Pmatrix M)3148 int PMrows(Pmatrix M)
3149 {
3150 return Mrows(M);
3151 }
PMcols(Pmatrix M)3152 int PMcols(Pmatrix M)
3153 {
3154 return Mcols(M);
3155 }
PMref1(Pmatrix M,int i,int j)3156 polynomial1 *PMref1(Pmatrix M, int i, int j)
3157 {
3158 return Mref1(M, i, j);
3159 }
3160
3161
3162
3163 /*
3164 ** Constructor/Destructors for Pmatrixes
3165 **
3166 ** Pmatrix Pmatrix_free(int r, int c);
3167 ** New Pmatrix cabable of holding r rows, and c collumbs.
3168 ** Pmatrix Pmatrix_new(Pmatrix V);
3169 */
3170
Pmatrix_new(int r,int c)3171 Pmatrix Pmatrix_new(int r, int c)
3172 {
3173 Pmatrix V;
3174 int i, j;
3175 V = (Pmatrix) mem_malloc(sizeof(struct Pmatrix_t));
3176 if (!V)
3177 bad_error("allocation failure in Pmatrix_new()");
3178 V->coords = (polynomial1 *) mem_malloc(r * c * sizeof(polynomial1));
3179 if (!V)
3180 bad_error("allocation failure 2 in Pmatrix_new()");
3181 Mstore(V) = r * c;
3182 Mrows(V) = r;
3183 Mcols(V) = c;
3184 MNcols(V) = c;
3185 for (i = 1; i <= r; i++) {
3186 for (j = 1; j <= c; j++)
3187 *Mref(V, i, j)=0;
3188 }
3189 return V;
3190 }
3191
Pmatrix_free(Pmatrix V)3192 void Pmatrix_free(Pmatrix V)
3193 {
3194 int i, j;
3195
3196 if (V != 0 && V->coords != 0) {
3197 for (i = 1; i <= Mrows(V); i++)
3198 for (j = 1; j <= Mcols(V); j++)
3199 freeP(*Mref(V, i, j));
3200 mem_free((char *) (V->coords));
3201 }
3202 if (V != 0)
3203 mem_free((char *) (V));
3204 }
3205
3206
3207 /*
3208 ** Pmatrix_resize(R,r,c)
3209 ** Reset R to hold an r,by c matrix.
3210 ** if R has enough storage to hold an rxc matrix resets
3211 ** row and columb entrees of r to r and c. otherwise
3212 ** frees R and reallocates an rxc matrix
3213 ** DOES NOT PRESERVE INDECIES OF EXISTING DATA
3214 */
Pmatrix_resize(Pmatrix R,int r,int c)3215 Pmatrix Pmatrix_resize(Pmatrix R, int r, int c)
3216 {
3217
3218 if (R == 0 || Mstore(R) < (r * c)) {
3219 if (R != 0) Pmatrix_free(R);
3220 R = Pmatrix_new(r, c);
3221 } else {
3222 Mrows(R) = r;
3223 Mcols(R) = c;
3224 MNcols(R) = c;
3225 }
3226 return R;
3227 }
3228
Pmatrix_submat(Pmatrix R,int r,int c)3229 Pmatrix Pmatrix_submat(Pmatrix R, int r, int c)
3230 {
3231
3232 if (R == 0 || c > Mcols(R) || r > Mrows(R) * MNcols(R)) {
3233 bad_error("bad subscripts or zero matrix in Pmatrix_submat()");
3234 } else {
3235 Mrows(R) = r;
3236 Mcols(R) = c;
3237 }
3238 return R;
3239 }
3240
3241 /*
3242 ** Pmatrix_print(M): print an Pmatrix
3243 ** if M is null print <<>> and return fail.
3244 ** otherwise print matrix and return true.
3245 */
Pmatrix_print(Pmatrix M)3246 Pmatrix Pmatrix_print(Pmatrix M)
3247 {
3248 int i, j;
3249
3250 if (M == 0) {
3251 printf("<<>>");
3252 return 0;
3253 }
3254 printf("<");
3255 for (i = 1; i <= Mrows(M); i++) {
3256 printf("<");
3257 for (j = 1; j <= Mcols(M); j++) {
3258 printP(*Mref(M, i, j));
3259 if (j < Mcols(M))
3260 printf(", ");
3261 }
3262 printf(">");
3263 if (i < Mrows(M))
3264 printf("\n");
3265 }
3266 printf(">");
3267 return M;
3268 }
3269
Pmatrix_copy(Pmatrix P)3270 Pmatrix Pmatrix_copy(Pmatrix P){
3271 Pmatrix C;
3272 int i,j;
3273 C=Pmatrix_new(Mrows(P),Mcols(P));
3274 for(i=1;i<=Mrows(P);i++){
3275 for(j=1;j<=Mcols(P);j++){
3276 *Mref(C,i,j)=copyP(*Mref(P,i,j));
3277 }
3278 }
3279 return C;
3280 }
3281
3282
3283 /* end Pmatrix.c */
3284
3285
3286 /**********************************************************************/
3287 /********* implementation code from the original Homotopies.c *********/
3288 /**********************************************************************/
3289
3290 /*------------------------------------------------------------------
3291 Homotopies.c created 9/15/1994 last modified 9/15/994
3292 Birk Huber (birk@math.cornell.edu
3293 ALL RIGHTS RESERVED
3294
3295 Store and evaluate Homotopies for Continuation module of Pelican.
3296 the functions rho and rhojac are defined to be usable by the
3297 hompack code (originaly in fortran translated with the help of f2c).
3298
3299 The representation of Homotopies and implementation is very similar
3300 to that used by hompack and described in Morgans book.
3301 -------------------------------------------------------------------*/
3302
3303 /*------------------------------------------------------------------
3304 globals representing the homotopy and macroes for accessing them
3305 The system of complex polynomial1s is represented in real form
3306 -----------------------------------------------------------------*/
3307
3308 /* parameters affecting the homotopy */
3309 int Hom_defd = 0; /* 0 no homotopy initialized , 1 else*/
3310 int Hom_use_proj = 1; /* 0 dont use proj trans, 1 else*/
3311 int Hom_num_vars = 0; /* number of complex vars in curr hom*/
3312
3313 /* private variables defining homotopy */
3314 static int NV,N,N1,M;
3315 static int *Starting_Monomial;
3316 static int *Number_of_Monomials;
3317 static int *Exponents;
3318 static int *Hdegree;
3319 static int *Edegree;
3320 static double *Coefitients;
3321 static int *Deformation;
3322 static double *Proj_Trans;
3323
3324 /* index in monomial list (starting at 0) of equation i*/
3325 #define monst(i,j) ((Starting_Monomial[(i)-1])+(j)-1)
3326
3327 /* number of monomials in the ith equation */
3328 #define Nmon(i) (Number_of_Monomials[(i)-1])
3329
3330 /* exponent of variable k in monomial j of equation i */
3331 #define Exp(i,j,k) Exponents[(monst(i,j))*NV+(k)-1]
3332
3333 /* get the exponent of the deformation parameter in the jth monomial
3334 of the ith equation */
3335 #define Def(i,j) Deformation[monst(i,j)]
3336
3337 /* get the real and imaginary part of jth monomial of ith equation */
3338 #define RCoef(i,j) Coefitients[2*(monst(i,j))]
3339 #define ICoef(i,j) Coefitients[2*(monst(i,j))+1]
3340
3341 /* real and imaginary parts of coordinates defining the
3342 projective transformation */
3343 #define RPtrans(j) Proj_Trans[2*(j)-2]
3344 #define IPtrans(j) Proj_Trans[2*(j)-1]
3345
3346 /* exponent of homogenizing variable in jthmonomial of ith equation*/
3347 #define Hdeg(i,j) Hdegree[monst(i,j)]
3348
3349 /* degree of equation i */
3350 #define Edeg(i) Edegree[i-1]
3351
3352 extern Pring Def_Ring;
psys_to_Pvec(psys sys)3353 Pvector psys_to_Pvec(psys sys){
3354 polynomial1 tmpm,tmpp;
3355 int j;
3356 Pvector PV;
3357 PV=Pvector_new(psys_d(sys));
3358 FORALL_POLY(sys,
3359 tmpp=0;
3360 FORALL_MONO(sys,
3361 tmpm=makeP(Def_Ring);
3362 *poly_coef(tmpm)=Complex(*psys_coef_real(sys),
3363 *psys_coef_imag(sys));
3364 *poly_def(tmpm)=*psys_def(sys);
3365 *poly_homog(tmpm)=*psys_homog(sys);
3366 for(j=1;j<=psys_d(sys);j++) *poly_exp(tmpm,j)=*psys_exp(sys,j);
3367 tmpp=addPPP(tmpp,tmpm,tmpp);
3368 freeP(tmpm);
3369 )
3370 *PVref(PV,psys_eqno(sys))=tmpp;
3371 )
3372 return PV;
3373 }
3374
3375 /*-------------------------------------------------------------------
3376 init_hom takes a Pvector and loads the above data structures
3377 to hold a representation of the system.
3378 -------------------------------------------------------------------*/
init_hom(psys PS)3379 int init_hom(psys PS){
3380 Pvector P;
3381 polynomial1 ptr;
3382 int i,j,k;
3383 double t;
3384 int seed=12;
3385
3386 P=psys_to_Pvec(PS);
3387 NV=poly_dim(*PMref(P,1,1)); /* number of complex variables*/
3388 N=2*NV; /* number of real coordinates */
3389 N1=N+1; /* number of realcoords incl lift coord*/
3390
3391
3392 /* check that P has n elts, and that Ring has NV vars */
3393 if ( NV != PMcols(P)) {
3394 printf("warning nonsquare homotopy in init_HPK\n");
3395 Hom_defd=0;
3396 }
3397
3398 /* Set number of vars for external use */
3399 Hom_num_vars=NV;
3400
3401 /* free all work space */
3402 Dfree(0); Ifree(0);
3403
3404 /* count monomials in equations to initialize monst and Nmon*/
3405 Starting_Monomial=Ires(NV+1);
3406 Number_of_Monomials=Ires(NV);
3407
3408 Starting_Monomial[0]=0;
3409 for(i=0;i<NV;i++){
3410 Number_of_Monomials[i]=0; ptr=*PMref(P,1,i+1);
3411 while(ptr!=0){Number_of_Monomials[i]++; ptr=poly_next(ptr);}
3412 Starting_Monomial[i+1]=
3413 Starting_Monomial[i]+Number_of_Monomials[i];
3414 }
3415 M=Starting_Monomial[NV];
3416
3417 Exponents=Ires(NV*M);
3418 Coefitients=Dres(M*2);
3419 Deformation=Ires(M);
3420 Hdegree=Ires(M);
3421 Edegree=Ires(NV);
3422 Proj_Trans=Dres(2*NV+2);
3423
3424 #if defined(HAVE_SRAND48)
3425 srand48(seed);
3426 #else
3427 srand(seed);
3428 #endif /* defined(HAVE_SRAND48) */
3429
3430 for(i=1;i<=NV;i++){
3431 j=1; ptr=*PMref(P,1,i); Edeg(i)=0;
3432 while(ptr!=0){
3433 Hdeg(i,j)=0;
3434 RCoef(i,j)=(*poly_coef(ptr)).r;
3435 ICoef(i,j)=(*poly_coef(ptr)).i;
3436 for(k=1;k<=NV;k++){
3437 Exp(i,j,k)=*poly_exp(ptr,k);
3438 Hdeg(i,j)+=Exp(i,j,k);
3439 }
3440 Edeg(i)=max(Edeg(i),Hdeg(i,j));
3441 Def(i,j)=*poly_def(ptr);
3442 j++; ptr=poly_next(ptr);
3443 }
3444 for(j=1;j<=Nmon(i);j++) Hdeg(i,j)=Edeg(i)-Hdeg(i,j);
3445 }
3446
3447 /*Define Projective transformation */
3448 for(j=1;j<=NV+1;j++){
3449 #if defined(HAVE_DRAND48)
3450 t=drand48()*2*PI;
3451 #else
3452 t=rand()*2*PI;
3453 #endif /* defined(HAVE_DRAND48) */
3454 RPtrans(j)=cos(t);
3455 IPtrans(j)=sin(t);
3456 }
3457 Hom_defd=1;
3458 Pvector_free(P);
3459 return Hom_defd;
3460 }
3461
print_homog(double * x,double * coord_r,double * coord_i)3462 void print_homog(double *x,double *coord_r,double *coord_i){
3463 int i;
3464 fcomplex PN;
3465 PN=Complex(RPtrans((NV+1)),IPtrans((NV+1)));
3466 for(i=1;i<=NV;i++)
3467 PN=Cadd(PN,Cmul(Complex(RPtrans(i),IPtrans(i)),
3468 Complex(x[2*i-2],x[2*i-1])));
3469 *coord_r=PN.r;
3470 *coord_i=PN.i;
3471 }
3472
print_proj_trans()3473 void print_proj_trans(){
3474 int i;
3475 #ifdef HOM_PRINT
3476 fprintf(Hom_LogFile,"T %d ", NV)
3477 #endif
3478 ;
3479 for(i=1;i<=NV+1;i++) {
3480 #ifdef HOM_PRINT
3481 fprintf(Hom_LogFile,"%10g %10g",RPtrans(i),IPtrans(i))
3482 #endif
3483 ;
3484 }
3485 #ifdef HOM_PRINT
3486 fprintf(Hom_LogFile,"\n")
3487 #endif
3488 ;
3489 }
3490
3491
3492 /*
3493 ** Htransform takes a vector in RP^(2*NV) and scales it so that
3494 ** it lies on the hyperplane defining the projective
3495 ** transformation. (Using the scaling factor
3496 ** L=pn1/(Z-(p1X1+...+pnXn)).)
3497 ** Huntransform takes a 2*NV+2 dvector whoose first 2*NV coords
3498 ** define a point in the affine chart defined by the
3499 ** projective transformation and fills in the last two
3500 ** coords with the value Z=pn1+p1X1+...+pnXn.
3501 **
3502 */
3503 #define X(i) (DVref(X,i))
Htransform(Dvector X)3504 void Htransform(Dvector X){
3505 int i;
3506 fcomplex C,L;
3507
3508 L=Complex(0.0,0.0);
3509 for(i=1;i<=NV;i++){
3510 L=Cadd(L,Cmul(Complex(RPtrans(i),IPtrans(i)),
3511 Complex(X(2*i+1),X(2*i+2))));
3512 }
3513 L=Cdiv(Complex(RPtrans(NV+1),IPtrans(NV+1)),
3514 Csub(Complex(X(1),X(2)),L));
3515 for(i=1;i<=NV+1;i++) {
3516 C=Cmul(L,Complex(X(2*i-1),X(2*i)));
3517 X(2*i-1)=C.r;
3518 X(2*i) =C.i;
3519 }
3520 }
3521
Huntransform(Dvector X)3522 void Huntransform(Dvector X){
3523 int i;
3524 fcomplex L;
3525 L=Complex(RPtrans(NV+1),IPtrans(NV+1));
3526 for(i=1;i<=NV;i++) {
3527 L=Cadd(L,Cmul(Complex(RPtrans(i),IPtrans(i)),
3528 Complex(X(2*i+1),X(2*i+2))));
3529 }
3530 X(1)=L.r;
3531 X(2)=L.i;
3532 }
3533 #undef X
3534
Hpath(double t)3535 fcomplex Hpath(double t){ return Complex(t,0.0);}
3536
DHpath(double t)3537 fcomplex DHpath(double t){return Complex(1.0,0.0);}
3538
rho_(double * a,double * lambda,double * x,double * v,double * par,int * ipar)3539 int rho_(double *a,
3540 double *lambda,
3541 double *x,
3542 double *v,
3543 double *par,
3544 int *ipar)
3545 {int i,j,h;
3546 fcomplex c,Hpath(double),PN;
3547 if (*lambda < 0.0) *lambda = 0.;
3548
3549 /* calculate projective coordinate */
3550 if (Hom_use_proj==1){
3551 PN=Complex(RPtrans(NV+1),IPtrans(NV+1));
3552 for(i=1;i<=NV;i++) PN=Cadd(PN,Cmul(Complex(RPtrans(i),IPtrans(i)),
3553 Complex(x[2*i-2],x[2*i-1])));
3554 }
3555 /* loop over equations */
3556 for(i=1;i<=NV;i++){
3557 v[2*i-2]=0.0;
3558 v[2*i-1]=0.0;
3559 for(j=1;j<=Nmon(i);j++){
3560 c=Cmul(Complex(RCoef(i,j),ICoef(i,j)),
3561 Cpow(Hpath(*lambda),Def(i,j)));
3562 for(h=1;h<=NV;h++)
3563 c=Cmul(c,Cpow(Complex(x[2*h-2],x[2*h-1]),Exp(i,j,h)));
3564 if (Hom_use_proj==1) c=Cmul(c,Cpow(PN,Hdeg(i,j)));
3565 v[2*i-2]+=c.r;
3566 v[2*i-1]+=c.i;
3567 }
3568 }
3569 return 0;
3570 }
3571
rhojac_(double * a,double * lambda,double * x,double * v,int * k,double * par,int * ipar)3572 int rhojac_(double *a,
3573 double *lambda,
3574 double *x,
3575 double *v,
3576 int *k,
3577 double *par,
3578 int *ipar)
3579 {
3580 int i,j,h,d;
3581 fcomplex c,Hpath(double), DHpath(double),PN;
3582 double t;
3583
3584 if (*lambda < 0.) *lambda = 0.;
3585 if(Hom_use_proj==1){
3586 PN=Complex(RPtrans((NV+1)),IPtrans((NV+1)));
3587 for(i=1;i<=NV;i++)
3588 PN=Cadd(PN,Cmul(Complex(RPtrans(i),IPtrans(i)),
3589 Complex(x[2*i-2],x[2*i-1])));
3590 }
3591 if (*k>1){
3592 d=(*k)/2; /* determine the complex var to diff by.*/
3593 for(i=1;i<=NV;i++){
3594 v[2*i-2]=0.0;
3595 v[2*i-1]=0.0;
3596 if (Hom_use_proj==1){
3597 for(j=1;j<=Nmon(i);j++){
3598 if (Hdeg(i,j)!=0){
3599 c=Cmul(Complex(RCoef(i,j),ICoef(i,j)),
3600 Cpow(Hpath(*lambda),Def(i,j)));
3601 for(h=1;h<=NV;h++)
3602 c=Cmul(c,Cpow(Complex(x[2*h-2],x[2*h-1]),Exp(i,j,h)));
3603 c=Cmul(c,Cpow(PN,Hdeg(i,j)-1));
3604 c=Cmul(c,Complex(t=Hdeg(i,j),0.0));
3605 c=Cmul(c,Complex(RPtrans(d),IPtrans(d)));
3606 v[2*i-2]+=c.r;
3607 v[2*i-1]+=c.i;
3608 }
3609 }
3610 }
3611
3612 for(j=1;j<=Nmon(i);j++){
3613 if (Exp(i,j,d)!=0){
3614 c=Cmul(Complex(RCoef(i,j),ICoef(i,j)),
3615 Cpow(Hpath(*lambda),Def(i,j)));
3616 for(h=1;h<d;h++)
3617 c=Cmul(c,Cpow(Complex(x[2*h-2],x[2*h-1]),Exp(i,j,h)));
3618 c=Cmul(c,Cpow(Complex(x[2*d-2],x[2*d-1]),Exp(i,j,d)-1));
3619 c=Cmul(c,Complex(t=Exp(i,j,d),0.0));
3620 for(h=d+1;h<=NV;h++)
3621 c=Cmul(c,Cpow(Complex(x[2*h-2],x[2*h-1]),Exp(i,j,h)));
3622 if(Hom_use_proj==1) c=Cmul(c,Cpow(PN,Hdeg(i,j)));
3623 v[2*i-2]+=c.r;
3624 v[2*i-1]+=c.i;
3625 }
3626 }
3627 if(((*k)%2)==1) {
3628 t=v[2*i-2];
3629 v[2*i-2]=-1.0*v[2*i-1];
3630 v[2*i-1]=t;
3631 }
3632 }
3633 }
3634 else {
3635 for(i=1;i<=NV;i++){
3636 v[2*i-2]=0.0;
3637 v[2*i-1]=0.0;
3638 for(j=1;j<=Nmon(i);j++){
3639 if(Def(i,j)!=0){
3640 c=Cmul(Complex(RCoef(i,j),ICoef(i,j)),
3641 Cpow(Hpath(*lambda),Def(i,j)-1));
3642 c=Cmul(c,Cmul(Complex(t=Def(i,j),0.0),DHpath(*lambda)));
3643 for(h=1;h<=NV;h++)
3644 c=Cmul(c,Cpow(Complex(x[2*h-2],x[2*h-1]),Exp(i,j,h)));
3645 if(Hom_use_proj==1) c=Cmul(c,Cpow(PN,Hdeg(i,j)));
3646 v[2*i-2]+=c.r;
3647 v[2*i-1]+=c.i;
3648 }
3649 }
3650 }
3651 }
3652 return 0;
3653 }
3654
3655 /* end, original Homotopies.c */
3656
3657
3658 /**************************************************************************/
3659 /********************* implementations from tangnf.c **********************/
3660 /**************************************************************************/
3661
3662 /* tangnf.f -- translated by f2c (version 19940305).
3663 You must link the resulting object file with the libraries:
3664 -lf2c -lm (in that order)
3665 */
3666
3667
3668 /* Table of constant values */
3669
tangnf_(double * rholen,double * y,double * yp,double * ypold,double * a,double * qr,double * alpha,double * tz,int * pivot,int * nfe,int * n,int * iflag,double * par,int * ipar)3670 /* Subroutine */ int tangnf_(double *rholen,
3671 double *y,
3672 double *yp,
3673 double *ypold,
3674 double *a,
3675 double *qr,
3676 double *alpha,
3677 double *tz,
3678 int *pivot,
3679 int *nfe,
3680 int *n,
3681 int *iflag,
3682 double *par,
3683 int *ipar)
3684 {
3685 /* System generated locals */
3686 int qr_dim1, qr_offset, i__1, i__2, i__3;
3687 double d__1;
3688
3689 /* Builtin functions */
3690 /* double sqrt(); CANT DECLARE BUILTINS UNDER C++ */
3691
3692 /* Local variables */
3693 /* extern */ /* Subroutine */ /* int fjac_(); */
3694 static double beta;
3695 static int jbar;
3696 /* extern double ddot_(); */
3697 static double qrkk;
3698 /* extern double dnrm2_(); */
3699 /* extern */ /* Subroutine */ /* int f_(); */
3700 static int i, j, k;
3701 static double sigma, lambda, alphak;
3702 /* extern */ /* Subroutine */ /* int rhojac_(); */
3703 static int kp1, np1, np2;
3704 static double ypnorm;
3705 /* extern */ /* Subroutine */ /* int rho_(); */
3706 static double sum;
3707
3708
3709 /* THIS SUBROUTINE BUILDS THE JACOBIAN MATRIX OF THE HOMOTOPY MAP, */
3710 /* COMPUTES A QR DECOMPOSITION OF THAT MATRIX, AND THEN CALCULATES THE */
3711 /* (UNIT) TANGENT VECTOR AND THE NEWTON STEP. */
3712
3713 /* ON INPUT: */
3714
3715 /* RHOLEN < 0 IF THE NORM OF THE HOMOTOPY MAP EVALUATED AT */
3716 /* (A, LAMBDA, X) IS TO BE COMPUTED. IF RHOLEN >= 0 THE NORM IS NOT
3717 */
3718 /* COMPUTED AND RHOLEN IS NOT CHANGED. */
3719
3720 /* Y(1:N+1) = CURRENT POINT (LAMBDA(S), X(S)). */
3721
3722 /* YPOLD(1:N+1) = UNIT TANGENT VECTOR AT PREVIOUS POINT ON THE ZERO */
3723 /* CURVE OF THE HOMOTOPY MAP. */
3724
3725 /* A(1:*) = PARAMETER VECTOR IN THE HOMOTOPY MAP. */
3726
3727 /* QR(1:N,1:N+2), ALPHA(1:N), TZ(1:N+1), PIVOT(1:N+1) ARE WORK ARRAYS */
3728 /* USED FOR THE QR FACTORIZATION. */
3729
3730 /* NFE = NUMBER OF JACOBIAN MATRIX EVALUATIONS = NUMBER OF HOMOTOPY */
3731 /* FUNCTION EVALUATIONS. */
3732
3733 /* N = DIMENSION OF X. */
3734
3735 /* IFLAG = -2, -1, OR 0, INDICATING THE PROBLEM TYPE. */
3736
3737 /* PAR(1:*) AND IPAR(1:*) ARE ARRAYS FOR (OPTIONAL) USER PARAMETERS, */
3738 /* WHICH ARE SIMPLY PASSED THROUGH TO THE USER WRITTEN SUBROUTINES */
3739 /* RHO, RHOJAC. */
3740
3741 /* ON OUTPUT: */
3742
3743 /* RHOLEN = ||RHO(A, LAMBDA(S), X(S)|| IF RHOLEN < 0 ON INPUT. */
3744 /* OTHERWISE RHOLEN IS UNCHANGED. */
3745
3746 /* Y, YPOLD, A, N ARE UNCHANGED. */
3747
3748 /* YP(1:N+1) = DY/DS = UNIT TANGENT VECTOR TO INTEGRAL CURVE OF */
3749 /* D(HOMOTOPY MAP)/DS = 0 AT Y(S) = (LAMBDA(S), X(S)) . */
3750
3751 /* TZ = THE NEWTON STEP = -(PSEUDO INVERSE OF (D RHO(A,Y(S))/D LAMBDA ,
3752 */
3753 /* D RHO(A,Y(S))/DX)) * RHO(A,Y(S)) . */
3754
3755 /* NFE HAS BEEN INCRMENTED BY 1. */
3756
3757 /* IFLAG IS UNCHANGED, UNLESS THE QR FACTORIZATION DETECTS A RANK < N, */
3758 /* IN WHICH CASE THE TANGENT AND NEWTON STEP VECTORS ARE NOT COMPUTED
3759 */
3760 /* AND TANGNF RETURNS WITH IFLAG = 4 . */
3761
3762
3763 /* CALLS DDOT , DNRM2 , F (OR RHO ), FJAC (OR RHOJAC ). */
3764
3765
3766 /* ***** ARRAY DECLARATIONS. ***** */
3767
3768
3769 /* ARRAYS FOR COMPUTING THE JACOBIAN MATRIX AND ITS KERNEL. */
3770
3771 /* ***** END OF DIMENSIONAL INFORMATION. ***** */
3772
3773
3774 /* Parameter adjustments */
3775 --pivot;
3776 --tz;
3777 --alpha;
3778 qr_dim1 = *n;
3779 qr_offset = qr_dim1 + 1;
3780 qr -= qr_offset;
3781 --a;
3782 --ypold;
3783 --yp;
3784 --y;
3785 --par;
3786 --ipar;
3787
3788 /* Function Body */
3789 lambda = y[1];
3790 np1 = *n + 1;
3791 np2 = *n + 2;
3792 ++(*nfe);
3793 /* NFE CONTAINS THE NUMBER OF JACOBIAN EVALUATIONS. */
3794 /* * * * * * * * * * * * * * * * * * */
3795
3796 /* COMPUTE THE JACOBIAN MATRIX, STORE IT AND HOMOTOPY MAP IN QR. */
3797
3798 if (*iflag == -2) {
3799
3800 /* QR = ( D RHO(A,LAMBDA,X)/D LAMBDA , D RHO(A,LAMBDA,X)/DX , */
3801 /* RHO(A,LAMBDA,X) ) .
3802 */
3803
3804 i__1 = np1;
3805 for (k = 1; k <= i__1; ++k) {
3806 rhojac_(&a[1], &lambda, &y[2], &qr[k * qr_dim1 + 1], &k, &par[1],
3807 &ipar[1]);
3808 /* L30: */
3809 }
3810 rho_(&a[1], &lambda, &y[2], &qr[np2 * qr_dim1 + 1], &par[1], &ipar[1])
3811 ;
3812 } else {
3813 f_(&y[2], &tz[1]);
3814 if (*iflag == 0) {
3815
3816 /* QR = ( A - F(X), I - LAMBDA*DF(X) , */
3817 /* X - A + LAMBDA*(A - F(X)) ) .
3818 */
3819
3820 i__1 = *n;
3821 for (j = 1; j <= i__1; ++j) {
3822 sigma = a[j];
3823 beta = sigma - tz[j];
3824 qr[j + qr_dim1] = beta;
3825 /* L100: */
3826 qr[j + np2 * qr_dim1] = y[j + 1] - sigma + lambda * beta;
3827 }
3828 i__1 = *n;
3829 for (k = 1; k <= i__1; ++k) {
3830 fjac_(&y[2], &tz[1], (integer *)&k);
3831 kp1 = k + 1;
3832 i__2 = *n;
3833 for (j = 1; j <= i__2; ++j) {
3834 /* L110: */
3835 qr[j + kp1 * qr_dim1] = -lambda * tz[j];
3836 }
3837 /* L120: */
3838 qr[k + kp1 * qr_dim1] += (float)1.;
3839 }
3840 } else {
3841
3842 /* QR = ( F(X) - X + A, LAMBDA*DF(X) + (1 - LAMBDA)*I , */
3843 /* X - A + LAMBDA*(F(X) - X + A)
3844 ) . */
3845
3846 /* L140: */
3847 i__1 = *n;
3848 for (j = 1; j <= i__1; ++j) {
3849 sigma = y[j + 1] - a[j];
3850 beta = tz[j] - sigma;
3851 qr[j + qr_dim1] = beta;
3852 /* L150: */
3853 qr[j + np2 * qr_dim1] = sigma + lambda * beta;
3854 }
3855 i__1 = *n;
3856 for (k = 1; k <= i__1; ++k) {
3857 fjac_(&y[2], &tz[1], (integer *)&k);
3858 kp1 = k + 1;
3859 i__2 = *n;
3860 for (j = 1; j <= i__2; ++j) {
3861 /* L160: */
3862 qr[j + kp1 * qr_dim1] = lambda * tz[j];
3863 }
3864 /* L170: */
3865 qr[k + kp1 * qr_dim1] = (float)1. - lambda + qr[k + kp1 *
3866 qr_dim1];
3867 }
3868 }
3869 }
3870
3871 /* * * * * * * * * * * * * * * * * * */
3872 /* COMPUTE THE NORM OF THE HOMOTOPY MAP IF IT WAS REQUESTED. */
3873 if (*rholen < (float)0.) {
3874 *rholen = dnrm2_((integer *)n, &qr[np2 * qr_dim1 + 1], &c__1);
3875 }
3876
3877 /* REDUCE THE JACOBIAN MATRIX TO UPPER TRIANGULAR FORM. */
3878
3879 /* THE FOLLOWING CODE IS A MODIFICATION OF THE ALGOL PROCEDURE */
3880 /* DECOMPOSE IN P. BUSINGER AND G. H. GOLUB, LINEAR LEAST */
3881 /* SQUARES SOLUTIONS BY HOUSEHOLDER TRANSFORMATIONS, */
3882 /* NUMER. MATH. 7 (1965) 269-276. */
3883
3884 i__1 = np1;
3885 for (j = 1; j <= i__1; ++j) {
3886 yp[j] = ddot_((integer *)n, &qr[j * qr_dim1 + 1],
3887 &c__1, &qr[j * qr_dim1 + 1], &c__1);
3888 /* L220: */
3889 pivot[j] = j;
3890 }
3891 i__1 = *n;
3892 for (k = 1; k <= i__1; ++k) {
3893 sigma = yp[k];
3894 jbar = k;
3895 kp1 = k + 1;
3896 i__2 = np1;
3897 for (j = kp1; j <= i__2; ++j) {
3898 if (sigma >= yp[j]) {
3899 goto L240;
3900 }
3901 sigma = yp[j];
3902 jbar = j;
3903 L240:
3904 ;
3905 }
3906 if (jbar == k) {
3907 goto L260;
3908 }
3909 i = pivot[k];
3910 pivot[k] = pivot[jbar];
3911 pivot[jbar] = i;
3912 yp[jbar] = yp[k];
3913 yp[k] = sigma;
3914 i__2 = *n;
3915 for (i = 1; i <= i__2; ++i) {
3916 sigma = qr[i + k * qr_dim1];
3917 qr[i + k * qr_dim1] = qr[i + jbar * qr_dim1];
3918 qr[i + jbar * qr_dim1] = sigma;
3919 /* L250: */
3920 }
3921 /* END OF COLUMN INTERCHANGE. */
3922 L260:
3923 i__2 = *n - k + 1;
3924 sigma = ddot_((integer *)&i__2, &qr[k + k * qr_dim1],
3925 &c__1, &qr[k + k * qr_dim1], &c__1);
3926 if (sigma == (float)0.) {
3927 *iflag = 4;
3928 return 0;
3929 }
3930 /* L270: */
3931 if (k == *n) {
3932 goto L300;
3933 }
3934 qrkk = qr[k + k * qr_dim1];
3935 alphak = -sqrt(sigma);
3936 if (qrkk < (float)0.) {
3937 alphak = -alphak;
3938 }
3939 alpha[k] = alphak;
3940 beta = (float)1. / (sigma - qrkk * alphak);
3941 qr[k + k * qr_dim1] = qrkk - alphak;
3942 i__2 = np2;
3943 for (j = kp1; j <= i__2; ++j) {
3944 i__3 = *n - k + 1;
3945 sigma = beta * ddot_((integer *)&i__3, &qr[k + k * qr_dim1],
3946 &c__1, &qr[k + j * qr_dim1], &c__1);
3947 i__3 = *n;
3948 for (i = k; i <= i__3; ++i) {
3949 qr[i + j * qr_dim1] -= qr[i + k * qr_dim1] * sigma;
3950 /* L280: */
3951 }
3952 if (j < np2) {
3953 /* Computing 2nd power */
3954 d__1 = qr[k + j * qr_dim1];
3955 yp[j] -= d__1 * d__1;
3956 }
3957 /* L290: */
3958 }
3959 L300:
3960 ;
3961 }
3962 alpha[*n] = qr[*n + *n * qr_dim1];
3963
3964
3965 /* COMPUTE KERNEL OF JACOBIAN, WHICH SPECIFIES YP=DY/DS. */
3966 tz[np1] = (float)1.;
3967 for (i = *n; i >= 1; --i) {
3968 sum = (float)0.;
3969 i__1 = np1;
3970 for (j = i + 1; j <= i__1; ++j) {
3971 /* L330: */
3972 sum += qr[i + j * qr_dim1] * tz[j];
3973 }
3974 /* L340: */
3975 tz[i] = -sum / alpha[i];
3976 }
3977 ypnorm = dnrm2_((integer *)&np1, &tz[1], &c__1);
3978 i__1 = np1;
3979 for (k = 1; k <= i__1; ++k) {
3980 /* L360: */
3981 yp[pivot[k]] = tz[k] / ypnorm;
3982 }
3983 if (ddot_((integer *)&np1, &yp[1], &c__1, &ypold[1], &c__1) >= (float)0.) {
3984 goto L380;
3985 }
3986 i__1 = np1;
3987 for (i = 1; i <= i__1; ++i) {
3988 /* L370: */
3989 yp[i] = -yp[i];
3990 }
3991 /* YP IS THE UNIT TANGENT VECTOR IN THE CORRECT DIRECTION. */
3992
3993 /* COMPUTE THE MINIMUM NORM SOLUTION OF [D RHO(Y(S))] V = -RHO(Y(S)). */
3994 /* V IS GIVEN BY P - (P,Q)Q , WHERE P IS ANY SOLUTION OF */
3995 /* [D RHO] V = -RHO AND Q IS A UNIT VECTOR IN THE KERNEL OF [D RHO]. */
3996
3997 L380:
3998 for (i = *n; i >= 1; --i) {
3999 sum = qr[i + np1 * qr_dim1] + qr[i + np2 * qr_dim1];
4000 i__1 = *n;
4001 for (j = i + 1; j <= i__1; ++j) {
4002 /* L430: */
4003 sum += qr[i + j * qr_dim1] * alpha[j];
4004 }
4005 /* L440: */
4006 alpha[i] = -sum / alpha[i];
4007 }
4008 i__1 = *n;
4009 for (k = 1; k <= i__1; ++k) {
4010 /* L450: */
4011 tz[pivot[k]] = alpha[k];
4012 }
4013 tz[pivot[np1]] = (float)1.;
4014 /* TZ NOW CONTAINS A PARTICULAR SOLUTION P, AND YP CONTAINS A VECTOR Q */
4015 /* IN THE KERNEL(THE TANGENT). */
4016 sigma = ddot_((integer *)&np1, &tz[1], &c__1, &yp[1], &c__1);
4017 i__1 = np1;
4018 for (j = 1; j <= i__1; ++j) {
4019 tz[j] -= sigma * yp[j];
4020 /* L470: */
4021 }
4022 /* TZ IS THE NEWTON STEP FROM THE CURRENT POINT Y(S) = (LAMBDA(S), X(S)).
4023 */
4024 return 0;
4025 } /* tangnf_ */
4026
4027 /* end tangnf.c */
4028