COPASI API  4.16.103
Public Member Functions | Private Attributes | List of all members
CTruncatedNewton Class Reference

#include <CTruncatedNewton.h>

Collaboration diagram for CTruncatedNewton:
Collaboration graph
[legend]

Public Member Functions

 CTruncatedNewton ()
 
int getptc_ (C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_INT *)
 
int gtims_ (C_FLOAT64 *v, C_FLOAT64 *gv, C_INT *n, C_FLOAT64 *x, C_FLOAT64 *g, C_FLOAT64 *w, C_INT *, FTruncatedNewton *sfun, C_INT *first, C_FLOAT64 *delta, C_FLOAT64 *accrcy, C_FLOAT64 *xnorm)
 
int initpc_ (C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_FLOAT64 *, C_INT *, C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *)
 
int linder_ (C_INT *, FTruncatedNewton *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_INT *, C_FLOAT64 *, C_INT *)
 
int lmqn_ (C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, FTruncatedNewton *, C_INT *, C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *)
 
int lmqnbc_ (C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, FTruncatedNewton *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_INT *, C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *)
 
int modlnp_ (C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_FLOAT64 *, C_INT *, C_INT *, C_INT *, C_INT *, C_INT *, C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, FTruncatedNewton *, C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *)
 
int msolve_ (C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_FLOAT64 *, C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_INT *)
 
int setpar_ (C_INT *)
 
int setucr_ (C_FLOAT64 *, C_INT *, C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, FTruncatedNewton *, C_FLOAT64 *, C_FLOAT64 *)
 
int tn_ (C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, FTruncatedNewton *)
 
int tnbc_ (C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, FTruncatedNewton *, C_FLOAT64 *, C_FLOAT64 *, C_INT *)
 
 ~CTruncatedNewton ()
 

Private Attributes

subscr_mpsubscr_
 

Detailed Description

Definition at line 77 of file CTruncatedNewton.h.

Constructor & Destructor Documentation

CTruncatedNewton::CTruncatedNewton ( )

Definition at line 88 of file CTruncatedNewton.cpp.

References mpsubscr_.

89 {
90  mpsubscr_ = new subscr_;
91 }
CTruncatedNewton::~CTruncatedNewton ( )

Definition at line 93 of file CTruncatedNewton.cpp.

References mpsubscr_.

94 {
95  if (mpsubscr_ != NULL) {delete mpsubscr_; mpsubscr_ = NULL;}
96 }

Member Function Documentation

int CTruncatedNewton::getptc_ ( C_FLOAT64 big,
C_FLOAT64 ,
C_FLOAT64 rtsmll,
C_FLOAT64 reltol,
C_FLOAT64 fabstol,
C_FLOAT64 tnytol,
C_FLOAT64 fpresn,
C_FLOAT64 eta,
C_FLOAT64 rmu,
C_FLOAT64 xbnd,
C_FLOAT64 u,
C_FLOAT64 fu,
C_FLOAT64 gu,
C_FLOAT64 xmin,
C_FLOAT64 fmin,
C_FLOAT64 gmin,
C_FLOAT64 xw,
C_FLOAT64 fw,
C_FLOAT64 gw,
C_FLOAT64 a,
C_FLOAT64 b,
C_FLOAT64 oldf,
C_FLOAT64 b1,
C_FLOAT64 scxbnd,
C_FLOAT64 e,
C_FLOAT64 step,
C_FLOAT64 factor,
C_INT braktd,
C_FLOAT64 gtest1,
C_FLOAT64 gtest2,
C_FLOAT64 tol,
C_INT ientry,
C_INT itest 
)

Definition at line 2841 of file CTruncatedNewton.cpp.

References C_FLOAT64, C_INT, FALSE_, and TRUE_.

Referenced by linder_().

2850 {
2851  /* System generated locals */
2852  C_FLOAT64 d__1, d__2;
2853 
2854  /* Local variables */
2855  C_FLOAT64 half, abgw, fabsr, five, zero, p, q, r__, s, scale,
2856  denom, three, a1, d1, d2, sumsq, point1, abgmin, chordm, eleven,
2857  chordu;
2858  C_INT convrg;
2859  C_FLOAT64 xmidpt, twotol, one;
2860 
2861  /* ************************************************************ */
2862  /* GETPTC, AN ALGORITHM FOR FINDING A STEPLENGTH, CALLED REPEATEDLY BY */
2863  /* ROUTINES WHICH REQUIRE A STEP LENGTH TO BE COMPUTED USING CUBIC */
2864  /* INTERPOLATION. THE PARAMETERS CONTAIN INFORMATION ABOUT THE INTERVAL */
2865  /* IN WHICH A LOWER POINT IS TO BE FOUND AND FROM THIS GETPTC COMPUTES A
2866  */
2867  /* POINT AT WHICH THE FUNCTION CAN BE EVALUATED BY THE CALLING PROGRAM. */
2868  /* THE VALUE OF THE C_INT PARAMETERS IENTRY DETERMINES THE PATH TAKEN */
2869  /* THROUGH THE CODE. */
2870  /* ************************************************************ */
2871 
2872  /* THE FOLLOWING STANDARD FUNCTIONS AND SYSTEM FUNCTIONS ARE CALLED */
2873  /* WITHIN GETPTC */
2874 
2875  zero = 0.;
2876  point1 = .1;
2877  half = .5;
2878  one = 1.;
2879  three = 3.;
2880  five = 5.;
2881  eleven = 11.;
2882 
2883  /* BRANCH TO APPROPRIATE SECTION OF CODE DEPENDING ON THE */
2884  /* VALUE OF IENTRY. */
2885 
2886  switch (*ientry)
2887  {
2888  case 1: goto L10;
2889  case 2: goto L20;
2890  }
2891 
2892  /* IENTRY=1 */
2893  /* CHECK INPUT PARAMETERS */
2894 
2895 L10:
2896  *itest = 2;
2897 
2898  if (*u <= zero || *xbnd <= *tnytol || *gu > zero)
2899  {
2900  return 0;
2901  }
2902 
2903  *itest = 1;
2904 
2905  if (*xbnd < *fabstol)
2906  {
2907  *fabstol = *xbnd;
2908  }
2909 
2910  *tol = *fabstol;
2911  twotol = *tol + *tol;
2912 
2913  /* A AND B DEFINE THE INTERVAL OF UNCERTAINTY, X AND XW ARE POINTS */
2914  /* WITH LOWEST AND SECOND LOWEST FUNCTION VALUES SO FAR OBTAINED. */
2915  /* INITIALIZE A,SMIN,XW AT ORIGIN AND CORRESPONDING VALUES OF */
2916  /* FUNCTION AND PROJECTION OF THE GRADIENT ALONG DIRECTION OF SEARCH */
2917  /* AT VALUES FOR LATEST ESTIMATE AT MINIMUM. */
2918 
2919  *a = zero;
2920  *xw = zero;
2921  *xmin = zero;
2922  *oldf = *fu;
2923  *fmin = *fu;
2924  *fw = *fu;
2925  *gw = *gu;
2926  *gmin = *gu;
2927  *step = *u;
2928  *factor = five;
2929 
2930  /* THE MINIMUM HAS NOT YET BEEN BRACKETED. */
2931 
2932  *braktd = FALSE_;
2933 
2934  /* SET UP XBND AS A BOUND ON THE STEP TO BE TAKEN. (XBND IS NOT COMPUTED
2935  */
2936  /* EXPLICITLY BUT SCXBND IS ITS SCALED VALUE.) SET THE UPPER BOUND */
2937  /* ON THE INTERVAL OF UNCERTAINTY INITIALLY TO XBND + TOL(XBND). */
2938 
2939  *scxbnd = *xbnd;
2940  *b = *scxbnd + *reltol * fabs(*scxbnd) + *fabstol;
2941  *e = *b + *b;
2942  *b1 = *b;
2943 
2944  /* COMPUTE THE CONSTANTS REQUIRED FOR THE TWO CONVERGENCE CRITERIA. */
2945 
2946  *gtest1 = -(*rmu) * *gu;
2947  *gtest2 = -(*eta) * *gu;
2948 
2949  /* SET IENTRY TO INDICATE THAT THIS IS THE FIRST ITERATION */
2950 
2951  *ientry = 2;
2952  goto L210;
2953 
2954  /* IENTRY = 2 */
2955 
2956  /* UPDATE A,B,XW, AND XMIN */
2957 
2958 L20:
2959 
2960  if (*fu > *fmin)
2961  {
2962  goto L60;
2963  }
2964 
2965  /* IF FUNCTION VALUE NOT INCREASED, NEW POINT BECOMES NEXT */
2966  /* ORIGIN AND OTHER POINTS ARE SCALED ACCORDINGLY. */
2967 
2968  chordu = *oldf - (*xmin + *u) * *gtest1;
2969 
2970  if (*fu <= chordu)
2971  {
2972  goto L30;
2973  }
2974 
2975  /* THE NEW FUNCTION VALUE DOES NOT SATISFY THE SUFFICIENT DECREASE */
2976  /* CRITERION. PREPARE TO MOVE THE UPPER BOUND TO THIS POINT AND */
2977  /* FORCE THE INTERPOLATION SCHEME TO EITHER BISECT THE INTERVAL OF */
2978  /* UNCERTAINTY OR TAKE THE LINEAR INTERPOLATION STEP WHICH ESTIMATES */
2979  /* THE ROOT OF F(ALPHA)=CHORD(ALPHA). */
2980 
2981  chordm = *oldf - *xmin * *gtest1;
2982  *gu = -(*gmin);
2983  denom = chordm - *fmin;
2984 
2985  if (fabs(denom) >= 1e-15)
2986  {
2987  goto L25;
2988  }
2989 
2990  denom = 1e-15;
2991 
2992  if (chordm - *fmin < 0.)
2993  {
2994  denom = -denom;
2995  }
2996 
2997 L25:
2998 
2999  if (*xmin != zero)
3000  {
3001  *gu = *gmin * (chordu - *fu) / denom;
3002  }
3003 
3004  *fu = half * *u * (*gmin + *gu) + *fmin;
3005 
3006  if (*fu < *fmin)
3007  {
3008  *fu = *fmin;
3009  }
3010 
3011  goto L60;
3012 L30:
3013  *fw = *fmin;
3014  *fmin = *fu;
3015  *gw = *gmin;
3016  *gmin = *gu;
3017  *xmin += *u;
3018  *a -= *u;
3019  *b -= *u;
3020  *xw = -(*u);
3021  *scxbnd -= *u;
3022 
3023  if (*gu <= zero)
3024  {
3025  goto L40;
3026  }
3027 
3028  *b = zero;
3029  *braktd = TRUE_;
3030  goto L50;
3031 L40:
3032  *a = zero;
3033 L50:
3034  *tol = fabs(*xmin) * *reltol + *fabstol;
3035  goto L90;
3036 
3037  /* IF FUNCTION VALUE INCREASED, ORIGIN REMAINS UNCHANGED */
3038  /* BUT NEW POINT MAY NOW QUALIFY AS W. */
3039 
3040 L60:
3041 
3042  if (*u < zero)
3043  {
3044  goto L70;
3045  }
3046 
3047  *b = *u;
3048  *braktd = TRUE_;
3049  goto L80;
3050 L70:
3051  *a = *u;
3052 L80:
3053  *xw = *u;
3054  *fw = *fu;
3055  *gw = *gu;
3056 L90:
3057  twotol = *tol + *tol;
3058  xmidpt = half * (*a + *b);
3059 
3060  /* CHECK TERMINATION CRITERIA */
3061 
3062  convrg = fabs(xmidpt) <= twotol - half * (*b - *a) ||
3063  (fabs(*gmin) <= *gtest2 &&
3064  *fmin < *oldf &&
3065  ((d__1 = *xmin - *xbnd, fabs(d__1)) > * tol ||
3066  !(*braktd)));
3067 
3068  if (! convrg)
3069  {
3070  goto L100;
3071  }
3072 
3073  *itest = 0;
3074 
3075  if (*xmin != zero)
3076  {
3077  return 0;
3078  }
3079 
3080  /* IF THE FUNCTION HAS NOT BEEN REDUCED, CHECK TO SEE THAT THE RELATIVE */
3081  /* CHANGE IN F(X) IS CONSISTENT WITH THE ESTIMATE OF THE DELTA- */
3082  /* UNIMODALITY CONSTANT, TOL. IF THE CHANGE IN F(X) IS LARGER THAN */
3083  /* EXPECTED, REDUCE THE VALUE OF TOL. */
3084 
3085  *itest = 3;
3086 
3087  if ((d__1 = *oldf - *fw, fabs(d__1)) <= *fpresn * (one + fabs(*oldf)))
3088  {
3089  return 0;
3090  }
3091 
3092  *tol = point1 * *tol;
3093 
3094  if (*tol < *tnytol)
3095  {
3096  return 0;
3097  }
3098 
3099  *reltol = point1 * *reltol;
3100  *fabstol = point1 * *fabstol;
3101  twotol = point1 * twotol;
3102 
3103  /* CONTINUE WITH THE COMPUTATION OF A TRIAL STEP LENGTH */
3104 
3105 L100:
3106  r__ = zero;
3107  q = zero;
3108  s = zero;
3109 
3110  if (fabs(*e) <= *tol)
3111  {
3112  goto L150;
3113  }
3114 
3115  /* FIT CUBIC THROUGH XMIN AND XW */
3116 
3117  r__ = three * (*fmin - *fw) / *xw + *gmin + *gw;
3118  fabsr = fabs(r__);
3119  q = fabsr;
3120 
3121  if (*gw == zero || *gmin == zero)
3122  {
3123  goto L140;
3124  }
3125 
3126  /* COMPUTE THE SQUARE ROOT OF (R*R - GMIN*GW) IN A WAY */
3127  /* WHICH AVOIDS UNDERFLOW AND OVERFLOW. */
3128 
3129  abgw = fabs(*gw);
3130  abgmin = fabs(*gmin);
3131  s = sqrt(abgmin) * sqrt(abgw);
3132 
3133  if (*gw / abgw * *gmin > zero)
3134  {
3135  goto L130;
3136  }
3137 
3138  /* COMPUTE THE SQUARE ROOT OF R*R + S*S. */
3139 
3140  sumsq = one;
3141  p = zero;
3142 
3143  if (fabsr >= s)
3144  {
3145  goto L110;
3146  }
3147 
3148  /* THERE IS A POSSIBILITY OF OVERFLOW. */
3149 
3150  if (s > *rtsmll)
3151  {
3152  p = s * *rtsmll;
3153  }
3154 
3155  if (fabsr >= p)
3156  {
3157  /* Computing 2nd power */
3158  d__1 = fabsr / s;
3159  sumsq = one + d__1 * d__1;
3160  }
3161 
3162  scale = s;
3163  goto L120;
3164 
3165  /* THERE IS A POSSIBILITY OF UNDERFLOW. */
3166 
3167 L110:
3168 
3169  if (fabsr > *rtsmll)
3170  {
3171  p = fabsr * *rtsmll;
3172  }
3173 
3174  if (s >= p)
3175  {
3176  /* Computing 2nd power */
3177  d__1 = s / fabsr;
3178  sumsq = one + d__1 * d__1;
3179  }
3180 
3181  scale = fabsr;
3182 L120:
3183  sumsq = sqrt(sumsq);
3184  q = *big;
3185 
3186  if (scale < *big / sumsq)
3187  {
3188  q = scale * sumsq;
3189  }
3190 
3191  goto L140;
3192 
3193  /* COMPUTE THE SQUARE ROOT OF R*R - S*S */
3194 
3195 L130:
3196  q = sqrt((d__1 = r__ + s, fabs(d__1))) * sqrt((d__2 = r__ - s, fabs(d__2)));
3197 
3198  if (r__ >= s || r__ <= -s)
3199  {
3200  goto L140;
3201  }
3202 
3203  r__ = zero;
3204  q = zero;
3205  goto L150;
3206 
3207  /* COMPUTE THE MINIMUM OF FITTED CUBIC */
3208 
3209 L140:
3210 
3211  if (*xw < zero)
3212  {
3213  q = -q;
3214  }
3215 
3216  s = *xw * (*gmin - r__ - q);
3217  q = *gw - *gmin + q + q;
3218 
3219  if (q > zero)
3220  {
3221  s = -s;
3222  }
3223 
3224  if (q <= zero)
3225  {
3226  q = -q;
3227  }
3228 
3229  r__ = *e;
3230 
3231  if (*b1 != *step || *braktd)
3232  {
3233  *e = *step;
3234  }
3235 
3236  /* CONSTRUCT AN ARTIFICIAL BOUND ON THE ESTIMATED STEPLENGTH */
3237 
3238 L150:
3239  a1 = *a;
3240  *b1 = *b;
3241  *step = xmidpt;
3242 
3243  if (*braktd)
3244  {
3245  goto L160;
3246  }
3247 
3248  *step = -(*factor) * *xw;
3249 
3250  if (*step > *scxbnd)
3251  {
3252  *step = *scxbnd;
3253  }
3254 
3255  if (*step != *scxbnd)
3256  {
3257  *factor = five * *factor;
3258  }
3259 
3260  goto L170;
3261 
3262  /* IF THE MINIMUM IS BRACKETED BY 0 AND XW THE STEP MUST LIE */
3263  /* WITHIN (A,B). */
3264 
3265 L160:
3266 
3267  if ((*a != zero || *xw >= zero) && (*b != zero || *xw <= zero))
3268  {
3269  goto L180;
3270  }
3271 
3272  /* IF THE MINIMUM IS NOT BRACKETED BY 0 AND XW THE STEP MUST LIE */
3273  /* WITHIN (A1,B1). */
3274 
3275  d1 = *xw;
3276  d2 = *a;
3277 
3278  if (*a == zero)
3279  {
3280  d2 = *b;
3281  }
3282 
3283  /* THIS LINE MIGHT BE */
3284  /* IF (A .EQ. ZERO) D2 = E */
3285  *u = -d1 / d2;
3286  *step = five * d2 * (point1 + one / *u) / eleven;
3287 
3288  if (*u < one)
3289  {
3290  *step = half * d2 * sqrt(*u);
3291  }
3292 
3293 L170:
3294 
3295  if (*step <= zero)
3296  {
3297  a1 = *step;
3298  }
3299 
3300  if (*step > zero)
3301  {
3302  *b1 = *step;
3303  }
3304 
3305  /* REJECT THE STEP OBTAINED BY INTERPOLATION IF IT LIES OUTSIDE THE */
3306  /* REQUIRED INTERVAL OR IT IS GREATER THAN HALF THE STEP OBTAINED */
3307  /* DURING THE LAST-BUT-ONE ITERATION. */
3308 
3309 L180:
3310 
3311  if (fabs(s) <= (d__1 = half * q * r__, fabs(d__1)) || s <= q * a1 || s >= q
3312  * *b1)
3313  {
3314  goto L200;
3315  }
3316 
3317  /* A CUBIC INTERPOLATION STEP */
3318 
3319  *step = s / q;
3320 
3321  /* THE FUNCTION MUST NOT BE EVALUTATED TOO CLOSE TO A OR B. */
3322 
3323  if (*step - *a >= twotol && *b - *step >= twotol)
3324  {
3325  goto L210;
3326  }
3327 
3328  if (xmidpt > zero)
3329  {
3330  goto L190;
3331  }
3332 
3333  *step = -(*tol);
3334  goto L210;
3335 L190:
3336  *step = *tol;
3337  goto L210;
3338 L200:
3339  *e = *b - *a;
3340 
3341  /* IF THE STEP IS TOO LARGE, REPLACE BY THE SCALED BOUND (SO AS TO */
3342  /* COMPUTE THE NEW POINT ON THE BOUNDARY). */
3343 
3344 L210:
3345 
3346  if (*step < *scxbnd)
3347  {
3348  goto L220;
3349  }
3350 
3351  *step = *scxbnd;
3352 
3353  /* MOVE SXBD TO THE LEFT SO THAT SBND + TOL(XBND) = XBND. */
3354 
3355  *scxbnd -= (*reltol * fabs(*xbnd) + *fabstol) / (one + *reltol);
3356 L220:
3357  *u = *step;
3358 
3359  if (fabs(*step) < *tol && *step < zero)
3360  {
3361  *u = -(*tol);
3362  }
3363 
3364  if (fabs(*step) < *tol && *step >= zero)
3365  {
3366  *u = *tol;
3367  }
3368 
3369  *itest = 1;
3370  return 0;
3371 } /* getptc_ */
#define C_INT
Definition: copasi.h:115
#define TRUE_
#define C_FLOAT64
Definition: copasi.h:92
#define FALSE_
int CTruncatedNewton::gtims_ ( C_FLOAT64 v,
C_FLOAT64 gv,
C_INT n,
C_FLOAT64 x,
C_FLOAT64 g,
C_FLOAT64 w,
C_INT ,
FTruncatedNewton sfun,
C_INT first,
C_FLOAT64 delta,
C_FLOAT64 accrcy,
C_FLOAT64 xnorm 
)

Definition at line 2266 of file CTruncatedNewton.cpp.

References C_FLOAT64, C_INT, FALSE_, and subscr_2.

Referenced by modlnp_().

2270 {
2271  /* System generated locals */
2272  C_INT i__1;
2273 
2274  /* Local variables */
2275  C_FLOAT64 dinv, f;
2276  C_INT i__, ihg;
2277 
2278  /* THIS ROUTINE COMPUTES THE PRODUCT OF THE MATRIX G TIMES THE VECTOR */
2279  /* V AND STORES THE RESULT IN THE VECTOR GV (FINITE-DIFFERENCE VERSION) */
2280 
2281  /* Parameter adjustments */
2282  --g;
2283  --x;
2284  --gv;
2285  --v;
2286  --w;
2287 
2288  /* Function Body */
2289  if (!(*first))
2290  {
2291  goto L20;
2292  }
2293 
2294  *delta = sqrt(*accrcy) * (*xnorm + 1.);
2295  *first = FALSE_;
2296 L20:
2297  dinv = 1. / *delta;
2298  ihg = subscr_2.lhg;
2299  i__1 = *n;
2300 
2301  for (i__ = 1; i__ <= i__1; ++i__)
2302  {
2303  w[ihg] = x[i__] + *delta * v[i__];
2304  ++ihg;
2305  /* L30: */
2306  }
2307 
2308  (*sfun)(n, &w[subscr_2.lhg], &f, &gv[1]);
2309  i__1 = *n;
2310 
2311  for (i__ = 1; i__ <= i__1; ++i__)
2312  {
2313  gv[i__] = (gv[i__] - g[i__]) * dinv;
2314  /* L40: */
2315  }
2316 
2317  return 0;
2318 } /* gtims_ */
#define C_INT
Definition: copasi.h:115
#define subscr_2
#define C_FLOAT64
Definition: copasi.h:92
#define FALSE_
int CTruncatedNewton::initpc_ ( C_FLOAT64 diagb,
C_FLOAT64 emat,
C_INT n,
C_FLOAT64 w,
C_INT ,
C_INT modet,
C_INT upd1,
C_FLOAT64 yksk,
C_FLOAT64 ,
C_FLOAT64 yrsr,
C_INT lreset 
)

Definition at line 2519 of file CTruncatedNewton.cpp.

References initp3_(), and subscr_2.

Referenced by modlnp_().

2522 {
2523  /* Parameter adjustments */
2524  --emat;
2525  --diagb;
2526  --w;
2527 
2528  /* Function Body */
2529  initp3_(&diagb[1], &emat[1], n, lreset, yksk, yrsr, &w[subscr_2.lhyk], &w[
2530  subscr_2.lsk], &w[subscr_2.lyk], &w[subscr_2.lsr], &w[
2531  subscr_2.lyr], modet, upd1);
2532  return 0;
2533 } /* initpc_ */
int initp3_(C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_INT *)
#define subscr_2
int CTruncatedNewton::linder_ ( C_INT n,
FTruncatedNewton sfun,
C_FLOAT64 small,
C_FLOAT64 epsmch,
C_FLOAT64 reltol,
C_FLOAT64 fabstol,
C_FLOAT64 tnytol,
C_FLOAT64 eta,
C_FLOAT64 ,
C_FLOAT64 xbnd,
C_FLOAT64 p,
C_FLOAT64 gtp,
C_FLOAT64 x,
C_FLOAT64 f,
C_FLOAT64 alpha,
C_FLOAT64 g,
C_INT nftotl,
C_INT iflag,
C_FLOAT64 w,
C_INT  
)

Definition at line 2689 of file CTruncatedNewton.cpp.

References c__1, C_FLOAT64, C_INT, dcopy_(), ddot_(), getptc_(), and lsout_().

Referenced by lmqn_(), and lmqnbc_().

2695 {
2696  /* System generated locals */
2697  C_INT i__1;
2698 
2699  /* Local variables */
2700  C_FLOAT64 oldf, fmin, gmin;
2701  C_INT numf;
2702  C_FLOAT64 step, xmin, a, b, e;
2703  C_INT i__, l;
2704  C_FLOAT64 u;
2705  C_INT itcnt;
2706  C_FLOAT64 b1;
2707  C_INT itest, nprnt;
2708  C_FLOAT64 gtest1, gtest2;
2709  C_INT lg;
2710  C_FLOAT64 fu, gu, fw, gw;
2711  C_INT lx;
2712  C_INT braktd;
2713  C_FLOAT64 ualpha, factor, scxbnd, xw;
2714  C_FLOAT64 fpresn;
2715  C_INT ientry;
2716  C_FLOAT64 rtsmll;
2717  C_INT lsprnt;
2718  C_FLOAT64 big, tol, rmu;
2719 
2720  /* THE FOLLOWING STANDARD FUNCTIONS AND SYSTEM FUNCTIONS ARE */
2721  /* CALLED WITHIN LINDER */
2722 
2723  /* ALLOCATE THE ADDRESSES FOR LOCAL WORKSPACE */
2724 
2725  /* Parameter adjustments */
2726  --g;
2727  --x;
2728  --p;
2729  --w;
2730 
2731  /* Function Body */
2732  lx = 1;
2733  lg = lx + *n;
2734  lsprnt = 0;
2735  nprnt = 10000;
2736  rtsmll = sqrt(*small);
2737  big = 1. / *small;
2738  itcnt = 0;
2739 
2740  /* SET THE ESTIMATED RELATIVE PRECISION IN F(X). */
2741 
2742  fpresn = *epsmch * 10.;
2743  numf = 0;
2744  u = *alpha;
2745  fu = *f;
2746  fmin = *f;
2747  gu = *gtp;
2748  rmu = 1e-4;
2749 
2750  /* FIRST ENTRY SETS UP THE INITIAL INTERVAL OF UNCERTAINTY. */
2751 
2752  ientry = 1;
2753 L10:
2754 
2755  /* TEST FOR TOO MANY ITERATIONS */
2756 
2757  ++itcnt;
2758  *iflag = 1;
2759 
2760  if (itcnt > 20)
2761  {
2762  goto L50;
2763  }
2764 
2765  *iflag = 0;
2766  CTruncatedNewton::getptc_(&big, small, &rtsmll, reltol, fabstol, tnytol, &fpresn, eta, &rmu,
2767  xbnd, &u, &fu, &gu, &xmin, &fmin, &gmin, &xw, &fw, &gw, &a, &b, &
2768  oldf, &b1, &scxbnd, &e, &step, &factor, &braktd, &gtest1, &gtest2,
2769  &tol, &ientry, &itest);
2770 
2771  /* LSOUT */
2772  if (lsprnt >= nprnt)
2773  {
2774  lsout_(&ientry, &itest, &xmin, &fmin, &gmin, &xw, &fw, &gw, &u, &a, &
2775  b, &tol, reltol, &scxbnd, xbnd);
2776  }
2777 
2778  /* IF ITEST=1, THE ALGORITHM REQUIRES THE FUNCTION VALUE TO BE */
2779  /* CALCULATED. */
2780 
2781  if (itest != 1)
2782  {
2783  goto L30;
2784  }
2785 
2786  ualpha = xmin + u;
2787  l = lx;
2788  i__1 = *n;
2789 
2790  for (i__ = 1; i__ <= i__1; ++i__)
2791  {
2792  w[l] = x[i__] + ualpha * p[i__];
2793  ++l;
2794  /* L20: */
2795  }
2796 
2797  (*sfun)(n, &w[lx], &fu, &w[lg]);
2798  ++numf;
2799  gu = ddot_(n, &w[lg], &c__1, &p[1], &c__1);
2800 
2801  /* THE GRADIENT VECTOR CORRESPONDING TO THE BEST POINT IS */
2802  /* OVERWRITTEN IF FU IS LESS THAN FMIN AND FU IS SUFFICIENTLY */
2803  /* LOWER THAN F AT THE ORIGIN. */
2804 
2805  if (fu <= fmin && fu <= oldf - ualpha * gtest1)
2806  {
2807  dcopy_(n, &w[lg], &c__1, &g[1], &c__1);
2808  }
2809 
2810  goto L10;
2811 
2812  /* IF ITEST=2 OR 3 A LOWER POINT COULD NOT BE FOUND */
2813 
2814 L30:
2815  *nftotl = numf;
2816  *iflag = 1;
2817 
2818  if (itest != 0)
2819  {
2820  goto L50;
2821  }
2822 
2823  /* IF ITEST=0 A SUCCESSFUL SEARCH HAS BEEN MADE */
2824 
2825  *iflag = 0;
2826  *f = fmin;
2827  *alpha = xmin;
2828  i__1 = *n;
2829 
2830  for (i__ = 1; i__ <= i__1; ++i__)
2831  {
2832  x[i__] += *alpha * p[i__];
2833  /* L40: */
2834  }
2835 
2836 L50:
2837  return 0;
2838 } /* linder_ */
#define C_INT
Definition: copasi.h:115
int lsout_(C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *)
int dcopy_(integer *n, doublereal *dx, integer *incx, doublereal *dy, integer *incy)
int getptc_(C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_INT *)
static C_INT c__1
doublereal ddot_(integer *n, doublereal *dx, integer *incx, doublereal *dy, integer *incy)
#define C_FLOAT64
Definition: copasi.h:92
int CTruncatedNewton::lmqn_ ( C_INT ifail,
C_INT n,
C_FLOAT64 x,
C_FLOAT64 f,
C_FLOAT64 g,
C_FLOAT64 w,
C_INT lw,
FTruncatedNewton sfun,
C_INT msglvl,
C_INT maxit,
C_INT maxfun,
C_FLOAT64 eta,
C_FLOAT64 stepmx,
C_FLOAT64 accrcy,
C_FLOAT64 xtol 
)

Definition at line 425 of file CTruncatedNewton.cpp.

References c__1, c_false, C_FLOAT64, C_INT, chkucp_(), dcopy_(), ddot_(), dnrm2_(), dxpy_(), FALSE_, linder_(), modlnp_(), setpar_(), setucr_(), step1_(), subscr_1, and TRUE_.

Referenced by tn_().

429 {
430  /* Format strings */
431  /* char fmt_800[] = "(//\002 NIT NF CG\002,9x,\002F\002,21x,\002"
432  "GTG\002,//)";
433  char fmt_810[] = "(\002 \002,i3,1x,i4,1x,i4,1x,1pd22.15,2x,1pd15."
434  "8)"; */
435 
436  /* System generated locals */
437  C_INT i__1;
438 
439  /* Builtin functions */
440  // C_INT s_wsfe(cilist *), e_wsfe(void), do_fio(C_INT *, char *, ftnlen);
441 
442  /* Local variables */
443  C_FLOAT64 fold, oldf;
444  C_FLOAT64 fnew;
445  C_INT numf;
446  C_FLOAT64 peps;
447  C_INT lhyr;
448  C_FLOAT64 zero, rtol, yksk, tiny;
449  C_INT nwhy;
450  C_FLOAT64 yrsr;
451  C_INT i__;
452  C_FLOAT64 alpha, fkeep;
453  C_INT ioldg;
454  C_FLOAT64 small;
455  C_INT modet;
456  C_INT niter;
457  C_FLOAT64 gnorm, ftest, fstop, pnorm, rteps, xnorm;
458  C_INT idiagb;
459  C_FLOAT64 fm, pe, difold;
460  C_INT icycle, nlincg, nfeval;
461  C_FLOAT64 difnew;
462  C_INT nmodif;
463  C_FLOAT64 epsmch;
464  C_FLOAT64 epsred, fabstol, oldgtp;
465  C_INT ireset;
466  C_INT lreset;
467  C_FLOAT64 reltol, gtpnew;
468  C_INT nftotl;
469  C_FLOAT64 toleps;
470  C_FLOAT64 rtleps;
471  C_INT ipivot, nm1;
472  C_FLOAT64 rtolsq, tnytol, gtg, one;
473  C_INT ipk;
474  C_FLOAT64 gsk, spe;
475  C_INT isk, iyk;
476  C_INT upd1;
477 
478  /* Fortran I/O blocks */
479  //remove the blocks
480  /* cilist io___27 = {0, 6, 0, fmt_800, 0 };
481  cilist io___56 = {0, 6, 0, fmt_810, 0 };
482  cilist io___81 = {0, 6, 0, fmt_810, 0 };
483  */
484 
485  /* THIS ROUTINE IS A TRUNCATED-NEWTON METHOD. */
486  /* THE TRUNCATED-NEWTON METHOD IS PRECONDITIONED BY A LIMITED-MEMORY */
487  /* QUASI-NEWTON METHOD (THIS PRECONDITIONING STRATEGY IS DEVELOPED */
488  /* IN THIS ROUTINE) WITH A FURTHER DIAGONAL SCALING (SEE ROUTINE NDIA3).
489  */
490  /* FOR FURTHER DETAILS ON THE PARAMETERS, SEE ROUTINE TN. */
491 
492  /* THE FOLLOWING IMSL AND STANDARD FUNCTIONS ARE USED */
493 
494  /* INITIALIZE PARAMETERS AND CONSTANTS */
495 
496  /* Parameter adjustments */
497  --g;
498  --x;
499  --w;
500 
501  /* Function Body */
502  /* if (*msglvl >= -2) {
503  s_wsfe(&io___27);
504  e_wsfe();
505  } */
506  setpar_(n);
507  upd1 = TRUE_;
508  ireset = 0;
509  nfeval = 0;
510  nmodif = 0;
511  nlincg = 0;
512  fstop = *f;
513  zero = 0.;
514  one = 1.;
515  nm1 = *n - 1;
516 
517  /* WITHIN THIS ROUTINE THE ARRAY W(LOLDG) IS SHARED BY W(LHYR) */
518 
519  lhyr = subscr_1.loldg;
520 
521  /* CHECK PARAMETERS AND SET CONSTANTS */
522 
523  chkucp_(&subscr_1.lwtest, maxfun, &nwhy, n, &alpha, &epsmch, eta, &peps, &
524  rteps, &rtol, &rtolsq, stepmx, &ftest, xtol, &xnorm, &x[1], lw, &
525  small, &tiny, accrcy);
526 
527  if (nwhy < 0)
528  {
529  goto L120;
530  }
531 
532  CTruncatedNewton::setucr_(&small, &nftotl, &niter, n, f, &fnew, &fm, &gtg, &oldf,
533  sfun, &g[1], &x[1]);
534  fold = fnew;
535 
536  if (*msglvl >= 1)
537  {
538  /* s_wsfe(&io___56);
539  do_fio(&c__1, (char *)&niter, (ftnlen)sizeof(C_INT));
540  do_fio(&c__1, (char *)&nftotl, (ftnlen)sizeof(C_INT));
541  do_fio(&c__1, (char *)&nlincg, (ftnlen)sizeof(C_INT));
542  do_fio(&c__1, (char *)&fnew, (ftnlen)sizeof(C_FLOAT64));
543  do_fio(&c__1, (char *)&gtg, (ftnlen)sizeof(C_FLOAT64));
544  e_wsfe();*/
545  }
546 
547  /* CHECK FOR SMALL GRADIENT AT THE STARTING POINT. */
548 
549  ftest = one + fabs(fnew);
550 
551  if (gtg < epsmch * 1e-4 * ftest * ftest)
552  {
553  goto L90;
554  }
555 
556  /* SET INITIAL VALUES TO OTHER PARAMETERS */
557 
558  icycle = nm1;
559  toleps = rtol + rteps;
560  rtleps = rtolsq + epsmch;
561  gnorm = sqrt(gtg);
562  difnew = zero;
563  epsred = .05;
564  fkeep = fnew;
565 
566  /* SET THE DIAGONAL OF THE APPROXIMATE HESSIAN TO UNITY. */
567 
568  idiagb = subscr_1.ldiagb;
569  i__1 = *n;
570 
571  for (i__ = 1; i__ <= i__1; ++i__)
572  {
573  w[idiagb] = one;
574  ++idiagb;
575  /* L10: */
576  }
577 
578  /* ..................START OF MAIN ITERATIVE LOOP.......... */
579 
580  /* COMPUTE THE NEW SEARCH DIRECTION */
581 
582  modet = *msglvl - 3;
583  CTruncatedNewton::modlnp_(&modet, &w[subscr_1.lpk], &w[subscr_1.lgv], &w[subscr_1.lz1], &w[
584  subscr_1.lv], &w[subscr_1.ldiagb], &w[subscr_1.lemat], &x[1], &g[
585  1], &w[subscr_1.lzk], n, &w[1], lw, &niter, maxit, &nfeval, &
586  nmodif, &nlincg, &upd1, &yksk, &gsk, &yrsr, &lreset, sfun, &
587  c_false, &ipivot, accrcy, &gtpnew, &gnorm, &xnorm);
588 L20:
589  dcopy_(n, &g[1], &c__1, &w[subscr_1.loldg], &c__1);
590  pnorm = dnrm2_(n, &w[subscr_1.lpk], &c__1);
591  oldf = fnew;
592  oldgtp = gtpnew;
593 
594  /* PREPARE TO COMPUTE THE STEP LENGTH */
595 
596  pe = pnorm + epsmch;
597 
598  /* COMPUTE THE fabsOLUTE AND RELATIVE TOLERANCES FOR THE LINEAR SEARCH */
599 
600  reltol = rteps * (xnorm + one) / pe;
601  fabstol = -epsmch * ftest / (oldgtp - epsmch);
602 
603  /* COMPUTE THE SMALLEST ALLOWABLE SPACING BETWEEN POINTS IN */
604  /* THE LINEAR SEARCH */
605 
606  tnytol = epsmch * (xnorm + one) / pe;
607  spe = *stepmx / pe;
608 
609  /* SET THE INITIAL STEP LENGTH. */
610 
611  alpha = step1_(&fnew, &fm, &oldgtp, &spe);
612 
613  /* PERFORM THE LINEAR SEARCH */
614 
615  CTruncatedNewton::linder_(n, sfun, &small, &epsmch, &reltol, &fabstol, &tnytol, eta, &
616  zero, &spe, &w[subscr_1.lpk], &oldgtp, &x[1], &fnew, &alpha, &g[1]
617  , &numf, &nwhy, &w[1], lw);
618 
619  fold = fnew;
620  ++niter;
621  nftotl += numf;
622  gtg = ddot_(n, &g[1], &c__1, &g[1], &c__1);
623 
624  if (*msglvl >= 1)
625  {
626  //remove the printing
627  /*
628  s_wsfe(&io___81);
629  do_fio(&c__1, (char *)&niter, (ftnlen)sizeof(C_INT));
630  do_fio(&c__1, (char *)&nftotl, (ftnlen)sizeof(C_INT));
631  do_fio(&c__1, (char *)&nlincg, (ftnlen)sizeof(C_INT));
632  do_fio(&c__1, (char *)&fnew, (ftnlen)sizeof(C_FLOAT64));
633  do_fio(&c__1, (char *)&gtg, (ftnlen)sizeof(C_FLOAT64));
634  e_wsfe();*/
635  }
636 
637  if (nwhy < 0)
638  {
639  goto L120;
640  }
641 
642  if (nwhy == 0 || nwhy == 2)
643  {
644  goto L30;
645  }
646 
647  /* THE LINEAR SEARCH HAS FAILED TO FIND A LOWER POINT */
648 
649  nwhy = 3;
650  goto L100;
651 L30:
652 
653  if (nwhy <= 1)
654  {
655  goto L40;
656  }
657 
658  (*sfun)(n, &x[1], &fnew, &g[1]);
659  ++nftotl;
660 
661  /* TERMINATE IF MORE THAN MAXFUN EVALUTATIONS HAVE BEEN MADE */
662 
663 L40:
664  nwhy = 2;
665 
666  if (nftotl > *maxfun)
667  {
668  goto L110;
669  }
670 
671  nwhy = 0;
672 
673  /* SET UP PARAMETERS USED IN CONVERGENCE AND RESETTING TESTS */
674 
675  difold = difnew;
676  difnew = oldf - fnew;
677 
678  /* IF THIS IS THE FIRST ITERATION OF A NEW CYCLE, COMPUTE THE */
679  /* PERCENTAGE REDUCTION FACTOR FOR THE RESETTING TEST. */
680 
681  if (icycle != 1)
682  {
683  goto L50;
684  }
685 
686  if (difnew > difold * 2.)
687  {
688  epsred += epsred;
689  }
690 
691  if (difnew < difold * .5)
692  {
693  epsred *= .5;
694  }
695 
696 L50:
697  gnorm = sqrt(gtg);
698  ftest = one + fabs(fnew);
699  xnorm = dnrm2_(n, &x[1], &c__1);
700 
701  /* TEST FOR CONVERGENCE */
702 
703  if ((alpha * pnorm < toleps * (one + xnorm) &&
704  fabs(difnew) < rtleps * ftest &&
705  gtg < peps * ftest * ftest) ||
706  gtg < *accrcy * 1e-4 * ftest * ftest)
707  {
708  goto L90;
709  }
710 
711  /* COMPUTE THE CHANGE IN THE ITERATES AND THE CORRESPONDING CHANGE */
712  /* IN THE GRADIENTS */
713 
714  isk = subscr_1.lsk;
715  ipk = subscr_1.lpk;
716  iyk = subscr_1.lyk;
717  ioldg = subscr_1.loldg;
718  i__1 = *n;
719 
720  for (i__ = 1; i__ <= i__1; ++i__)
721  {
722  w[iyk] = g[i__] - w[ioldg];
723  w[isk] = alpha * w[ipk];
724  ++ipk;
725  ++isk;
726  ++iyk;
727  ++ioldg;
728  /* L60: */
729  }
730 
731  /* SET UP PARAMETERS USED IN UPDATING THE DIRECTION OF SEARCH. */
732 
733  yksk = ddot_(n, &w[subscr_1.lyk], &c__1, &w[subscr_1.lsk], &c__1);
734  lreset = FALSE_;
735 
736  if (icycle == nm1 || difnew < epsred * (fkeep - fnew))
737  {
738  lreset = TRUE_;
739  }
740 
741  if (lreset)
742  {
743  goto L70;
744  }
745 
746  yrsr = ddot_(n, &w[subscr_1.lyr], &c__1, &w[subscr_1.lsr], &c__1);
747 
748  if (yrsr <= zero)
749  {
750  lreset = TRUE_;
751  }
752 
753 L70:
754  upd1 = FALSE_;
755 
756  /* COMPUTE THE NEW SEARCH DIRECTION */
757 
758  modet = *msglvl - 3;
759  CTruncatedNewton::modlnp_(&modet, &w[subscr_1.lpk], &w[subscr_1.lgv], &w[subscr_1.lz1], &w[
760  subscr_1.lv], &w[subscr_1.ldiagb], &w[subscr_1.lemat], &x[1], &g[
761  1], &w[subscr_1.lzk], n, &w[1], lw, &niter, maxit, &nfeval, &
762  nmodif, &nlincg, &upd1, &yksk, &gsk, &yrsr, &lreset, sfun, &
763  c_false, &ipivot, accrcy, &gtpnew, &gnorm, &xnorm);
764 
765  if (lreset)
766  {
767  goto L80;
768  }
769 
770  /* STORE THE ACCUMULATED CHANGE IN THE POINT AND GRADIENT AS AN */
771  /* "AVERAGE" DIRECTION FOR PRECONDITIONING. */
772 
773  dxpy_(n, &w[subscr_1.lsk], &c__1, &w[subscr_1.lsr], &c__1);
774  dxpy_(n, &w[subscr_1.lyk], &c__1, &w[subscr_1.lyr], &c__1);
775  ++icycle;
776  goto L20;
777 
778  /* RESET */
779 
780 L80:
781  ++ireset;
782 
783  /* INITIALIZE THE SUM OF ALL THE CHANGES IN X. */
784 
785  dcopy_(n, &w[subscr_1.lsk], &c__1, &w[subscr_1.lsr], &c__1);
786  dcopy_(n, &w[subscr_1.lyk], &c__1, &w[subscr_1.lyr], &c__1);
787  fkeep = fnew;
788  icycle = 1;
789  goto L20;
790 
791  /* ...............END OF MAIN ITERATION....................... */
792 
793 L90:
794  *ifail = 0;
795  *f = fnew;
796  return 0;
797 L100:
798  oldf = fnew;
799 
800  /* LOCAL SEARCH HERE COULD BE INSTALLED HERE */
801 
802 L110:
803  *f = oldf;
804 
805  /* SET IFAIL */
806 
807 L120:
808  *ifail = nwhy;
809  return 0;
810 } /* lmqn_ */
#define C_INT
Definition: copasi.h:115
C_FLOAT64 step1_(C_FLOAT64 *fnew, C_FLOAT64 *fm, C_FLOAT64 *gtp, C_FLOAT64 *smax)
int linder_(C_INT *, FTruncatedNewton *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_INT *, C_FLOAT64 *, C_INT *)
int dcopy_(integer *n, doublereal *dx, integer *incx, doublereal *dy, integer *incy)
static C_INT c__1
#define TRUE_
doublereal dnrm2_(integer *n, doublereal *x, integer *incx)
doublereal ddot_(integer *n, doublereal *dx, integer *incx, doublereal *dy, integer *incy)
int dxpy_(C_INT *, C_FLOAT64 *, C_INT *, C_FLOAT64 *, C_INT *)
#define subscr_1
int modlnp_(C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_FLOAT64 *, C_INT *, C_INT *, C_INT *, C_INT *, C_INT *, C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, FTruncatedNewton *, C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *)
static C_INT c_false
int setucr_(C_FLOAT64 *, C_INT *, C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, FTruncatedNewton *, C_FLOAT64 *, C_FLOAT64 *)
int chkucp_(C_INT *, C_INT *, C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *)
#define C_FLOAT64
Definition: copasi.h:92
#define FALSE_
int CTruncatedNewton::lmqnbc_ ( C_INT ifail,
C_INT n,
C_FLOAT64 x,
C_FLOAT64 f,
C_FLOAT64 g,
C_FLOAT64 w,
C_INT lw,
FTruncatedNewton sfun,
C_FLOAT64 low,
C_FLOAT64 up,
C_INT ipivot,
C_INT msglvl,
C_INT maxit,
C_INT maxfun,
C_FLOAT64 eta,
C_FLOAT64 stepmx,
C_FLOAT64 accrcy,
C_FLOAT64 xtol 
)

Definition at line 812 of file CTruncatedNewton.cpp.

References c__1, C_FLOAT64, C_INT, c_true, chkucp_(), cnvtst_(), crash_(), dcopy_(), ddot_(), dnrm2_(), dxpy_(), FALSE_, linder_(), modlnp_(), modz_(), monit_(), setpar_(), setucr_(), step1_(), stpmax_(), subscr_1, TRUE_, and ztime_().

Referenced by tnbc_().

817 {
818  /* Format strings */
819  //remove
820  /*
821  char fmt_800[] = "(\002 THERE IS NO FEASIBLE POINT; TERMINATING A"
822  "LGORITHM\002)";
823  char fmt_810[] = "(//\002 NIT NF CG\002,9x,\002F\002,21x,"
824  "\002GTG\002,//)";
825  char fmt_820[] = "(\002 LINESEARCH RESULTS: ALPHA,PNOR"
826  "M\002,2(1pd12.4))";
827  char fmt_830[] = "(\002 UPD1 IS TRUE - TRIVIAL PRECONDITIONING"
828  "\002)";
829  char fmt_840[] = "(\002 NEWCON IS TRUE - CONSTRAINT ADDED IN LINE"
830  "SEARCH\002)";
831  */
832 
833  /* System generated locals */
834  C_INT i__1;
835  C_FLOAT64 d__1;
836 
837  /* Builtin functions */
838  // C_INT s_wsfe(cilist *), e_wsfe(void);
839  // C_INT do_fio(C_INT *, char *, ftnlen);
840 
841  /* Local variables */
842  C_FLOAT64 fold, oldf;
843  C_FLOAT64 fnew;
844  C_INT numf;
845  C_INT conv;
846  C_FLOAT64 peps;
847  C_INT lhyr;
848  C_FLOAT64 zero, rtol, yksk, tiny;
849  C_INT nwhy;
850  C_FLOAT64 yrsr;
851  C_INT i__;
852  C_FLOAT64 alpha, fkeep;
853  C_INT ioldg;
854  C_FLOAT64 small, flast;
855  C_INT modet;
856  C_INT niter;
857  C_FLOAT64 gnorm, ftest;
858  C_FLOAT64 fstop, pnorm, rteps, xnorm;
859  C_INT idiagb;
860  C_FLOAT64 fm, pe, difold;
861  C_INT icycle, nlincg, nfeval;
862  C_FLOAT64 difnew;
863  C_INT nmodif;
864  C_FLOAT64 epsmch;
865  C_FLOAT64 epsred, fabstol, oldgtp;
866  C_INT newcon;
867  C_INT ireset;
868  C_INT lreset;
869  C_FLOAT64 reltol, gtpnew;
870  C_INT nftotl;
871  C_FLOAT64 toleps;
872  C_FLOAT64 rtleps;
873  C_INT nm1;
874  C_FLOAT64 rtolsq, tnytol;
875  C_INT ier;
876  C_FLOAT64 gtg, one;
877  C_INT ipk;
878  C_FLOAT64 gsk, spe;
879  C_INT isk, iyk;
880  C_INT upd1;
881 
882  /* Fortran I/O blocks */
883  /*
884  cilist io___88 = {0, 6, 0, fmt_800, 0 };
885  cilist io___89 = {0, 6, 0, fmt_810, 0 };
886  cilist io___144 = {0, 6, 0, fmt_820, 0 };
887  cilist io___150 = {0, 6, 0, fmt_830, 0 };
888  cilist io___151 = {0, 6, 0, fmt_840, 0 }; */
889 
890  /* THIS ROUTINE IS A BOUNDS-CONSTRAINED TRUNCATED-NEWTON METHOD. */
891  /* THE TRUNCATED-NEWTON METHOD IS PRECONDITIONED BY A LIMITED-MEMORY */
892  /* QUASI-NEWTON METHOD (THIS PRECONDITIONING STRATEGY IS DEVELOPED */
893  /* IN THIS ROUTINE) WITH A FURTHER DIAGONAL SCALING (SEE ROUTINE NDIA3).
894  */
895  /* FOR FURTHER DETAILS ON THE PARAMETERS, SEE ROUTINE TNBC. */
896 
897  /* THE FOLLOWING STANDARD FUNCTIONS AND SYSTEM FUNCTIONS ARE USED */
898 
899  /* CHECK THAT INITIAL X IS FEASIBLE AND THAT THE BOUNDS ARE CONSISTENT */
900 
901  /* Parameter adjustments */
902  --ipivot;
903  --up;
904  --low;
905  --g;
906  --x;
907  --w;
908 
909  /* Function Body */
910  crash_(n, &x[1], &ipivot[1], &low[1], &up[1], &ier);
911  /*if (ier != 0) {
912  s_wsfe(&io___88);
913  e_wsfe();
914  }
915  if (ier != 0) {
916  return 0;
917  }
918  if (*msglvl >= 1) {
919  s_wsfe(&io___89);
920  e_wsfe();
921  } */
922 
923  /* INITIALIZE VARIABLES */
924 
925  setpar_(n);
926  upd1 = TRUE_;
927  ireset = 0;
928  nfeval = 0;
929  nmodif = 0;
930  nlincg = 0;
931  fstop = *f;
932  conv = FALSE_;
933  zero = 0.;
934  one = 1.;
935  nm1 = *n - 1;
936 
937  /* WITHIN THIS ROUTINE THE ARRAY W(LOLDG) IS SHARED BY W(LHYR) */
938 
939  lhyr = subscr_1.loldg;
940 
941  /* CHECK PARAMETERS AND SET CONSTANTS */
942 
943  chkucp_(&subscr_1.lwtest, maxfun, &nwhy, n, &alpha, &epsmch, eta, &peps, &
944  rteps, &rtol, &rtolsq, stepmx, &ftest, xtol, &xnorm, &x[1], lw, &
945  small, &tiny, accrcy);
946 
947  if (nwhy < 0)
948  {
949  goto L160;
950  }
951 
952  CTruncatedNewton::setucr_(&small, &nftotl, &niter, n, f, &fnew, &fm, &gtg, &oldf,
953  sfun, &g[1], &x[1]);
954  fold = fnew;
955  flast = fnew;
956 
957  /* TEST THE LAGRANGE MULTIPLIERS TO SEE IF THEY ARE NON-NEGATIVE. */
958  /* BECAUSE THE CONSTRAINTS ARE ONLY LOWER BOUNDS, THE COMPONENTS */
959  /* OF THE GRADIENT CORRESPONDING TO THE ACTIVE CONSTRAINTS ARE THE */
960  /* LAGRANGE MULTIPLIERS. AFTERWORDS, THE PROJECTED GRADIENT IS FORMED. */
961 
962  i__1 = *n;
963 
964  for (i__ = 1; i__ <= i__1; ++i__)
965  {
966  if (ipivot[i__] == 2)
967  {
968  goto L10;
969  }
970 
971  if (-ipivot[i__] * g[i__] >= 0.)
972  {
973  goto L10;
974  }
975 
976  ipivot[i__] = 0;
977 L10:
978  ;
979  }
980 
981  ztime_(n, &g[1], &ipivot[1]);
982  gtg = ddot_(n, &g[1], &c__1, &g[1], &c__1);
983 
984  if (*msglvl >= 1)
985  {
986  monit_(n, &x[1], &fnew, &g[1], &niter, &nftotl, &nfeval, &lreset, &
987  ipivot[1]);
988  }
989 
990  /* CHECK IF THE INITIAL POINT IS A LOCAL MINIMUM. */
991 
992  ftest = one + fabs(fnew);
993 
994  if (gtg < epsmch * 1e-4 * ftest * ftest)
995  {
996  goto L130;
997  }
998 
999  /* SET INITIAL VALUES TO OTHER PARAMETERS */
1000 
1001  icycle = nm1;
1002  toleps = rtol + rteps;
1003  rtleps = rtolsq + epsmch;
1004  gnorm = sqrt(gtg);
1005  difnew = zero;
1006  epsred = .05;
1007  fkeep = fnew;
1008 
1009  /* SET THE DIAGONAL OF THE APPROXIMATE HESSIAN TO UNITY. */
1010 
1011  idiagb = subscr_1.ldiagb;
1012  i__1 = *n;
1013 
1014  for (i__ = 1; i__ <= i__1; ++i__)
1015  {
1016  w[idiagb] = one;
1017  ++idiagb;
1018  /* L15: */
1019  }
1020 
1021  /* ..................START OF MAIN ITERATIVE LOOP.......... */
1022 
1023  /* COMPUTE THE NEW SEARCH DIRECTION */
1024 
1025  modet = *msglvl - 3;
1026  CTruncatedNewton::modlnp_(&modet, &w[subscr_1.lpk], &w[subscr_1.lgv], &w[subscr_1.lz1], &w[
1027  subscr_1.lv], &w[subscr_1.ldiagb], &w[subscr_1.lemat], &x[1], &g[
1028  1], &w[subscr_1.lzk], n, &w[1], lw, &niter, maxit, &nfeval, &
1029  nmodif, &nlincg, &upd1, &yksk, &gsk, &yrsr, &lreset, sfun, &
1030  c_true, &ipivot[1], accrcy, &gtpnew, &gnorm, &xnorm);
1031 L20:
1032 
1033  /* added manually by Pedro Mendes 12/2/1998 */
1034  // if(callback != 0/*NULL*/) callback(fnew);
1035 
1036  dcopy_(n, &g[1], &c__1, &w[subscr_1.loldg], &c__1);
1037  pnorm = dnrm2_(n, &w[subscr_1.lpk], &c__1);
1038  oldf = fnew;
1039  oldgtp = gtpnew;
1040 
1041  /* PREPARE TO COMPUTE THE STEP LENGTH */
1042 
1043  pe = pnorm + epsmch;
1044 
1045  /* COMPUTE THE fabsOLUTE AND RELATIVE TOLERANCES FOR THE LINEAR SEARCH */
1046 
1047  reltol = rteps * (xnorm + one) / pe;
1048  fabstol = -epsmch * ftest / (oldgtp - epsmch);
1049 
1050  /* COMPUTE THE SMALLEST ALLOWABLE SPACING BETWEEN POINTS IN */
1051  /* THE LINEAR SEARCH */
1052 
1053  tnytol = epsmch * (xnorm + one) / pe;
1054  stpmax_(stepmx, &pe, &spe, n, &x[1], &w[subscr_1.lpk], &ipivot[1], &low[1]
1055  , &up[1]);
1056 
1057  /* SET THE INITIAL STEP LENGTH. */
1058 
1059  alpha = step1_(&fnew, &fm, &oldgtp, &spe);
1060 
1061  /* PERFORM THE LINEAR SEARCH */
1062 
1063  CTruncatedNewton::linder_(n, sfun, &small, &epsmch, &reltol, &fabstol, &tnytol, eta, &
1064  zero, &spe, &w[subscr_1.lpk], &oldgtp, &x[1], &fnew, &alpha, &g[1]
1065  , &numf, &nwhy, &w[1], lw);
1066  newcon = FALSE_;
1067 
1068  if ((d__1 = alpha - spe, fabs(d__1)) > epsmch * 10.)
1069  {
1070  goto L30;
1071  }
1072 
1073  newcon = TRUE_;
1074  nwhy = 0;
1075  modz_(n, &x[1], &w[subscr_1.lpk], &ipivot[1], &epsmch, &low[1], &up[1], &
1076  flast, &fnew);
1077  flast = fnew;
1078 
1079 L30:
1080 
1081  if (*msglvl >= 3)
1082  {
1083  /* s_wsfe(&io___144);
1084  do_fio(&c__1, (char *)&alpha, (ftnlen)sizeof(C_FLOAT64));
1085  do_fio(&c__1, (char *)&pnorm, (ftnlen)sizeof(C_FLOAT64));
1086  e_wsfe();*/
1087  }
1088 
1089  fold = fnew;
1090  ++niter;
1091  nftotl += numf;
1092 
1093  /* IF REQUIRED, PRINT THE DETAILS OF THIS ITERATION */
1094 
1095  /* ############################ */
1096  if (*msglvl >= 1)
1097  {
1098  monit_(n, &x[1], &fnew, &g[1], &niter, &nftotl, &nfeval, &lreset, &
1099  ipivot[1]);
1100  }
1101 
1102  if (nwhy < 0)
1103  {
1104  goto L160;
1105  }
1106 
1107  if (nwhy == 0 || nwhy == 2)
1108  {
1109  goto L40;
1110  }
1111 
1112  /* THE LINEAR SEARCH HAS FAILED TO FIND A LOWER POINT */
1113 
1114  nwhy = 3;
1115  goto L140;
1116 L40:
1117 
1118  if (nwhy <= 1)
1119  {
1120  goto L50;
1121  }
1122 
1123  (*sfun)(n, &x[1], &fnew, &g[1]);
1124  ++nftotl;
1125 
1126  /* TERMINATE IF MORE THAN MAXFUN EVALUATIONS HAVE BEEN MADE */
1127 
1128 L50:
1129  nwhy = 2;
1130 
1131  if (nftotl > *maxfun)
1132  {
1133  goto L150;
1134  }
1135 
1136  nwhy = 0;
1137 
1138  /* SET UP PARAMETERS USED IN CONVERGENCE AND RESETTING TESTS */
1139 
1140  difold = difnew;
1141  difnew = oldf - fnew;
1142 
1143  /* IF THIS IS THE FIRST ITERATION OF A NEW CYCLE, COMPUTE THE */
1144  /* PERCENTAGE REDUCTION FACTOR FOR THE RESETTING TEST. */
1145 
1146  if (icycle != 1)
1147  {
1148  goto L60;
1149  }
1150 
1151  if (difnew > difold * 2.)
1152  {
1153  epsred += epsred;
1154  }
1155 
1156  if (difnew < difold * .5)
1157  {
1158  epsred *= .5;
1159  }
1160 
1161 L60:
1162  dcopy_(n, &g[1], &c__1, &w[subscr_1.lgv], &c__1);
1163  ztime_(n, &w[subscr_1.lgv], &ipivot[1]);
1164  gtg = ddot_(n, &w[subscr_1.lgv], &c__1, &w[subscr_1.lgv], &c__1);
1165  gnorm = sqrt(gtg);
1166  ftest = one + fabs(fnew);
1167  xnorm = dnrm2_(n, &x[1], &c__1);
1168 
1169  /* TEST FOR CONVERGENCE */
1170 
1171  cnvtst_(&conv, &alpha, &pnorm, &toleps, &xnorm, &difnew, &rtleps, &ftest,
1172  &gtg, &peps, &epsmch, &gtpnew, &fnew, &flast, &g[1], &ipivot[1],
1173  n, accrcy);
1174 
1175  if (conv)
1176  {
1177  goto L130;
1178  }
1179 
1180  ztime_(n, &g[1], &ipivot[1]);
1181 
1182  /* COMPUTE THE CHANGE IN THE ITERATES AND THE CORRESPONDING CHANGE */
1183  /* IN THE GRADIENTS */
1184 
1185  if (newcon)
1186  {
1187  goto L90;
1188  }
1189 
1190  isk = subscr_1.lsk;
1191  ipk = subscr_1.lpk;
1192  iyk = subscr_1.lyk;
1193  ioldg = subscr_1.loldg;
1194  i__1 = *n;
1195 
1196  for (i__ = 1; i__ <= i__1; ++i__)
1197  {
1198  w[iyk] = g[i__] - w[ioldg];
1199  w[isk] = alpha * w[ipk];
1200  ++ipk;
1201  ++isk;
1202  ++iyk;
1203  ++ioldg;
1204  /* L70: */
1205  }
1206 
1207  /* SET UP PARAMETERS USED IN UPDATING THE PRECONDITIONING STRATEGY. */
1208 
1209  yksk = ddot_(n, &w[subscr_1.lyk], &c__1, &w[subscr_1.lsk], &c__1);
1210  lreset = FALSE_;
1211 
1212  if (icycle == nm1 || difnew < epsred * (fkeep - fnew))
1213  {
1214  lreset = TRUE_;
1215  }
1216 
1217  if (lreset)
1218  {
1219  goto L80;
1220  }
1221 
1222  yrsr = ddot_(n, &w[subscr_1.lyr], &c__1, &w[subscr_1.lsr], &c__1);
1223 
1224  if (yrsr <= zero)
1225  {
1226  lreset = TRUE_;
1227  }
1228 
1229 L80:
1230  upd1 = FALSE_;
1231 
1232  /* COMPUTE THE NEW SEARCH DIRECTION */
1233 
1234 L90:
1235 
1236  if (upd1 && *msglvl >= 3)
1237  {
1238  /*
1239  s_wsfe(&io___150);
1240  e_wsfe();*/
1241  }
1242 
1243  if (newcon && *msglvl >= 3)
1244  {
1245  /*
1246  s_wsfe(&io___151);
1247  e_wsfe();*/
1248  }
1249 
1250  modet = *msglvl - 3;
1251  CTruncatedNewton::modlnp_(&modet, &w[subscr_1.lpk], &w[subscr_1.lgv], &w[subscr_1.lz1], &w[
1252  subscr_1.lv], &w[subscr_1.ldiagb], &w[subscr_1.lemat], &x[1], &g[
1253  1], &w[subscr_1.lzk], n, &w[1], lw, &niter, maxit, &nfeval, &
1254  nmodif, &nlincg, &upd1, &yksk, &gsk, &yrsr, &lreset, sfun, &
1255  c_true, &ipivot[1], accrcy, &gtpnew, &gnorm, &xnorm);
1256 
1257  if (newcon)
1258  {
1259  goto L20;
1260  }
1261 
1262  if (lreset)
1263  {
1264  goto L110;
1265  }
1266 
1267  /* COMPUTE THE ACCUMULATED STEP AND ITS CORRESPONDING */
1268  /* GRADIENT DIFFERENCE. */
1269 
1270  dxpy_(n, &w[subscr_1.lsk], &c__1, &w[subscr_1.lsr], &c__1);
1271  dxpy_(n, &w[subscr_1.lyk], &c__1, &w[subscr_1.lyr], &c__1);
1272  ++icycle;
1273  goto L20;
1274 
1275  /* RESET */
1276 
1277 L110:
1278  ++ireset;
1279 
1280  /* INITIALIZE THE SUM OF ALL THE CHANGES IN X. */
1281 
1282  dcopy_(n, &w[subscr_1.lsk], &c__1, &w[subscr_1.lsr], &c__1);
1283  dcopy_(n, &w[subscr_1.lyk], &c__1, &w[subscr_1.lyr], &c__1);
1284  fkeep = fnew;
1285  icycle = 1;
1286  goto L20;
1287 
1288  /* ...............END OF MAIN ITERATION....................... */
1289 
1290 L130:
1291  *ifail = 0;
1292  *f = fnew;
1293  return 0;
1294 L140:
1295  oldf = fnew;
1296 
1297  /* LOCAL SEARCH COULD BE INSTALLED HERE */
1298 
1299 L150:
1300  *f = oldf;
1301 
1302  if (*msglvl >= 1)
1303  {
1304  monit_(n, &x[1], f, &g[1], &niter, &nftotl, &nfeval, &ireset, &ipivot[
1305  1]);
1306  }
1307 
1308  /* SET IFAIL */
1309 
1310 L160:
1311  *ifail = nwhy;
1312  return 0;
1313 } /* lmqnbc_ */
#define C_INT
Definition: copasi.h:115
C_FLOAT64 step1_(C_FLOAT64 *fnew, C_FLOAT64 *fm, C_FLOAT64 *gtp, C_FLOAT64 *smax)
int linder_(C_INT *, FTruncatedNewton *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_INT *, C_FLOAT64 *, C_INT *)
static C_INT c_true
int dcopy_(integer *n, doublereal *dx, integer *incx, doublereal *dy, integer *incy)
int modz_(C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *)
static C_INT c__1
#define TRUE_
int crash_(C_INT *, C_FLOAT64 *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_INT *)
int stpmax_(C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_FLOAT64 *, C_FLOAT64 *)
doublereal dnrm2_(integer *n, doublereal *x, integer *incx)
doublereal ddot_(integer *n, doublereal *dx, integer *incx, doublereal *dy, integer *incy)
int dxpy_(C_INT *, C_FLOAT64 *, C_INT *, C_FLOAT64 *, C_INT *)
#define subscr_1
int modlnp_(C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_FLOAT64 *, C_INT *, C_INT *, C_INT *, C_INT *, C_INT *, C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, FTruncatedNewton *, C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *)
int monit_(C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_INT *, C_INT *, C_INT *, C_INT *)
int setucr_(C_FLOAT64 *, C_INT *, C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, FTruncatedNewton *, C_FLOAT64 *, C_FLOAT64 *)
int chkucp_(C_INT *, C_INT *, C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *)
#define C_FLOAT64
Definition: copasi.h:92
int ztime_(C_INT *, C_FLOAT64 *, C_INT *)
#define FALSE_
int cnvtst_(C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_INT *, C_FLOAT64 *)
int CTruncatedNewton::modlnp_ ( C_INT modet,
C_FLOAT64 zsol,
C_FLOAT64 gv,
C_FLOAT64 r__,
C_FLOAT64 v,
C_FLOAT64 diagb,
C_FLOAT64 emat,
C_FLOAT64 x,
C_FLOAT64 g,
C_FLOAT64 zk,
C_INT n,
C_FLOAT64 w,
C_INT lw,
C_INT ,
C_INT maxit,
C_INT nfeval,
C_INT ,
C_INT nlincg,
C_INT upd1,
C_FLOAT64 yksk,
C_FLOAT64 gsk,
C_FLOAT64 yrsr,
C_INT lreset,
FTruncatedNewton sfun,
C_INT bounds,
C_INT ipivot,
C_FLOAT64 accrcy,
C_FLOAT64 gtp,
C_FLOAT64 gnorm,
C_FLOAT64 xnorm 
)

Definition at line 1666 of file CTruncatedNewton.cpp.

References c__1, C_FLOAT64, C_INT, daxpy_(), dcopy_(), ddot_(), gtims_(), initpc_(), msolve_(), ndia3_(), negvec_(), TRUE_, and ztime_().

Referenced by lmqn_(), and lmqnbc_().

1674 {
1675  /* Format strings */
1676  //remove the format
1677  /* char fmt_800[] = "(\002 \002,//,\002 ENTERING MODLNP\002)";
1678  char fmt_810[] = "(\002 \002,//,\002 ### ITERATION \002,i2,\002 #"
1679  "##\002)";
1680  char fmt_820[] = "(\002 ALPHA\002,1pd16.8)";
1681  char fmt_830[] = "(\002 G(T)Z POSITIVE AT ITERATION \002,i2,\002 "
1682  "- TRUNCATING METHOD\002,/)";
1683  char fmt_840[] = "(\002 \002,10x,\002HESSIAN NOT POSITIVE-DEFIN"
1684  "ITE\002)";
1685  char fmt_850[] = "(\002 \002,/,8x,\002MODLAN TRUNCATED AFTER \002"
1686  ",i3,\002 ITERATIONS\002,\002 RNORM = \002,1pd14.6)";
1687  char fmt_860[] = "(\002 PRECONDITIONING NOT POSITIVE-DEFINITE\002)"
1688  ;*/
1689 
1690  /* System generated locals */
1691  C_INT i__1, i__2;
1692  C_FLOAT64 d__1;
1693 
1694  /* Builtin functions */
1695  // C_INT s_wsfe(cilist *), e_wsfe(void), do_fio(C_INT *, char *, ftnlen);
1696 
1697  /* Local variables */
1698  C_FLOAT64 beta = 0.0;
1699  C_FLOAT64 qold, qnew;
1700  C_INT i__, k;
1701  C_FLOAT64 alpha, delta;
1702  C_INT first;
1703  C_FLOAT64 rzold = 0.0, qtest, pr;
1704  C_FLOAT64 rz;
1705  C_FLOAT64 rhsnrm, tol, vgv;
1706 
1707  /* THIS ROUTINE PERFORMS A PRECONDITIONED CONJUGATE-GRADIENT */
1708  /* ITERATION IN ORDER TO SOLVE THE NEWTON EQUATIONS FOR A SEARCH */
1709  /* DIRECTION FOR A TRUNCATED-NEWTON ALGORITHM. WHEN THE VALUE OF THE */
1710  /* QUADRATIC MODEL IS SUFFICIENTLY REDUCED, */
1711  /* THE ITERATION IS TERMINATED. */
1712 
1713  /* PARAMETERS */
1714 
1715  /* MODET - C_INT WHICH CONTROLS AMOUNT OF OUTPUT */
1716  /* ZSOL - COMPUTED SEARCH DIRECTION */
1717  /* G - CURRENT GRADIENT */
1718  /* GV,GZ1,V - SCRATCH VECTORS */
1719  /* R - RESIDUAL */
1720  /* DIAGB,EMAT - DIAGONAL PRECONDITONING MATRIX */
1721  /* NITER - NONLINEAR ITERATION # */
1722  /* FEVAL - VALUE OF QUADRATIC FUNCTION */
1723 
1724  /* ************************************************************* */
1725  /* INITIALIZATION */
1726  /* ************************************************************* */
1727 
1728  /* GENERAL INITIALIZATION */
1729 
1730  /* Parameter adjustments */
1731  --zk;
1732  --g;
1733  --x;
1734  --emat;
1735  --diagb;
1736  --v;
1737  --r__;
1738  --gv;
1739  --zsol;
1740  --w;
1741  --ipivot;
1742 
1743  /* Function Body */
1744  if (*modet > 0)
1745  {
1746  //remove
1747  /*
1748  s_wsfe(&io___167);
1749  e_wsfe();*/
1750  }
1751 
1752  if (*maxit == 0)
1753  {
1754  return 0;
1755  }
1756 
1757  first = TRUE_;
1758  rhsnrm = *gnorm;
1759  tol = 1e-12;
1760  qold = 0.;
1761 
1762  /* INITIALIZATION FOR PRECONDITIONED CONJUGATE-GRADIENT ALGORITHM */
1763 
1764  CTruncatedNewton::initpc_(&diagb[1], &emat[1], n, &w[1], lw, modet, upd1, yksk, gsk, yrsr,
1765  lreset);
1766  i__1 = *n;
1767 
1768  for (i__ = 1; i__ <= i__1; ++i__)
1769  {
1770  r__[i__] = -g[i__];
1771  v[i__] = 0.;
1772  zsol[i__] = 0.;
1773  /* L10: */
1774  }
1775 
1776  /* ************************************************************ */
1777  /* MAIN ITERATION */
1778  /* ************************************************************ */
1779 
1780  i__1 = *maxit;
1781 
1782  for (k = 1; k <= i__1; ++k)
1783  {
1784  ++(*nlincg);
1785 
1786  if (*modet > 1)
1787  {
1788  /*
1789  s_wsfe(&io___174);
1790  do_fio(&c__1, (char *)&k, (ftnlen)sizeof(C_INT));
1791  e_wsfe();*/
1792  }
1793 
1794  /* CG ITERATION TO SOLVE SYSTEM OF EQUATIONS */
1795 
1796  if (*bounds)
1797  {
1798  ztime_(n, &r__[1], &ipivot[1]);
1799  }
1800 
1801  CTruncatedNewton::msolve_(&r__[1], &zk[1], n, &w[1], lw, upd1, yksk, gsk, yrsr, lreset,
1802  &first);
1803 
1804  if (*bounds)
1805  {
1806  ztime_(n, &zk[1], &ipivot[1]);
1807  }
1808 
1809  rz = ddot_(n, &r__[1], &c__1, &zk[1], &c__1);
1810 
1811  if (rz / rhsnrm < tol)
1812  {
1813  goto L80;
1814  }
1815 
1816  if (k == 1)
1817  {
1818  beta = 0.;
1819  }
1820 
1821  if (k > 1)
1822  {
1823  beta = rz / rzold;
1824  }
1825 
1826  i__2 = *n;
1827 
1828  for (i__ = 1; i__ <= i__2; ++i__)
1829  {
1830  v[i__] = zk[i__] + beta * v[i__];
1831  /* L20: */
1832  }
1833 
1834  if (*bounds)
1835  {
1836  ztime_(n, &v[1], &ipivot[1]);
1837  }
1838 
1839  CTruncatedNewton::gtims_(&v[1], &gv[1], n, &x[1], &g[1], &w[1], lw, sfun, &first,
1840  &delta, accrcy, xnorm);
1841 
1842  if (*bounds)
1843  {
1844  ztime_(n, &gv[1], &ipivot[1]);
1845  }
1846 
1847  ++(*nfeval);
1848  vgv = ddot_(n, &v[1], &c__1, &gv[1], &c__1);
1849 
1850  if (vgv / rhsnrm < tol)
1851  {
1852  goto L50;
1853  }
1854 
1855  ndia3_(n, &emat[1], &v[1], &gv[1], &r__[1], &vgv, modet);
1856 
1857  /* COMPUTE LINEAR STEP LENGTH */
1858 
1859  alpha = rz / vgv;
1860 
1861  if (*modet >= 1)
1862  {
1863  /*
1864  s_wsfe(&io___181);
1865  do_fio(&c__1, (char *)&alpha, (ftnlen)sizeof(C_FLOAT64));
1866  e_wsfe();*/
1867  }
1868 
1869  /* COMPUTE CURRENT SOLUTION AND RELATED VECTORS */
1870 
1871  daxpy_(n, &alpha, &v[1], &c__1, &zsol[1], &c__1);
1872  d__1 = -alpha;
1873  daxpy_(n, &d__1, &gv[1], &c__1, &r__[1], &c__1);
1874 
1875  /* TEST FOR CONVERGENCE */
1876 
1877  *gtp = ddot_(n, &zsol[1], &c__1, &g[1], &c__1);
1878  pr = ddot_(n, &r__[1], &c__1, &zsol[1], &c__1);
1879  qnew = (*gtp + pr) * .5;
1880  qtest = k * (1. - qold / qnew);
1881 
1882  if (qtest < 0.)
1883  {
1884  goto L70;
1885  }
1886 
1887  qold = qnew;
1888 
1889  if (qtest <= .5)
1890  {
1891  goto L70;
1892  }
1893 
1894  /* PERFORM CAUTIONARY TEST */
1895 
1896  if (*gtp > 0.)
1897  {
1898  goto L40;
1899  }
1900 
1901  rzold = rz;
1902  /* L30: */
1903  }
1904 
1905  /* TERMINATE ALGORITHM */
1906 
1907  --k;
1908  goto L70;
1909 
1910  /* TRUNCATE ALGORITHM IN CASE OF AN EMERGENCY */
1911 
1912 L40:
1913 
1914  if (*modet >= -1)
1915  {
1916  /*
1917  s_wsfe(&io___185);
1918  do_fio(&c__1, (char *)&k, (ftnlen)sizeof(C_INT));
1919  e_wsfe();*/
1920  }
1921 
1922  d__1 = -alpha;
1923  daxpy_(n, &d__1, &v[1], &c__1, &zsol[1], &c__1);
1924  *gtp = ddot_(n, &zsol[1], &c__1, &g[1], &c__1);
1925  goto L90;
1926 L50:
1927 
1928  if (*modet > -2)
1929  {
1930  /*
1931  s_wsfe(&io___186);
1932  e_wsfe();*/
1933  }
1934 
1935  /* L60: */
1936  if (k > 1)
1937  {
1938  goto L70;
1939  }
1940 
1941  CTruncatedNewton::msolve_(&g[1], &zsol[1], n, &w[1], lw, upd1, yksk, gsk, yrsr, lreset, &
1942  first);
1943  negvec_(n, &zsol[1]);
1944 
1945  if (*bounds)
1946  {
1947  ztime_(n, &zsol[1], &ipivot[1]);
1948  }
1949 
1950  *gtp = ddot_(n, &zsol[1], &c__1, &g[1], &c__1);
1951 L70:
1952 
1953  if (*modet >= -1)
1954  {
1955  /*
1956  s_wsfe(&io___187);
1957  do_fio(&c__1, (char *)&k, (ftnlen)sizeof(C_INT));
1958  do_fio(&c__1, (char *)&rnorm, (ftnlen)sizeof(C_FLOAT64));
1959  e_wsfe();*/
1960  }
1961 
1962  goto L90;
1963 L80:
1964 
1965  if (*modet >= -1)
1966  {
1967  /*
1968  s_wsfe(&io___189);
1969  e_wsfe();*/
1970  }
1971 
1972  if (k > 1)
1973  {
1974  goto L70;
1975  }
1976 
1977  dcopy_(n, &g[1], &c__1, &zsol[1], &c__1);
1978  negvec_(n, &zsol[1]);
1979 
1980  if (*bounds)
1981  {
1982  ztime_(n, &zsol[1], &ipivot[1]);
1983  }
1984 
1985  *gtp = ddot_(n, &zsol[1], &c__1, &g[1], &c__1);
1986  goto L70;
1987 
1988  /* STORE (OR RESTORE) DIAGONAL PRECONDITIONING */
1989 
1990 L90:
1991  dcopy_(n, &emat[1], &c__1, &diagb[1], &c__1);
1992  return 0;
1993 } /* modlnp_ */
#define C_INT
Definition: copasi.h:115
int daxpy_(integer *n, doublereal *da, doublereal *dx, integer *incx, doublereal *dy, integer *incy)
int dcopy_(integer *n, doublereal *dx, integer *incx, doublereal *dy, integer *incy)
int negvec_(C_INT *, C_FLOAT64 *)
int ndia3_(C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *)
static C_INT c__1
int initpc_(C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_FLOAT64 *, C_INT *, C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *)
#define TRUE_
doublereal ddot_(integer *n, doublereal *dx, integer *incx, doublereal *dy, integer *incy)
int gtims_(C_FLOAT64 *v, C_FLOAT64 *gv, C_INT *n, C_FLOAT64 *x, C_FLOAT64 *g, C_FLOAT64 *w, C_INT *, FTruncatedNewton *sfun, C_INT *first, C_FLOAT64 *delta, C_FLOAT64 *accrcy, C_FLOAT64 *xnorm)
#define C_FLOAT64
Definition: copasi.h:92
int ztime_(C_INT *, C_FLOAT64 *, C_INT *)
int msolve_(C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_FLOAT64 *, C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_INT *)
int CTruncatedNewton::msolve_ ( C_FLOAT64 g,
C_FLOAT64 y,
C_INT n,
C_FLOAT64 w,
C_INT ,
C_INT upd1,
C_FLOAT64 yksk,
C_FLOAT64 gsk,
C_FLOAT64 yrsr,
C_INT lreset,
C_INT first 
)

Definition at line 2320 of file CTruncatedNewton.cpp.

References mslv_(), and subscr_2.

Referenced by modlnp_().

2323 {
2324  /* THIS ROUTINE SETS UPT THE ARRAYS FOR MSLV */
2325 
2326  /* Parameter adjustments */
2327  --y;
2328  --g;
2329  --w;
2330 
2331  /* Function Body */
2332  mslv_(&g[1], &y[1], n, &w[subscr_2.lsk], &w[subscr_2.lyk], &w[
2333  subscr_2.ldiagb], &w[subscr_2.lsr], &w[subscr_2.lyr], &w[
2334  subscr_2.lhyr], &w[subscr_2.lhg], &w[subscr_2.lhyk], upd1, yksk,
2335  gsk, yrsr, lreset, first);
2336  return 0;
2337 } /* msolve_ */
int mslv_(C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_INT *)
#define subscr_2
int CTruncatedNewton::setpar_ ( C_INT n)

Definition at line 2671 of file CTruncatedNewton.cpp.

References C_INT, and subscr_3.

Referenced by lmqn_(), and lmqnbc_().

2672 {
2673  C_INT i__;
2674 
2675  /* SET UP PARAMETERS FOR THE OPTIMIZATION ROUTINE */
2676 
2677  for (i__ = 1; i__ <= 14; ++i__)
2678  {
2679  subscr_3.lsub[i__ - 1] = (i__ - 1) * *n + 1;
2680  /* L10: */
2681  }
2682 
2683  subscr_3.lwtest = subscr_3.lsub[13] + *n - 1;
2684  return 0;
2685 } /* setpar_ */
#define C_INT
Definition: copasi.h:115
#define subscr_3
int CTruncatedNewton::setucr_ ( C_FLOAT64 ,
C_INT nftotl,
C_INT niter,
C_INT n,
C_FLOAT64 f,
C_FLOAT64 fnew,
C_FLOAT64 fm,
C_FLOAT64 gtg,
C_FLOAT64 oldf,
FTruncatedNewton sfun,
C_FLOAT64 g,
C_FLOAT64 x 
)

Definition at line 2238 of file CTruncatedNewton.cpp.

References c__1, and ddot_().

Referenced by lmqn_(), and lmqnbc_().

2242 {
2243  /* CHECK INPUT PARAMETERS, COMPUTE THE INITIAL FUNCTION VALUE, SET */
2244  /* CONSTANTS FOR THE SUBSEQUENT MINIMIZATION */
2245 
2246  /* Parameter adjustments */
2247  --x;
2248  --g;
2249 
2250  /* Function Body */
2251  *fm = *f;
2252 
2253  /* COMPUTE THE INITIAL FUNCTION VALUE */
2254 
2255  (*sfun)(n, &x[1], fnew, &g[1]);
2256  *nftotl = 1;
2257 
2258  /* SET CONSTANTS FOR LATER */
2259 
2260  *niter = 0;
2261  *oldf = *fnew;
2262  *gtg = ddot_(n, &g[1], &c__1, &g[1], &c__1);
2263  return 0;
2264 } /* setucr_ */
static C_INT c__1
doublereal ddot_(integer *n, doublereal *dx, integer *incx, doublereal *dy, integer *incy)
int CTruncatedNewton::tn_ ( C_INT ierror,
C_INT n,
C_FLOAT64 x,
C_FLOAT64 f,
C_FLOAT64 g,
C_FLOAT64 w,
C_INT lw,
FTruncatedNewton sfun 
)

Definition at line 105 of file CTruncatedNewton.cpp.

References C_FLOAT64, C_INT, lmqn_(), and mchpr1_().

107 {
108  /* Format strings */
109  /* static char fmt_800[] = "(//,\002 ERROR CODE =\002,i3)";
110  static char fmt_810[] = "(//,\002 OPTIMAL FUNCTION VALUE = \002,1pd22.15)"
111  ;
112  static char fmt_820[] = "(10x,\002CURRENT SOLUTION IS (AT MOST 10 COMPON"
113  "ENTS)\002,/,14x,\002I\002,11x,\002X(I)\002)";
114  static char fmt_830[] = "(10x,i5,2x,1pd22.15)";*/
115 
116  /* Local variables */
117  C_FLOAT64 xtol;
118  C_INT maxit;
119  C_FLOAT64 accrcy;
120  C_INT maxfun, msglvl;
121  C_FLOAT64 stepmx, eta;
122 
123  /* Fortran I/O blocks */ /*
124  cilist io___8 = {0, 6, 0, fmt_800, 0 };
125  cilist io___9 = {0, 6, 0, fmt_810, 0 };
126  cilist io___10 = {0, 6, 0, fmt_820, 0 };
127  cilist io___12 = {0, 6, 0, fmt_830, 0 };*/
128 
129  /* THIS ROUTINE SOLVES THE OPTIMIZATION PROBLEM */
130 
131  /* MINIMIZE F(X) */
132  /* X */
133 
134  /* WHERE X IS A VECTOR OF N REAL VARIABLES. THE METHOD USED IS */
135  /* A TRUNCATED-NEWTON ALGORITHM (SEE "NEWTON-TYPE MINIMIZATION VIA */
136  /* THE LANCZOS METHOD" BY S.G. NASH (SIAM J. NUMER. ANAL. 21 (1984), */
137  /* PP. 770-778). THIS ALGORITHM FINDS A LOCAL MINIMUM OF F(X). IT DOES
138  */
139  /* NOT ASSUME THAT THE FUNCTION F IS CONVEX (AND SO CANNOT GUARANTEE A */
140  /* GLOBAL SOLUTION), BUT DOES ASSUME THAT THE FUNCTION IS BOUNDED BELOW.
141  */
142  /* IT CAN SOLVE PROBLEMS HAVING ANY NUMBER OF VARIABLES, BUT IT IS */
143  /* ESPECIALLY USEFUL WHEN THE NUMBER OF VARIABLES (N) IS LARGE. */
144 
145  /* SUBROUTINE PARAMETERS: */
146 
147  /* IERROR - (C_INT) ERROR CODE */
148  /* (0 => NORMAL RETURN) */
149  /* (2 => MORE THAN MAXFUN EVALUATIONS) */
150  /* (3 => LINE SEARCH FAILED TO FIND */
151  /* (LOWER POINT (MAY NOT BE SERIOUS) */
152  /* (-1 => ERROR IN INPUT PARAMETERS) */
153  /* N - (C_INT) NUMBER OF VARIABLES */
154  /* X - (REAL*8) VECTOR OF LENGTH AT LEAST N; ON INPUT, AN INITIAL */
155  /* ESTIMATE OF THE SOLUTION; ON OUTPUT, THE COMPUTED SOLUTION. */
156  /* G - (REAL*8) VECTOR OF LENGTH AT LEAST N; ON OUTPUT, THE FINAL */
157  /* VALUE OF THE GRADIENT */
158  /* F - (REAL*8) ON INPUT, A ROUGH ESTIMATE OF THE VALUE OF THE */
159  /* OBJECTIVE FUNCTION AT THE SOLUTION; ON OUTPUT, THE VALUE */
160  /* OF THE OBJECTIVE FUNCTION AT THE SOLUTION */
161  /* W - (REAL*8) WORK VECTOR OF LENGTH AT LEAST 14*N */
162  /* LW - (C_INT) THE DECLARED DIMENSION OF W */
163  /* SFUN - A USER-SPECIFIED SUBROUTINE THAT COMPUTES THE FUNCTION */
164  /* AND GRADIENT OF THE OBJECTIVE FUNCTION. IT MUST HAVE */
165  /* THE CALLING SEQUENCE */
166  /* SUBROUTINE SFUN (N, X, F, G) */
167  /* C_INT N */
168  /* DOUBLE PRECISION X(N), G(N), F */
169 
170  /* THIS IS AN EASY-TO-USE DRIVER FOR THE MAIN OPTIMIZATION ROUTINE */
171  /* LMQN. MORE EXPERIENCED USERS WHO WISH TO CUSTOMIZE PERFORMANCE */
172  /* OF THIS ALGORITHM SHOULD CALL LMQN DIRECTLY. */
173 
174  /* ----------------------------------------------------------------------
175  */
176  /* THIS ROUTINE SETS UP ALL THE PARAMETERS FOR THE TRUNCATED-NEWTON */
177  /* ALGORITHM. THE PARAMETERS ARE: */
178 
179  /* ETA - SEVERITY OF THE LINESEARCH */
180  /* MAXFUN - MAXIMUM ALLOWABLE NUMBER OF FUNCTION EVALUATIONS */
181  /* XTOL - DESIRED ACCURACY FOR THE SOLUTION X* */
182  /* STEPMX - MAXIMUM ALLOWABLE STEP IN THE LINESEARCH */
183  /* ACCRCY - ACCURACY OF COMPUTED FUNCTION VALUES */
184  /* MSGLVL - DETERMINES QUANTITY OF PRINTED OUTPUT */
185  /* 0 = NONE, 1 = ONE LINE PER MAJOR ITERATION. */
186  /* MAXIT - MAXIMUM NUMBER OF INNER ITERATIONS PER STEP */
187 
188  /* SET UP PARAMETERS FOR THE OPTIMIZATION ROUTINE */
189 
190  /* Parameter adjustments */
191  --g;
192  --x;
193  --w;
194 
195  /* Function Body */
196  maxit = *n / 2;
197 
198  if (maxit > 50)
199  {
200  maxit = 50;
201  }
202 
203  if (maxit <= 0)
204  {
205  maxit = 1;
206  }
207 
208  msglvl = 0;
209  maxfun = *n * 150;
210  eta = .25;
211  stepmx = 10.;
212  accrcy = mchpr1_() * 100.;
213  xtol = sqrt(accrcy);
214 
215  /* MINIMIZE THE FUNCTION */
216 
217  CTruncatedNewton::lmqn_(ierror, n, &x[1], f, &g[1], &w[1], lw, sfun, &msglvl, &maxit,
218  &maxfun, &eta, &stepmx, &accrcy, &xtol);
219 
220  /* PRINT THE RESULTS */
221  //remove the printing
222  /* if (*ierror != 0) {
223  s_wsfe(&io___8);
224  do_fio(&c__1, (char *)&(*ierror), (ftnlen)sizeof(C_INT));
225  e_wsfe();
226  }
227  s_wsfe(&io___9);
228  do_fio(&c__1, (char *)&(*f), (ftnlen)sizeof(C_FLOAT64));
229  e_wsfe();
230  if (msglvl < 1) {
231  return 0;
232  }
233  s_wsfe(&io___10);
234  e_wsfe();
235  nmax = 10;
236  if (*n < nmax) {
237  nmax = *n;
238  }
239  s_wsfe(&io___12);
240  i__1 = nmax;
241  for (i__ = 1; i__ <= i__1; ++i__) {
242  do_fio(&c__1, (char *)&i__, (ftnlen)sizeof(C_INT));
243  do_fio(&c__1, (char *)&x[i__], (ftnlen)sizeof(C_FLOAT64));
244  }
245  e_wsfe();*/
246  return 0;
247 } /* tn_ */
#define C_INT
Definition: copasi.h:115
C_FLOAT64 mchpr1_(void)
#define C_FLOAT64
Definition: copasi.h:92
int lmqn_(C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, FTruncatedNewton *, C_INT *, C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *)
int CTruncatedNewton::tnbc_ ( C_INT ierror,
C_INT n,
C_FLOAT64 x,
C_FLOAT64 f,
C_FLOAT64 g,
C_FLOAT64 w,
C_INT lw,
FTruncatedNewton sfun,
C_FLOAT64 low,
C_FLOAT64 up,
C_INT ipivot 
)

Definition at line 249 of file CTruncatedNewton.cpp.

References C_FLOAT64, C_INT, lmqnbc_(), and mchpr1_().

Referenced by COptMethodTruncatedNewton::optimise().

252 {
253  /* Format strings */
254  /* char fmt_800[] = "(//,\002 ERROR CODE =\002,i3)";
255  char fmt_810[] = "(//,\002 OPTIMAL FUNCTION VALUE = \002,1pd22.15)"
256  ;
257  char fmt_820[] = "(10x,\002CURRENT SOLUTION IS (AT MOST 10 COMPON"
258  "ENTS)\002,/,14x,\002I\002,11x,\002X(I)\002)";
259  char fmt_830[] = "(10x,i5,2x,1pd22.15)";*/
260 
261  /* System generated locals */
262  C_INT i__1;
263 
264  // C_INT s_wsfe(cilist *), do_fio(C_INT *, char *, ftnlen), e_wsfe(void);
265 
266  /* Local variables */
267  C_INT nmax;
268  C_FLOAT64 xtol;
269  C_INT i__, maxit;
270  C_FLOAT64 accrcy;
271  C_INT maxfun, msglvl;
272  C_FLOAT64 stepmx, eta;
273 
274  /* Fortran I/O blocks */
275  //remove
276  /*
277  cilist io___21 = {0, 6, 0, fmt_800, 0 };
278  cilist io___22 = {0, 6, 0, fmt_810, 0 };
279  cilist io___23 = {0, 6, 0, fmt_820, 0 };
280  cilist io___25 = {0, 6, 0, fmt_830, 0 };*/
281 
282  /* THIS ROUTINE SOLVES THE OPTIMIZATION PROBLEM */
283 
284  /* MINIMIZE F(X) */
285  /* X */
286  /* SUBJECT TO LOW <= X <= UP */
287 
288  /* WHERE X IS A VECTOR OF N REAL VARIABLES. THE METHOD USED IS */
289  /* A TRUNCATED-NEWTON ALGORITHM (SEE "NEWTON-TYPE MINIMIZATION VIA */
290  /* THE LANCZOS ALGORITHM" BY S.G. NASH (TECHNICAL REPORT 378, MATH. */
291  /* THE LANCZOS METHOD" BY S.G. NASH (SIAM J. NUMER. ANAL. 21 (1984), */
292  /* PP. 770-778). THIS ALGORITHM FINDS A LOCAL MINIMUM OF F(X). IT DOES
293  */
294  /* NOT ASSUME THAT THE FUNCTION F IS CONVEX (AND SO CANNOT GUARANTEE A */
295  /* GLOBAL SOLUTION), BUT DOES ASSUME THAT THE FUNCTION IS BOUNDED BELOW.
296  */
297  /* IT CAN SOLVE PROBLEMS HAVING ANY NUMBER OF VARIABLES, BUT IT IS */
298  /* ESPECIALLY USEFUL WHEN THE NUMBER OF VARIABLES (N) IS LARGE. */
299 
300  /* SUBROUTINE PARAMETERS: */
301 
302  /* IERROR - (C_INT) ERROR CODE */
303  /* (0 => NORMAL RETURN */
304  /* (2 => MORE THAN MAXFUN EVALUATIONS */
305  /* (3 => LINE SEARCH FAILED TO FIND LOWER */
306  /* (POINT (MAY NOT BE SERIOUS) */
307  /* (-1 => ERROR IN INPUT PARAMETERS */
308  /* N - (C_INT) NUMBER OF VARIABLES */
309  /* X - (REAL*8) VECTOR OF LENGTH AT LEAST N; ON INPUT, AN INITIAL */
310  /* ESTIMATE OF THE SOLUTION; ON OUTPUT, THE COMPUTED SOLUTION.
311  */
312  /* G - (REAL*8) VECTOR OF LENGTH AT LEAST N; ON OUTPUT, THE FINAL */
313  /* VALUE OF THE GRADIENT */
314  /* F - (REAL*8) ON INPUT, A ROUGH ESTIMATE OF THE VALUE OF THE */
315  /* OBJECTIVE FUNCTION AT THE SOLUTION; ON OUTPUT, THE VALUE */
316  /* OF THE OBJECTIVE FUNCTION AT THE SOLUTION */
317  /* W - (REAL*8) WORK VECTOR OF LENGTH AT LEAST 14*N */
318  /* LW - (C_INT) THE DECLARED DIMENSION OF W */
319  /* SFUN - A USER-SPECIFIED SUBROUTINE THAT COMPUTES THE FUNCTION */
320  /* AND GRADIENT OF THE OBJECTIVE FUNCTION. IT MUST HAVE */
321  /* THE CALLING SEQUENCE */
322  /* SUBROUTINE SFUN (N, X, F, G) */
323  /* C_INT N */
324  /* DOUBLE PRECISION X(N), G(N), F */
325  /* LOW, UP - (REAL*8) VECTORS OF LENGTH AT LEAST N CONTAINING */
326  /* THE LOWER AND UPPER BOUNDS ON THE VARIABLES. IF */
327  /* THERE ARE NO BOUNDS ON A PARTICULAR VARIABLE, SET */
328  /* THE BOUNDS TO -1.D38 AND 1.D38, RESPECTIVELY. */
329  /* IPIVOT - (C_INT) WORK VECTOR OF LENGTH AT LEAST N, USED */
330  /* TO RECORD WHICH VARIABLES ARE AT THEIR BOUNDS. */
331 
332  /* THIS IS AN EASY-TO-USE DRIVER FOR THE MAIN OPTIMIZATION ROUTINE */
333  /* LMQNBC. MORE EXPERIENCED USERS WHO WISH TO CUSTOMIZE PERFORMANCE */
334  /* OF THIS ALGORITHM SHOULD CALL LMQBC DIRECTLY. */
335 
336  /* ----------------------------------------------------------------------
337  */
338  /* THIS ROUTINE SETS UP ALL THE PARAMETERS FOR THE TRUNCATED-NEWTON */
339  /* ALGORITHM. THE PARAMETERS ARE: */
340 
341  /* ETA - SEVERITY OF THE LINESEARCH */
342  /* MAXFUN - MAXIMUM ALLOWABLE NUMBER OF FUNCTION EVALUATIONS */
343  /* XTOL - DESIRED ACCURACY FOR THE SOLUTION X* */
344  /* STEPMX - MAXIMUM ALLOWABLE STEP IN THE LINESEARCH */
345  /* ACCRCY - ACCURACY OF COMPUTED FUNCTION VALUES */
346  /* MSGLVL - CONTROLS QUANTITY OF PRINTED OUTPUT */
347  /* 0 = NONE, 1 = ONE LINE PER MAJOR ITERATION. */
348  /* MAXIT - MAXIMUM NUMBER OF INNER ITERATIONS PER STEP */
349 
350  /* SET PARAMETERS FOR THE OPTIMIZATION ROUTINE */
351 
352  /* Parameter adjustments */
353  --ipivot;
354  --up;
355  --low;
356  --g;
357  --x;
358  --w;
359 
360  /* Function Body */
361  maxit = *n / 2;
362 
363  if (maxit > 50)
364  {
365  maxit = 50;
366  }
367 
368  if (maxit <= 0)
369  {
370  maxit = 1;
371  }
372 
373  msglvl = 0;
374  maxfun = *n * 150;
375  eta = .25;
376  stepmx = 10.;
377  accrcy = mchpr1_() * 100.;
378  xtol = sqrt(accrcy);
379 
380  /* MINIMIZE FUNCTION */
381 
382  CTruncatedNewton::lmqnbc_(ierror, n, &x[1], f, &g[1], &w[1], lw, sfun, &low[1], &up[1]
383  , &ipivot[1], &msglvl, &maxit, &maxfun, &eta, &stepmx, &accrcy, &
384  xtol);
385 
386  /* PRINT RESULTS */
387 
388  if (*ierror != 0)
389  {
390  // s_wsfe(&io___21);
391  // do_fio(&c__1, (char *)&(*ierror), (ftnlen)sizeof(C_INT));
392  // e_wsfe();
393  }
394 
395  // s_wsfe(&io___22);
396  // do_fio(&c__1, (char *)&(*f), (ftnlen)sizeof(C_FLOAT64));
397  // e_wsfe();
398  if (msglvl < 1)
399  {
400  return 0;
401  }
402 
403  // s_wsfe(&io___23);
404  // e_wsfe();
405  nmax = 10;
406 
407  if (*n < nmax)
408  {
409  nmax = *n;
410  }
411 
412  // s_wsfe(&io___25);
413  i__1 = nmax;
414 
415  for (i__ = 1; i__ <= i__1; ++i__)
416  {
417  // do_fio(&c__1, (char *)&i__, (ftnlen)sizeof(C_INT));
418  // do_fio(&c__1, (char *)&x[i__], (ftnlen)sizeof(C_FLOAT64));
419  }
420 
421  // e_wsfe();
422  return 0;
423 } /* tnbc_ */
#define C_INT
Definition: copasi.h:115
C_FLOAT64 mchpr1_(void)
int lmqnbc_(C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, FTruncatedNewton *, C_FLOAT64 *, C_FLOAT64 *, C_INT *, C_INT *, C_INT *, C_INT *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *, C_FLOAT64 *)
#define C_FLOAT64
Definition: copasi.h:92

Member Data Documentation

subscr_* CTruncatedNewton::mpsubscr_
private

Definition at line 160 of file CTruncatedNewton.h.

Referenced by CTruncatedNewton(), and ~CTruncatedNewton().


The documentation for this class was generated from the following files: