cloudy trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 #include "cddefines.h" 00004 #include "optimize.h" 00005 00006 STATIC void calcc(long,realnum*,long,long,int,realnum[]); 00007 STATIC double cdasum(long,realnum[],long); 00008 STATIC void cdaxpy(long,double,realnum[],long,realnum[],long); 00009 STATIC void cdcopy(long,realnum[],long,realnum[],long); 00010 STATIC void csscal(long,double,realnum[],long); 00011 STATIC double dist(long,realnum[],realnum[]); 00012 STATIC void evalf(long,long[],realnum[],long,realnum[],realnum*,long*); 00013 STATIC void fstats(double,long,int); 00014 STATIC void newpt(long,double,realnum[],realnum[],int,realnum[],int*); 00015 STATIC void order(long,realnum[],long*,long*,long*); 00016 STATIC void partx(long,long[],realnum[],long*,long[]); 00017 STATIC void setstp(long,long,realnum[],realnum[]); 00018 STATIC void simplx(long,realnum[],long,long[],long,int,realnum[], 00019 realnum*,long*,realnum*,realnum[],long*); 00020 STATIC void sortd(long,realnum[],long[]); 00021 STATIC void start(long,realnum[],realnum[],long,long[],realnum*,int*); 00022 STATIC void subopt(long); 00023 00024 /*lint -e801 goto is deprecated */ 00025 /*lint -e64 type mismatch */ 00026 /********************************************************************* 00027 * 00028 * NB this is much the original coding and does not use strong typing 00029 *gjf mod: to avoid conflict with exemplar parallel version of lapack, 00030 *dasum changed to cdasum 00031 *daxpy changed to cdaxpy 00032 *dcopy changed to cdcopy 00033 *sscal changed to csscal 00034 * 00035 *********************************************************************** */ 00036 struct t_usubc { 00037 realnum alpha, 00038 beta, 00039 gamma, 00040 delta, 00041 psi, 00042 omega; 00043 long int nsmin, 00044 nsmax, 00045 irepl, 00046 ifxsw; 00047 realnum bonus, 00048 fstop; 00049 long int nfstop, 00050 nfxe; 00051 realnum fxstat[4], 00052 ftest; 00053 int minf, 00054 initx, 00055 newx; 00056 } usubc; 00057 struct t_isubc { 00058 realnum fbonus, 00059 sfstop, 00060 sfbest; 00061 int IntNew; 00062 } isubc; 00063 00064 void optimize_subplex(long int n, 00065 double tol, 00066 long int maxnfe, 00067 long int mode, 00068 realnum scale[], 00069 realnum x[], 00070 realnum *fx, 00071 long int *nfe, 00072 realnum work[], 00073 long int iwork[], 00074 long int *iflag) 00075 { 00076 static int cmode; 00077 static long int i, 00078 i_, 00079 ifsptr, 00080 ins, 00081 insfnl, 00082 insptr, 00083 ipptr, 00084 isptr, 00085 istep, 00086 istptr, 00087 nnn, 00088 ns, 00089 nsubs; 00090 static realnum dum, 00091 scl[1], 00092 sfx, 00093 xpscl, 00094 xstop, 00095 xstop2; 00096 00097 static const realnum bnsfac[2][3] = { {-1.0f,-2.0f,0.0f}, {1.0f,0.0f,2.0f} }; 00098 00099 DEBUG_ENTRY( "optimize_subplex()" ); 00100 00101 /* Coded by Tom Rowan 00102 * Department of Computer Sciences 00103 * University of Texas at Austin 00104 * 00105 * Jason Ferguson: 00106 * Changed on 8/4/94, in order to incorparate into cloudy 00107 * changes made are: double precision to real, 00108 * a 2-D function f(n,x) into 1-D f(x), 00109 * and the termination check. 00110 * 00111 * optimize_subplex uses the subplex method to solve unconstrained 00112 * optimization problems. The method is well suited for 00113 * optimizing objective functions that are noisy or are 00114 * discontinuous at the solution. 00115 * 00116 * optimize_subplex sets default optimization options by calling the 00117 * subroutine subopt. The user can override these defaults 00118 * by calling subopt prior to calling optimize_subplex, changing the 00119 * appropriate common variables, and setting the value of 00120 * mode as indicated below. 00121 * 00122 * By default, optimize_subplex performs minimization. 00123 * 00124 * input 00125 * 00126 * ffff - user supplied function f(n,x) to be optimized, 00127 * declared external in calling routine 00128 * always uses optimize_func - change 97 dec 8 00129 * 00130 * n - problem dimension 00131 * 00132 * tol - relative error tolerance for x (tol .ge. 0.) 00133 * 00134 * maxnfe - maximum number of function evaluations 00135 * 00136 * mode - integer mode switch with binary expansion 00137 * (bit 1) (bit 0) : 00138 * bit 0 = 0 : first call to optimize_subplex 00139 * = 1 : continuation of previous call 00140 * bit 1 = 0 : use default options 00141 * = 1 : user set options 00142 * 00143 * scale - scale and initial stepsizes for corresponding 00144 * components of x 00145 * (If scale(1) .lt. 0., 00146 * ABS(scale(1)) is used for all components of x, 00147 * and scale(2),...,scale(n) are not referenced.) 00148 * 00149 * x - starting guess for optimum 00150 * 00151 * work - real work array of dimension .ge. 00152 * 2*n + nsmax*(nsmax+4) + 1 00153 * (nsmax is set in subroutine subopt. 00154 * default: nsmax = MIN(5,n)) 00155 * 00156 * iwork - integer work array of dimension .ge. 00157 * n + INT(n/nsmin) 00158 * (nsmin is set in subroutine subopt. 00159 * default: nsmin = MIN(2,n)) 00160 * 00161 * output 00162 * 00163 * x - computed optimum 00164 * 00165 * fx - value of f at x 00166 * 00167 * nfe - number of function evaluations 00168 * 00169 * iflag - error flag 00170 * = -2 : invalid input 00171 * = -1 : maxnfe exceeded 00172 * = 0 : tol satisfied 00173 * = 1 : limit of machine precision 00174 * = 2 : fstop reached (fstop usage is determined 00175 * by values of options minf, nfstop, and 00176 * irepl. default: f(x) not tested against 00177 * fstop) 00178 * iflag should not be reset between calls to 00179 * optimize_subplex. 00180 * 00181 * common 00182 * */ 00183 00184 if( (mode%2) == 0 ) 00185 { 00186 00187 /* first call, check input 00188 * */ 00189 if( n < 1 ) 00190 goto L_120; 00191 if( tol < 0.0f ) 00192 goto L_120; 00193 if( maxnfe < 1 ) 00194 goto L_120; 00195 if( scale[0] >= 0.0f ) 00196 { 00197 for( i=1; i <= n; i++ ) 00198 { 00199 i_ = i - 1; 00200 xpscl = x[i_] + scale[i_]; 00201 if( fp_equal( xpscl, x[i_] ) ) 00202 goto L_120; 00203 } 00204 } 00205 else 00206 { 00207 scl[0] = (realnum)fabs(scale[0]); 00208 for( i=1; i <= n; i++ ) 00209 { 00210 i_ = i - 1; 00211 xpscl = x[i_] + scl[0]; 00212 if( fp_equal( xpscl, x[i_] ) ) 00213 goto L_120; 00214 } 00215 } 00216 if( ((mode/2)%2) == 0 ) 00217 { 00218 subopt(n); 00219 } 00220 else if( usubc.alpha <= 0.0f ) 00221 { 00222 goto L_120; 00223 } 00224 else if( usubc.beta <= 0.0f || usubc.beta >= 1.0f ) 00225 { 00226 goto L_120; 00227 } 00228 else if( usubc.gamma <= 1.0f ) 00229 { 00230 goto L_120; 00231 } 00232 else if( usubc.delta <= 0.0f || usubc.delta >= 1.0f ) 00233 { 00234 goto L_120; 00235 } 00236 else if( usubc.psi <= 0.0f || usubc.psi >= 1.0f ) 00237 { 00238 goto L_120; 00239 } 00240 else if( usubc.omega <= 0.0f || usubc.omega >= 1.0f ) 00241 { 00242 goto L_120; 00243 } 00244 else if( (usubc.nsmin < 1 || usubc.nsmax < usubc.nsmin) || 00245 n < usubc.nsmax ) 00246 { 00247 goto L_120; 00248 } 00249 else if( n < ((n - 1)/usubc.nsmax + 1)*usubc.nsmin ) 00250 { 00251 goto L_120; 00252 } 00253 else if( usubc.irepl < 0 || usubc.irepl > 2 ) 00254 { 00255 goto L_120; 00256 } 00257 else if( usubc.ifxsw < 1 || usubc.ifxsw > 3 ) 00258 { 00259 goto L_120; 00260 } 00261 else if( usubc.bonus < 0.0f ) 00262 { 00263 goto L_120; 00264 } 00265 else if( usubc.nfstop < 0 ) 00266 { 00267 goto L_120; 00268 } 00269 00270 /* initialization 00271 * */ 00272 istptr = n + 1; 00273 isptr = istptr + n; 00274 ifsptr = isptr + usubc.nsmax*(usubc.nsmax + 3); 00275 insptr = n + 1; 00276 if( scale[0] > 0.0f ) 00277 { 00278 cdcopy(n,scale,1,work,1); 00279 cdcopy(n,scale,1,&work[istptr-1],1); 00280 } 00281 else 00282 { 00283 cdcopy(n,scl,0,work,1); 00284 cdcopy(n,scl,0,&work[istptr-1],1); 00285 } 00286 for( i=1; i <= n; i++ ) 00287 { 00288 i_ = i - 1; 00289 iwork[i_] = i; 00290 } 00291 *nfe = 0; 00292 usubc.nfxe = 1; 00293 if( usubc.irepl == 0 ) 00294 { 00295 isubc.fbonus = 0.0f; 00296 } 00297 else if( usubc.minf ) 00298 { 00299 isubc.fbonus = bnsfac[0][usubc.ifxsw-1]*usubc.bonus; 00300 } 00301 else 00302 { 00303 isubc.fbonus = bnsfac[1][usubc.ifxsw-1]*usubc.bonus; 00304 } 00305 if( usubc.nfstop == 0 ) 00306 { 00307 isubc.sfstop = 0.0f; 00308 } 00309 else if( usubc.minf ) 00310 { 00311 isubc.sfstop = usubc.fstop; 00312 } 00313 else 00314 { 00315 isubc.sfstop = -usubc.fstop; 00316 } 00317 usubc.ftest = 0.0f; 00318 cmode = false; 00319 isubc.IntNew = true; 00320 usubc.initx = true; 00321 nnn = 0; 00322 evalf(nnn,iwork,(realnum*)&dum,n,x,&sfx,nfe); 00323 usubc.initx = false; 00324 00325 /* continuation of previous call 00326 * */ 00327 } 00328 else if( *iflag == 2 ) 00329 { 00330 if( usubc.minf ) 00331 { 00332 isubc.sfstop = usubc.fstop; 00333 } 00334 else 00335 { 00336 isubc.sfstop = -usubc.fstop; 00337 } 00338 cmode = true; 00339 goto L_70; 00340 } 00341 else if( *iflag == -1 ) 00342 { 00343 cmode = true; 00344 goto L_70; 00345 } 00346 else if( *iflag == 0 ) 00347 { 00348 cmode = false; 00349 goto L_90; 00350 } 00351 else 00352 { 00353 return; 00354 } 00355 00356 /* subplex loop 00357 * */ 00358 L_40: 00359 00360 for( i=0; i < n; i++ ) 00361 { 00362 work[i] = (realnum)fabs(work[i]); 00363 } 00364 00365 sortd(n,work,iwork); 00366 partx(n,iwork,work,&nsubs,&iwork[insptr-1]); 00367 cdcopy(n,x,1,work,1); 00368 ins = insptr; 00369 insfnl = insptr + nsubs - 1; 00370 ipptr = 1; 00371 00372 /* simplex loop 00373 * */ 00374 00375 L_60: 00376 ns = iwork[ins-1]; 00377 00378 L_70: 00379 simplx(n,&work[istptr-1],ns,&iwork[ipptr-1],maxnfe,cmode,x,&sfx, 00380 nfe,&work[isptr-1],&work[ifsptr-1],iflag); 00381 00382 cmode = false; 00383 if( *iflag != 0 ) 00384 goto L_121; 00385 00386 if( ins < insfnl ) 00387 { 00388 ins += 1; 00389 ipptr += ns; 00390 goto L_60; 00391 } 00392 00393 /* end simplex loop 00394 * */ 00395 for( i=0; i < n; i++ ) 00396 { 00397 work[i] = x[i] - work[i]; 00398 } 00399 00400 /* check termination 00401 * */ 00402 L_90: 00403 /* new stop criterion: calculate distance between 00404 * previous and new best model, if distance is 00405 * smaller than tol return. */ 00406 istep = istptr-1; 00407 xstop = 0.f; 00408 xstop2 = 0.f; 00409 for( i=0; i < n; i++ ) 00410 { 00411 realnum ttemp; 00412 xstop += POW2(work[i]); 00413 ttemp = (realnum)fabs(work[istep]); 00414 /* chng to avoid side effects 00415 * xstop2 = MAX2(xstop2,(realnum)fabs(work[istep]));*/ 00416 xstop2 = MAX2( xstop2 , ttemp ); 00417 istep++; 00418 } 00419 00420 if( sqrt(xstop) > tol || xstop2*usubc.psi > (realnum)tol ) 00421 { 00422 setstp(nsubs,n,work,&work[istptr-1]); 00423 goto L_40; 00424 } 00425 00426 /* end subplex loop 00427 * */ 00428 00429 *iflag = 0; 00430 L_121: 00431 if( usubc.minf ) 00432 { 00433 *fx = sfx; 00434 } 00435 else 00436 { 00437 *fx = -sfx; 00438 } 00439 return; 00440 00441 /* invalid input 00442 * */ 00443 L_120: 00444 00445 *iflag = -2; 00446 return; 00447 } 00448 /********************************************************************* */ 00449 STATIC void subopt(long int n) 00450 { 00451 00452 DEBUG_ENTRY( "subopt()" ); 00453 00454 00455 /* Coded by Tom Rowan 00456 * Department of Computer Sciences 00457 * University of Texas at Austin 00458 * 00459 * subopt sets options for optimize_subplex. 00460 * 00461 * input 00462 * 00463 * n - problem dimension 00464 * 00465 * common 00466 * */ 00467 00468 00469 00470 /* subroutines and functions 00471 * 00472 * fortran 00473 * 00474 *----------------------------------------------------------- 00475 * 00476 ************************************************************ 00477 * simplex method strategy parameters 00478 ************************************************************ 00479 * 00480 * alpha - reflection coefficient 00481 * alpha .gt. 0 00482 * */ 00483 usubc.alpha = 1.0f; 00484 00485 /* beta - contraction coefficient 00486 * 0 .lt. beta .lt. 1 00487 * */ 00488 usubc.beta = .50f; 00489 00490 /* gamma - expansion coefficient 00491 * gamma .gt. 1 00492 * */ 00493 usubc.gamma = 2.0f; 00494 00495 /* delta - shrinkage (massive contraction) coefficient 00496 * 0 .lt. delta .lt. 1 00497 * */ 00498 usubc.delta = .50f; 00499 00500 /************************************************************ 00501 * subplex method strategy parameters 00502 ************************************************************ 00503 * 00504 * psi - simplex reduction coefficient 00505 * 0 .lt. psi .lt. 1 00506 * */ 00507 usubc.psi = .250f; 00508 00509 /* omega - step reduction coefficient 00510 * 0 .lt. omega .lt. 1 00511 * */ 00512 usubc.omega = .10f; 00513 00514 /* nsmin and nsmax specify a range of subspace dimensions. 00515 * In addition to satisfying 1 .le. nsmin .le. nsmax .le. n, 00516 * nsmin and nsmax must be chosen so that n can be expressed 00517 * as a sum of positive integers where each of these integers 00518 * ns(i) satisfies nsmin .le. ns(i) .ge. nsmax. 00519 * Specifically, 00520 * nsmin*ceil(n/nsmax) .le. n must be true. 00521 * 00522 * nsmin - subspace dimension minimum 00523 * */ 00524 usubc.nsmin = MIN2(2,n); 00525 00526 /* nsmax - subspace dimension maximum 00527 * */ 00528 usubc.nsmax = MIN2(5,n); 00529 00530 /************************************************************ 00531 * subplex method special cases 00532 ************************************************************ 00533 * nelder-mead simplex method with periodic restarts 00534 * nsmin = nsmax = n 00535 ************************************************************ 00536 * nelder-mead simplex method 00537 * nsmin = nsmax = n, psi = small positive 00538 ************************************************************ 00539 * 00540 * irepl, ifxsw, and bonus deal with measurement replication. 00541 * Objective functions subject to large amounts of noise can 00542 * cause an optimization method to halt at a false optimum. 00543 * An expensive solution to this problem is to evaluate f 00544 * several times at each point and return the average (or max 00545 * or min) of these trials as the function value. optimize_subplex 00546 * performs measurement replication only at the current best 00547 * point. The longer a point is retained as best, the more 00548 * accurate its function value becomes. 00549 * 00550 * The common variable nfxe contains the number of function 00551 * evaluations at the current best point. fxstat contains the 00552 * mean, max, min, and standard deviation of these trials. 00553 * 00554 * irepl - measurement replication switch 00555 * irepl = 0, 1, or 2 00556 * = 0 : no measurement replication 00557 * = 1 : optimize_subplex performs measurement replication 00558 * = 2 : user performs measurement replication 00559 * (This is useful when optimizing on the mean, 00560 * max, or min of trials is insufficient. Common 00561 * variable initx is true for first function 00562 * evaluation. newx is true for first trial at 00563 * this point. The user uses subroutine fstats 00564 * within his objective function to maintain 00565 * fxstat. By monitoring newx, the user can tell 00566 * whether to return the function evaluation 00567 * (newx = .true.) or to use the new function 00568 * evaluation to refine the function evaluation 00569 * of the current best point (newx = .false.). 00570 * The common variable ftest gives the function 00571 * value that a new point must beat to be 00572 * considered the new best point.) 00573 * */ 00574 usubc.irepl = 0; 00575 00576 /* ifxsw - measurement replication optimization switch 00577 * ifxsw = 1, 2, or 3 00578 * = 1 : retain mean of trials as best function value 00579 * = 2 : retain max 00580 * = 3 : retain min 00581 * */ 00582 usubc.ifxsw = 1; 00583 00584 /* Since the current best point will also be the most 00585 * accurately evaluated point whenever irepl .gt. 0, a bonus 00586 * should be added to the function value of the best point 00587 * so that the best point is not replaced by a new point 00588 * that only appears better because of noise. 00589 * optimize_subplex uses bonus to determine how many multiples of 00590 * fxstat(4) should be added as a bonus to the function 00591 * evaluation. (The bonus is adjusted automatically by 00592 * optimize_subplex when ifxsw or minf is changed.) 00593 * 00594 * bonus - measurement replication bonus coefficient 00595 * bonus .ge. 0 (normally, bonus = 0 or 1) 00596 * = 0 : bonus not used 00597 * = 1 : bonus used 00598 * */ 00599 usubc.bonus = 1.0f; 00600 00601 /* nfstop = 0 : f(x) is not tested against fstop 00602 * = 1 : if f(x) has reached fstop, optimize_subplex returns 00603 * iflag = 2 00604 * = 2 : (only valid when irepl .gt. 0) 00605 * if f(x) has reached fstop and 00606 * nfxe .gt. nfstop, optimize_subplex returns iflag = 2 00607 * */ 00608 usubc.nfstop = 0; 00609 00610 /* fstop - f target value 00611 * Its usage is determined by the value of nfstop. 00612 * 00613 * minf - logical switch 00614 * = .true. : optimize_subplex performs minimization 00615 * = .false. : optimize_subplex performs maximization 00616 * */ 00617 usubc.minf = true; 00618 return; 00619 } 00620 /**********************************************************************/ 00621 /* >>chng 01 jan 03, cleaned up -1 and formatting in this routine */ 00622 STATIC void cdcopy(long int n, 00623 realnum dx[], 00624 long int incx, 00625 realnum dy[], 00626 long int incy) 00627 { 00628 long int i, 00629 ix, 00630 iy, 00631 m; 00632 00633 DEBUG_ENTRY( "cdcopy()" ); 00634 00635 /* 00636 * copies a vector, x, to a vector, y. 00637 * uses unrolled loops for increments equal to one. 00638 * Jack Dongarra, lapack, 3/11/78. 00639 */ 00640 00641 if( n > 0 ) 00642 { 00643 if( incx == 1 && incy == 1 ) 00644 { 00645 00646 /* code for both increments equal to 1 */ 00647 00648 /* first the clean-up loop */ 00649 m = n%7; 00650 if( m != 0 ) 00651 { 00652 for( i=0; i < m; i++ ) 00653 { 00654 dy[i] = dx[i]; 00655 } 00656 if( n < 7 ) 00657 { 00658 return; 00659 } 00660 } 00661 00662 for( i=m; i < n; i += 7 ) 00663 { 00664 dy[i] = dx[i]; 00665 dy[i+1] = dx[i+1]; 00666 dy[i+2] = dx[i+2]; 00667 dy[i+3] = dx[i+3]; 00668 dy[i+4] = dx[i+4]; 00669 dy[i+5] = dx[i+5]; 00670 dy[i+6] = dx[i+6]; 00671 } 00672 } 00673 else 00674 { 00675 00676 /* code for unequal increments or equal increments 00677 * not equal to 1 */ 00678 00679 ix = 1; 00680 iy = 1; 00681 if( incx < 0 ) 00682 ix = (-n + 1)*incx + 1; 00683 if( incy < 0 ) 00684 iy = (-n + 1)*incy + 1; 00685 for( i=0; i < n; i++ ) 00686 { 00687 dy[iy-1] = dx[ix-1]; 00688 ix += incx; 00689 iy += incy; 00690 } 00691 } 00692 } 00693 return; 00694 } 00695 /********************************************************************* */ 00696 STATIC void evalf(long int ns, 00697 long int ips[], 00698 realnum xs[], 00699 long int n, 00700 realnum x[], 00701 realnum *sfx, 00702 long int *nfe) 00703 { 00704 static int newbst; 00705 static long int i, 00706 i_; 00707 static realnum fx; 00708 /* gary delete since in header */ 00709 /*double optimize_func();*/ 00710 00711 DEBUG_ENTRY( "evalf()" ); 00712 /* gary add, not used, so trick compiler notes */ 00713 i_ = n; 00714 00715 /*>>chng 97 dec 07, implicit nonte, rid of first function argument 00716 * 00717 * Coded by Tom Rowan 00718 * Department of Computer Sciences 00719 * University of Texas at Austin 00720 * 00721 * evalf evaluates the function f at a point defined by x 00722 * with ns of its components replaced by those in xs. 00723 * 00724 * input 00725 * 00726 * f - user supplied function f(n,x) to be optimized 00727 * changed to optimize_func - cannot specify arbitrary function now 00728 * 00729 * ns - subspace dimension 00730 * 00731 * ips - permutation vector 00732 * 00733 * xs - real ns-vector to be mapped to x 00734 * 00735 * n - problem dimension 00736 * 00737 * x - real n-vector 00738 * 00739 * nfe - number of function evaluations 00740 * 00741 * output 00742 * 00743 * sfx - signed value of f evaluated at x 00744 * 00745 * nfe - incremented number of function evaluations 00746 * 00747 * common 00748 * */ 00749 00750 00751 00752 00753 /* local variables 00754 * */ 00755 00756 00757 /* subroutines and functions 00758 * */ 00759 00760 /*----------------------------------------------------------- 00761 * */ 00762 for( i=1; i <= ns; i++ ) 00763 { 00764 i_ = i - 1; 00765 x[ips[i_]-1] = xs[i_]; 00766 } 00767 usubc.newx = isubc.IntNew || usubc.irepl != 2; 00768 fx = (realnum)optimize_func(x); 00769 /* fx = f(n,x) */ 00770 if( usubc.irepl == 0 ) 00771 { 00772 if( usubc.minf ) 00773 { 00774 *sfx = fx; 00775 } 00776 else 00777 { 00778 *sfx = -fx; 00779 } 00780 } 00781 else if( isubc.IntNew ) 00782 { 00783 if( usubc.minf ) 00784 { 00785 *sfx = fx; 00786 newbst = fx < usubc.ftest; 00787 } 00788 else 00789 { 00790 *sfx = -fx; 00791 newbst = fx > usubc.ftest; 00792 } 00793 if( usubc.initx || newbst ) 00794 { 00795 if( usubc.irepl == 1 ) 00796 fstats(fx,1,true); 00797 usubc.ftest = fx; 00798 isubc.sfbest = *sfx; 00799 } 00800 } 00801 else 00802 { 00803 if( usubc.irepl == 1 ) 00804 { 00805 fstats(fx,1,false); 00806 fx = usubc.fxstat[usubc.ifxsw-1]; 00807 } 00808 usubc.ftest = fx + isubc.fbonus*usubc.fxstat[3]; 00809 if( usubc.minf ) 00810 { 00811 *sfx = usubc.ftest; 00812 isubc.sfbest = fx; 00813 } 00814 else 00815 { 00816 *sfx = -usubc.ftest; 00817 isubc.sfbest = -fx; 00818 } 00819 } 00820 *nfe += 1; 00821 return; 00822 } 00823 00824 /********************************************************************* */ 00825 STATIC void setstp(long int nsubs, 00826 long int n, 00827 realnum deltax[], 00828 realnum step[]) 00829 { 00830 static long int i, 00831 i_; 00832 static realnum stpfac; 00833 00834 DEBUG_ENTRY( "setstp()" ); 00835 00836 00837 /* Coded by Tom Rowan 00838 * Department of Computer Sciences 00839 * University of Texas at Austin 00840 * 00841 * setstp sets the stepsizes for the corresponding components 00842 * of the solution vector. 00843 * 00844 * input 00845 * 00846 * nsubs - number of subspaces 00847 * 00848 * n - number of components (problem dimension) 00849 * 00850 * deltax - vector of change in solution vector 00851 * 00852 * step - stepsizes for corresponding components of 00853 * solution vector 00854 * 00855 * output 00856 * 00857 * step - new stepsizes 00858 * 00859 * common 00860 * */ 00861 00862 00863 /* local variables 00864 * */ 00865 00866 00867 /* subroutines and functions 00868 * 00869 * blas */ 00870 /* fortran 00871 * 00872 *----------------------------------------------------------- 00873 * 00874 * set new step 00875 * */ 00876 if( nsubs > 1 ) 00877 { 00878 double a , b , c; 00879 a = cdasum(n,deltax,1); 00880 b = cdasum(n,step,1); 00881 c = MAX2(a/b ,usubc.omega); 00882 stpfac = (realnum)MIN2(c , 1.0f/usubc.omega); 00883 } 00884 else 00885 { 00886 stpfac = usubc.psi; 00887 } 00888 csscal(n,stpfac,step,1); 00889 00890 /* reorient simplex 00891 * */ 00892 for( i=1; i <= n; i++ ) 00893 { 00894 i_ = i - 1; 00895 if( deltax[i_] == 0.f ) 00896 { 00897 step[i_] = -step[i_]; 00898 } 00899 else 00900 { 00901 step[i_] = (realnum)sign(step[i_],deltax[i_]); 00902 } 00903 } 00904 return; 00905 } 00906 00907 /********************************************************************* */ 00908 STATIC void sortd(long int n, 00909 realnum xkey[], 00910 long int ix[]) 00911 { 00912 long int i, 00913 i_, 00914 ifirst, 00915 ilast, 00916 iswap, 00917 ixi, 00918 ixip1; 00919 00920 DEBUG_ENTRY( "sortd()" ); 00921 00922 00923 /* Coded by Tom Rowan 00924 * Department of Computer Sciences 00925 * University of Texas at Austin 00926 * 00927 * sortd uses the shakersort method to sort an array of keys 00928 * in decreasing order. The sort is performed implicitly by 00929 * modifying a vector of indices. 00930 * 00931 * For nearly sorted arrays, sortd requires O(n) comparisons. 00932 * for completely unsorted arrays, sortd requires O(n**2) 00933 * comparisons and will be inefficient unless n is small. 00934 * 00935 * input 00936 * 00937 * n - number of components 00938 * 00939 * xkey - real vector of keys 00940 * 00941 * ix - integer vector of indices 00942 * 00943 * output 00944 * 00945 * ix - indices satisfy xkey(ix(i)) .ge. xkey(ix(i+1)) 00946 * for i = 1,...,n-1 00947 * 00948 * local variables 00949 * */ 00950 00951 /*----------------------------------------------------------- 00952 * */ 00953 ifirst = 1; 00954 iswap = 1; 00955 ilast = n - 1; 00956 while( ifirst <= ilast ) 00957 { 00958 for( i=ifirst; i <= ilast; i++ ) 00959 { 00960 i_ = i - 1; 00961 ixi = ix[i_]; 00962 ixip1 = ix[i_+1]; 00963 if( xkey[ixi-1] < xkey[ixip1-1] ) 00964 { 00965 ix[i_] = ixip1; 00966 ix[i_+1] = ixi; 00967 iswap = i; 00968 } 00969 } 00970 ilast = iswap - 1; 00971 for( i=ilast; i >= ifirst; i-- ) 00972 { 00973 i_ = i - 1; 00974 ixi = ix[i_]; 00975 ixip1 = ix[i_+1]; 00976 if( xkey[ixi-1] < xkey[ixip1-1] ) 00977 { 00978 ix[i_] = ixip1; 00979 ix[i_+1] = ixi; 00980 iswap = i; 00981 } 00982 } 00983 ifirst = iswap + 1; 00984 } 00985 return; 00986 } 00987 00988 /********************************************************************* */ 00989 STATIC void partx(long int n, 00990 long int ip[], 00991 realnum absdx[], 00992 long int *nsubs, 00993 long int nsvals[]) 00994 { 00995 static long int i, 00996 limit, 00997 nleft, 00998 ns1, 00999 ns1_, 01000 ns2, 01001 nused; 01002 static realnum as1, 01003 as1max, 01004 as2, 01005 asleft, 01006 gap, 01007 gapmax; 01008 01009 DEBUG_ENTRY( "partx()" ); 01010 01011 01012 /* Coded by Tom Rowan 01013 * Department of Computer Sciences 01014 * University of Texas at Austin 01015 * 01016 * partx partitions the vector x by grouping components of 01017 * similar magnitude of change. 01018 * 01019 * input 01020 * 01021 * n - number of components (problem dimension) 01022 * 01023 * ip - permutation vector 01024 * 01025 * absdx - vector of magnitude of change in x 01026 * 01027 * nsvals - integer array dimensioned .ge. INT(n/nsmin) 01028 * 01029 * output 01030 * 01031 * nsubs - number of subspaces 01032 * 01033 * nsvals - integer array of subspace dimensions 01034 * 01035 * common 01036 * */ 01037 01038 01039 /* local variables 01040 * */ 01041 01042 01043 /* subroutines and functions 01044 * 01045 * 01046 *----------------------------------------------------------- 01047 * */ 01048 *nsubs = 0; 01049 nused = 0; 01050 nleft = n; 01051 asleft = absdx[0]; 01052 for( i=1; i < n; i++ ) 01053 { 01054 asleft += absdx[i]; 01055 } 01056 01057 while( nused < n ) 01058 { 01059 *nsubs += 1; 01060 as1 = 0.0f; 01061 for( i=0; i < (usubc.nsmin - 1); i++ ) 01062 { 01063 as1 += absdx[ip[nused+i]-1]; 01064 } 01065 01066 gapmax = -1.0f; 01067 limit = MIN2(usubc.nsmax,nleft); 01068 for( ns1=usubc.nsmin; ns1 <= limit; ns1++ ) 01069 { 01070 ns1_ = ns1 - 1; 01071 as1 += absdx[ip[nused+ns1_]-1]; 01072 ns2 = nleft - ns1; 01073 if( ns2 > 0 ) 01074 { 01075 if( ns2 >= ((ns2 - 1)/usubc.nsmax + 1)*usubc.nsmin ) 01076 { 01077 as2 = asleft - as1; 01078 gap = as1/ns1 - as2/ns2; 01079 if( gap > gapmax ) 01080 { 01081 gapmax = gap; 01082 nsvals[*nsubs-1] = ns1; 01083 as1max = as1; 01084 } 01085 } 01086 } 01087 else if( as1/ns1 > gapmax ) 01088 { 01089 goto L_21; 01090 } 01091 } 01092 nused += nsvals[*nsubs-1]; 01093 nleft = n - nused; 01094 asleft -= as1max; 01095 } 01096 return; 01097 L_21: 01098 nsvals[*nsubs-1] = ns1; 01099 return; 01100 } 01101 01102 /********************************************************************* */ 01103 01104 #undef S 01105 #define S(I_,J_) (*(s+(I_)*(ns)+(J_))) 01106 01107 STATIC void simplx(long int n, 01108 realnum step[], 01109 long int ns, 01110 long int ips[], 01111 long int maxnfe, 01112 int cmode, 01113 realnum x[], 01114 realnum *fx, 01115 long int *nfe, 01116 realnum *s, 01117 realnum fs[], 01118 long int *iflag) 01119 { 01120 static int small, 01121 updatc; 01122 01123 static long int i, 01124 icent, 01125 ih, 01126 il, 01127 inew = 0, 01128 is, 01129 itemp, 01130 j, 01131 j_, 01132 npts; 01133 01134 static realnum dum, 01135 fc, 01136 fe, 01137 fr, 01138 tol; 01139 01140 DEBUG_ENTRY( "simplx()" ); 01141 01142 01143 /* Coded by Tom Rowan 01144 * Department of Computer Sciences 01145 * University of Texas at Austin 01146 * 01147 * simplx uses the Nelder-Mead simplex method to minimize the 01148 * function f on a subspace. 01149 * 01150 * input 01151 * 01152 * ffff - function to be minimized, declared external in 01153 * calling routine 01154 * removed - always calls optimize_func 01155 * 01156 * n - problem dimension 01157 * 01158 * step - stepsizes for corresponding components of x 01159 * 01160 * ns - subspace dimension 01161 * 01162 * ips - permutation vector 01163 * 01164 * maxnfe - maximum number of function evaluations 01165 * 01166 * cmode - logical switch 01167 * = .true. : continuation of previous call 01168 * = .false. : first call 01169 * 01170 * x - starting guess for minimum 01171 * 01172 * fx - value of f at x 01173 * 01174 * nfe - number of function evaluations 01175 * 01176 * s - real work array of dimension .ge. 01177 * ns*(ns+3) used to store simplex 01178 * 01179 * fs - real work array of dimension .ge. 01180 * ns+1 used to store function values of simplex 01181 * vertices 01182 * 01183 * output 01184 * 01185 * x - computed minimum 01186 * 01187 * fx - value of f at x 01188 * 01189 * nfe - incremented number of function evaluations 01190 * 01191 * iflag - error flag 01192 * = -1 : maxnfe exceeded 01193 * = 0 : simplex reduced by factor of psi 01194 * = 1 : limit of machine precision 01195 * = 2 : reached fstop 01196 * 01197 * common 01198 * */ 01199 01200 01201 01202 01203 /* local variables 01204 * */ 01205 01206 01207 /* subroutines and functions 01208 * 01209 * external optimize_func,calcc,dist,evalf,newpt,order,start */ 01210 /* blas */ 01211 01212 /*----------------------------------------------------------- 01213 * */ 01214 if( cmode ) 01215 goto L_50; 01216 npts = ns + 1; 01217 icent = ns + 2; 01218 itemp = ns + 3; 01219 updatc = false; 01220 start(n,x,step,ns,ips,s,&small); 01221 if( small ) 01222 { 01223 *iflag = 1; 01224 return; 01225 } 01226 else 01227 { 01228 if( usubc.irepl > 0 ) 01229 { 01230 isubc.IntNew = false; 01231 evalf(ns,ips,&S(0,0),n,x,&fs[0],nfe); 01232 } 01233 else 01234 { 01235 fs[0] = *fx; 01236 } 01237 isubc.IntNew = true; 01238 for( j=2; j <= npts; j++ ) 01239 { 01240 j_ = j - 1; 01241 evalf(ns,ips,&S(j_,0),n,x,&fs[j_],nfe); 01242 } 01243 il = 1; 01244 order(npts,fs,&il,&is,&ih); 01245 tol = (realnum)(usubc.psi*dist(ns,&S(ih-1,0),&S(il-1,0))); 01246 } 01247 01248 /* main loop 01249 * */ 01250 L_20: 01251 calcc(ns,s,ih,inew,updatc,&S(icent-1,0)); 01252 updatc = true; 01253 inew = ih; 01254 01255 /* reflect 01256 * */ 01257 newpt(ns,usubc.alpha,&S(icent-1,0),&S(ih-1,0),true,&S(itemp-1,0), 01258 &small); 01259 if( !small ) 01260 { 01261 evalf(ns,ips,&S(itemp-1,0),n,x,&fr,nfe); 01262 if( fr < fs[il-1] ) 01263 { 01264 01265 /* expand 01266 * */ 01267 newpt(ns,-usubc.gamma,&S(icent-1,0),&S(itemp-1,0),true, 01268 &S(ih-1,0),&small); 01269 if( small ) 01270 goto L_40; 01271 evalf(ns,ips,&S(ih-1,0),n,x,&fe,nfe); 01272 if( fe < fr ) 01273 { 01274 fs[ih-1] = fe; 01275 } 01276 else 01277 { 01278 cdcopy(ns,&S(itemp-1,0),1,&S(ih-1,0),1); 01279 fs[ih-1] = fr; 01280 } 01281 } 01282 else if( fr < fs[is-1] ) 01283 { 01284 01285 /* accept reflected point 01286 * */ 01287 cdcopy(ns,&S(itemp-1,0),1,&S(ih-1,0),1); 01288 fs[ih-1] = fr; 01289 } 01290 else 01291 { 01292 01293 /* contract 01294 * */ 01295 if( fr > fs[ih-1] ) 01296 { 01297 newpt(ns,-usubc.beta,&S(icent-1,0),&S(ih-1,0),true, 01298 &S(itemp-1,0),&small); 01299 } 01300 else 01301 { 01302 newpt(ns,-usubc.beta,&S(icent-1,0),&S(itemp-1,0),false, 01303 (realnum*)&dum,&small); 01304 } 01305 if( small ) 01306 goto L_40; 01307 evalf(ns,ips,&S(itemp-1,0),n,x,&fc,nfe); 01308 if( fc < (realnum)MIN2(fr,fs[ih-1]) ) 01309 { 01310 cdcopy(ns,&S(itemp-1,0),1,&S(ih-1,0),1); 01311 fs[ih-1] = fc; 01312 } 01313 else 01314 { 01315 01316 /* shrink simplex 01317 * */ 01318 for( j=1; j <= npts; j++ ) 01319 { 01320 j_ = j - 1; 01321 if( j != il ) 01322 { 01323 newpt(ns,-usubc.delta,&S(il-1,0),&S(j_,0), 01324 false,(realnum*)&dum,&small); 01325 if( small ) 01326 goto L_40; 01327 evalf(ns,ips,&S(j_,0),n,x,&fs[j_],nfe); 01328 } 01329 } 01330 } 01331 updatc = false; 01332 } 01333 order(npts,fs,&il,&is,&ih); 01334 } 01335 01336 /* check termination 01337 * */ 01338 01339 L_40: 01340 if( usubc.irepl == 0 ) 01341 { 01342 *fx = fs[il-1]; 01343 } 01344 else 01345 { 01346 *fx = isubc.sfbest; 01347 } 01348 01349 L_50: 01350 if( (usubc.nfstop > 0 && *fx <= isubc.sfstop) && usubc.nfxe >= 01351 usubc.nfstop ) 01352 goto L_51; 01353 if( *nfe >= maxnfe ) 01354 goto L_52; 01355 if( !(dist(ns,&S(ih-1,0),&S(il-1,0)) <= tol || small) ) 01356 goto L_20; 01357 *iflag = 0; 01358 goto L_53; 01359 01360 L_52: 01361 *iflag = -1; 01362 goto L_53; 01363 01364 L_51: 01365 *iflag = 2; 01366 01367 /* end main loop, return best point 01368 * */ 01369 01370 L_53: 01371 for( i=0; i < ns; i++ ) 01372 { 01373 x[ips[i]-1] = S(il-1,i); 01374 } 01375 return; 01376 } 01377 01378 /********************************************************************* */ 01379 STATIC void fstats(double fx, 01380 long int ifxwt, 01381 int reset) 01382 { 01383 static long int nsv; 01384 static realnum f1sv, 01385 fscale; 01386 01387 DEBUG_ENTRY( "fstats()" ); 01388 01389 01390 /* Coded by Tom Rowan 01391 * Department of Computer Sciences 01392 * University of Texas at Austin 01393 * 01394 * fstats modifies the common /usubc/ variables nfxe,fxstat. 01395 * 01396 * input 01397 * 01398 * fx - most recent evaluation of f at best x 01399 * 01400 * ifxwt - integer weight for fx 01401 * 01402 * reset - logical switch 01403 * = .true. : initialize nfxe,fxstat 01404 * = .false. : update nfxe,fxstat 01405 * 01406 * common 01407 * */ 01408 01409 01410 /* local variables 01411 * */ 01412 01413 01414 /* subroutines and functions 01415 * 01416 * 01417 *----------------------------------------------------------- 01418 * */ 01419 if( reset ) 01420 { 01421 usubc.nfxe = ifxwt; 01422 usubc.fxstat[0] = (realnum)fx; 01423 usubc.fxstat[1] = (realnum)fx; 01424 usubc.fxstat[2] = (realnum)fx; 01425 usubc.fxstat[3] = 0.0f; 01426 } 01427 else 01428 { 01429 nsv = usubc.nfxe; 01430 f1sv = usubc.fxstat[0]; 01431 usubc.nfxe += ifxwt; 01432 usubc.fxstat[0] += (realnum)(ifxwt*(fx - usubc.fxstat[0])/usubc.nfxe); 01433 usubc.fxstat[1] = MAX2(usubc.fxstat[1],(realnum)fx); 01434 usubc.fxstat[2] = MIN2(usubc.fxstat[2],(realnum)fx); 01435 fscale = (realnum)MAX3(fabs(usubc.fxstat[1]),fabs(usubc.fxstat[2]),1.); 01436 usubc.fxstat[3] = (realnum)(fscale*sqrt(((nsv-1)*POW2(usubc.fxstat[3]/ 01437 fscale)+nsv*POW2((usubc.fxstat[0]-f1sv)/fscale)+ifxwt* 01438 POW2((fx-usubc.fxstat[0])/fscale))/(usubc.nfxe-1))); 01439 } 01440 return; 01441 } 01442 01443 STATIC double cdasum(long int n, 01444 realnum dx[], 01445 long int incx) 01446 { 01447 /* 01448 * 01449 * takes the sum of the absolute values. 01450 * uses unrolled loops for increment equal to one. 01451 * jack dongarra, lapack, 3/11/78. 01452 * modified to correct problem with negative increment, 8/21/90. 01453 * 01454 */ 01455 long int i, 01456 ix, 01457 m; 01458 realnum cdasum_v, 01459 dtemp; 01460 01461 DEBUG_ENTRY( "cdasum()" ); 01462 01463 cdasum_v = 0.00f; 01464 dtemp = 0.00f; 01465 if( n > 0 ) 01466 { 01467 if( incx == 1 ) 01468 { 01469 01470 /* code for increment equal to 1 01471 * 01472 * 01473 * clean-up loop 01474 * */ 01475 m = n%6; 01476 if( m != 0 ) 01477 { 01478 for( i=0; i < m; i++ ) 01479 { 01480 dtemp += (realnum)fabs(dx[i]); 01481 } 01482 if( n < 6 ) 01483 goto L_60; 01484 } 01485 01486 for( i=m; i < n; i += 6 ) 01487 { 01488 dtemp += (realnum)(fabs(dx[i]) + fabs(dx[i+1]) + fabs(dx[i+2]) + 01489 fabs(dx[i+3]) + fabs(dx[i+4]) + fabs(dx[i+5])); 01490 } 01491 L_60: 01492 cdasum_v = dtemp; 01493 } 01494 else 01495 { 01496 01497 /* code for increment not equal to 1 01498 * */ 01499 ix = 1; 01500 01501 if( incx < 0 ) 01502 ix = (-n + 1)*incx + 1; 01503 01504 for( i=0; i < n; i++ ) 01505 { 01506 dtemp += (realnum)fabs(dx[ix-1]); 01507 ix += incx; 01508 } 01509 cdasum_v = dtemp; 01510 } 01511 } 01512 return( cdasum_v ); 01513 } 01514 01515 STATIC void csscal(long int n, 01516 double da, 01517 realnum dx[], 01518 long int incx) 01519 { 01520 long int 01521 i, 01522 i_, 01523 m, 01524 mp1, 01525 nincx; 01526 01527 DEBUG_ENTRY( "csscal()" ); 01528 01529 /* scales a vector by a constant. 01530 * uses unrolled loops for increment equal to one. 01531 * jack dongarra, lapack, 3/11/78. 01532 * modified 3/93 to return if incx .le. 0. 01533 * modified 12/3/93, array(1) declarations changed to array(*) 01534 * changed to single precisions 01535 * */ 01536 01537 if( !(n <= 0 || incx <= 0) ) 01538 { 01539 if( incx == 1 ) 01540 { 01541 01542 /* code for increment equal to 1 01543 * 01544 * 01545 * clean-up loop 01546 * */ 01547 m = n%5; 01548 if( m != 0 ) 01549 { 01550 for( i=1; i <= m; i++ ) 01551 { 01552 i_ = i - 1; 01553 dx[i_] *= (realnum)da; 01554 } 01555 if( n < 5 ) 01556 { 01557 return; 01558 } 01559 } 01560 mp1 = m + 1; 01561 for( i=mp1; i <= n; i += 5 ) 01562 { 01563 i_ = i - 1; 01564 dx[i_] *= (realnum)da; 01565 dx[i_+1] *= (realnum)da; 01566 dx[i_+2] *= (realnum)da; 01567 dx[i_+3] *= (realnum)da; 01568 dx[i_+4] *= (realnum)da; 01569 } 01570 } 01571 else 01572 { 01573 01574 /* code for increment not equal to 1 01575 * */ 01576 nincx = n*incx; 01577 /*for( i=1, _do0=DOCNT(i,nincx,_do1 = incx); _do0 > 0; i += _do1, _do0-- )*/ 01578 /* gary change forc */ 01579 for( i=0; i<nincx; i=i+incx) 01580 /*for( i=1, _do0=DOCNT(i,nincx,_do1 = incx); _do0 > 0; i += _do1, _do0-- )*/ 01581 { 01582 dx[i] *= (realnum)da; 01583 } 01584 } 01585 } 01586 return; 01587 } 01588 01589 /********************************************************************* */ 01590 01591 #undef S 01592 #define S(I_,J_) (*(s+(I_)*(ns)+(J_))) 01593 01594 STATIC void start(long int n, 01595 realnum x[], 01596 realnum step[], 01597 long int ns, 01598 long int ips[], 01599 realnum *s, 01600 int *small) 01601 { 01602 long int i, 01603 i_, 01604 j, 01605 j_; 01606 01607 DEBUG_ENTRY( "start()" ); 01608 01609 /* gary add, not used so set to i_ to mute compiler note */ 01610 i_ = n; 01611 01612 01613 /* Coded by Tom Rowan 01614 * Department of Computer Sciences 01615 * University of Texas at Austin 01616 * 01617 * start creates the initial simplex for simplx minimization. 01618 * 01619 * input 01620 * 01621 * n - problem dimension 01622 * 01623 * x - current best point 01624 * 01625 * step - stepsizes for corresponding components of x 01626 * 01627 * ns - subspace dimension 01628 * 01629 * ips - permutation vector 01630 * 01631 * 01632 * output 01633 * 01634 * s - first ns+1 columns contain initial simplex 01635 * 01636 * small - logical flag 01637 * = .true. : coincident points 01638 * = .false. : otherwise 01639 * 01640 * local variables 01641 * */ 01642 01643 /* subroutines and functions 01644 * 01645 * blas */ 01646 /* fortran */ 01647 01648 /*----------------------------------------------------------- 01649 * */ 01650 for( i=1; i <= ns; i++ ) 01651 { 01652 i_ = i - 1; 01653 S(0,i_) = x[ips[i_]-1]; 01654 } 01655 01656 for( j=2; j <= (ns + 1); j++ ) 01657 { 01658 j_ = j - 1; 01659 cdcopy(ns,&S(0,0),1,&S(j_,0),1); 01660 S(j_,j_-1) = S(0,j_-1) + step[ips[j_-1]-1]; 01661 } 01662 01663 /* check for coincident points 01664 * */ 01665 for( j=2; j <= (ns + 1); j++ ) 01666 { 01667 j_ = j - 1; 01668 if( (double)(S(j_,j_-1)) == (double)(S(0,j_-1)) ) 01669 goto L_40; 01670 } 01671 *small = false; 01672 return; 01673 01674 /* coincident points 01675 * */ 01676 L_40: 01677 *small = true; 01678 return; 01679 } 01680 01681 /********************************************************************* */ 01682 STATIC void order(long int npts, 01683 realnum fs[], 01684 long int *il, 01685 long int *is, 01686 long int *ih) 01687 { 01688 long int i, 01689 il0, 01690 j; 01691 01692 DEBUG_ENTRY( "order()" ); 01693 01694 01695 /* Coded by Tom Rowan 01696 * Department of Computer Sciences 01697 * University of Texas at Austin 01698 * 01699 * order determines the indices of the vertices with the 01700 * lowest, second highest, and highest function values. 01701 * 01702 * input 01703 * 01704 * npts - number of points in simplex 01705 * 01706 * fs - real vector of function values of 01707 * simplex 01708 * 01709 * il - index to vertex with lowest function value 01710 * 01711 * output 01712 * 01713 * il - new index to vertex with lowest function value 01714 * 01715 * is - new index to vertex with second highest 01716 * function value 01717 * 01718 * ih - new index to vertex with highest function value 01719 * 01720 * local variables 01721 * */ 01722 01723 /* subroutines and functions 01724 * 01725 * 01726 *----------------------------------------------------------- 01727 * */ 01728 il0 = *il; 01729 j = (il0%npts) + 1; 01730 01731 if( fs[j-1] >= fs[*il-1] ) 01732 { 01733 *ih = j; 01734 *is = il0; 01735 } 01736 else 01737 { 01738 *ih = il0; 01739 *is = j; 01740 *il = j; 01741 } 01742 01743 for( i=il0 + 1; i <= (il0 + npts - 2); i++ ) 01744 { 01745 j = (i%npts) + 1; 01746 if( fs[j-1] >= fs[*ih-1] ) 01747 { 01748 *is = *ih; 01749 *ih = j; 01750 } 01751 else if( fs[j-1] > fs[*is-1] ) 01752 { 01753 *is = j; 01754 } 01755 else if( fs[j-1] < fs[*il-1] ) 01756 { 01757 *il = j; 01758 } 01759 } 01760 return; 01761 } 01762 01763 STATIC double dist(long int n, 01764 realnum x[], 01765 realnum y[]) 01766 { 01767 /* 01768 * 01769 */ 01770 long int i, 01771 i_; 01772 realnum absxmy, 01773 dist_v, 01774 scale, 01775 sum; 01776 01777 DEBUG_ENTRY( "dist()" ); 01778 01779 /* Coded by Tom Rowan 01780 * Department of Computer Sciences 01781 * University of Texas at Austin 01782 * 01783 * dist calculates the distance between the points x,y. 01784 * 01785 * input 01786 * 01787 * n - number of components 01788 * 01789 * x - point in n-space 01790 * 01791 * y - point in n-space 01792 * 01793 * local variables 01794 * */ 01795 01796 /* subroutines and functions 01797 * 01798 * fortran 01799 * 01800 *----------------------------------------------------------- 01801 * */ 01802 absxmy = (realnum)fabs(x[0]-y[0]); 01803 if( absxmy <= 1.0f ) 01804 { 01805 sum = absxmy*absxmy; 01806 scale = 1.0f; 01807 } 01808 else 01809 { 01810 sum = 1.0f; 01811 scale = absxmy; 01812 } 01813 01814 for( i=2; i <= n; i++ ) 01815 { 01816 i_ = i - 1; 01817 absxmy = (realnum)fabs(x[i_]-y[i_]); 01818 if( absxmy <= scale ) 01819 { 01820 sum += POW2(absxmy/scale); 01821 } 01822 else 01823 { 01824 sum = 1.0f + sum*POW2(scale/absxmy); 01825 scale = absxmy; 01826 } 01827 } 01828 dist_v = (realnum)(scale*sqrt(sum)); 01829 return( dist_v ); 01830 } 01831 01832 /********************************************************************* */ 01833 01834 #undef S 01835 #define S(I_,J_) (*(s+(I_)*(ns)+(J_))) 01836 01837 STATIC void calcc(long int ns, 01838 realnum *s, 01839 long int ih, 01840 long int inew, 01841 int updatc, 01842 realnum c[]) 01843 { 01844 long int i, 01845 j, 01846 j_; 01847 /* >>chng 99 apr 21, following was not initialized, caught by Peter van Hoof */ 01848 realnum xNothing[1] = { 0.0f }; 01849 01850 DEBUG_ENTRY( "calcc()" ); 01851 01852 01853 /* Coded by Tom Rowan 01854 * Department of Computer Sciences 01855 * University of Texas at Austin 01856 * 01857 * calcc calculates the centroid of the simplex without the 01858 * vertex with highest function value. 01859 * 01860 * input 01861 * 01862 * ns - subspace dimension 01863 * 01864 * s - real work space of dimension .ge. 01865 * ns*(ns+3) used to store simplex 01866 * 01867 * ih - index to vertex with highest function value 01868 * 01869 * inew - index to new point 01870 * 01871 * updatc - logical switch 01872 * = .true. : update centroid 01873 * = .false. : calculate centroid from scratch 01874 * 01875 * c - centroid of the simplex without vertex with 01876 * highest function value 01877 * 01878 * output 01879 * 01880 * c - new centroid 01881 * 01882 * local variables 01883 * */ 01884 /* added to get prototypes to work */ 01885 01886 /* subroutines and functions 01887 * 01888 * blas */ 01889 01890 /*----------------------------------------------------------- 01891 * */ 01892 if( !updatc ) 01893 { 01894 /* call cdcopy (ns,0.0,0,c,1) 01895 * xNothing will not be used since 0 is increment */ 01896 cdcopy(ns,xNothing,0,c,1); 01897 for( j=1; j <= (ns + 1); j++ ) 01898 { 01899 j_ = j - 1; 01900 if( j != ih ) 01901 cdaxpy(ns,1.0f,&S(j_,0),1,c,1); 01902 } 01903 csscal(ns,1.0f/ns,c,1); 01904 } 01905 else if( ih != inew ) 01906 { 01907 for( i=0; i < ns; i++ ) 01908 { 01909 c[i] += (S(inew-1,i) - S(ih-1,i))/ns; 01910 } 01911 } 01912 return; 01913 } 01914 01915 /********************************************************************* */ 01916 STATIC void newpt(long int ns, 01917 double coef, 01918 realnum xbase[], 01919 realnum xold[], 01920 int IntNew, 01921 realnum xnew[], 01922 int *small) 01923 { 01924 int eqbase, 01925 eqold; 01926 long int i; 01927 realnum xoldi; 01928 01929 DEBUG_ENTRY( "newpt()" ); 01930 01931 01932 /* Coded by Tom Rowan 01933 * Department of Computer Sciences 01934 * University of Texas at Austin 01935 * 01936 * newpt performs reflections, expansions, contractions, and 01937 * shrinkages (massive contractions) by computing: 01938 * 01939 * xbase + coef * (xbase - xold) 01940 * 01941 * The result is stored in xnew if IntNew .eq. .true., 01942 * in xold otherwise. 01943 * 01944 * use : coef .gt. 0 to reflect 01945 * coef .lt. 0 to expand, contract, or shrink 01946 * 01947 * input 01948 * 01949 * ns - number of components (subspace dimension) 01950 * 01951 * coef - one of four simplex method coefficients 01952 * 01953 * xbase - real ns-vector representing base 01954 * point 01955 * 01956 * xold - real ns-vector representing old 01957 * point 01958 * 01959 * IntNew - logical switch 01960 * = .true. : store result in xnew 01961 * = .false. : store result in xold, xnew is not 01962 * referenced 01963 * 01964 * output 01965 * 01966 * xold - unchanged if IntNew .eq. .true., contains new 01967 * point otherwise 01968 * 01969 * xnew - real ns-vector representing new 01970 * point if new .eq. .true., not referenced 01971 * otherwise 01972 * 01973 * small - logical flag 01974 * = .true. : coincident points 01975 * = .false. : otherwise 01976 * 01977 * local variables 01978 * */ 01979 01980 /* subroutines and functions 01981 * 01982 * fortran */ 01983 01984 /*----------------------------------------------------------- 01985 * */ 01986 eqbase = true; 01987 eqold = true; 01988 if( IntNew ) 01989 { 01990 for( i=0; i < ns; i++ ) 01991 { 01992 xnew[i] = (realnum)(xbase[i] + coef*(xbase[i] - xold[i])); 01993 eqbase = eqbase && ((double)(xnew[i]) == (double)(xbase[i])); 01994 eqold = eqold && ((double)(xnew[i]) == (double)(xold[i])); 01995 } 01996 } 01997 else 01998 { 01999 for( i=0; i < ns; i++ ) 02000 { 02001 xoldi = xold[i]; 02002 xold[i] = (realnum)(xbase[i] + coef*(xbase[i] - xold[i])); 02003 eqbase = eqbase && ((double)(xold[i]) == (double)(xbase[i])); 02004 eqold = eqold && ((double)(xold[i]) == (double)(xoldi)); 02005 } 02006 } 02007 *small = eqbase || eqold; 02008 return; 02009 } 02010 /********************************************************************* */ 02011 STATIC void cdaxpy(long int n, 02012 double da, 02013 realnum dx[], 02014 long int incx, 02015 realnum dy[], 02016 long int incy) 02017 { 02018 long int i, 02019 i_, 02020 ix, 02021 iy, 02022 m; 02023 02024 DEBUG_ENTRY( "cdaxpy()" ); 02025 02026 /* constant times a vector plus a vector. 02027 * uses unrolled loops for increments equal to one. 02028 * jack dongarra, lapack, 3/11/78. 02029 * */ 02030 02031 if( n > 0 ) 02032 { 02033 if( da != 0.00f ) 02034 { 02035 if( incx == 1 && incy == 1 ) 02036 { 02037 02038 /* code for both increments equal to 1 02039 * 02040 * 02041 * clean-up loop 02042 * */ 02043 m = n%4; 02044 if( m != 0 ) 02045 { 02046 for( i=1; i <= m; i++ ) 02047 { 02048 i_ = i - 1; 02049 dy[i_] += (realnum)(da*dx[i_]); 02050 } 02051 if( n < 4 ) 02052 { 02053 return; 02054 } 02055 } 02056 02057 for( i=m; i < n; i += 4 ) 02058 { 02059 dy[i] += (realnum)(da*dx[i]); 02060 dy[i+1] += (realnum)(da*dx[i+1]); 02061 dy[i+2] += (realnum)(da*dx[i+2]); 02062 dy[i+3] += (realnum)(da*dx[i+3]); 02063 } 02064 } 02065 else 02066 { 02067 02068 /* code for unequal increments or equal increments 02069 * not equal to 1 02070 * */ 02071 ix = 1; 02072 iy = 1; 02073 if( incx < 0 ) 02074 ix = (-n + 1)*incx + 1; 02075 if( incy < 0 ) 02076 iy = (-n + 1)*incy + 1; 02077 for( i=0; i < n; i++ ) 02078 { 02079 dy[iy-1] += (realnum)(da*dx[ix-1]); 02080 ix += incx; 02081 iy += incy; 02082 } 02083 } 02084 } 02085 } 02086 return; 02087 } 02088 /*lint +e801 goto is deprecated */ 02089 /*lint +e64 type mismatch */