00001
00002
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
00025
00026
00027
00028
00029
00030
00031
00032
00033
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
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184 if( (mode%2) == 0 )
00185 {
00186
00187
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
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
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
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
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
00394
00395 for( i=0; i < n; i++ )
00396 {
00397 work[i] = x[i] - work[i];
00398 }
00399
00400
00401
00402 L_90:
00403
00404
00405
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
00415
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
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
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
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483 usubc.alpha = 1.0f;
00484
00485
00486
00487
00488 usubc.beta = .50f;
00489
00490
00491
00492
00493 usubc.gamma = 2.0f;
00494
00495
00496
00497
00498 usubc.delta = .50f;
00499
00500
00501
00502
00503
00504
00505
00506
00507 usubc.psi = .250f;
00508
00509
00510
00511
00512 usubc.omega = .10f;
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524 usubc.nsmin = MIN2(2,n);
00525
00526
00527
00528 usubc.nsmax = MIN2(5,n);
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574 usubc.irepl = 0;
00575
00576
00577
00578
00579
00580
00581
00582 usubc.ifxsw = 1;
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599 usubc.bonus = 1.0f;
00600
00601
00602
00603
00604
00605
00606
00607
00608 usubc.nfstop = 0;
00609
00610
00611
00612
00613
00614
00615
00616
00617 usubc.minf = true;
00618 return;
00619 }
00620
00621
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
00637
00638
00639
00640
00641 if( n > 0 )
00642 {
00643 if( incx == 1 && incy == 1 )
00644 {
00645
00646
00647
00648
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
00677
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
00709
00710
00711 DEBUG_ENTRY( "evalf()" );
00712
00713 i_ = n;
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751
00752
00753
00754
00755
00756
00757
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
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
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874
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
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
00924
00925
00926
00927
00928
00929
00930
00931
00932
00933
00934
00935
00936
00937
00938
00939
00940
00941
00942
00943
00944
00945
00946
00947
00948
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
01013
01014
01015
01016
01017
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027
01028
01029
01030
01031
01032
01033
01034
01035
01036
01037
01038
01039
01040
01041
01042
01043
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
01144
01145
01146
01147
01148
01149
01150
01151
01152
01153
01154
01155
01156
01157
01158
01159
01160
01161
01162
01163
01164
01165
01166
01167
01168
01169
01170
01171
01172
01173
01174
01175
01176
01177
01178
01179
01180
01181
01182
01183
01184
01185
01186
01187
01188
01189
01190
01191
01192
01193
01194
01195
01196
01197
01198
01199
01200
01201
01202
01203
01204
01205
01206
01207
01208
01209
01210
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
01249
01250 L_20:
01251 calcc(ns,s,ih,inew,updatc,&S(icent-1,0));
01252 updatc = true;
01253 inew = ih;
01254
01255
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
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
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
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
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
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
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
01391
01392
01393
01394
01395
01396
01397
01398
01399
01400
01401
01402
01403
01404
01405
01406
01407
01408
01409
01410
01411
01412
01413
01414
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
01450
01451
01452
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
01471
01472
01473
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
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
01530
01531
01532
01533
01534
01535
01536
01537 if( !(n <= 0 || incx <= 0) )
01538 {
01539 if( incx == 1 )
01540 {
01541
01542
01543
01544
01545
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
01575
01576 nincx = n*incx;
01577
01578
01579 for( i=0; i<nincx; i=i+incx)
01580
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
01610 i_ = n;
01611
01612
01613
01614
01615
01616
01617
01618
01619
01620
01621
01622
01623
01624
01625
01626
01627
01628
01629
01630
01631
01632
01633
01634
01635
01636
01637
01638
01639
01640
01641
01642
01643
01644
01645
01646
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
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
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
01696
01697
01698
01699
01700
01701
01702
01703
01704
01705
01706
01707
01708
01709
01710
01711
01712
01713
01714
01715
01716
01717
01718
01719
01720
01721
01722
01723
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
01780
01781
01782
01783
01784
01785
01786
01787
01788
01789
01790
01791
01792
01793
01794
01795
01796
01797
01798
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
01848 realnum xNothing[1] = { 0.0f };
01849
01850 DEBUG_ENTRY( "calcc()" );
01851
01852
01853
01854
01855
01856
01857
01858
01859
01860
01861
01862
01863
01864
01865
01866
01867
01868
01869
01870
01871
01872
01873
01874
01875
01876
01877
01878
01879
01880
01881
01882
01883
01884
01885
01886
01887
01888
01889
01890
01891
01892 if( !updatc )
01893 {
01894
01895
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
01933
01934
01935
01936
01937
01938
01939
01940
01941
01942
01943
01944
01945
01946
01947
01948
01949
01950
01951
01952
01953
01954
01955
01956
01957
01958
01959
01960
01961
01962
01963
01964
01965
01966
01967
01968
01969
01970
01971
01972
01973
01974
01975
01976
01977
01978
01979
01980
01981
01982
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
02027
02028
02029
02030
02031 if( n > 0 )
02032 {
02033 if( da != 0.00f )
02034 {
02035 if( incx == 1 && incy == 1 )
02036 {
02037
02038
02039
02040
02041
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
02069
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
02089