/* Homotopies.c */ /* This is the implementation code from the various files that used to reside in the Cont subdirectory of Pelican0.80/source. */ #include "pelhomot.h" /**********************************************************************/ /*************** implementation code from Hom_params.c ****************/ /**********************************************************************/ double Hom_tol=.000001; /* tolerence for path tracking */ void print_Hom_params(FILE *outfile){ #ifdef HOM_PRINT fprintf(outfile," Hominuation Parameters:\n"); fprintf(outfile," Hom_tol=%g\n",Hom_tol) #endif ; } /* end Hom_params.c */ /************************************************************************/ /***************************** code from fixpnf.c ***********************/ /************************************************************************/ #define TRUE_ (1) #define FALSE_ (0) double get_abs_homog(); double coord_r, coord_i; /* Table of constant values */ static integer c__4 = 4; static integer c__1 = 1; /* Subroutine */ int 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) /* added to adjust step size */ { /* System generated locals */ int qr_dim1, qr_offset, i__1, i__2; double d__1 /* ,abx UNUSED */; /* Builtin functions */ /* double sqrt(double); CANT DECLARE BUILTINS UNDER C++ */ integer s_wsfe(), do_fio(), e_wsfe(); /* in Hom_params.c these are int's */ /* Local variables */ static int nfec; static double hold; static int iter; extern double dnrm2_(integer *n, doublereal *dx, integer *incx); static double h, s; static long int crash; static int limit; extern double d1mach_(integer *); static long int start; static int nc, iflagc, jw; static double abserr, relerr; extern /* Subroutine */ int 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); static double curtol; extern /* Subroutine */ int 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); static int np1; static long int polsys; /* SUBROUTINE FIXPNF FINDS A FIXED POINT OR ZERO OF THE */ /* N-DIMENSIONAL VECTOR FUNCTION F(X), OR TRACKS A ZERO CURVE */ /* OF A GENERAL HOMOTOPY MAP RHO(A,LAMBDA,X). FOR THE FIXED */ /* POINT PROBLEM F(X) IS ASSUMED TO BE A C2 MAP OF SOME BALL */ /* INTO ITSELF. THE EQUATION X = F(X) IS SOLVED BY */ /* FOLLOWING THE ZERO CURVE OF THE HOMOTOPY MAP */ /* LAMBDA*(X - F(X)) + (1 - LAMBDA)*(X - A) , */ /* STARTING FROM LAMBDA = 0, X = A. THE CURVE IS PARAMETERIZED */ /* BY ARC LENGTH S, AND IS FOLLOWED BY SOLVING THE ORDINARY */ /* DIFFERENTIAL EQUATION D(HOMOTOPY MAP)/DS = 0 FOR */ /* Y(S) = (LAMBDA(S), X(S)) USING A HERMITE CUBIC PREDICTOR AND A */ /* CORRECTOR WHICH RETURNS TO THE ZERO CURVE ALONG THE FLOW NORMAL */ /* TO THE DAVIDENKO FLOW (WHICH CONSISTS OF THE INTEGRAL CURVES OF */ /* D(HOMOTOPY MAP)/DS ). */ /* FOR THE ZERO FINDING PROBLEM F(X) IS ASSUMED TO BE A C2 MAP */ /* SUCH THAT FOR SOME R > 0, X*F(X) >= 0 WHENEVER NORM(X) = R. */ /* THE EQUATION F(X) = 0 IS SOLVED BY FOLLOWING THE ZERO CURVE */ /* OF THE HOMOTOPY MAP */ /* LAMBDA*F(X) + (1 - LAMBDA)*(X - A) */ /* EMANATING FROM LAMBDA = 0, X = A. */ /* A MUST BE AN INTERIOR POINT OF THE ABOVE MENTIONED BALLS. */ /* FOR THE CURVE TRACKING PROBLEM RHO(A,LAMBDA,X) IS ASSUMED TO */ /* BE A C2 MAP FROM E**M X [0,1) X E**N INTO E**N, WHICH FOR */ /* ALMOST ALL PARAMETER VECTORS A IN SOME NONEMPTY OPEN SUBSET */ /* OF E**M SATISFIES */ /* RANK [D RHO(A,LAMBDA,X)/D LAMBDA , D RHO(A,LAMBDA,X)/DX] = N */ /* FOR ALL POINTS (LAMBDA,X) SUCH THAT RHO(A,LAMBDA,X)=0. IT IS */ /* FURTHER ASSUMED THAT */ /* RANK [ D RHO(A,0,X0)/DX ] = N . */ /* WITH A FIXED, THE ZERO CURVE OF RHO(A,LAMBDA,X) EMANATING */ /* FROM LAMBDA = 0, X = X0 IS TRACKED UNTIL LAMBDA = 1 BY */ /* SOLVING THE ORDINARY DIFFERENTIAL EQUATION */ /* D RHO(A,LAMBDA(S),X(S))/DS = 0 FOR Y(S) = (LAMBDA(S), X(S)), */ /* WHERE S IS ARC LENGTH ALONG THE ZERO CURVE. ALSO THE HOMOTOPY */ /* MAP RHO(A,LAMBDA,X) IS ASSUMED TO BE CONSTRUCTED SUCH THAT */ /* D LAMBDA(0)/DS > 0 . */ /* FOR THE FIXED POINT AND ZERO FINDING PROBLEMS, THE USER MUST SUPPLY */ /* A SUBROUTINE F(X,V) WHICH EVALUATES F(X) AT X AND RETURNS THE */ /* VECTOR F(X) IN V, AND A SUBROUTINE FJAC(X,V,K) WHICH RETURNS IN V */ /* THE KTH COLUMN OF THE JACOBIAN MATRIX OF F(X) EVALUATED AT X. FOR */ /* THE CURVE TRACKING PROBLEM, THE USER MUST SUPPLY A SUBROUTINE */ /* RHO(A,LAMBDA,X,V,PAR,IPAR) WHICH EVALUATES THE HOMOTOPY MAP RHO AT */ /* (A,LAMBDA,X) AND RETURNS THE VECTOR RHO(A,LAMBDA,X) IN V, AND A */ /* SUBROUTINE RHOJAC(A,LAMBDA,X,V,K,PAR,IPAR) WHICH RETURNS IN V THE KTH */ /* COLUMN OF THE N X (N+1) JACOBIAN MATRIX [D RHO/D LAMBDA, D RHO/DX] */ /* EVALUATED AT (A,LAMBDA,X). FIXPNF DIRECTLY OR INDIRECTLY USES */ /* THE SUBROUTINES STEPNF , TANGNF , ROOTNF , ROOT , F (OR RHO ), */ /* FJAC (OR RHOJAC ), D1MACH , AND THE BLAS FUNCTIONS DDOT AND */ /* DNRM2 . ONLY D1MACH CONTAINS MACHINE DEPENDENT CONSTANTS. */ /* NO OTHER MODIFICATIONS BY THE USER ARE REQUIRED. */ /* ON INPUT: */ /* N IS THE DIMENSION OF X, F(X), AND RHO(A,LAMBDA,X). */ /* Y IS AN ARRRAY OF LENGTH N + 1. (Y(2),...,Y(N+1)) = A IS THE */ /* STARTING POINT FOR THE ZERO CURVE FOR THE FIXED POINT AND */ /* ZERO FINDING PROBLEMS. (Y(2),...,Y(N+1)) = X0 FOR THE CURVE */ /* TRACKING PROBLEM. */ /* IFLAG CAN BE -2, -1, 0, 2, OR 3. IFLAG SHOULD BE 0 ON THE */ /* FIRST CALL TO FIXPNF FOR THE PROBLEM X=F(X), -1 FOR THE */ /* PROBLEM F(X)=0, AND -2 FOR THE PROBLEM RHO(A,LAMBDA,X)=0. */ /* IN CERTAIN SITUATIONS IFLAG IS SET TO 2 OR 3 BY FIXPNF, */ /* AND FIXPNF CAN BE CALLED AGAIN WITHOUT CHANGING IFLAG. */ /* ARCRE , ARCAE ARE THE RELATIVE AND ABSOLUTE ERRORS, RESPECTIVELY, */ /* ALLOWED THE NORMAL FLOW ITERATION ALONG THE ZERO CURVE. IF */ /* ARC?E .LE. 0.0 ON INPUT IT IS RESET TO .5*SQRT(ANS?E) . */ /* NORMALLY ARC?E SHOULD BE CONSIDERABLY LARGER THAN ANS?E . */ /* ANSRE , ANSAE ARE THE RELATIVE AND ABSOLUTE ERROR VALUES USED FOR */ /* THE ANSWER AT LAMBDA = 1. THE ACCEPTED ANSWER Y = (LAMBDA, X) */ /* SATISFIES */ /* |Y(1) - 1| .LE. ANSRE + ANSAE .AND. */ /* ||Z|| .LE. ANSRE*||X|| + ANSAE WHERE */ /* (.,Z) IS THE NEWTON STEP TO Y. */ /* TRACE IS AN INTEGER SPECIFYING THE LOGICAL I/O UNIT FOR */ /* INTERMEDIATE OUTPUT. IF TRACE .GT. 0 THE POINTS COMPUTED ON */ /* THE ZERO CURVE ARE WRITTEN TO I/O UNIT TRACE . */ /* A(1:*) CONTAINS THE PARAMETER VECTOR A . FOR THE FIXED POINT */ /* AND ZERO FINDING PROBLEMS, A NEED NOT BE INITIALIZED BY THE */ /* USER, AND IS ASSUMED TO HAVE LENGTH N. FOR THE CURVE */ /* TRACKING PROBLEM, A MUST BE INITIALIZED BY THE USER. */ /* YP(1:N+1) IS A WORK ARRAY CONTAINING THE TANGENT VECTOR TO */ /* THE ZERO CURVE AT THE CURRENT POINT Y . */ /* YOLD(1:N+1) IS A WORK ARRAY CONTAINING THE PREVIOUS POINT FOUND */ /* ON THE ZERO CURVE. */ /* YPOLD(1:N+1) IS A WORK ARRAY CONTAINING THE TANGENT VECTOR TO */ /* THE ZERO CURVE AT YOLD . */ /* QR(1:N,1:N+2), ALPHA(1:N), TZ(1:N+1), PIVOT(1:N+1) , W(1:N+1) , */ /* WP(1:N+1) , Z0(1:N+1) , Z1(1:N+1) ARE ALL WORK ARRAYS USED BY */ /* STEPNF TO CALCULATE THE TANGENT VECTORS AND NEWTON STEPS. */ /* SSPAR(1:8) = (LIDEAL, RIDEAL, DIDEAL, HMIN, HMAX, BMIN, BMAX, P) IS */ /* A VECTOR OF PARAMETERS USED FOR THE OPTIMAL STEP SIZE ESTIMATION. */ /* IF SSPAR(J) .LE. 0.0 ON INPUT, IT IS RESET TO A DEFAULT VALUE */ /* BY FIXPNF . OTHERWISE THE INPUT VALUE OF SSPAR(J) IS USED. */ /* SEE THE COMMENTS BELOW AND IN STEPNF FOR MORE INFORMATION ABOUT */ /* THESE CONSTANTS. */ /* PAR(1:*) AND IPAR(1:*) ARE ARRAYS FOR (OPTIONAL) USER PARAMETERS, */ /* WHICH ARE SIMPLY PASSED THROUGH TO THE USER WRITTEN SUBROUTINES */ /* RHO, RHOJAC. */ /* ON OUTPUT: */ /* N , TRACE , A ARE UNCHANGED. */ /* Y(1) = LAMBDA, (Y(2),...,Y(N+1)) = X, AND Y IS AN APPROXIMATE */ /* ZERO OF THE HOMOTOPY MAP. NORMALLY LAMBDA = 1 AND X IS A */ /* FIXED POINT(ZERO) OF F(X). IN ABNORMAL SITUATIONS LAMBDA */ /* MAY ONLY BE NEAR 1 AND X IS NEAR A FIXED POINT(ZERO). */ /* IFLAG = */ /* -2 CAUSES FIXPNF TO INITIALIZE EVERYTHING FOR THE PROBLEM */ /* RHO(A,LAMBDA,X) = 0 (USE ON FIRST CALL). */ /* -1 CAUSES FIXPNF TO INITIALIZE EVERYTHING FOR THE PROBLEM */ /* F(X) = 0 (USE ON FIRST CALL). */ /* 0 CAUSES FIXPNF TO INITIALIZE EVERYTHING FOR THE PROBLEM */ /* X = F(X) (USE ON FIRST CALL). */ /* 1 NORMAL RETURN. */ /* 2 SPECIFIED ERROR TOLERANCE CANNOT BE MET. SOME OR ALL OF */ /* ARCRE , ARCAE , ANSRE , ANSAE HAVE BEEN INCREASED TO */ /* SUITABLE VALUES. TO CONTINUE, JUST CALL FIXPNF AGAIN */ /* WITHOUT CHANGING ANY PARAMETERS. */ /* 3 STEPNF HAS BEEN CALLED 1000 TIMES. TO CONTINUE, CALL */ /* FIXPNF AGAIN WITHOUT CHANGING ANY PARAMETERS. */ /* 4 JACOBIAN MATRIX DOES NOT HAVE FULL RANK. THE ALGORITHM */ /* HAS FAILED (THE ZERO CURVE OF THE HOMOTOPY MAP CANNOT BE */ /* FOLLOWED ANY FURTHER). */ /* 5 THE TRACKING ALGORITHM HAS LOST THE ZERO CURVE OF THE */ /* HOMOTOPY MAP AND IS NOT MAKING PROGRESS. THE ERROR TOLERANCES */ /* ARC?E AND ANS?E WERE TOO LENIENT. THE PROBLEM SHOULD BE */ /* RESTARTED BY CALLING FIXPNF WITH SMALLER ERROR TOLERANCES */ /* AND IFLAG = 0 (-1, -2). */ /* 6 THE NORMAL FLOW NEWTON ITERATION IN STEPNF OR ROOTNF */ /* FAILED TO CONVERGE. THE ERROR TOLERANCES ANS?E MAY BE TOO */ /* STRINGENT. */ /* 7 ILLEGAL INPUT PARAMETERS, A FATAL ERROR. */ /* ARCRE , ARCAE , ANSRE , ANSAE ARE UNCHANGED AFTER A NORMAL RETURN */ /* (IFLAG = 1). THEY ARE INCREASED TO APPROPRIATE VALUES ON THE */ /* RETURN IFLAG = 2 . */ /* NFE IS THE NUMBER OF FUNCTION EVALUATIONS (= NUMBER OF */ /* JACOBIAN EVALUATIONS). */ /* ARCLEN IS THE LENGTH OF THE PATH FOLLOWED. */ /* ***** ARRAY DECLARATIONS. ***** */ /* ***** END OF DIMENSIONAL INFORMATION. ***** */ /* LIMITD IS AN UPPER BOUND ON THE NUMBER OF STEPS. IT MAY BE */ /* CHANGED BY CHANGING THE FOLLOWING PARAMETER STATEMENT: */ /* SWITCH FROM THE TOLERANCE ARC?E TO THE (FINER) TOLERANCE ANS?E IF */ /* THE CURVATURE OF ANY COMPONENT OF Y EXCEEDS CURSW. */ /* : : : : : : : : : : : : : : : : : : : : : : : : */ /* SET LOGICAL SWITCH TO REFLECT ENTRY POINT. */ /* Parameter adjustments */ --y; --a; --yp; --yold; --ypold; qr_dim1 = *n; qr_offset = qr_dim1 + 1; qr -= qr_offset; --alpha; --tz; --pivot; --w; --wp; --z0; --z1; --sspar; --par; --ipar; /* Function Body */ switch(n__) { case 1: goto L_polynf; } polsys = FALSE_; goto L11; L_polynf: polsys = TRUE_; L11: if (*n <= 0 || *ansre <= (float)0. || *ansae < (float)0.) { *iflag = 7; } if (*iflag >= -2 && *iflag <= 0) { goto L20; } if (*iflag == 2) { goto L120; } if (*iflag == 3) { goto L90; } /* ONLY VALID INPUT FOR IFLAG IS -2, -1, 0, 2, 3. */ *iflag = 7; return 0; /* ***** INITIALIZATION BLOCK. ***** */ L20: *arclen = (float)0.; if (*arcre <= (float)0.) { *arcre = sqrt(*ansre) * (float).5; } if (*arcae <= (float)0.) { *arcae = sqrt(*ansae) * (float).5; } nc = *n; nfec = 0; iflagc = *iflag; np1 = *n + 1; /* SET INITIAL CONDITIONS FOR FIRST CALL TO STEPNF . */ start = TRUE_; crash = FALSE_; hold = (float)1.; h = (float).1; s = (float)0.; ypold[1] = (float)1.; yp[1] = (float)1.; y[1] = (float)0.; i__1 = np1; for (jw = 2; jw <= i__1; ++jw) { ypold[jw] = (float)0.; yp[jw] = (float)0.; /* L40: */ } /* SET OPTIMAL STEP SIZE ESTIMATION PARAMETERS. */ /* LET Z[K] DENOTE THE NEWTON ITERATES ALONG THE FLOW NORMAL TO THE */ /* DAVIDENKO FLOW AND Y THEIR LIMIT. */ /* IDEAL CONTRACTION FACTOR: ||Z[2] - Z[1]|| / ||Z[1] - Z[0]|| */ if (sspar[1] <= (float)0.) { sspar[1] = (float).5; } /* IDEAL RESIDUAL FACTOR: ||RHO(A, Z[1])|| / ||RHO(A, Z[0])|| */ if (sspar[2] <= (float)0.) { sspar[2] = (float).01; } /* IDEAL DISTANCE FACTOR: ||Z[1] - Y|| / ||Z[0] - Y|| */ if (sspar[3] <= (float)0.) { sspar[3] = (float).5; } /* MINIMUM STEP SIZE HMIN . */ if (sspar[4] <= (float)0.) { sspar[4] = (sqrt(*n + (float)1.) + (float)4.) * d1mach_(&c__4); } /* MAXIMUM STEP SIZE HMAX . */ if (sspar[5] <= (float)0.) { sspar[5] = (float)1.; } /* MINIMUM STEP SIZE REDUCTION FACTOR BMIN . */ if (sspar[6] <= (float)0.) { sspar[6] = (float).1; } /* MAXIMUM STEP SIZE EXPANSION FACTOR BMAX . */ if (sspar[7] <= (float)0.) { sspar[7] = (float)3.; } /* ASSUMED OPERATING ORDER P . */ if (sspar[8] <= (float)0.) { sspar[8] = (float)2.; } /* Adjustments of step size */ if (tweak == 1) sspar[5] = (float).1; else if (tweak == 2) sspar[5] = (float).01; else if (tweak == 3) sspar[5] = (float).001; /* LOAD A FOR THE FIXED POINT AND ZERO FINDING PROBLEMS. */ if (iflagc >= -1) { i__1 = np1; for (jw = 2; jw <= i__1; ++jw) { a[jw - 1] = y[jw]; /* L60: */ } } L90: limit = 1000; /* Adjustment of maximum number of steps for smaller step size */ if (tweak == 1) limit = 10000; else if (tweak == 2) limit = 100000; else if (tweak == 3) limit = 1000000; /* ***** END OF INITIALIZATION BLOCK. ***** */ /* ***** MAIN LOOP. ***** */ L120: i__1 = limit; for (iter = 1; iter <= i__1; ++iter) { if (y[1] < (float)0.) { *arclen = s; *iflag = 5; return 0; } /* SET DIFFERENT ERROR TOLERANCE IF THE TRAJECTORY Y(S) HAS ANY HIGH */ /* CURVATURE COMPONENTS. */ /* L140: */ curtol = hold * 10.; relerr = *arcre; abserr = *arcae; i__2 = np1; for (jw = 1; jw <= i__2; ++jw) { if ((d__1 = yp[jw] - ypold[jw], abs((int)d__1)) > curtol) { relerr = *ansre; abserr = *ansae; goto L200; } /* L160: */ } /* TAKE A STEP ALONG THE CURVE. */ L200: stepnf_((integer *)&nc, (integer *)&nfec, (integer *)&iflagc, &start, &crash, &hold, &h, &relerr, &abserr, &s, &y[1], &yp[1], &yold[1], &ypold[1], &a[1], &qr[ qr_offset], &alpha[1], &tz[1], (integer *)&pivot[1], &w[1], &wp[1], &z0[ 1], &z1[1], &sspar[1], &par[1], (integer *)&ipar[1]); /* PRINT LATEST POINT ON CURVE IF REQUESTED. */ if (*trace > 0) { print_homog(y+2,&coord_r, &coord_i); #ifdef HOM_PRINT fprintf(Hom_LogFile,"C %g %g",coord_r,coord_i) #endif ; for(jw=2; jw<=np1; jw++){ #ifdef HOM_PRINT fprintf(Hom_LogFile," %d = %g",jw,y[jw]) #endif ; } #ifdef HOM_PRINT fprintf(Hom_LogFile," %g",y[1]) #endif ; } #ifdef HOM_PRINT fprintf(Hom_LogFile," 4 0 %d %d %f\n",iter,nfec,s) #endif ; *nfe = nfec; /* CHECK IF THE STEP WAS SUCCESSFUL. */ if (iflagc > 0) { *arclen = s; *iflag = iflagc; return 0; } if (crash) { /* RETURN CODE FOR ERROR TOLERANCE TOO SMALL. */ *iflag = 2; /* CHANGE ERROR TOLERANCES. */ if (*arcre < relerr) { *arcre = relerr; } if (*ansre < relerr) { *ansre = relerr; } if (*arcae < abserr) { *arcae = abserr; } if (*ansae < abserr) { *ansae = abserr; } /* CHANGE LIMIT ON NUMBER OF ITERATIONS. */ limit -= iter; return 0; } if (y[1] >= (float)1.) { /* USE HERMITE CUBIC INTERPOLATION AND NEWTON ITERATION TO GET THE */ /* ANSWER AT LAMBDA = 1.0 . */ /* SAVE YOLD FOR ARC LENGTH CALCULATION LATER. */ i__2 = np1; for (jw = 1; jw <= i__2; ++jw) { z0[jw] = yold[jw]; /* L260: */ } rootnf_(&nc, &nfec, &iflagc, ansre, ansae, &y[1], &yp[1], &yold[1], &ypold[1], &a[1], &qr[qr_offset], &alpha[1], &tz[1], &pivot[1], &w[1], &wp[1], &par[1], &ipar[1]); *nfe = nfec; *iflag = 1; /* SET ERROR FLAG IF ROOTNF COULD NOT GET THE POINT ON THE ZERO */ /* CURVE AT LAMBDA = 1.0 . */ if (iflagc > 0) { *iflag = iflagc; } /* CALCULATE FINAL ARC LENGTH. */ i__2 = np1; for (jw = 1; jw <= i__2; ++jw) { w[jw] = y[jw] - z0[jw]; /* L290: */ } *arclen = s - hold + dnrm2_((integer *)&np1, &w[1], &c__1); return 0; } /* FOR POLYNOMIAL SYSTEMS AND THE POLSYS HOMOTOPY MAP, */ /* D LAMBDA/DS .GE. 0 NECESSARILY. THIS CONDITION IS FORCED HERE IF */ /* THE ENTRY POINT WAS POLYNF . */ /* if (polsys) { */ if (yp[1] < (float)0.) { /* REVERSE TANGENT DIRECTION SO D LAMBDA/DS = YP(1) > 0 . */ i__2 = np1; for (jw = 1; jw <= i__2; ++jw) { yp[jw] = -yp[jw]; ypold[jw] = yp[jw]; /* L310: */ /* } */ /* FORCE STEPNF TO USE THE LINEAR PREDICTOR FOR THE NEXT STE P ONLY. */ start = TRUE_; } } /* L400: */ } /* ***** END OF MAIN LOOP. ***** */ /* LAMBDA HAS NOT REACHED 1 IN 1000 STEPS. */ *iflag = 3; *arclen = s; return 0; } /* fixpnf_ */ /* Subroutine */ int 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) { return fixpnf_0_(0, n, y, iflag, arcre, arcae, ansre, ansae, trace, a, nfe, arclen, yp, yold, ypold, qr, alpha, tz, pivot, w, wp, z0, z1, sspar, par, ipar, tweak); } /* Subroutine */ int 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) { return fixpnf_0_(1, n, y, iflag, arcre, arcae, ansre, ansae, trace, a, nfe, arclen, yp, yold, ypold, qr, alpha, tz, pivot, w, wp, z0, z1, sspar, par, ipar, tweak); } /* end fixpnf.c */ /**********************************************************************/ /***************** implementation code from Hom_Mem.c *****************/ /**********************************************************************/ /*-------------------------------------------------------------------- Workspace and Workspace management macros -------------------------------------------------------------------*/ #define TOPN 20 /*Assume no more than 20 variables*/ #define TOPM 500 /* Assume no more than 500 monomials total*/ static int Istore[TOPM*(TOPN+1)+3*TOPN+2+TOPM]; static double Dstore[TOPM*2+12*(TOPN+1)+8+TOPN*(TOPN+2)+1]; static int didx=0,iidx=0; /* initialization for storage currently static should be made dynamic*/ /* access funtions to double storage */ double *Dres(int sz){ int v=didx; didx+=sz; return Dstore+v;} int Dtop(){return didx;} int Dfree(int ntop){ didx=ntop; return 0;} /* access funtions to int storage */ int *Ires(int sz){ int v=iidx; iidx+=sz; return Istore+v;} int Itop(){return iidx;} int Ifree(int ntop){ iidx=ntop; return 0;} /* end Hom_Mem.c */ /**********************************************************************/ /***************** implementation code from d1mach.c ******************/ /**********************************************************************/ /* Table of constant values */ doublereal d1mach_(integer *i) { /* Initialized data */ static integer sc = 987; static struct { integer e_1[10]; doublereal e_2; } equiv_4 = { { 0, 1048576, -1, 2146435071, 0, 1017118720, 0, 1018167296, 1352628735, 1070810131 }, 0.0 }; /* System generated locals */ doublereal ret_val = 0.0; /* Builtin functions */ /* Subroutine */ int s_stop(); integer s_wsfe(), do_fio(), e_wsfe(); /* Local variables */ #define log10 ((integer *)&equiv_4 + 8) #define dmach ((doublereal *)&equiv_4) #define large ((integer *)&equiv_4 + 2) #define small ((integer *)&equiv_4) #define diver ((integer *)&equiv_4 + 6) #define right ((integer *)&equiv_4 + 4) /* DOUBLE-PRECISION MACHINE CONSTANTS */ /* D1MACH( 1) = B**(EMIN-1), THE SMALLEST POSITIVE MAGNITUDE. */ /* D1MACH( 2) = B**EMAX*(1 - B**(-T)), THE LARGEST MAGNITUDE. */ /* D1MACH( 3) = B**(-T), THE SMALLEST RELATIVE SPACING. */ /* D1MACH( 4) = B**(1-T), THE LARGEST RELATIVE SPACING. */ /* D1MACH( 5) = LOG10(B) */ /* TO ALTER THIS FUNCTION FOR A PARTICULAR ENVIRONMENT, */ /* THE DESIRED SET OF DATA STATEMENTS SHOULD BE ACTIVATED BY */ /* REMOVING THE C FROM COLUMN 1. */ /* ON RARE MACHINES A STATIC STATEMENT MAY NEED TO BE ADDED. */ /* (BUT PROBABLY MORE SYSTEMS PROHIBIT IT THAN REQUIRE IT.) */ /* FOR IEEE-ARITHMETIC MACHINES (BINARY STANDARD), ONE OF THE FIRST */ /* TWO SETS OF CONSTANTS BELOW SHOULD BE APPROPRIATE. IF YOU DO NOT */ /* KNOW WHICH SET TO USE, TRY BOTH AND SEE WHICH GIVES PLAUSIBLE */ /* VALUES. */ /* WHERE POSSIBLE, DECIMAL, OCTAL OR HEXADECIMAL CONSTANTS ARE USED */ /* TO SPECIFY THE CONSTANTS EXACTLY. SOMETIMES THIS REQUIRES USING */ /* EQUIVALENT INTEGER ARRAYS. IF YOUR COMPILER USES HALF-WORD */ /* INTEGERS BY DEFAULT (SOMETIMES CALLED INTEGER*2), YOU MAY NEED TO */ /* CHANGE INTEGER TO INTEGER*4 OR OTHERWISE INSTRUCT YOUR COMPILER */ /* TO USE FULL-WORD INTEGERS IN THE NEXT 5 DECLARATIONS. */ /* COMMENTS JUST BEFORE THE END STATEMENT (LINES STARTING WITH *) */ /* GIVE C SOURCE FOR D1MACH. */ /* MACHINE CONSTANTS FOR BIG-ENDIAN IEEE ARITHMETIC (BINARY FORMAT) */ /* MACHINES IN WHICH THE MOST SIGNIFICANT BYTE IS STORED FIRST, */ /* SUCH AS THE AT&T 3B SERIES, MOTOROLA 68000 BASED MACHINES (E.G. */ /* SUN 3), AND MACHINES THAT USE SPARC, HP, OR IBM RISC CHIPS. */ /* DATA SMALL(1),SMALL(2) / 1048576, 0 / */ /* DATA LARGE(1),LARGE(2) / 2146435071, -1 / */ /* DATA RIGHT(1),RIGHT(2) / 1017118720, 0 / */ /* DATA DIVER(1),DIVER(2) / 1018167296, 0 / */ /* DATA LOG10(1),LOG10(2) / 1070810131, 1352628735 /, SC/987/ */ /* MACHINE CONSTANTS FOR LITTLE-ENDIAN (BINARY) IEEE ARITHMETIC */ /* MACHINES IN WHICH THE LEAST SIGNIFICANT BYTE IS STORED FIRST, */ /* E.G. IBM PCS AND OTHER MACHINES THAT USE INTEL 80X87 OR DEC */ /* ALPHA CHIPS. */ /* MACHINE CONSTANTS FOR AMDAHL MACHINES. */ /* DATA SMALL(1),SMALL(2) / 1048576, 0 / */ /* DATA LARGE(1),LARGE(2) / 2147483647, -1 / */ /* DATA RIGHT(1),RIGHT(2) / 856686592, 0 / */ /* DATA DIVER(1),DIVER(2) / 873463808, 0 / */ /* DATA LOG10(1),LOG10(2) / 1091781651, 1352628735 /, SC/987/ */ /* MACHINE CONSTANTS FOR THE BURROUGHS 1700 SYSTEM. */ /* DATA SMALL(1) / ZC00800000 / */ /* DATA SMALL(2) / Z000000000 / */ /* DATA LARGE(1) / ZDFFFFFFFF / */ /* DATA LARGE(2) / ZFFFFFFFFF / */ /* DATA RIGHT(1) / ZCC5800000 / */ /* DATA RIGHT(2) / Z000000000 / */ /* DATA DIVER(1) / ZCC6800000 / */ /* DATA DIVER(2) / Z000000000 / */ /* DATA LOG10(1) / ZD00E730E7 / */ /* DATA LOG10(2) / ZC77800DC0 /, SC/987/ */ /* MACHINE CONSTANTS FOR THE BURROUGHS 5700 SYSTEM. */ /* DATA SMALL(1) / O1771000000000000 / */ /* DATA SMALL(2) / O0000000000000000 / */ /* DATA LARGE(1) / O0777777777777777 / */ /* DATA LARGE(2) / O0007777777777777 / */ /* DATA RIGHT(1) / O1461000000000000 / */ /* DATA RIGHT(2) / O0000000000000000 / */ /* DATA DIVER(1) / O1451000000000000 / */ /* DATA DIVER(2) / O0000000000000000 / */ /* DATA LOG10(1) / O1157163034761674 / */ /* DATA LOG10(2) / O0006677466732724 /, SC/987/ */ /* MACHINE CONSTANTS FOR THE BURROUGHS 6700/7700 SYSTEMS. */ /* DATA SMALL(1) / O1771000000000000 / */ /* DATA SMALL(2) / O7770000000000000 / */ /* DATA LARGE(1) / O0777777777777777 / */ /* DATA LARGE(2) / O7777777777777777 / */ /* DATA RIGHT(1) / O1461000000000000 / */ /* DATA RIGHT(2) / O0000000000000000 / */ /* DATA DIVER(1) / O1451000000000000 / */ /* DATA DIVER(2) / O0000000000000000 / */ /* DATA LOG10(1) / O1157163034761674 / */ /* DATA LOG10(2) / O0006677466732724 /, SC/987/ */ /* MACHINE CONSTANTS FOR FTN4 ON THE CDC 6000/7000 SERIES. */ /* DATA SMALL(1) / 00564000000000000000B / */ /* DATA SMALL(2) / 00000000000000000000B / */ /* DATA LARGE(1) / 37757777777777777777B / */ /* DATA LARGE(2) / 37157777777777777774B / */ /* DATA RIGHT(1) / 15624000000000000000B / */ /* DATA RIGHT(2) / 00000000000000000000B / */ /* DATA DIVER(1) / 15634000000000000000B / */ /* DATA DIVER(2) / 00000000000000000000B / */ /* DATA LOG10(1) / 17164642023241175717B / */ /* DATA LOG10(2) / 16367571421742254654B /, SC/987/ */ /* MACHINE CONSTANTS FOR FTN5 ON THE CDC 6000/7000 SERIES. */ /* DATA SMALL(1) / O"00564000000000000000" / */ /* DATA SMALL(2) / O"00000000000000000000" / */ /* DATA LARGE(1) / O"37757777777777777777" / */ /* DATA LARGE(2) / O"37157777777777777774" / */ /* DATA RIGHT(1) / O"15624000000000000000" / */ /* DATA RIGHT(2) / O"00000000000000000000" / */ /* DATA DIVER(1) / O"15634000000000000000" / */ /* DATA DIVER(2) / O"00000000000000000000" / */ /* DATA LOG10(1) / O"17164642023241175717" / */ /* DATA LOG10(2) / O"16367571421742254654" /, SC/987/ */ /* MACHINE CONSTANTS FOR CONVEX C-1 */ /* DATA SMALL(1),SMALL(2) / '00100000'X, '00000000'X / */ /* DATA LARGE(1),LARGE(2) / '7FFFFFFF'X, 'FFFFFFFF'X / */ /* DATA RIGHT(1),RIGHT(2) / '3CC00000'X, '00000000'X / */ /* DATA DIVER(1),DIVER(2) / '3CD00000'X, '00000000'X / */ /* DATA LOG10(1),LOG10(2) / '3FF34413'X, '509F79FF'X /, SC/987/ */ /* MACHINE CONSTANTS FOR THE CRAY 1, XMP, 2, AND 3. */ /* DATA SMALL(1) / 201354000000000000000B / */ /* DATA SMALL(2) / 000000000000000000000B / */ /* DATA LARGE(1) / 577767777777777777777B / */ /* DATA LARGE(2) / 000007777777777777776B / */ /* DATA RIGHT(1) / 376434000000000000000B / */ /* DATA RIGHT(2) / 000000000000000000000B / */ /* DATA DIVER(1) / 376444000000000000000B / */ /* DATA DIVER(2) / 000000000000000000000B / */ /* DATA LOG10(1) / 377774642023241175717B / */ /* DATA LOG10(2) / 000007571421742254654B /, SC/987/ */ /* MACHINE CONSTANTS FOR THE DATA GENERAL ECLIPSE S/200 */ /* SMALL, LARGE, RIGHT, DIVER, LOG10 SHOULD BE DECLARED */ /* INTEGER SMALL(4), LARGE(4), RIGHT(4), DIVER(4), LOG10(4) */ /* NOTE - IT MAY BE APPROPRIATE TO INCLUDE THE FOLLOWING LINE - */ /* STATIC DMACH(5) */ /* DATA SMALL/20K,3*0/,LARGE/77777K,3*177777K/ */ /* DATA RIGHT/31420K,3*0/,DIVER/32020K,3*0/ */ /* DATA LOG10/40423K,42023K,50237K,74776K/, SC/987/ */ /* MACHINE CONSTANTS FOR THE HARRIS SLASH 6 AND SLASH 7 */ /* DATA SMALL(1),SMALL(2) / '20000000, '00000201 / */ /* DATA LARGE(1),LARGE(2) / '37777777, '37777577 / */ /* DATA RIGHT(1),RIGHT(2) / '20000000, '00000333 / */ /* DATA DIVER(1),DIVER(2) / '20000000, '00000334 / */ /* DATA LOG10(1),LOG10(2) / '23210115, '10237777 /, SC/987/ */ /* MACHINE CONSTANTS FOR THE HONEYWELL DPS 8/70 SERIES. */ /* DATA SMALL(1),SMALL(2) / O402400000000, O000000000000 / */ /* DATA LARGE(1),LARGE(2) / O376777777777, O777777777777 / */ /* DATA RIGHT(1),RIGHT(2) / O604400000000, O000000000000 / */ /* DATA DIVER(1),DIVER(2) / O606400000000, O000000000000 / */ /* DATA LOG10(1),LOG10(2) / O776464202324, O117571775714 /, SC/987/ */ /* MACHINE CONSTANTS FOR THE IBM 360/370 SERIES, */ /* THE XEROX SIGMA 5/7/9 AND THE SEL SYSTEMS 85/86. */ /* DATA SMALL(1),SMALL(2) / Z00100000, Z00000000 / */ /* DATA LARGE(1),LARGE(2) / Z7FFFFFFF, ZFFFFFFFF / */ /* DATA RIGHT(1),RIGHT(2) / Z33100000, Z00000000 / */ /* DATA DIVER(1),DIVER(2) / Z34100000, Z00000000 / */ /* DATA LOG10(1),LOG10(2) / Z41134413, Z509F79FF /, SC/987/ */ /* MACHINE CONSTANTS FOR THE INTERDATA 8/32 */ /* WITH THE UNIX SYSTEM FORTRAN 77 COMPILER. */ /* FOR THE INTERDATA FORTRAN VII COMPILER REPLACE */ /* THE Z'S SPECIFYING HEX CONSTANTS WITH Y'S. */ /* DATA SMALL(1),SMALL(2) / Z'00100000', Z'00000000' / */ /* DATA LARGE(1),LARGE(2) / Z'7EFFFFFF', Z'FFFFFFFF' / */ /* DATA RIGHT(1),RIGHT(2) / Z'33100000', Z'00000000' / */ /* DATA DIVER(1),DIVER(2) / Z'34100000', Z'00000000' / */ /* DATA LOG10(1),LOG10(2) / Z'41134413', Z'509F79FF' /, SC/987/ */ /* MACHINE CONSTANTS FOR THE PDP-10 (KA PROCESSOR). */ /* DATA SMALL(1),SMALL(2) / "033400000000, "000000000000 / */ /* DATA LARGE(1),LARGE(2) / "377777777777, "344777777777 / */ /* DATA RIGHT(1),RIGHT(2) / "113400000000, "000000000000 / */ /* DATA DIVER(1),DIVER(2) / "114400000000, "000000000000 / */ /* DATA LOG10(1),LOG10(2) / "177464202324, "144117571776 /, SC/987/ */ /* MACHINE CONSTANTS FOR THE PDP-10 (KI PROCESSOR). */ /* DATA SMALL(1),SMALL(2) / "000400000000, "000000000000 / */ /* DATA LARGE(1),LARGE(2) / "377777777777, "377777777777 / */ /* DATA RIGHT(1),RIGHT(2) / "103400000000, "000000000000 / */ /* DATA DIVER(1),DIVER(2) / "104400000000, "000000000000 / */ /* DATA LOG10(1),LOG10(2) / "177464202324, "047674776746 /, SC/987/ */ /* MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING */ /* 32-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). */ /* DATA SMALL(1),SMALL(2) / 8388608, 0 / */ /* DATA LARGE(1),LARGE(2) / 2147483647, -1 / */ /* DATA RIGHT(1),RIGHT(2) / 612368384, 0 / */ /* DATA DIVER(1),DIVER(2) / 620756992, 0 / */ /* DATA LOG10(1),LOG10(2) / 1067065498, -2063872008 /, SC/987/ */ /* DATA SMALL(1),SMALL(2) / O00040000000, O00000000000 / */ /* DATA LARGE(1),LARGE(2) / O17777777777, O37777777777 / */ /* DATA RIGHT(1),RIGHT(2) / O04440000000, O00000000000 / */ /* DATA DIVER(1),DIVER(2) / O04500000000, O00000000000 / */ /* DATA LOG10(1),LOG10(2) / O07746420232, O20476747770 /, SC/987/ */ /* MACHINE CONSTANTS FOR PDP-11 FORTRANS SUPPORTING */ /* 16-BIT INTEGERS (EXPRESSED IN INTEGER AND OCTAL). */ /* SMALL, LARGE, RIGHT, DIVER, LOG10 SHOULD BE DECLARED */ /* INTEGER SMALL(4), LARGE(4), RIGHT(4), DIVER(4), LOG10(4) */ /* DATA SMALL(1),SMALL(2) / 128, 0 / */ /* DATA SMALL(3),SMALL(4) / 0, 0 / */ /* DATA LARGE(1),LARGE(2) / 32767, -1 / */ /* DATA LARGE(3),LARGE(4) / -1, -1 / */ /* DATA RIGHT(1),RIGHT(2) / 9344, 0 / */ /* DATA RIGHT(3),RIGHT(4) / 0, 0 / */ /* DATA DIVER(1),DIVER(2) / 9472, 0 / */ /* DATA DIVER(3),DIVER(4) / 0, 0 / */ /* DATA LOG10(1),LOG10(2) / 16282, 8346 / */ /* DATA LOG10(3),LOG10(4) / -31493, -12296 /, SC/987/ */ /* DATA SMALL(1),SMALL(2) / O000200, O000000 / */ /* DATA SMALL(3),SMALL(4) / O000000, O000000 / */ /* DATA LARGE(1),LARGE(2) / O077777, O177777 / */ /* DATA LARGE(3),LARGE(4) / O177777, O177777 / */ /* DATA RIGHT(1),RIGHT(2) / O022200, O000000 / */ /* DATA RIGHT(3),RIGHT(4) / O000000, O000000 / */ /* DATA DIVER(1),DIVER(2) / O022400, O000000 / */ /* DATA DIVER(3),DIVER(4) / O000000, O000000 / */ /* DATA LOG10(1),LOG10(2) / O037632, O020232 / */ /* DATA LOG10(3),LOG10(4) / O102373, O147770 /, SC/987/ */ /* MACHINE CONSTANTS FOR THE PRIME 50 SERIES SYSTEMS */ /* WITH 32-BIT INTEGERS AND 64V MODE INSTRUCTIONS, */ /* SUPPLIED BY IGOR BRAY. */ /* DATA SMALL(1),SMALL(2) / :10000000000, :00000100001 / */ /* DATA LARGE(1),LARGE(2) / :17777777777, :37777677775 / */ /* DATA RIGHT(1),RIGHT(2) / :10000000000, :00000000122 / */ /* DATA DIVER(1),DIVER(2) / :10000000000, :00000000123 / */ /* DATA LOG10(1),LOG10(2) / :11504046501, :07674600177 /, SC/987/ */ /* MACHINE CONSTANTS FOR THE SEQUENT BALANCE 8000 */ /* DATA SMALL(1),SMALL(2) / $00000000, $00100000 / */ /* DATA LARGE(1),LARGE(2) / $FFFFFFFF, $7FEFFFFF / */ /* DATA RIGHT(1),RIGHT(2) / $00000000, $3CA00000 / */ /* DATA DIVER(1),DIVER(2) / $00000000, $3CB00000 / */ /* DATA LOG10(1),LOG10(2) / $509F79FF, $3FD34413 /, SC/987/ */ /* MACHINE CONSTANTS FOR THE UNIVAC 1100 SERIES. */ /* DATA SMALL(1),SMALL(2) / O000040000000, O000000000000 / */ /* DATA LARGE(1),LARGE(2) / O377777777777, O777777777777 / */ /* DATA RIGHT(1),RIGHT(2) / O170540000000, O000000000000 / */ /* DATA DIVER(1),DIVER(2) / O170640000000, O000000000000 / */ /* DATA LOG10(1),LOG10(2) / O177746420232, O411757177572 /, SC/987/ */ /* MACHINE CONSTANTS FOR THE VAX UNIX F77 COMPILER */ /* DATA SMALL(1),SMALL(2) / 128, 0 / */ /* DATA LARGE(1),LARGE(2) / -32769, -1 / */ /* DATA RIGHT(1),RIGHT(2) / 9344, 0 / */ /* DATA DIVER(1),DIVER(2) / 9472, 0 / */ /* DATA LOG10(1),LOG10(2) / 546979738, -805796613 /, SC/987/ */ /* MACHINE CONSTANTS FOR THE VAX-11 WITH */ /* FORTRAN IV-PLUS COMPILER */ /* DATA SMALL(1),SMALL(2) / Z00000080, Z00000000 / */ /* DATA LARGE(1),LARGE(2) / ZFFFF7FFF, ZFFFFFFFF / */ /* DATA RIGHT(1),RIGHT(2) / Z00002480, Z00000000 / */ /* DATA DIVER(1),DIVER(2) / Z00002500, Z00000000 / */ /* DATA LOG10(1),LOG10(2) / Z209A3F9A, ZCFF884FB /, SC/987/ */ /* MACHINE CONSTANTS FOR VAX/VMS VERSION 2.2 */ /* DATA SMALL(1),SMALL(2) / '80'X, '0'X / */ /* DATA LARGE(1),LARGE(2) / 'FFFF7FFF'X, 'FFFFFFFF'X / */ /* DATA RIGHT(1),RIGHT(2) / '2480'X, '0'X / */ /* DATA DIVER(1),DIVER(2) / '2500'X, '0'X / */ /* DATA LOG10(1),LOG10(2) / '209A3F9A'X, 'CFF884FB'X /, SC/987/ */ /* *** ISSUE STOP 779 IF ALL DATA STATEMENTS ARE COMMENTED... */ if (sc != 987) { /* s_stop("779", 3L); */exit(0); } /* *** ISSUE STOP 778 IF ALL DATA STATEMENTS ARE OBVIOUSLY WRONG... */ if (dmach[3] >= 1.) { /* s_stop("778", 3L); */ exit(0); } if (*i < 1 || *i > 5) { goto L999; } ret_val = dmach[*i - 1]; return ret_val; L999: /* s_wsfe(&io___8); do_fio(&c__1, (char *)&(*i), (ftnlen)sizeof(integer)); e_wsfe(); s_stop("", 0L); */ return ret_val; } /* d1mach_ */ #undef right #undef diver #undef small #undef large #undef dmach #undef log10 /* end d1mach.c */ /**********************************************************************/ /***************** implementation code from daxpy.c *******************/ /**********************************************************************/ /* Subroutine */ int daxpy_(integer *n, doublereal *da, doublereal *dx, integer *incx, doublereal *dy, integer *incy) { /* System generated locals */ integer i__1; /* Local variables */ static integer i, m, ix, iy, mp1; /* constant times a vector plus a vector. */ /* uses unrolled loops for increments equal to one. */ /* jack dongarra, linpack, 3/11/78. */ /* Parameter adjustments */ --dy; --dx; /* Function Body */ if (*n <= 0) { return 0; } if (*da == 0.) { return 0; } if (*incx == 1 && *incy == 1) { goto L20; } /* code for unequal increments or equal increments */ /* not equal to 1 */ ix = 1; iy = 1; if (*incx < 0) { ix = (-(*n) + 1) * *incx + 1; } if (*incy < 0) { iy = (-(*n) + 1) * *incy + 1; } i__1 = *n; for (i = 1; i <= i__1; ++i) { dy[iy] += *da * dx[ix]; ix += *incx; iy += *incy; /* L10: */ } return 0; /* code for both increments equal to 1 */ /* clean-up loop */ L20: m = *n % 4; if (m == 0) { goto L40; } i__1 = m; for (i = 1; i <= i__1; ++i) { dy[i] += *da * dx[i]; /* L30: */ } if (*n < 4) { return 0; } L40: mp1 = m + 1; i__1 = *n; for (i = mp1; i <= i__1; i += 4) { dy[i] += *da * dx[i]; dy[i + 1] += *da * dx[i + 1]; dy[i + 2] += *da * dx[i + 2]; dy[i + 3] += *da * dx[i + 3]; /* L50: */ } return 0; } /* daxpy_ */ /* end daxpy.c */ /**********************************************************************/ /***************** implementation code from dcopy.c *******************/ /**********************************************************************/ /* Subroutine */ int dcopy_(integer *n, doublereal *dx, integer *incx, doublereal *dy, integer *incy) { /* System generated locals */ integer i__1; /* Local variables */ static integer i, m, ix, iy, mp1; /* copies a vector, x, to a vector, y. */ /* uses unrolled loops for increments equal to one. */ /* jack dongarra, linpack, 3/11/78. */ /* Parameter adjustments */ --dy; --dx; /* Function Body */ if (*n <= 0) { return 0; } if (*incx == 1 && *incy == 1) { goto L20; } /* code for unequal increments or equal increments */ /* not equal to 1 */ ix = 1; iy = 1; if (*incx < 0) { ix = (-(*n) + 1) * *incx + 1; } if (*incy < 0) { iy = (-(*n) + 1) * *incy + 1; } i__1 = *n; for (i = 1; i <= i__1; ++i) { dy[iy] = dx[ix]; ix += *incx; iy += *incy; /* L10: */ } return 0; /* code for both increments equal to 1 */ /* clean-up loop */ L20: m = *n % 7; if (m == 0) { goto L40; } i__1 = m; for (i = 1; i <= i__1; ++i) { dy[i] = dx[i]; /* L30: */ } if (*n < 7) { return 0; } L40: mp1 = m + 1; i__1 = *n; for (i = mp1; i <= i__1; i += 7) { dy[i] = dx[i]; dy[i + 1] = dx[i + 1]; dy[i + 2] = dx[i + 2]; dy[i + 3] = dx[i + 3]; dy[i + 4] = dx[i + 4]; dy[i + 5] = dx[i + 5]; dy[i + 6] = dx[i + 6]; /* L50: */ } return 0; } /* dcopy_ */ /* end dcopy.c */ /**********************************************************************/ /***************** implementation code from dcpose.c ******************/ /**********************************************************************/ /* Table of constant values */ /* static integer c__1 = 1; NOW REDUNDANT */ /* Subroutine */ int dcpose_(integer *ndim, integer *n, doublereal *qr, doublereal *alpha, integer *pivot, integer *ierr, doublereal *y, doublereal *sum) { /* System generated locals */ integer qr_dim1, qr_offset, i__1, i__2, i__3; doublereal d__1; /* Builtin functions */ /* double sqrt(double); CANT DECLARE BUILTINS UNDER C++ */ /* Local variables */ static doublereal beta; static integer jbar; extern doublereal ddot_(integer *n, doublereal *dx, integer *incx, doublereal *dy, integer *incy); static doublereal qrkk; static integer i, j, k; static doublereal sigma, alphak; static integer kp1, np1; /* SUBROUTINE DCPOSE IS A MODIFICATION OF THE ALGOL PROCEDURE */ /* DECOMPOSE IN P. BUSINGER AND G. H. GOLUB, LINEAR LEAST */ /* SQUARES SOLUTIONS BY HOUSEHOLDER TRANSFORMATIONS, */ /* NUMER. MATH. 7 (1965) 269-276. */ /* Parameter adjustments */ qr_dim1 = *ndim; qr_offset = qr_dim1 + 1; qr -= qr_offset; --alpha; --pivot; --y; --sum; /* Function Body */ *ierr = 0; np1 = *n + 1; i__1 = np1; for (j = 1; j <= i__1; ++j) { sum[j] = ddot_(n, &qr[j * qr_dim1 + 1], &c__1, &qr[j * qr_dim1 + 1], & c__1); /* L20: */ pivot[j] = j; } i__1 = *n; for (k = 1; k <= i__1; ++k) { sigma = sum[k]; jbar = k; kp1 = k + 1; i__2 = np1; for (j = kp1; j <= i__2; ++j) { if (sigma >= sum[j]) { goto L40; } sigma = sum[j]; jbar = j; L40: ; } if (jbar == k) { goto L70; } i = pivot[k]; pivot[k] = pivot[jbar]; pivot[jbar] = i; sum[jbar] = sum[k]; sum[k] = sigma; i__2 = *n; for (i = 1; i <= i__2; ++i) { sigma = qr[i + k * qr_dim1]; qr[i + k * qr_dim1] = qr[i + jbar * qr_dim1]; qr[i + jbar * qr_dim1] = sigma; /* L50: */ } /* END OF COLUMN INTERCHANGE. */ L70: i__2 = *n - k + 1; sigma = ddot_(&i__2, &qr[k + k * qr_dim1], &c__1, &qr[k + k * qr_dim1] , &c__1); if (sigma != (float)0.) { goto L60; } *ierr = 1; return 0; L60: if (k == *n) { goto L500; } qrkk = qr[k + k * qr_dim1]; alphak = -sqrt(sigma); if (qrkk < (float)0.) { alphak = -alphak; } alpha[k] = alphak; beta = (float)1. / (sigma - qrkk * alphak); qr[k + k * qr_dim1] = qrkk - alphak; i__2 = np1; for (j = kp1; j <= i__2; ++j) { /* L80: */ i__3 = *n - k + 1; y[j] = beta * ddot_(&i__3, &qr[k + k * qr_dim1], &c__1, &qr[k + j * qr_dim1], &c__1); } i__3 = np1; for (j = kp1; j <= i__3; ++j) { i__2 = *n; for (i = k; i <= i__2; ++i) { qr[i + j * qr_dim1] -= qr[i + k * qr_dim1] * y[j]; /* L90: */ } /* Computing 2nd power */ d__1 = qr[k + j * qr_dim1]; sum[j] -= d__1 * d__1; /* L100: */ } L500: ; } alpha[*n] = qr[*n + *n * qr_dim1]; return 0; } /* dcpose_ */ /* end dcpose.c */ /**********************************************************************/ /****************** implementation code from ddot.c *******************/ /**********************************************************************/ doublereal ddot_(integer *n, doublereal *dx, integer *incx, doublereal *dy, integer *incy) { /* System generated locals */ integer i__1; doublereal ret_val; /* Local variables */ static integer i, m; static doublereal dtemp; static integer ix, iy, mp1; /* forms the dot product of two vectors. */ /* uses unrolled loops for increments equal to one. */ /* jack dongarra, linpack, 3/11/78. */ /* Parameter adjustments */ --dy; --dx; /* Function Body */ ret_val = 0.; dtemp = 0.; if (*n <= 0) { return ret_val; } if (*incx == 1 && *incy == 1) { goto L20; } /* code for unequal increments or equal increments */ /* not equal to 1 */ ix = 1; iy = 1; if (*incx < 0) { ix = (-(*n) + 1) * *incx + 1; } if (*incy < 0) { iy = (-(*n) + 1) * *incy + 1; } i__1 = *n; for (i = 1; i <= i__1; ++i) { dtemp += dx[ix] * dy[iy]; ix += *incx; iy += *incy; /* L10: */ } ret_val = dtemp; return ret_val; /* code for both increments equal to 1 */ /* clean-up loop */ L20: m = *n % 5; if (m == 0) { goto L40; } i__1 = m; for (i = 1; i <= i__1; ++i) { dtemp += dx[i] * dy[i]; /* L30: */ } if (*n < 5) { goto L60; } L40: mp1 = m + 1; i__1 = *n; for (i = mp1; i <= i__1; i += 5) { dtemp = dtemp + dx[i] * dy[i] + dx[i + 1] * dy[i + 1] + dx[i + 2] * dy[i + 2] + dx[i + 3] * dy[i + 3] + dx[i + 4] * dy[i + 4]; /* L50: */ } L60: ret_val = dtemp; return ret_val; } /* ddot_ */ /* end ddot.c */ /**********************************************************************/ /******************* implementation code from divp.c *********************/ /**********************************************************************/ /* Table of constant values */ static integer c__2 = 2; /* Subroutine */ int divp_(doublereal *xxxx, doublereal *yyyy, doublereal *zzzz, integer *ierr) { static doublereal xnum, denom; extern doublereal d1mach_(integer *i); /* THIS SUBROUTINE PERFORMS DIVISION OF COMPLEX NUMBERS: */ /* ZZZZ = XXXX/YYYY */ /* ON INPUT: */ /* XXXX IS AN ARRAY OF LENGTH TWO REPRESENTING THE FIRST COMPLEX */ /* NUMBER, WHERE XXXX(1) = REAL PART OF XXXX AND XXXX(2) = */ /* IMAGINARY PART OF XXXX. */ /* YYYY IS AN ARRAY OF LENGTH TWO REPRESENTING THE SECOND COMPLEX */ /* NUMBER, WHERE YYYY(1) = REAL PART OF YYYY AND YYYY(2) = */ /* IMAGINARY PART OF YYYY. */ /* ON OUTPUT: */ /* ZZZZ IS AN ARRAY OF LENGTH TWO REPRESENTING THE RESULT OF */ /* THE DIVISION, ZZZZ = XXXX/YYYY, WHERE ZZZZ(1) = */ /* REAL PART OF ZZZZ AND ZZZZ(2) = IMAGINARY PART OF ZZZZ. */ /* IERR = */ /* 1 IF DIVISION WOULD HAVE CAUSED OVERFLOW. IN THIS CASE, THE */ /* APPROPRIATE PARTS OF ZZZZ ARE SET EQUAL TO THE LARGEST */ /* FLOATING POINT NUMBER, AS GIVEN BY FUNCTION D1MACH . */ /* 0 IF DIVISION DOES NOT CAUSE OVERFLOW. */ /* DECLARATION OF INPUT */ /* DECLARATION OF OUTPUT */ /* DECLARATION OF VARIABLES */ /* Parameter adjustments */ --zzzz; --yyyy; --xxxx; /* Function Body */ *ierr = 0; denom = yyyy[1] * yyyy[1] + yyyy[2] * yyyy[2]; xnum = xxxx[1] * yyyy[1] + xxxx[2] * yyyy[2]; if (dblabs(denom) >= (float)1. || (dblabs(denom) < (float)1. && dblabs(xnum) / d1mach_(&c__2) < dblabs(denom)) ) { zzzz[1] = xnum / denom; } else { zzzz[1] = d1mach_(&c__2); *ierr = 1; } xnum = xxxx[2] * yyyy[1] - xxxx[1] * yyyy[2]; if (dblabs(denom) >= (float)1. || (dblabs(denom) < (float)1. && dblabs(xnum) / d1mach_(&c__2) < dblabs(denom))) { zzzz[2] = xnum / denom; } else { zzzz[2] = d1mach_(&c__2); *ierr = 1; } return 0; } /* divp_ */ /* end divp.c */ /**********************************************************************/ /****************** implementation code from dnrm2.c ******************/ /**********************************************************************/ doublereal dnrm2_(integer *n, doublereal *dx, integer *incx) { /* Initialized data */ static doublereal zero = 0.; static doublereal one = 1.; static doublereal cutlo = 8.232e-11; static doublereal cuthi = 1.304e19; /* Format strings */ static char fmt_30[] = ""; static char fmt_50[] = ""; static char fmt_70[] = ""; static char fmt_110[] = ""; /* System generated locals */ integer i__1; doublereal ret_val, d__1; /* Builtin functions */ /* double sqrt(double); CANT DECLARE BUILTINS UNDER C++ */ /* Local variables */ static doublereal xmax; static integer next, i, j, ix; static doublereal hitest, sum; /* Assigned format variables */ char *next_fmt; /* Parameter adjustments */ --dx; /* Function Body */ /* euclidean norm of the n-vector stored in dx() with storage */ /* increment incx . */ /* if n .le. 0 return with result = 0. */ /* if n .ge. 1 then incx must be .ge. 1 */ /* c.l.lawson, 1978 jan 08 */ /* modified to correct failure to update ix, 1/25/92. */ /* modified 3/93 to return if incx .le. 0. */ /* four phase method using two built-in constants that are */ /* hopefully applicable to all machines. */ /* cutlo = maximum of dsqrt(u/eps) over all known machines. */ /* cuthi = minimum of dsqrt(v) over all known machines. */ /* where */ /* eps = smallest no. such that eps + 1. .gt. 1. */ /* u = smallest positive no. (underflow limit) */ /* v = largest no. (overflow limit) */ /* brief outline of algorithm.. */ /* phase 1 scans zero components. */ /* move to phase 2 when a component is nonzero and .le. cutlo */ /* move to phase 3 when a component is .gt. cutlo */ /* move to phase 4 when a component is .ge. cuthi/m */ /* where m = n for x() real and m = 2*n for complex. */ /* values for cutlo and cuthi.. */ /* from the environmental parameters listed in the imsl converter */ /* document the limiting values are as follows.. */ /* cutlo, s.p. u/eps = 2**(-102) for honeywell. close seconds are */ /* univac and dec at 2**(-103) */ /* thus cutlo = 2**(-51) = 4.44089e-16 */ /* cuthi, s.p. v = 2**127 for univac, honeywell, and dec. */ /* thus cuthi = 2**(63.5) = 1.30438e19 */ /* cutlo, d.p. u/eps = 2**(-67) for honeywell and dec. */ /* thus cutlo = 2**(-33.5) = 8.23181d-11 */ /* cuthi, d.p. same as s.p. cuthi = 1.30438d19 */ /* data cutlo, cuthi / 8.232d-11, 1.304d19 / */ /* data cutlo, cuthi / 4.441e-16, 1.304e19 / */ if (*n > 0 && *incx > 0) { goto L10; } ret_val = zero; goto L300; L10: next = 0; next_fmt = fmt_30; sum = zero; i = 1; ix = 1; /* begin main loop */ L20: switch ((int)next) { case 0: goto L30; case 1: goto L50; case 2: goto L70; case 3: goto L110; } L30: if ((d__1 = dx[i], dblabs(d__1)) > cutlo) { goto L85; } next = 1; next_fmt = fmt_50; xmax = zero; /* phase 1. sum is zero */ L50: if (dx[i] == zero) { goto L200; } if ((d__1 = dx[i], dblabs(d__1)) > cutlo) { goto L85; } /* prepare for phase 2. */ next = 2; next_fmt = fmt_70; goto L105; /* prepare for phase 4. */ L100: ix = j; next = 3; next_fmt = fmt_110; sum = sum / dx[i] / dx[i]; L105: xmax = (d__1 = dx[i], dblabs(d__1)); goto L115; /* phase 2. sum is small. */ /* scale to avoid destructive underflow. */ L70: if ((d__1 = dx[i], dblabs(d__1)) > cutlo) { goto L75; } /* common code for phases 2 and 4. */ /* in phase 4 sum is large. scale to avoid overflow. */ L110: if ((d__1 = dx[i], dblabs(d__1)) <= xmax) { goto L115; } /* Computing 2nd power */ d__1 = xmax / dx[i]; sum = one + sum * (d__1 * d__1); xmax = (d__1 = dx[i], dblabs(d__1)); goto L200; L115: /* Computing 2nd power */ d__1 = dx[i] / xmax; sum += d__1 * d__1; goto L200; /* prepare for phase 3. */ L75: sum = sum * xmax * xmax; /* for real or d.p. set hitest = cuthi/n */ /* for complex set hitest = cuthi/(2*n) */ L85: hitest = cuthi / (real) (*n); /* phase 3. sum is mid-range. no scaling. */ i__1 = *n; for (j = ix; j <= i__1; ++j) { if ((d__1 = dx[i], dblabs(d__1)) >= hitest) { goto L100; } /* Computing 2nd power */ d__1 = dx[i]; sum += d__1 * d__1; i += *incx; /* L95: */ } ret_val = sqrt(sum); goto L300; L200: ++ix; i += *incx; if (ix <= *n) { goto L20; } /* end of main loop. */ /* compute square root and adjust for scaling. */ ret_val = xmax * sqrt(sum); L300: return ret_val; } /* dnrm2_ */ /* end dnrm2.c */ /**********************************************************************/ /***************** implementation code from dscal.c *******************/ /**********************************************************************/ /* Subroutine */ int dscal_(integer *n, doublereal *da, doublereal *dx, integer *incx) { /* System generated locals */ integer i__1, i__2; /* Local variables */ static integer i, m, nincx, mp1; /* scales a vector by a constant. */ /* uses unrolled loops for increment equal to one. */ /* jack dongarra, linpack, 3/11/78. */ /* modified 3/93 to return if incx .le. 0. */ /* Parameter adjustments */ --dx; /* Function Body */ if (*n <= 0 || *incx <= 0) { return 0; } if (*incx == 1) { goto L20; } /* code for increment not equal to 1 */ nincx = *n * *incx; i__1 = nincx; i__2 = *incx; for (i = 1; i__2 < 0 ? i >= i__1 : i <= i__1; i += i__2) { dx[i] = *da * dx[i]; /* L10: */ } return 0; /* code for increment equal to 1 */ /* clean-up loop */ L20: m = *n % 5; if (m == 0) { goto L40; } i__2 = m; for (i = 1; i <= i__2; ++i) { dx[i] = *da * dx[i]; /* L30: */ } if (*n < 5) { return 0; } L40: mp1 = m + 1; i__2 = *n; for (i = mp1; i <= i__2; i += 5) { dx[i] = *da * dx[i]; dx[i + 1] = *da * dx[i + 1]; dx[i + 2] = *da * dx[i + 2]; dx[i + 3] = *da * dx[i + 3]; dx[i + 4] = *da * dx[i + 4]; /* L50: */ } return 0; } /* dscal_ */ /* end dscal.c */ /**********************************************************************/ /******************* implementation code from f.c *********************/ /**********************************************************************/ /* Subroutine */ int f_(doublereal *x, doublereal *v) { /* EVALUATE F(X) AND RETURN IN THE VECTOR V . */ /* Parameter adjustments */ --v; --x; /* Function Body */ return 0; } /* f_ */ /* end f.c */ /**********************************************************************/ /***************** implementation code from fjac.c ********************/ /**********************************************************************/ /* Subroutine */ int fjac_(doublereal *x, doublereal *v, integer *k) { /* RETURN IN V THE KTH COLUMN OF THE JACOBIAN MATRIX OF */ /* F(X) EVALUATED AT X . */ /* Parameter adjustments */ --v; --x; /* Function Body */ return 0; } /* fjac_ */ /* end fjac.c */ /**********************************************************************/ /***************** implementation code from idamax.c ******************/ /**********************************************************************/ integer idamax_(integer *n, doublereal *dx, integer *incx) { /* System generated locals */ integer ret_val, i__1; doublereal d__1; /* Local variables */ static doublereal dmax_; static integer i, ix; /* finds the index of element having max. absolute value. */ /* jack dongarra, linpack, 3/11/78. */ /* modified 3/93 to return if incx .le. 0. */ /* Parameter adjustments */ --dx; /* Function Body */ ret_val = 0; if (*n < 1 || *incx <= 0) { return ret_val; } ret_val = 1; if (*n == 1) { return ret_val; } if (*incx == 1) { goto L20; } /* code for increment not equal to 1 */ ix = 1; dmax_ = dblabs(dx[1]); ix += *incx; i__1 = *n; for (i = 2; i <= i__1; ++i) { if ((d__1 = dx[ix], dblabs(d__1)) <= dmax_) { goto L5; } ret_val = i; dmax_ = (d__1 = dx[ix], dblabs(d__1)); L5: ix += *incx; /* L10: */ } return ret_val; /* code for increment equal to 1 */ L20: dmax_ = dblabs(dx[1]); i__1 = *n; for (i = 2; i <= i__1; ++i) { if ((d__1 = dx[i], dblabs(d__1)) <= dmax_) { goto L30; } ret_val = i; dmax_ = (d__1 = dx[i], dblabs(d__1)); L30: ; } return ret_val; } /* idamax_ */ /* end idamax.c */ /**********************************************************************/ /****************** implementation code from root.c *******************/ /**********************************************************************/ double d_sign(double *arg1, double *arg2){ if (*arg2>=0.0) return dblabs(*arg1); else return -1.0*dblabs(*arg1); } /* Table of constant values */ /* static integer c__4 = 4; NOW REDUNDANT */ static doublereal c_b17 = 1.; /* Subroutine */ int root_(doublereal *t, doublereal *ft, doublereal *b, doublereal *c, doublereal *relerr, doublereal *abserr, integer *iflag) { /* System generated locals */ doublereal d__1, d__2; /* Builtin functions */ double d_sign(double *arg1, double *arg2); /* Local variables */ static doublereal acmb, acbs, a, p, q, u; extern doublereal d1mach_(integer *i); static integer kount; static doublereal ae, fa, fb, fc; static integer ic; static doublereal re, fx, cmb, tol; /* ROOT COMPUTES A ROOT OF THE NONLINEAR EQUATION F(X)=0 */ /* WHERE F(X) IS A CONTINOUS REAL FUNCTION OF A SINGLE REAL */ /* VARIABLE X. THE METHOD USED IS A COMBINATION OF BISECTION */ /* AND THE SECANT RULE. */ /* NORMAL INPUT CONSISTS OF A CONTINUOS FUNCTION F AND AN */ /* INTERVAL (B,C) SUCH THAT F(B)*F(C).LE.0.0. EACH ITERATION */ /* FINDS NEW VALUES OF B AND C SUCH THAT THE INTERVAL(B,C) IS */ /* SHRUNK AND F(B)*F(C).LE.0.0. THE STOPPING CRITERION IS */ /* DABS(B-C).LE.2.0*(RELERR*DABS(B)+ABSERR) */ /* WHERE RELERR=RELATIVE ERROR AND ABSERR=ABSOLUTE ERROR ARE */ /* INPUT QUANTITIES. SET THE FLAG, IFLAG, POSITIVE TO INITIALIZE */ /* THE COMPUTATION. AS B,C AND IFLAG ARE USED FOR BOTH INPUT AND */ /* OUTPUT, THEY MUST BE VARIABLES IN THE CALLING PROGRAM. */ /* IF 0 IS A POSSIBLE ROOT, ONE SHOULD NOT CHOOSE ABSERR=0.0. */ /* THE OUTPUT VALUE OF B IS THE BETTER APPROXIMATION TO A ROOT */ /* AS B AND C ARE ALWAYS REDEFINED SO THAT DABS(F(B)).LE.DABS(F(C)). */ /* TO SOLVE THE EQUATION, ROOT MUST EVALUATE F(X) REPEATEDLY. THIS */ /* IS DONE IN THE CALLING PROGRAM. WHEN AN EVALUATION OF F IS */ /* NEEDED AT T, ROOT RETURNS WITH IFLAG NEGATIVE. EVALUATE FT=F(T) */ /* AND CALL ROOT AGAIN. DO NOT ALTER IFLAG. */ /* WHEN THE COMPUTATION IS COMPLETE, ROOT RETURNS TO THE CALLING */ /* PROGRAM WITH IFLAG POSITIVE= */ /* IFLAG=1 IF F(B)*F(C).LT.0 AND THE STOPPING CRITERION IS MET. */ /* =2 IF A VALUE B IS FOUND SUCH THAT THE COMPUTED VALUE */ /* F(B) IS EXACTLY ZERO. THE INTERVAL (B,C) MAY NOT */ /* SATISFY THE STOPPING CRITERION. */ /* =3 IF DABS(F(B)) EXCEEDS THE INPUT VALUES DABS(F(B)), */ /* DABS(F(C)). IN THIS CASE IT IS LIKELY THAT B IS CLOSE */ /* TO A POLE OF F. */ /* =4 IF NO ODD ORDER ROOT WAS FOUND IN THE INTERVAL. A */ /* LOCAL MINIMUM MAY HAVE BEEN OBTAINED. */ /* =5 IF TOO MANY FUNCTION EVALUATIONS WERE MADE. */ /* (AS PROGRAMMED, 500 ARE ALLOWED.) */ /* THIS CODE IS A MODIFICATION OF THE CODE ZEROIN WHICH IS COMPLETELY */ /* EXPLAINED AND DOCUMENTED IN THE TEXT NUMERICAL COMPUTING: AN */ /* INTRODUCTION, BY L. F. SHAMPINE AND R. C. ALLEN. */ /* CALLS D1MACH . */ if (*iflag >= 0) { goto L100; } *iflag = dblabs(*iflag); switch ((int)*iflag) { case 1: goto L200; case 2: goto L300; case 3: goto L400; } L100: u = d1mach_(&c__4); re = max(*relerr,u); ae = max(*abserr,0.); ic = 0; acbs = (d__1 = *b - *c, dblabs(d__1)); a = *c; *t = a; *iflag = -1; return 0; L200: fa = *ft; *t = *b; *iflag = -2; return 0; L300: fb = *ft; fc = fa; kount = 2; /* Computing MAX */ d__1 = dblabs(fb), d__2 = dblabs(fc); fx = max(d__1,d__2); L1: if (dblabs(fc) >= dblabs(fb)) { goto L2; } /* INTERCHANGE B AND C SO THAT ABS(F(B)).LE.ABS(F(C)). */ a = *b; fa = fb; *b = *c; fb = fc; *c = a; fc = fa; L2: cmb = (*c - *b) * (float).5; acmb = dblabs(cmb); tol = re * dblabs(*b) + ae; /* TEST STOPPING CRITERION AND FUNCTION COUNT. */ if (acmb <= tol) { goto L8; } if (kount >= 500) { goto L12; } /* CALCULATE NEW ITERATE EXPLICITLY AS B+P/Q */ /* WHERE WE ARRANGE P.GE.0. THE IMPLICIT */ /* FORM IS USED TO PREVENT OVERFLOW. */ p = (*b - a) * fb; q = fa - fb; if (p >= (float)0.) { goto L3; } p = -p; q = -q; /* UPDATE A, CHECK IF REDUCTION IN THE SIZE OF BRACKETING */ /* INTERVAL IS SATISFACTORY. IF NOT BISECT UNTIL IT IS. */ L3: a = *b; fa = fb; ++ic; if (ic < 4) { goto L4; } if (acmb * (float)8. >= acbs) { goto L6; } ic = 0; acbs = acmb; /* TEST FOR TOO SMALL A CHANGE. */ L4: if (p > dblabs(q) * tol) { goto L5; } /* INCREMENT BY TOLERANCE */ /* sign gives sign(arg2)*|arg1| */ *b += d_sign(&tol, &cmb); goto L7; /* ROOT OUGHT TO BE BETWEEN B AND (C+B)/2 */ L5: if (p >= cmb * q) { goto L6; } /* USE SECANT RULE. */ *b += p / q; goto L7; /* USE BISECTION. */ L6: *b = (*c + *b) * (float).5; /* HAVE COMPLETED COMPUTATION FOR NEW ITERATE B. */ L7: *t = *b; *iflag = -3; return 0; L400: fb = *ft; if (fb == (float)0.) { goto L9; } ++kount; if (d_sign(&c_b17, &fb) != d_sign(&c_b17, &fc)) { goto L1; } *c = a; fc = fa; goto L1; /* FINISHED. SET IFLAG. */ L8: if (d_sign(&c_b17, &fb) == d_sign(&c_b17, &fc)) { goto L11; } if (dblabs(fb) > fx) { goto L10; } *iflag = 1; return 0; L9: *iflag = 2; return 0; L10: *iflag = 3; return 0; L11: *iflag = 4; return 0; L12: *iflag = 5; return 0; } /* root_ */ /* end root.c */ /**********************************************************************/ /***************** implementation code from rootnf.c ******************/ /**********************************************************************/ /* Table of constant values */ /* static int c__4 = 4; static int c__1 = 1; */ /* Subroutine */ int 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) { /* System generated locals */ int qr_dim1, qr_offset, i__1; double d__1; /* Local variables */ static double dels, aerr, rerr; static int judy; extern /* Subroutine */ int root_(doublereal *t, doublereal *ft, doublereal *b, doublereal *c, doublereal *relerr, doublereal *abserr, integer *iflag); static double sout; extern double dnrm2_(integer *n, doublereal *dx, integer *incx); static double u; static int lcode; extern double d1mach_(integer *); static double qsout, sa, sb; static int jw; /* extern */ /* Subroutine */ /* int tangnf_(); IN pelutils.h */ static int np1; /* ROOTNF FINDS THE POINT YBAR = (1, XBAR) ON THE ZERO CURVE OF THE */ /* HOMOTOPY MAP. IT STARTS WITH TWO POINTS YOLD=(LAMBDAOLD,XOLD) AND */ /* Y=(LAMBDA,X) SUCH THAT LAMBDAOLD < 1 <= LAMBDA , AND ALTERNATES */ /* BETWEEN HERMITE CUBIC INTERPOLATION AND NEWTON ITERATION UNTIL */ /* CONVERGENCE. */ /* ON INPUT: */ /* N = DIMENSION OF X AND THE HOMOTOPY MAP. */ /* NFE = NUMBER OF JACOBIAN MATRIX EVALUATIONS. */ /* IFLAG = -2, -1, OR 0, INDICATING THE PROBLEM TYPE. */ /* RELERR, ABSERR = RELATIVE AND ABSOLUTE ERROR VALUES. THE ITERATION IS */ /* CONSIDERED TO HAVE CONVERGED WHEN A POINT Y=(LAMBDA,X) IS FOUND */ /* SUCH THAT */ /* |Y(1) - 1| <= RELERR + ABSERR AND */ /* ||Z|| <= RELERR*||X|| + ABSERR , WHERE */ /* (?,Z) IS THE NEWTON STEP TO Y=(LAMBDA,X). */ /* Y(1:N+1) = POINT (LAMBDA(S), X(S)) ON ZERO CURVE OF HOMOTOPY MAP. */ /* YP(1:N+1) = UNIT TANGENT VECTOR TO THE ZERO CURVE OF THE HOMOTOPY MAP */ /* AT Y . */ /* YOLD(1:N+1) = A POINT DIFFERENT FROM Y ON THE ZERO CURVE. */ /* YPOLD(1:N+1) = UNIT TANGENT VECTOR TO THE ZERO CURVE OF THE HOMOTOPY */ /* MAP AT YOLD . */ /* A(1:*) = PARAMETER VECTOR IN THE HOMOTOPY MAP. */ /* QR(1:N,1:N+2), ALPHA(1:N), TZ(1:N+1), PIVOT(1:N+1), W(1:N+1), */ /* WP(1:N+1) ARE WORK ARRAYS USED FOR THE QR FACTORIZATION (IN THE */ /* NEWTON STEP CALCULATION) AND THE INTERPOLATION. */ /* PAR(1:*) AND IPAR(1:*) ARE ARRAYS FOR (OPTIONAL) USER PARAMETERS, */ /* WHICH ARE SIMPLY PASSED THROUGH TO THE USER WRITTEN SUBROUTINES */ /* RHO, RHOJAC. */ /* ON OUTPUT: */ /* N , RELERR , ABSERR , A ARE UNCHANGED. */ /* NFE HAS BEEN UPDATED. */ /* IFLAG */ /* = -2, -1, OR 0 (UNCHANGED) ON A NORMAL RETURN. */ /* = 4 IF A JACOBIAN MATRIX WITH RANK < N HAS OCCURRED. THE */ /* ITERATION WAS NOT COMPLETED. */ /* = 6 IF THE ITERATION FAILED TO CONVERGE. Y AND YOLD CONTAIN */ /* THE LAST TWO POINTS FOUND ON THE ZERO CURVE. */ /* Y IS THE POINT ON THE ZERO CURVE OF THE HOMOTOPY MAP AT LAMBDA = 1 . */ /* CALLS D1MACH , DNRM2 , ROOT , TANGNF . */ /* ***** ARRAY DECLARATIONS. ***** */ /* ***** END OF DIMENSIONAL INFORMATION. ***** */ /* THE LIMIT ON THE NUMBER OF ITERATIONS ALLOWED MAY BE CHANGED BY */ /* CHANGING THE FOLLOWING PARAMETER STATEMENT: */ /* DEFINITION OF HERMITE CUBIC INTERPOLANT VIA DIVIDED DIFFERENCES. */ /* Parameter adjustments */ --wp; --w; --pivot; --tz; --alpha; qr_dim1 = *n; qr_offset = qr_dim1 + 1; qr -= qr_offset; --a; --ypold; --yold; --yp; --y; --par; --ipar; /* Function Body */ u = d1mach_(&c__4); rerr = max(*relerr,u); aerr = max(*abserr,0.); np1 = *n + 1; /* ***** MAIN LOOP. ***** */ /* L100: */ for (judy = 1; judy <= 20; ++judy) { i__1 = np1; for (jw = 1; jw <= i__1; ++jw) { tz[jw] = y[jw] - yold[jw]; /* L110: */ } dels = dnrm2_((integer *)&np1, &tz[1], &c__1); /* USING TWO POINTS AND TANGENTS ON THE HOMOTOPY ZERO CURVE, CONSTRUCT */ /* THE HERMITE CUBIC INTERPOLANT Q(S). THEN USE ROOT TO FIND THE S */ /* CORRESPONDING TO LAMBDA = 1 . THE TWO POINTS ON THE ZERO CURVE AR E */ /* ALWAYS CHOSEN TO BRACKET LAMBDA=1, WITH THE BRACKETING INTERVAL */ /* ALWAYS BEING [0, DELS]. */ sa = (float)0.; sb = dels; lcode = 1; L130: root_(&sout, &qsout, &sa, &sb, &rerr, &aerr, (integer *)&lcode); if (lcode > 0) { goto L140; } qsout = ((((yp[1] - (y[1] - yold[1]) / dels) / dels - ((y[1] - yold[1] ) / dels - ypold[1]) / dels) / dels * (sout - dels) + ((y[1] - yold[1]) / dels - ypold[1]) / dels) * sout + ypold[1]) * sout + yold[1] - (float)1.; goto L130; /* IF LAMBDA = 1 WERE BRACKETED, ROOT CANNOT FAIL. */ L140: if (lcode > 2) { *iflag = 6; return 0; } /* CALCULATE Q(SA) AS THE INITIAL POINT FOR A NEWTON ITERATION. */ i__1 = np1; for (jw = 1; jw <= i__1; ++jw) { w[jw] = ((((yp[jw] - (y[jw] - yold[jw]) / dels) / dels - ((y[jw] - yold[jw]) / dels - ypold[jw]) / dels) / dels * (sa - dels) + ((y[jw] - yold[jw]) / dels - ypold[jw]) / dels) * sa + ypold[jw]) * sa + yold[jw]; /* L150: */ } /* CALCULATE NEWTON STEP AT Q(SA). */ tangnf_(&sa, &w[1], &wp[1], &ypold[1], &a[1], &qr[qr_offset], &alpha[ 1], &tz[1], &pivot[1], nfe, n, iflag, &par[1], &ipar[1]); if (*iflag > 0) { return 0; } /* NEXT POINT = CURRENT POINT + NEWTON STEP. */ i__1 = np1; for (jw = 1; jw <= i__1; ++jw) { w[jw] += tz[jw]; /* L160: */ } /* GET THE TANGENT WP AT W AND THE NEXT NEWTON STEP IN TZ . */ tangnf_(&sa, &w[1], &wp[1], &ypold[1], &a[1], &qr[qr_offset], &alpha[ 1], &tz[1], &pivot[1], nfe, n, iflag, &par[1], &ipar[1]); if (*iflag > 0) { return 0; } /* TAKE NEWTON STEP AND CHECK CONVERGENCE. */ i__1 = np1; for (jw = 1; jw <= i__1; ++jw) { w[jw] += tz[jw]; /* L170: */ } if ((d__1 = w[1] - (float)1., dblabs(d__1)) <= rerr + aerr && dnrm2_((integer *)n, & tz[2], &c__1) <= rerr * dnrm2_((integer*)n, &w[2], &c__1) + aerr) { i__1 = np1; for (jw = 1; jw <= i__1; ++jw) { y[jw] = w[jw]; /* L180: */ } return 0; } /* IF THE ITERATION HAS NOT CONVERGED, DISCARD ONE OF THE OLD POINTS */ /* SUCH THAT LAMBDA = 1 IS STILL BRACKETED. */ if ((yold[1] - (float)1.) * (w[1] - (float)1.) > (float)0.) { i__1 = np1; for (jw = 1; jw <= i__1; ++jw) { yold[jw] = w[jw]; ypold[jw] = wp[jw]; /* L200: */ } } else { i__1 = np1; for (jw = 1; jw <= i__1; ++jw) { y[jw] = w[jw]; yp[jw] = wp[jw]; /* L210: */ } } /* L300: */ } /* ***** END OF MAIN LOOP. ***** */ /* THE ALTERNATING OSCULATORY CUBIC INTERPOLATION AND NEWTON ITERATION */ /* HAS NOT CONVERGED IN LIMIT STEPS. ERROR RETURN. */ *iflag = 6; return 0; } /* rootnf_ */ /* end rootnf.c */ /**********************************************************************/ /******************* implementation code from stepnf.c *********************/ /**********************************************************************/ /* Table of constant values */ /* static integer c__4 = 4; static integer c__1 = 1; */ /* Subroutine */ int 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) { /* System generated locals */ integer qr_dim1, qr_offset, i__1; doublereal d__1, d__2, d__3, d__4; /* Builtin functions */ /* double sqrt(double), pow(double,double); CANT DECLARE BUILTINS UNDER C++ */ /* Local variables */ static logical fail; static doublereal temp; static integer judy; static doublereal twou; extern doublereal dnrm2_(integer *n, doublereal *dx, integer *incx); static doublereal dcalc; static integer j; static doublereal lcalc, hfail, rcalc; static integer itnum; extern doublereal d1mach_(integer *i); static doublereal fouru, ht; /*extern */ /* Subroutine */ /* int tangnf_(); IN pelutils.h */ static doublereal rholen; static integer np1; /* STEPNF TAKES ONE STEP ALONG THE ZERO CURVE OF THE HOMOTOPY MAP */ /* USING A PREDICTOR-CORRECTOR ALGORITHM. THE PREDICTOR USES A HERMITE */ /* CUBIC INTERPOLANT, AND THE CORRECTOR RETURNS TO THE ZERO CURVE ALONG */ /* THE FLOW NORMAL TO THE DAVIDENKO FLOW. STEPNF ALSO ESTIMATES A */ /* STEP SIZE H FOR THE NEXT STEP ALONG THE ZERO CURVE. NORMALLY */ /* STEPNF IS USED INDIRECTLY THROUGH FIXPNF , AND SHOULD BE CALLED */ /* DIRECTLY ONLY IF IT IS NECESSARY TO MODIFY THE STEPPING ALGORITHM'S */ /* PARAMETERS. */ /* ON INPUT: */ /* N = DIMENSION OF X AND THE HOMOTOPY MAP. */ /* NFE = NUMBER OF JACOBIAN MATRIX EVALUATIONS. */ /* IFLAG = -2, -1, OR 0, INDICATING THE PROBLEM TYPE. */ /* START = .TRUE. ON FIRST CALL TO STEPNF , .FALSE. OTHERWISE. */ /* HOLD = ||Y - YOLD||; SHOULD NOT BE MODIFIED BY THE USER. */ /* H = UPPER LIMIT ON LENGTH OF STEP THAT WILL BE ATTEMPTED. H MUST BE */ /* SET TO A POSITIVE NUMBER ON THE FIRST CALL TO STEPNF . */ /* THEREAFTER STEPNF CALCULATES AN OPTIMAL VALUE FOR H , AND H */ /* SHOULD NOT BE MODIFIED BY THE USER. */ /* RELERR, ABSERR = RELATIVE AND ABSOLUTE ERROR VALUES. THE ITERATION IS */ /* CONSIDERED TO HAVE CONVERGED WHEN A POINT W=(LAMBDA,X) IS FOUND */ /* SUCH THAT */ /* ||Z|| <= RELERR*||W|| + ABSERR , WHERE */ /* Z IS THE NEWTON STEP TO W=(LAMBDA,X). */ /* S = (APPROXIMATE) ARC LENGTH ALONG THE HOMOTOPY ZERO CURVE UP TO */ /* Y(S) = (LAMBDA(S), X(S)). */ /* Y(1:N+1) = PREVIOUS POINT (LAMBDA(S), X(S)) FOUND ON THE ZERO CURVE OF */ /* THE HOMOTOPY MAP. */ /* YP(1:N+1) = UNIT TANGENT VECTOR TO THE ZERO CURVE OF THE HOMOTOPY MAP */ /* AT Y . */ /* YOLD(1:N+1) = A POINT BEFORE Y ON THE ZERO CURVE OF THE HOMOTOPY MAP. */ /* YPOLD(1:N+1) = UNIT TANGENT VECTOR TO THE ZERO CURVE OF THE HOMOTOPY */ /* MAP AT YOLD . */ /* A(1:*) = PARAMETER VECTOR IN THE HOMOTOPY MAP. */ /* QR(1:N,1:N+2), ALPHA(1:N), TZ(1:N+1), PIVOT(1:N+1), W(1:N+1), */ /* WP(1:N+1) ARE WORK ARRAYS USED FOR THE QR FACTORIZATION (IN THE */ /* NEWTON STEP CALCULATION) AND THE INTERPOLATION. */ /* Z0(1:N+1), Z1(1:N+1) ARE WORK ARRAYS USED FOR THE ESTIMATION OF THE */ /* NEXT STEP SIZE H . */ /* SSPAR(1:8) = (LIDEAL, RIDEAL, DIDEAL, HMIN, HMAX, BMIN, BMAX, P) IS */ /* A VECTOR OF PARAMETERS USED FOR THE OPTIMAL STEP SIZE ESTIMATION. */ /* PAR(1:*) AND IPAR(1:*) ARE ARRAYS FOR (OPTIONAL) USER PARAMETERS, */ /* WHICH ARE SIMPLY PASSED THROUGH TO THE USER WRITTEN SUBROUTINES */ /* RHO, RHOJAC. */ /* ON OUTPUT: */ /* N , A , SSPAR ARE UNCHANGED. */ /* NFE HAS BEEN UPDATED. */ /* IFLAG */ /* = -2, -1, OR 0 (UNCHANGED) ON A NORMAL RETURN. */ /* = 4 IF A JACOBIAN MATRIX WITH RANK < N HAS OCCURRED. THE */ /* ITERATION WAS NOT COMPLETED. */ /* = 6 IF THE ITERATION FAILED TO CONVERGE. W CONTAINS THE LAST */ /* NEWTON ITERATE. */ /* START = .FALSE. ON A NORMAL RETURN. */ /* CRASH */ /* = .FALSE. ON A NORMAL RETURN. */ /* = .TRUE. IF THE STEP SIZE H WAS TOO SMALL. H HAS BEEN */ /* INCREASED TO AN ACCEPTABLE VALUE, WITH WHICH STEPNF MAY BE */ /* CALLED AGAIN. */ /* = .TRUE. IF RELERR AND/OR ABSERR WERE TOO SMALL. THEY HAVE */ /* BEEN INCREASED TO ACCEPTABLE VALUES, WITH WHICH STEPNF MAY */ /* BE CALLED AGAIN. */ /* HOLD = ||Y - YOLD||. */ /* H = OPTIMAL VALUE FOR NEXT STEP TO BE ATTEMPTED. NORMALLY H SHOULD */ /* NOT BE MODIFIED BY THE USER. */ /* RELERR, ABSERR ARE UNCHANGED ON A NORMAL RETURN. */ /* S = (APPROXIMATE) ARC LENGTH ALONG THE ZERO CURVE OF THE HOMOTOPY MAP */ /* UP TO THE LATEST POINT FOUND, WHICH IS RETURNED IN Y . */ /* Y, YP, YOLD, YPOLD CONTAIN THE TWO MOST RECENT POINTS AND TANGENT */ /* VECTORS FOUND ON THE ZERO CURVE OF THE HOMOTOPY MAP. */ /* CALLS D1MACH , DNRM2 , TANGNF . */ /* ***** ARRAY DECLARATIONS. ***** */ /* ***** END OF DIMENSIONAL INFORMATION. ***** */ /* THE LIMIT ON THE NUMBER OF NEWTON ITERATIONS ALLOWED BEFORE REDUCING */ /* THE STEP SIZE H MAY BE CHANGED BY CHANGING THE FOLLOWING PARAMETER */ /* STATEMENT: */ /* DEFINITION OF HERMITE CUBIC INTERPOLANT VIA DIVIDED DIFFERENCES. */ /* Parameter adjustments */ --z1; --z0; --wp; --w; --pivot; --tz; --alpha; qr_dim1 = *n; qr_offset = qr_dim1 + 1; qr -= qr_offset; --a; --ypold; --yold; --yp; --y; --sspar; --par; --ipar; /* Function Body */ twou = d1mach_(&c__4) * (float)2.; fouru = twou + twou; np1 = *n + 1; *crash = TRUE_; /* THE ARCLENGTH S MUST BE NONNEGATIVE. */ if (*s < (float)0.) { return 0; } /* IF STEP SIZE IS TOO SMALL, DETERMINE AN ACCEPTABLE ONE. */ if (*h < fouru * (*s + (float)1.)) { *h = fouru * (*s + (float)1.); return 0; } /* IF ERROR TOLERANCES ARE TOO SMALL, INCREASE THEM TO ACCEPTABLE VALUES. */ temp = dnrm2_(&np1, &y[1], &c__1); if ((*relerr * temp + *abserr) * (float).5 >= twou * temp) { goto L40; } if (*relerr != (float)0.) { *relerr = fouru * (fouru + (float)1.); *abserr = max(*abserr,0.); } else { *abserr = fouru * temp; } return 0; L40: *crash = FALSE_; if (! (*start)) { goto L300; } /* ***** STARTUP SECTION(FIRST STEP ALONG ZERO CURVE. ***** */ fail = FALSE_; *start = FALSE_; /* DETERMINE SUITABLE INITIAL STEP SIZE. */ /* Computing MIN */ d__1 = min(*h,.1), d__2 = sqrt(sqrt(*relerr * temp + *abserr)); *h = min(d__1,d__2); /* USE LINEAR PREDICTOR ALONG TANGENT DIRECTION TO START NEWTON ITERATION. */ ypold[1] = (float)1.; i__1 = np1; for (j = 2; j <= i__1; ++j) { ypold[j] = (float)0.; /* L50: */ } tangnf_(s, &y[1], &yp[1], &ypold[1], &a[1], &qr[qr_offset], &alpha[1], &tz[1], (int *)&pivot[1], (int *)nfe, (int *)n, (int *)iflag, &par[1], (int *)&ipar[1]); if (*iflag > 0) { return 0; } L70: i__1 = np1; for (j = 1; j <= i__1; ++j) { temp = y[j] + *h * yp[j]; w[j] = temp; z0[j] = temp; /* L80: */ } for (judy = 1; judy <= 4; ++judy) { rholen = (float)-1.; /* CALCULATE THE NEWTON STEP TZ AT THE CURRENT POINT W . */ tangnf_(&rholen, &w[1], &wp[1], &ypold[1], &a[1], &qr[qr_offset], &alpha[1], &tz[1], (int *)&pivot[1], (int *)nfe, (int *)n, (int *)iflag, &par[1], (int *)&ipar[1]) ; if (*iflag > 0) { return 0; } /* TAKE NEWTON STEP AND CHECK CONVERGENCE. */ i__1 = np1; for (j = 1; j <= i__1; ++j) { w[j] += tz[j]; /* L90: */ } itnum = judy; /* COMPUTE QUANTITIES USED FOR OPTIMAL STEP SIZE ESTIMATION. */ if (judy == 1) { lcalc = dnrm2_(&np1, &tz[1], &c__1); rcalc = rholen; i__1 = np1; for (j = 1; j <= i__1; ++j) { z1[j] = w[j]; /* L110: */ } } else if (judy == 2) { lcalc = dnrm2_(&np1, &tz[1], &c__1) / lcalc; rcalc = rholen / rcalc; } /* GO TO MOP-UP SECTION AFTER CONVERGENCE. */ if (dnrm2_(&np1, &tz[1], &c__1) <= *relerr * dnrm2_(&np1, &w[1], & c__1) + *abserr) { goto L600; } /* L200: */ } /* NO CONVERGENCE IN LITFH ITERATIONS. REDUCE H AND TRY AGAIN. */ if (*h <= fouru * (*s + (float)1.)) { *iflag = 6; return 0; } *h *= (float).5; goto L70; /* ***** END OF STARTUP SECTION. ***** */ /* ***** PREDICTOR SECTION. ***** */ L300: fail = FALSE_; /* COMPUTE POINT PREDICTED BY HERMITE INTERPOLANT. USE STEP SIZE H */ /* COMPUTED ON LAST CALL TO STEPNF . */ L320: i__1 = np1; for (j = 1; j <= i__1; ++j) { d__1 = *hold + *h; temp = ((((yp[j] - (y[j] - yold[j]) / *hold) / *hold - ((y[j] - yold[ j]) / *hold - ypold[j]) / *hold) / *hold * (d__1 - *hold) + (( y[j] - yold[j]) / *hold - ypold[j]) / *hold) * d__1 + ypold[j] ) * d__1 + yold[j]; w[j] = temp; z0[j] = temp; /* L330: */ } /* ***** END OF PREDICTOR SECTION. ***** */ /* ***** CORRECTOR SECTION. ***** */ for (judy = 1; judy <= 4; ++judy) { rholen = (float)-1.; /* CALCULATE THE NEWTON STEP TZ AT THE CURRENT POINT W . */ tangnf_(&rholen, &w[1], &wp[1], &yp[1], &a[1], &qr[qr_offset], &alpha[1], &tz[1], (int *)&pivot[1], (int *)nfe, (int *)n, (int *)iflag, &par[1], (int *)&ipar[1]); if (*iflag > 0) { return 0; } /* TAKE NEWTON STEP AND CHECK CONVERGENCE. */ i__1 = np1; for (j = 1; j <= i__1; ++j) { w[j] += tz[j]; /* L420: */ } itnum = judy; /* COMPUTE QUANTITIES USED FOR OPTIMAL STEP SIZE ESTIMATION. */ if (judy == 1) { lcalc = dnrm2_(&np1, &tz[1], &c__1); rcalc = rholen; i__1 = np1; for (j = 1; j <= i__1; ++j) { z1[j] = w[j]; /* L440: */ } } else if (judy == 2) { lcalc = dnrm2_(&np1, &tz[1], &c__1) / lcalc; rcalc = rholen / rcalc; } /* GO TO MOP-UP SECTION AFTER CONVERGENCE. */ if (dnrm2_(&np1, &tz[1], &c__1) <= *relerr * dnrm2_(&np1, &w[1], & c__1) + *abserr) { goto L600; } /* L500: */ } /*NO CONVERGENCE IN LITFH ITERATIONS. RECORD FAILURE AT CALCULATED H , */ /* SAVE THIS STEP SIZE, REDUCE H AND TRY AGAIN. */ fail = TRUE_; hfail = *h; if (*h <= fouru * (*s + (float)1.)) { *iflag = 6; return 0; } *h *= (float).5; goto L320; /* ***** END OF CORRECTOR SECTION. ***** */ /* ***** MOP-UP SECTION. ***** */ /* YOLD AND Y ALWAYS CONTAIN THE LAST TWO POINTS FOUND ON THE ZERO */ /* CURVE OF THE HOMOTOPY MAP. YPOLD AND YP CONTAIN THE TANGENT */ /* VECTORS TO THE ZERO CURVE AT YOLD AND Y , RESPECTIVELY. */ L600: i__1 = np1; for (j = 1; j <= i__1; ++j) { yold[j] = y[j]; ypold[j] = yp[j]; y[j] = w[j]; yp[j] = wp[j]; w[j] = y[j] - yold[j]; /* L620: */ } /* UPDATE ARC LENGTH. */ *hold = dnrm2_(&np1, &w[1], &c__1); *s += *hold; /* ***** END OF MOP-UP SECTION. ***** */ /* ***** OPTIMAL STEP SIZE ESTIMATION SECTION. ***** */ /* CALCULATE THE DISTANCE FACTOR DCALC . */ /* L700: */ i__1 = np1; for (j = 1; j <= i__1; ++j) { tz[j] = z0[j] - y[j]; w[j] = z1[j] - y[j]; /* L710: */ } dcalc = dnrm2_(&np1, &tz[1], &c__1); if (dcalc != (float)0.) { dcalc = dnrm2_(&np1, &w[1], &c__1) / dcalc; } /* THE OPTIMAL STEP SIZE HBAR IS DEFINED BY */ /* HT=HOLD * [MIN(LIDEAL/LCALC, RIDEAL/RCALC, DIDEAL/DCALC)]**(1/P) */ /* HBAR = MIN [ MAX(HT, BMIN*HOLD, HMIN), BMAX*HOLD, HMAX ] */ /* IF CONVERGENCE HAD OCCURRED AFTER 1 ITERATION, SET THE CONTRACTION */ /* FACTOR LCALC TO ZERO. */ if (itnum == 1) { lcalc = (float)0.; } /* FORMULA FOR OPTIMAL STEP SIZE. */ if (lcalc + rcalc + dcalc == (float)0.) { ht = sspar[7] * *hold; } else { /* Computing MAX */ d__2 = lcalc / sspar[1], d__3 = rcalc / sspar[2], d__2 = max(d__2, d__3), d__3 = dcalc / sspar[3]; d__1 = (float)1. / max(d__2,d__3); d__4 = (float)1. / sspar[8]; /* ht = pow_dd(&d__1, &d__4) * *hold; */ ht=pow(d__1,d__4) * *hold; } /* HT CONTAINS THE ESTIMATED OPTIMAL STEP SIZE. NOW PUT IT WITHIN */ /* REASONABLE BOUNDS. */ /* Computing MIN */ /* Computing MAX */ d__3 = ht, d__4 = sspar[6] * *hold, d__3 = max(d__3,d__4); d__1 = max(d__3,sspar[4]), d__2 = sspar[7] * *hold, d__1 = min(d__1,d__2); *h = min(d__1,sspar[5]); if (itnum == 1) { /* IF CONVERGENCE HAD OCCURRED AFTER 1 ITERATION, DON'T DECREASE H . */ *h = max(*h,*hold); } else if (itnum == 4) { /* IF CONVERGENCE REQUIRED THE MAXIMUM LITFH ITERATIONS, DON'T */ /* INCREASE H . */ *h = min(*h,*hold); } /* IF CONVERGENCE DID NOT OCCUR IN LITFH ITERATIONS FOR A PARTICULAR */ /* H = HFAIL , DON'T CHOOSE THE NEW STEP SIZE LARGER THAN HFAIL . */ if (fail) { *h = min(*h,hfail); } return 0; } /* stepnf_ */ /* end stepnf.c */ /**********************************************************************/ /***************** implementation code from Pmatrix.c *****************/ /**********************************************************************/ #define min(i,j) ((i) < (j) ? (i): (j)) #define FALSE 0 #define TRUE 1 /*--------------------------------------------------------------- vector/matrix type a linear array of int, whith auxilary info. *) the number of elements that can be stored is in elt[0] *) the current number of rows is in elt[1] *) the current number of collumbs is in elt[2] The actual data are then stored in row major order from elt[3] on ---------------------------------------------------------------*/ struct Pmatrix_t { int store; /* maximum number of helts reserved */ int topc; /* effective number of columbs */ int topr; /* effective number of rows */ int ncols; /* number of elts between first elts of rows */ polynomial1 *coords; }; /* ** a couple of functions which are defined in files that ** I don't want to include right now, when things stableize ** bad_error -- used mainly to stop things when a function ** is passed bad input, as things stablizes ** less catastrophic reactions may be introduced. ** mem_alloc, ** mem_free -- call malloc and free with some extra book-keeping. */ void bad_error(char *); /* WEIRD - The first pair of declarations below are the global one, which are now (6/13/99) included from pelutils.h. They replace the second pair, which had been declared at this point, specifically for this file, EVEN THOUGH THE ARGUMENT LIST IS DIFFERENT. I do not understand. - AMM */ /* void *mem_malloc(int); void mem_free(void *); */ /* char *mem_malloc(size_t); void mem_free(char *); */ /*------------------------------------------------------------- vector access macroes (which ignore any rows except for first) -------------------------------------------------------------*/ #define Vstore(V) ((V->store)) /* maximum #elts available */ #define Vlength(V) ((V->topr)) /* actual #elts stored */ #define Vref1(V,i) (&(((V->coords)[i-1]))) /*acces ith elt (starting at 1) */ #define Vref0(V,i) (&(((V->coords)[i]))) /*acces ith elt (starting at 0) */ #define Vref(V,i) Vref1(V,i) /*------------------------------------------------------------ matrix access macroes -------------------------------------------------------------*/ #define Mstore(V) ((V->store)) /* maximum #elts available */ #define MMrows(V) ((V->store/V->ncols)) /* maximum #rows */ #define Mrows(V) ((V->topr)) /* number rows stored */ #define Mcols(V) ((V->topc)) /* number cols stored */ #define MNcols(V) ((V->ncols)) /* number cols stored */ #define Mref1(V,i,j) (&(((V->coords)[(i-1)*(V->ncols)+j-1]))) #define Mref0(V,i,j) (&(((V->coords)[i*(V->ncols)+j]))) #define Mref(V,i,j) Mref1((V),i,j) int PMstore(Pmatrix M) { return Mstore(M); } int PMMrows(Pmatrix M) { return MMrows(M); } int PMrows(Pmatrix M) { return Mrows(M); } int PMcols(Pmatrix M) { return Mcols(M); } polynomial1 *PMref1(Pmatrix M, int i, int j) { return Mref1(M, i, j); } /* ** Constructor/Destructors for Pmatrixes ** ** Pmatrix Pmatrix_free(int r, int c); ** New Pmatrix cabable of holding r rows, and c collumbs. ** Pmatrix Pmatrix_new(Pmatrix V); */ Pmatrix Pmatrix_new(int r, int c) { Pmatrix V; int i, j; V = (Pmatrix) mem_malloc(sizeof(struct Pmatrix_t)); if (!V) bad_error("allocation failure in Pmatrix_new()"); V->coords = (polynomial1 *) mem_malloc(r * c * sizeof(polynomial1)); if (!V) bad_error("allocation failure 2 in Pmatrix_new()"); Mstore(V) = r * c; Mrows(V) = r; Mcols(V) = c; MNcols(V) = c; for (i = 1; i <= r; i++) { for (j = 1; j <= c; j++) *Mref(V, i, j)=0; } return V; } void Pmatrix_free(Pmatrix V) { int i, j; if (V != 0 && V->coords != 0) { for (i = 1; i <= Mrows(V); i++) for (j = 1; j <= Mcols(V); j++) freeP(*Mref(V, i, j)); mem_free((char *) (V->coords)); } if (V != 0) mem_free((char *) (V)); } /* ** Pmatrix_resize(R,r,c) ** Reset R to hold an r,by c matrix. ** if R has enough storage to hold an rxc matrix resets ** row and columb entrees of r to r and c. otherwise ** frees R and reallocates an rxc matrix ** DOES NOT PRESERVE INDECIES OF EXISTING DATA */ Pmatrix Pmatrix_resize(Pmatrix R, int r, int c) { if (R == 0 || Mstore(R) < (r * c)) { if (R != 0) Pmatrix_free(R); R = Pmatrix_new(r, c); } else { Mrows(R) = r; Mcols(R) = c; MNcols(R) = c; } return R; } Pmatrix Pmatrix_submat(Pmatrix R, int r, int c) { if (R == 0 || c > Mcols(R) || r > Mrows(R) * MNcols(R)) { bad_error("bad subscripts or zero matrix in Pmatrix_submat()"); } else { Mrows(R) = r; Mcols(R) = c; } return R; } /* ** Pmatrix_print(M): print an Pmatrix ** if M is null print <<>> and return fail. ** otherwise print matrix and return true. */ Pmatrix Pmatrix_print(Pmatrix M) { int i, j; if (M == 0) { printf("<<>>"); return 0; } printf("<"); for (i = 1; i <= Mrows(M); i++) { printf("<"); for (j = 1; j <= Mcols(M); j++) { printP(*Mref(M, i, j)); if (j < Mcols(M)) printf(", "); } printf(">"); if (i < Mrows(M)) printf("\n"); } printf(">"); return M; } Pmatrix Pmatrix_copy(Pmatrix P){ Pmatrix C; int i,j; C=Pmatrix_new(Mrows(P),Mcols(P)); for(i=1;i<=Mrows(P);i++){ for(j=1;j<=Mcols(P);j++){ *Mref(C,i,j)=copyP(*Mref(P,i,j)); } } return C; } /* end Pmatrix.c */ /**********************************************************************/ /********* implementation code from the original Homotopies.c *********/ /**********************************************************************/ /*------------------------------------------------------------------ Homotopies.c created 9/15/1994 last modified 9/15/994 Birk Huber (birk@math.cornell.edu ALL RIGHTS RESERVED Store and evaluate Homotopies for Continuation module of Pelican. the functions rho and rhojac are defined to be usable by the hompack code (originaly in fortran translated with the help of f2c). The representation of Homotopies and implementation is very similar to that used by hompack and described in Morgans book. -------------------------------------------------------------------*/ /*------------------------------------------------------------------ globals representing the homotopy and macroes for accessing them The system of complex polynomial1s is represented in real form -----------------------------------------------------------------*/ /* parameters affecting the homotopy */ int Hom_defd = 0; /* 0 no homotopy initialized , 1 else*/ int Hom_use_proj = 1; /* 0 dont use proj trans, 1 else*/ int Hom_num_vars = 0; /* number of complex vars in curr hom*/ /* private variables defining homotopy */ static int NV,N,N1,M; static int *Starting_Monomial; static int *Number_of_Monomials; static int *Exponents; static int *Hdegree; static int *Edegree; static double *Coefitients; static int *Deformation; static double *Proj_Trans; /* index in monomial list (starting at 0) of equation i*/ #define monst(i,j) ((Starting_Monomial[(i)-1])+(j)-1) /* number of monomials in the ith equation */ #define Nmon(i) (Number_of_Monomials[(i)-1]) /* exponent of variable k in monomial j of equation i */ #define Exp(i,j,k) Exponents[(monst(i,j))*NV+(k)-1] /* get the exponent of the deformation parameter in the jth monomial of the ith equation */ #define Def(i,j) Deformation[monst(i,j)] /* get the real and imaginary part of jth monomial of ith equation */ #define RCoef(i,j) Coefitients[2*(monst(i,j))] #define ICoef(i,j) Coefitients[2*(monst(i,j))+1] /* real and imaginary parts of coordinates defining the projective transformation */ #define RPtrans(j) Proj_Trans[2*(j)-2] #define IPtrans(j) Proj_Trans[2*(j)-1] /* exponent of homogenizing variable in jthmonomial of ith equation*/ #define Hdeg(i,j) Hdegree[monst(i,j)] /* degree of equation i */ #define Edeg(i) Edegree[i-1] extern Pring Def_Ring; Pvector psys_to_Pvec(psys sys){ polynomial1 tmpm,tmpp; int j; Pvector PV; PV=Pvector_new(psys_d(sys)); FORALL_POLY(sys, tmpp=0; FORALL_MONO(sys, tmpm=makeP(Def_Ring); *poly_coef(tmpm)=Complex(*psys_coef_real(sys), *psys_coef_imag(sys)); *poly_def(tmpm)=*psys_def(sys); *poly_homog(tmpm)=*psys_homog(sys); for(j=1;j<=psys_d(sys);j++) *poly_exp(tmpm,j)=*psys_exp(sys,j); tmpp=addPPP(tmpp,tmpm,tmpp); freeP(tmpm); ) *PVref(PV,psys_eqno(sys))=tmpp; ) return PV; } // WARNING: RDM added the following just to get to compile under BCC // I have no idea if this is correct!!! #if !defined(HAVE_DRAND48) double drand48() { return rand(); } #endif // HAVE_DRAND48 /*------------------------------------------------------------------- init_hom takes a Pvector and loads the above data structures to hold a representation of the system. -------------------------------------------------------------------*/ int init_hom(psys PS){ Pvector P; polynomial1 ptr; int i,j,k; double t; /* double drand48(); CANT DECLARE BUILTINS UNDER C++ */ int seed=12; P=psys_to_Pvec(PS); NV=poly_dim(*PMref(P,1,1)); /* number of complex variables*/ N=2*NV; /* number of real coordinates */ N1=N+1; /* number of realcoords incl lift coord*/ /* check that P has n elts, and that Ring has NV vars */ if ( NV != PMcols(P)) { printf("warning nonsquare homotopy in init_HPK\n"); Hom_defd=0; } /* Set number of vars for external use */ Hom_num_vars=NV; /* free all work space */ Dfree(0); Ifree(0); /* count monomials in equations to initialize monst and Nmon*/ Starting_Monomial=Ires(NV+1); Number_of_Monomials=Ires(NV); Starting_Monomial[0]=0; for(i=0;i1){ d=(*k)/2; /* determine the complex var to diff by.*/ for(i=1;i<=NV;i++){ v[2*i-2]=0.0; v[2*i-1]=0.0; if (Hom_use_proj==1){ for(j=1;j<=Nmon(i);j++){ if (Hdeg(i,j)!=0){ c=Cmul(Complex(RCoef(i,j),ICoef(i,j)), Cpow(Hpath(*lambda),Def(i,j))); for(h=1;h<=NV;h++) c=Cmul(c,Cpow(Complex(x[2*h-2],x[2*h-1]),Exp(i,j,h))); c=Cmul(c,Cpow(PN,Hdeg(i,j)-1)); c=Cmul(c,Complex(t=Hdeg(i,j),0.0)); c=Cmul(c,Complex(RPtrans(d),IPtrans(d))); v[2*i-2]+=c.r; v[2*i-1]+=c.i; } } } for(j=1;j<=Nmon(i);j++){ if (Exp(i,j,d)!=0){ c=Cmul(Complex(RCoef(i,j),ICoef(i,j)), Cpow(Hpath(*lambda),Def(i,j))); for(h=1;h= 0 THE NORM IS NOT */ /* COMPUTED AND RHOLEN IS NOT CHANGED. */ /* Y(1:N+1) = CURRENT POINT (LAMBDA(S), X(S)). */ /* YPOLD(1:N+1) = UNIT TANGENT VECTOR AT PREVIOUS POINT ON THE ZERO */ /* CURVE OF THE HOMOTOPY MAP. */ /* A(1:*) = PARAMETER VECTOR IN THE HOMOTOPY MAP. */ /* QR(1:N,1:N+2), ALPHA(1:N), TZ(1:N+1), PIVOT(1:N+1) ARE WORK ARRAYS */ /* USED FOR THE QR FACTORIZATION. */ /* NFE = NUMBER OF JACOBIAN MATRIX EVALUATIONS = NUMBER OF HOMOTOPY */ /* FUNCTION EVALUATIONS. */ /* N = DIMENSION OF X. */ /* IFLAG = -2, -1, OR 0, INDICATING THE PROBLEM TYPE. */ /* PAR(1:*) AND IPAR(1:*) ARE ARRAYS FOR (OPTIONAL) USER PARAMETERS, */ /* WHICH ARE SIMPLY PASSED THROUGH TO THE USER WRITTEN SUBROUTINES */ /* RHO, RHOJAC. */ /* ON OUTPUT: */ /* RHOLEN = ||RHO(A, LAMBDA(S), X(S)|| IF RHOLEN < 0 ON INPUT. */ /* OTHERWISE RHOLEN IS UNCHANGED. */ /* Y, YPOLD, A, N ARE UNCHANGED. */ /* YP(1:N+1) = DY/DS = UNIT TANGENT VECTOR TO INTEGRAL CURVE OF */ /* D(HOMOTOPY MAP)/DS = 0 AT Y(S) = (LAMBDA(S), X(S)) . */ /* TZ = THE NEWTON STEP = -(PSEUDO INVERSE OF (D RHO(A,Y(S))/D LAMBDA , */ /* D RHO(A,Y(S))/DX)) * RHO(A,Y(S)) . */ /* NFE HAS BEEN INCRMENTED BY 1. */ /* IFLAG IS UNCHANGED, UNLESS THE QR FACTORIZATION DETECTS A RANK < N, */ /* IN WHICH CASE THE TANGENT AND NEWTON STEP VECTORS ARE NOT COMPUTED */ /* AND TANGNF RETURNS WITH IFLAG = 4 . */ /* CALLS DDOT , DNRM2 , F (OR RHO ), FJAC (OR RHOJAC ). */ /* ***** ARRAY DECLARATIONS. ***** */ /* ARRAYS FOR COMPUTING THE JACOBIAN MATRIX AND ITS KERNEL. */ /* ***** END OF DIMENSIONAL INFORMATION. ***** */ /* Parameter adjustments */ --pivot; --tz; --alpha; qr_dim1 = *n; qr_offset = qr_dim1 + 1; qr -= qr_offset; --a; --ypold; --yp; --y; --par; --ipar; /* Function Body */ lambda = y[1]; np1 = *n + 1; np2 = *n + 2; ++(*nfe); /* NFE CONTAINS THE NUMBER OF JACOBIAN EVALUATIONS. */ /* * * * * * * * * * * * * * * * * * */ /* COMPUTE THE JACOBIAN MATRIX, STORE IT AND HOMOTOPY MAP IN QR. */ if (*iflag == -2) { /* QR = ( D RHO(A,LAMBDA,X)/D LAMBDA , D RHO(A,LAMBDA,X)/DX , */ /* RHO(A,LAMBDA,X) ) . */ i__1 = np1; for (k = 1; k <= i__1; ++k) { rhojac_(&a[1], &lambda, &y[2], &qr[k * qr_dim1 + 1], &k, &par[1], &ipar[1]); /* L30: */ } rho_(&a[1], &lambda, &y[2], &qr[np2 * qr_dim1 + 1], &par[1], &ipar[1]) ; } else { f_(&y[2], &tz[1]); if (*iflag == 0) { /* QR = ( A - F(X), I - LAMBDA*DF(X) , */ /* X - A + LAMBDA*(A - F(X)) ) . */ i__1 = *n; for (j = 1; j <= i__1; ++j) { sigma = a[j]; beta = sigma - tz[j]; qr[j + qr_dim1] = beta; /* L100: */ qr[j + np2 * qr_dim1] = y[j + 1] - sigma + lambda * beta; } i__1 = *n; for (k = 1; k <= i__1; ++k) { fjac_(&y[2], &tz[1], (integer *)&k); kp1 = k + 1; i__2 = *n; for (j = 1; j <= i__2; ++j) { /* L110: */ qr[j + kp1 * qr_dim1] = -lambda * tz[j]; } /* L120: */ qr[k + kp1 * qr_dim1] += (float)1.; } } else { /* QR = ( F(X) - X + A, LAMBDA*DF(X) + (1 - LAMBDA)*I , */ /* X - A + LAMBDA*(F(X) - X + A) ) . */ /* L140: */ i__1 = *n; for (j = 1; j <= i__1; ++j) { sigma = y[j + 1] - a[j]; beta = tz[j] - sigma; qr[j + qr_dim1] = beta; /* L150: */ qr[j + np2 * qr_dim1] = sigma + lambda * beta; } i__1 = *n; for (k = 1; k <= i__1; ++k) { fjac_(&y[2], &tz[1], (integer *)&k); kp1 = k + 1; i__2 = *n; for (j = 1; j <= i__2; ++j) { /* L160: */ qr[j + kp1 * qr_dim1] = lambda * tz[j]; } /* L170: */ qr[k + kp1 * qr_dim1] = (float)1. - lambda + qr[k + kp1 * qr_dim1]; } } } /* * * * * * * * * * * * * * * * * * */ /* COMPUTE THE NORM OF THE HOMOTOPY MAP IF IT WAS REQUESTED. */ if (*rholen < (float)0.) { *rholen = dnrm2_((integer *)n, &qr[np2 * qr_dim1 + 1], &c__1); } /* REDUCE THE JACOBIAN MATRIX TO UPPER TRIANGULAR FORM. */ /* THE FOLLOWING CODE IS A MODIFICATION OF THE ALGOL PROCEDURE */ /* DECOMPOSE IN P. BUSINGER AND G. H. GOLUB, LINEAR LEAST */ /* SQUARES SOLUTIONS BY HOUSEHOLDER TRANSFORMATIONS, */ /* NUMER. MATH. 7 (1965) 269-276. */ i__1 = np1; for (j = 1; j <= i__1; ++j) { yp[j] = ddot_((integer *)n, &qr[j * qr_dim1 + 1], &c__1, &qr[j * qr_dim1 + 1], &c__1); /* L220: */ pivot[j] = j; } i__1 = *n; for (k = 1; k <= i__1; ++k) { sigma = yp[k]; jbar = k; kp1 = k + 1; i__2 = np1; for (j = kp1; j <= i__2; ++j) { if (sigma >= yp[j]) { goto L240; } sigma = yp[j]; jbar = j; L240: ; } if (jbar == k) { goto L260; } i = pivot[k]; pivot[k] = pivot[jbar]; pivot[jbar] = i; yp[jbar] = yp[k]; yp[k] = sigma; i__2 = *n; for (i = 1; i <= i__2; ++i) { sigma = qr[i + k * qr_dim1]; qr[i + k * qr_dim1] = qr[i + jbar * qr_dim1]; qr[i + jbar * qr_dim1] = sigma; /* L250: */ } /* END OF COLUMN INTERCHANGE. */ L260: i__2 = *n - k + 1; sigma = ddot_((integer *)&i__2, &qr[k + k * qr_dim1], &c__1, &qr[k + k * qr_dim1], &c__1); if (sigma == (float)0.) { *iflag = 4; return 0; } /* L270: */ if (k == *n) { goto L300; } qrkk = qr[k + k * qr_dim1]; alphak = -sqrt(sigma); if (qrkk < (float)0.) { alphak = -alphak; } alpha[k] = alphak; beta = (float)1. / (sigma - qrkk * alphak); qr[k + k * qr_dim1] = qrkk - alphak; i__2 = np2; for (j = kp1; j <= i__2; ++j) { i__3 = *n - k + 1; sigma = beta * ddot_((integer *)&i__3, &qr[k + k * qr_dim1], &c__1, &qr[k + j * qr_dim1], &c__1); i__3 = *n; for (i = k; i <= i__3; ++i) { qr[i + j * qr_dim1] -= qr[i + k * qr_dim1] * sigma; /* L280: */ } if (j < np2) { /* Computing 2nd power */ d__1 = qr[k + j * qr_dim1]; yp[j] -= d__1 * d__1; } /* L290: */ } L300: ; } alpha[*n] = qr[*n + *n * qr_dim1]; /* COMPUTE KERNEL OF JACOBIAN, WHICH SPECIFIES YP=DY/DS. */ tz[np1] = (float)1.; for (i = *n; i >= 1; --i) { sum = (float)0.; i__1 = np1; for (j = i + 1; j <= i__1; ++j) { /* L330: */ sum += qr[i + j * qr_dim1] * tz[j]; } /* L340: */ tz[i] = -sum / alpha[i]; } ypnorm = dnrm2_((integer *)&np1, &tz[1], &c__1); i__1 = np1; for (k = 1; k <= i__1; ++k) { /* L360: */ yp[pivot[k]] = tz[k] / ypnorm; } if (ddot_((integer *)&np1, &yp[1], &c__1, &ypold[1], &c__1) >= (float)0.) { goto L380; } i__1 = np1; for (i = 1; i <= i__1; ++i) { /* L370: */ yp[i] = -yp[i]; } /* YP IS THE UNIT TANGENT VECTOR IN THE CORRECT DIRECTION. */ /* COMPUTE THE MINIMUM NORM SOLUTION OF [D RHO(Y(S))] V = -RHO(Y(S)). */ /* V IS GIVEN BY P - (P,Q)Q , WHERE P IS ANY SOLUTION OF */ /* [D RHO] V = -RHO AND Q IS A UNIT VECTOR IN THE KERNEL OF [D RHO]. */ L380: for (i = *n; i >= 1; --i) { sum = qr[i + np1 * qr_dim1] + qr[i + np2 * qr_dim1]; i__1 = *n; for (j = i + 1; j <= i__1; ++j) { /* L430: */ sum += qr[i + j * qr_dim1] * alpha[j]; } /* L440: */ alpha[i] = -sum / alpha[i]; } i__1 = *n; for (k = 1; k <= i__1; ++k) { /* L450: */ tz[pivot[k]] = alpha[k]; } tz[pivot[np1]] = (float)1.; /* TZ NOW CONTAINS A PARTICULAR SOLUTION P, AND YP CONTAINS A VECTOR Q */ /* IN THE KERNEL(THE TANGENT). */ sigma = ddot_((integer *)&np1, &tz[1], &c__1, &yp[1], &c__1); i__1 = np1; for (j = 1; j <= i__1; ++j) { tz[j] -= sigma * yp[j]; /* L470: */ } /* TZ IS THE NEWTON STEP FROM THE CURRENT POINT Y(S) = (LAMBDA(S), X(S)). */ return 0; } /* tangnf_ */ /* end tangnf.c */