1 /*-------------------------------------------------------------------
2 Call_Hompack.c created 9/15/1994 last modified 9/15/1994
3 Birk Huber (birk@math.cornell.edu
4 ALL RIGHTS RESERVED
5
6 This File represents the interface between Pelican and Hompacks
7 FIXPNF path tracker. The HOMPACK routines it calls have actually
8 been translated from fortran into c with f2c and then modified a
9 little so as not to require the f2c library.
10
11 The two routines the user needs to be aware of are init_HPK
12 which takes a pelican Pvector, converts it to "tableau" format,
13 and initialies all the nescessary global variables to represent
14 the homotopy. Call_HPK_Hom wich takes a double vector and uses
15 it as a starting point for path continuation.
16 --------------------------------------------------------------------*/
17
18 #include "pelclhpk.h"
19
20 #define X(i) (DVref(X,i))
21
HPK_cont(Dvector X,int tweak)22 int HPK_cont(Dvector X, int tweak)
23 {
24 int i, ist, dst,N,N1;
25 static int iflag,trace,nfe,*pivot,*ipar;
26 static double *yold,*a,arcae;
27 static double *qr, arclen, *wp, *yp, *tz,*par,*z0,*z1;
28 static double ansre, *ypold, *sspar,*alpha,ansae,*w,*y,arcre;
29 /* extern int fixpnf_(); IN Homotopies.h */
30
31 if (Hom_defd!=1) return 0;
32
33 /* save start of free storage space */
34 ist=Itop(); dst=Dtop();
35
36 /* init numbers of real eqs with and without hom param*/
37 N=2*Hom_num_vars;
38 N1=N+1;
39
40 /* get space for local arrays from store */
41 ipar=Ires(1); pivot=Ires(N1); yold=Dres(N1); a=Dres(N1);
42 alpha=Dres(N1); w=Dres(N1); y=Dres(N1); ypold=Dres(N1);
43 sspar=Dres(8); z0=Dres(N1);
44 z1=Dres(N1); qr=Dres(N*(N1+1)); wp=Dres(N1); yp=Dres(N1);
45 tz=Dres(N1); par=Dres(1);
46
47 /* initialize parameters */
48 iflag = -2; /* should not be changed switch to tell hompack to do path tracking*/
49 arcre = -1.; arcae = -1.;
50 /* errors allowed during normal flow iteration will be
51 automatically reset to appropriat values later*/
52 ansre = Hom_tol;
53 ansae = Hom_tol;
54 trace = 1; /* 1 keep a log , 0 dont */
55 nfe = 10; /* I am not sure what this one is for */
56 for(i=0;i<8;i++) sspar[i] = -1.; /* sspar holds a number of flags used to determine
57 optimal step size, set to -1 will cause hompack
58 to choose them by its own heuristics */
59 Htransform(X);
60 /* print Starting point to Log File*/
61 print_proj_trans();
62 for(i=1;i<=N+3;i++)
63 #ifdef CHP_PRINT
64 fprintf(Hom_LogFile,"S %g", X(i));
65 fprintf(Hom_LogFile," 0 \n")
66 #endif
67 ;
68 y[0]=X(N+3);/* y0 holds the initial deformation parameter */
69 for(i=3;i<=N+2;i++){ /* y1 and on hold the starting coordinates */
70 y[i-2]=X(i);
71 }
72 fixpnf_(&N, y, &iflag, &arcre, &arcae, &ansre, &ansae, &trace,
73 a, &nfe, &arclen, yp, yold, ypold, qr, alpha, tz,
74 pivot, w, wp, z0, z1, sspar, par, ipar,
75 tweak); /* tweak is used to refine step size */
76 /*
77 #ifdef CHP_PRINT
78 fprintf(Hom_OutFile,"Done arclen=%f, LAMBDA=%f, return flag %d\n", arclen,y[0],iflag)
79 #endif
80
81 */
82 ;
83 for(i=1;i<=N;i++) X(i+2)=y[i];
84 X(N+3)=y[0];
85 Huntransform(X);
86 /* print ending point to log file */
87 #ifdef CHP_PRINT
88 fprintf(Hom_LogFile,"E")
89 #endif
90 ;
91 for(i=1;i<=N+3;i++)
92 #ifdef CHP_PRINT
93 fprintf(Hom_LogFile," %g", X(i))
94 #endif
95 ;
96 #ifdef CHP_PRINT
97 fprintf(Hom_LogFile," 4 1 %d %f %f",iflag,arclen,y[0]);
98 fprintf(Hom_LogFile,"\n")
99 #endif
100 ;
101
102 /*free space*/
103 Ifree(ist), Dfree(dst);
104 return 0;
105 }
106
107 #undef X
108