00001
00002
00003
00004
00005 #if defined(unix) || defined(__unix) || defined(__unix__)
00006
00007
00008
00009
00010
00011 #define _INCLUDE_POSIX_SOURCE
00012 #endif
00013 #include "cddefines.h"
00014 #include "version.h"
00015 #include "optimize.h"
00016 #ifdef __unix
00017 #include <cstddef>
00018 #include <sys/types.h>
00019 #include <sys/wait.h>
00020 #include <unistd.h>
00021 #else
00022 #define pid_t int
00023 #define fork() TotalInsanityAsStub<pid_t>()
00024 #define wait(X) TotalInsanityAsStub<pid_t>()
00025 #endif
00026 #include "mpi_utilities.h"
00027
00028
00031 const char* STATEFILE = "continue.pmr";
00032 const char* STATEFILE_BACKUP = "continue.bak";
00033
00034
00035
00036
00037 inline void wr_block(const void*,size_t,const char*);
00038 inline void rd_block(void*,size_t,const char*);
00039
00040
00041 template<class X, class Y, int NP, int NSTR>
00042 void phymir_state<X,Y,NP,NSTR>::p_clear1()
00043 {
00044 DEBUG_ENTRY( "p_clear1()" );
00045
00046
00047
00048 memset( this, 0, sizeof(*this) );
00049
00050 p_xmax = numeric_limits<X>::max();
00051 p_ymax = numeric_limits<Y>::max() / Y(2.);
00052 for( int i=0; i < 2*NP+1; ++i )
00053 {
00054 for( int j=0; j < NP; ++j )
00055 p_xp[i][j] = -p_xmax;
00056 p_yp[i] = -p_ymax;
00057 }
00058 for( int i=0; i < NP; ++i )
00059 {
00060 p_absmin[i] = -p_xmax;
00061 p_absmax[i] = p_xmax;
00062 p_varmin[i] = p_xmax;
00063 p_varmax[i] = -p_xmax;
00064 for( int j=0; j < NP; ++j )
00065 p_a2[i][j] = -p_xmax;
00066 p_c1[i] = -p_xmax;
00067 p_c2[i] = -p_xmax;
00068 p_xc[i] = -p_xmax;
00069 p_xcold[i] = -p_xmax;
00070 }
00071 p_vers = VRSNEW;
00072 p_toler = -p_xmax;
00073 p_dmax = X(0.);
00074 p_dold = X(0.);
00075 p_ymin = p_ymax;
00076 p_dim = int32(NP);
00077 p_sdim = int32(NSTR);
00078 p_nvar = int32(0);
00079 p_noptim = int32(0);
00080 p_maxiter = int32(0);
00081 p_jmin = int32(-1);
00082 p_maxcpu = int32(1);
00083 p_curcpu = int32(0);
00084 p_mode = PHYMIR_ILL;
00085
00086
00087
00088
00089
00090
00091 p_size = static_cast<uint32>(reinterpret_cast<size_t>(&p_size) - reinterpret_cast<size_t>(this));
00092 p_func = NULL;
00093 }
00094
00095 template<class X, class Y, int NP, int NSTR>
00096 void phymir_state<X,Y,NP,NSTR>::p_wr_state(const char *fnam) const
00097 {
00098 DEBUG_ENTRY( "p_wr_state()" );
00099
00100 if( cpu.lgMaster() && strlen(fnam) > 0 )
00101 {
00102 FILE *fdes = open_data( fnam, "wb", AS_LOCAL_ONLY_TRY );
00103 bool lgErr = ( fdes == NULL );
00104 lgErr = lgErr || ( fwrite( &p_size, sizeof(uint32), 1, fdes ) != 1 );
00105 lgErr = lgErr || ( fwrite( this, static_cast<size_t>(p_size), 1, fdes ) != 1 );
00106 lgErr = lgErr || ( fclose(fdes) != 0 );
00107 if( lgErr )
00108 {
00109 printf( "p_wr_state: error writing file: %s\n", fnam );
00110 remove(fnam);
00111 }
00112 }
00113 return;
00114 }
00115
00116 template<class X, class Y, int NP, int NSTR>
00117 void phymir_state<X,Y,NP,NSTR>::p_rd_state(const char *fnam)
00118 {
00119 DEBUG_ENTRY( "p_rd_state()" );
00120
00121 FILE *fdes = open_data( fnam, "rb", AS_LOCAL_ONLY );
00122 uint32 wrsize;
00123 bool lgErr = ( fread( &wrsize, sizeof(uint32), 1, fdes ) != 1 );
00124 lgErr = lgErr || ( p_size != wrsize );
00125 lgErr = lgErr || ( fread( this, static_cast<size_t>(p_size), 1, fdes ) != 1 );
00126 lgErr = lgErr || ( fclose(fdes) != 0 );
00127 if( lgErr )
00128 {
00129 printf( "p_rd_state: error reading file: %s\n", fnam );
00130 cdEXIT(EXIT_FAILURE);
00131 }
00132 return;
00133 }
00134
00135 template<class X, class Y, int NP, int NSTR>
00136 Y phymir_state<X,Y,NP,NSTR>::p_execute_job(const X x[],
00137 int jj,
00138 int runNr)
00139 {
00140 DEBUG_ENTRY( "p_execute_job()" );
00141
00142 pid_t pid;
00143 switch( p_mode )
00144 {
00145 case PHYMIR_SEQ:
00146 return p_lgLimitExceeded(x) ? p_ymax : p_func(x,runNr);
00147 case PHYMIR_FORK:
00148 p_curcpu++;
00149 if( p_curcpu > p_maxcpu )
00150 {
00151
00152 (void)wait(NULL);
00153 p_curcpu--;
00154 }
00155
00156 fflush(NULL);
00157 pid = fork();
00158 if( pid < 0 )
00159 {
00160 fprintf( ioQQQ, "creating the child process failed\n" );
00161 cdEXIT(EXIT_FAILURE);
00162 }
00163 else if( pid == 0 )
00164 {
00165
00166 p_execute_job_parallel( x, jj, runNr );
00167
00168
00169 ioQQQ = NULL;
00170 cdEXIT(EXIT_SUCCESS);
00171 }
00172
00173 return p_ymax;
00174 case PHYMIR_MPI:
00175 if( (jj%cpu.nCPU()) == cpu.nRANK() )
00176 p_execute_job_parallel( x, jj, runNr );
00177
00178 return p_ymax;
00179 default:
00180 TotalInsanity();
00181 }
00182 }
00183
00184
00185 template<class X, class Y, int NP, int NSTR>
00186 void phymir_state<X,Y,NP,NSTR>::p_execute_job_parallel(const X x[],
00187 int jj,
00188 int runNr) const
00189 {
00190 DEBUG_ENTRY( "p_execute_job_parallel()" );
00191
00192 char fnam1[20], fnam2[20];
00193 sprintf(fnam1,"yval_%d",jj);
00194 sprintf(fnam2,"output_%d",jj);
00195
00196 FILE *ioQQQ_old = ioQQQ;
00197 ioQQQ = open_data( fnam2, "w", AS_LOCAL_ONLY );
00198
00199 Y yval = p_ymax;
00200 wr_block(&yval,sizeof(yval),fnam1);
00201 if( !p_lgLimitExceeded(x) )
00202 {
00203 try
00204 {
00205 yval = p_func(x,runNr);
00206 }
00207 catch( ... )
00208 {
00209
00210 fflush(NULL);
00211 yval = p_ymax;
00212 }
00213 wr_block(&yval,sizeof(yval),fnam1);
00214 }
00215 fclose( ioQQQ );
00216 ioQQQ = ioQQQ_old;
00217 }
00218
00219 template<class X, class Y, int NP, int NSTR>
00220 void phymir_state<X,Y,NP,NSTR>::p_barrier(int jlo,
00221 int jhi)
00222 {
00223 DEBUG_ENTRY( "p_barrier()" );
00224
00225 switch( p_mode )
00226 {
00227 case PHYMIR_SEQ:
00228
00229 break;
00230 case PHYMIR_FORK:
00231
00232 while( p_curcpu > 0 )
00233 {
00234 (void)wait(NULL);
00235 p_curcpu--;
00236 }
00237 p_process_output( jlo, jhi );
00238 break;
00239 case PHYMIR_MPI:
00240
00241 MPI::COMM_WORLD.Barrier();
00242 p_process_output( jlo, jhi );
00243
00244 MPI::COMM_WORLD.Bcast( &p_yp[jlo], jhi-jlo+1, MPI::type(p_yp), 0 );
00245 break;
00246 default:
00247 TotalInsanity();
00248 }
00249 }
00250
00251 template<class X, class Y, int NP, int NSTR>
00252 void phymir_state<X,Y,NP,NSTR>::p_process_output(int jlo,
00253 int jhi)
00254 {
00255 DEBUG_ENTRY( "p_process_output()" );
00256
00257 if( cpu.lgMaster() )
00258 {
00259 char fnam[20];
00260 for( int jj=jlo; jj <= jhi; jj++ )
00261 {
00262 sprintf(fnam,"yval_%d",jj);
00263 rd_block(&p_yp[jj],sizeof(p_yp[jj]),fnam);
00264 remove(fnam);
00265 sprintf(fnam,"output_%d",jj);
00266 append_file(ioQQQ,fnam);
00267 remove(fnam);
00268 }
00269 }
00270 }
00271
00272 template<class X, class Y, int NP, int NSTR>
00273 void phymir_state<X,Y,NP,NSTR>::p_evaluate_hyperblock()
00274 {
00275 DEBUG_ENTRY( "p_evaluate_hyperblock()" );
00276
00277 int jlo = 1, jhi = 0;
00278 for( int j=0; j < p_nvar; j++ )
00279 {
00280 X sgn = X(1.);
00281 for( int jj=2*j+1; jj <= 2*j+2; jj++ )
00282 {
00283 sgn = -sgn;
00284 for( int i=0; i < p_nvar; i++ )
00285 {
00286 p_xp[jj][i] = p_xc[i] + sgn*p_dmax*p_c2[j]*p_a2[j][i];
00287 p_varmax[i] = max(p_varmax[i],p_xp[jj][i]);
00288 p_varmin[i] = min(p_varmin[i],p_xp[jj][i]);
00289 }
00290 if( !lgMaxIterExceeded() )
00291 {
00292
00293
00294 p_yp[jj] = p_execute_job( p_xp[jj], jj, p_noptim++ );
00295 jhi = jj;
00296 }
00297 }
00298 }
00299
00300 p_barrier( jlo, jhi );
00301 }
00302
00303 template<class X, class Y, int NP, int NSTR>
00304 void phymir_state<X,Y,NP,NSTR>::p_setup_next_hyperblock()
00305 {
00306 DEBUG_ENTRY( "p_setup_next_hyperblock()" );
00307
00308 const Y F0 = Y(1.4142136);
00309 const X F1 = X(0.7071068);
00310 const X F2 = X(0.1);
00311
00312
00313 for( int jj=1; jj <= 2*p_nvar; jj++ )
00314 {
00315 if( p_yp[jj] < p_ymin )
00316 {
00317 p_ymin = p_yp[jj];
00318 p_jmin = jj;
00319 }
00320 }
00321 bool lgNewCnt = p_jmin > 0;
00322
00323
00324 bool lgNegd2 = false;
00325 X xnrm = X(0.);
00326 X xhlp[NP];
00327 for( int i=0; i < p_nvar; i++ )
00328 {
00329 Y d1 = p_yp[2*i+2] - p_yp[2*i+1];
00330 Y d2 = Y(0.5)*p_yp[2*i+2] - p_yp[0] + Y(0.5)*p_yp[2*i+1];
00331 if( d2 <= Y(0.) )
00332 lgNegd2 = true;
00333 xhlp[i] = -p_dmax*p_c2[i]*(static_cast<X>(max(min((Y(0.25)*d1)/max(d2,Y(1.e-10)),F0),-F0)) -
00334 p_delta(2*i+1,p_jmin) + p_delta(2*i+2,p_jmin));
00335 xnrm += pow2(xhlp[i]);
00336 }
00337 xnrm = static_cast<X>(sqrt(xnrm));
00338
00339 int imax = 0;
00340 X amax = X(0.);
00341 for( int j=0; j < p_nvar; j++ )
00342 {
00343 for( int i=0; i < p_nvar; i++ )
00344 {
00345 if( xnrm > X(0.) )
00346 {
00347 if( j == 0 )
00348 {
00349 p_a2[0][i] *= xhlp[0];
00350 }
00351 else
00352 {
00353 p_a2[0][i] += xhlp[j]*p_a2[j][i];
00354 p_a2[j][i] = p_delta(i,j);
00355 if( j == p_nvar-1 && abs(p_a2[0][i]) > amax )
00356 {
00357 imax = i;
00358 amax = abs(p_a2[0][i]);
00359 }
00360 }
00361 }
00362 else
00363 {
00364 p_a2[j][i] = p_delta(i,j);
00365 }
00366 }
00367 }
00368
00369 if( imax > 0 )
00370 {
00371 p_a2[imax][0] = X(1.);
00372 p_a2[imax][imax] = X(0.);
00373 }
00374
00375 p_phygrm( p_a2, p_nvar );
00376
00377
00378 for( int i=0; i < p_nvar; i++ )
00379 {
00380 p_c2[i] = X(0.);
00381 for( int j=0; j < p_nvar; j++ )
00382 {
00383 p_c2[i] += pow2(p_a2[i][j]/p_c1[j]);
00384 }
00385 p_c2[i] = static_cast<X>(1./sqrt(p_c2[i]));
00386 p_xc[i] = p_xp[p_jmin][i];
00387 p_xp[0][i] = p_xc[i];
00388 }
00389 p_yp[0] = p_yp[p_jmin];
00390 p_jmin = 0;
00391
00392 X r1, r2;
00393 if( lgNegd2 )
00394 {
00395
00396
00397
00398 r1 = F1;
00399 r2 = F1;
00400 }
00401 else
00402 {
00403 r1 = F2;
00404 if( lgNewCnt )
00405 {
00406
00407 r2 = sqrt(1.f/F1);
00408 }
00409 else
00410 {
00411
00412
00413 r2 = sqrt(F1);
00414 }
00415 }
00416 p_dmax = min(max(xnrm/p_c2[0],p_dmax*r1),p_dmax*r2);
00417
00418 p_dmax = MIN2(p_dmax,p_dold);
00419 }
00420
00421 template<class X, class Y, int NP, int NSTR>
00422 void phymir_state<X,Y,NP,NSTR>::p_reset_hyperblock()
00423 {
00424 DEBUG_ENTRY( "p_reset_hyperblock()" );
00425
00426 if( !lgConvergedRestart() )
00427 {
00428
00429 for( int i=0; i < p_nvar; i++ )
00430 {
00431 p_xcold[i] = p_xc[i];
00432 p_c2[i] = p_c1[i];
00433 }
00434 p_dmax = p_dold;
00435 p_reset_transformation_matrix();
00436 }
00437 }
00438
00439 template<class X, class Y, int NP, int NSTR>
00440 void phymir_state<X,Y,NP,NSTR>::p_phygrm(X a[][NP],
00441 int n)
00442 {
00443 DEBUG_ENTRY( "p_phygrm()" );
00444
00445
00446 for( int i=0; i < n; i++ )
00447 {
00448 X ip = X(0.);
00449 for( int k=0; k < n; k++ )
00450 ip += pow2(a[i][k]);
00451 ip = sqrt(ip);
00452 for( int k=0; k < n; k++ )
00453 a[i][k] /= ip;
00454 for( int j=i+1; j < n; j++ )
00455 {
00456 X ip = X(0.);
00457 for( int k=0; k < n; k++ )
00458 ip += a[i][k]*a[j][k];
00459 for( int k=0; k < n; k++ )
00460 a[j][k] -= ip*a[i][k];
00461 }
00462 }
00463 return;
00464 }
00465
00466 template<class X, class Y, int NP, int NSTR>
00467 bool phymir_state<X,Y,NP,NSTR>::p_lgLimitExceeded(const X x[]) const
00468 {
00469 DEBUG_ENTRY( "p_lgLimitExceeded()" );
00470
00471 for( int i=0; i < p_nvar; i++ )
00472 {
00473 if( x[i] < p_absmin[i] )
00474 return true;
00475 if( x[i] > p_absmax[i] )
00476 return true;
00477 }
00478 return false;
00479 }
00480
00481 template<class X, class Y, int NP, int NSTR>
00482 void phymir_state<X,Y,NP,NSTR>::p_reset_transformation_matrix()
00483 {
00484 DEBUG_ENTRY( "p_reset_transformation_matrix()" );
00485
00486
00487 for( int i=0; i < p_nvar; i++ )
00488 for( int j=0; j < p_nvar; j++ )
00489 p_a2[j][i] = p_delta(i,j);
00490 }
00491
00492 template<class X, class Y, int NP, int NSTR>
00493 void phymir_state<X,Y,NP,NSTR>::init_minmax(const X pmin[],
00494 const X pmax[],
00495 int nvar)
00496 {
00497 DEBUG_ENTRY( "init_minmax()" );
00498
00499 ASSERT( !lgInitialized() );
00500
00501 for( int i=0; i < nvar; i++ )
00502 {
00503 p_absmin[i] = pmin[i];
00504 p_absmax[i] = pmax[i];
00505 }
00506 }
00507
00508 template<class X, class Y, int NP, int NSTR>
00509 void phymir_state<X,Y,NP,NSTR>::init_strings(const char* date,
00510 const char* version,
00511 const char* host_name)
00512 {
00513 DEBUG_ENTRY( "init_strings()" );
00514
00515
00516 if( date != NULL )
00517 strncpy( p_chStr1, date, NSTR-1 );
00518 if( version != NULL )
00519 strncpy( p_chStr2, version, NSTR-1 );
00520 if( host_name != NULL )
00521 strncpy( p_chStr3, host_name, NSTR-1 );
00522 }
00523
00524 template<class X, class Y, int NP, int NSTR>
00525 void phymir_state<X,Y,NP,NSTR>::init_state_file_name(const char* fnam)
00526 {
00527 DEBUG_ENTRY( "init_state_file_name()" );
00528
00529
00530 strncpy( p_chState, fnam, NSTR-1 );
00531 }
00532
00533 template<class X, class Y, int NP, int NSTR>
00534 void phymir_state<X,Y,NP,NSTR>::initial_run(Y (*func)(const X[],int),
00535 int nvar,
00536 const X start[],
00537 const X del[],
00538 X toler,
00539 int maxiter,
00540 phymir_mode mode,
00541 int maxcpu)
00542 {
00543 DEBUG_ENTRY( "initial_run()" );
00544
00545 ASSERT( nvar > 0 && nvar <= NP );
00546
00547 p_func = func;
00548 p_nvar = nvar;
00549 p_toler = toler;
00550 p_maxiter = maxiter;
00551 p_mode = mode;
00552 p_maxcpu = maxcpu;
00553 p_noptim = 0;
00554
00555
00556 p_dmax = X(0.);
00557 for( int i=0; i < p_nvar; i++ )
00558 p_dmax = max(p_dmax,abs(del[i]));
00559
00560 p_dold = p_dmax;
00561 for( int i=0; i < p_nvar; i++ )
00562 {
00563 p_xc[i] = start[i];
00564 p_xcold[i] = p_xc[i] + X(10.)*p_toler;
00565 p_c1[i] = abs(del[i])/p_dmax;
00566 p_c2[i] = p_c1[i];
00567 p_xp[0][i] = p_xc[i];
00568 p_varmax[i] = max(p_varmax[i],p_xc[i]);
00569 p_varmin[i] = min(p_varmin[i],p_xc[i]);
00570 }
00571
00572
00573
00574 p_yp[0] = p_execute_job( p_xc, 0, p_noptim++ );
00575 p_barrier( 0, 0 );
00576
00577 p_ymin = p_yp[0];
00578 p_jmin = 0;
00579
00580 p_reset_transformation_matrix();
00581
00582 p_wr_state( p_chState );
00583 }
00584
00585 template<class X, class Y, int NP, int NSTR>
00586 void phymir_state<X,Y,NP,NSTR>::continue_from_state(Y (*func)(const X[],int),
00587 int nvar,
00588 const char* fnam,
00589 X toler,
00590 int maxiter,
00591 phymir_mode mode,
00592 int maxcpu)
00593 {
00594 DEBUG_ENTRY( "continue_from_state()" );
00595
00596 p_rd_state( fnam );
00597
00598 if( !fp_equal( p_vers, VRSNEW ) )
00599 {
00600 printf( "optimize continue - file has incompatible version, sorry\n" );
00601 cdEXIT(EXIT_FAILURE);
00602 }
00603 if( p_dim != NP )
00604 {
00605 printf( "optimize continue - arrays have wrong dimension, sorry\n" );
00606 cdEXIT(EXIT_FAILURE);
00607 }
00608 if( p_sdim != NSTR )
00609 {
00610 printf( "optimize continue - strings have wrong length, sorry\n" );
00611 cdEXIT(EXIT_FAILURE);
00612 }
00613 if( p_nvar != nvar )
00614 {
00615 printf( "optimize continue - wrong number of free parameters, sorry\n" );
00616 cdEXIT(EXIT_FAILURE);
00617 }
00618
00619 p_func = func;
00620
00621
00622 p_toler = toler;
00623
00624 p_maxiter = maxiter;
00625
00626 p_mode = mode;
00627 p_maxcpu = maxcpu;
00628 }
00629
00630 template<class X, class Y, int NP, int NSTR>
00631 void phymir_state<X,Y,NP,NSTR>::optimize()
00632 {
00633 DEBUG_ENTRY( "optimize()" );
00634
00635 ASSERT( lgInitialized() );
00636
00637 while( !lgConverged() )
00638 {
00639 p_evaluate_hyperblock();
00640 if( lgMaxIterExceeded() )
00641 break;
00642 p_setup_next_hyperblock();
00643 p_wr_state( p_chState );
00644 }
00645 }
00646
00647 template<class X, class Y, int NP, int NSTR>
00648 void phymir_state<X,Y,NP,NSTR>::optimize_with_restart()
00649 {
00650 DEBUG_ENTRY( "optimize_with_restart()" );
00651
00652 ASSERT( lgInitialized() );
00653
00654 while( !lgConvergedRestart() )
00655 {
00656 optimize();
00657 if( lgMaxIterExceeded() )
00658 break;
00659 p_reset_hyperblock();
00660 }
00661 }
00662
00663 template<class X, class Y, int NP, int NSTR>
00664 bool phymir_state<X,Y,NP,NSTR>::lgConvergedRestart() const
00665 {
00666 DEBUG_ENTRY( "lgConvergedRestart()" );
00667
00668 if( lgConverged() )
00669 {
00670 X dist = X(0.);
00671 for( int i=0; i < p_nvar; i++ )
00672 dist += pow2(p_xc[i] - p_xcold[i]);
00673 dist = static_cast<X>(sqrt(dist));
00674 return ( dist <= p_toler );
00675 }
00676 return false;
00677 }
00678
00679 void optimize_phymir(realnum xc[],
00680 const realnum del[],
00681 long int nvarPhymir,
00682 chi2_type *ymin,
00683 realnum toler)
00684 {
00685 DEBUG_ENTRY( "optimize_phymir()" );
00686
00687 if( nvarPhymir > LIMPAR )
00688 {
00689 fprintf( ioQQQ, "optimize_phymir: too many parameters are varied, increase LIMPAR\n" );
00690 cdEXIT(EXIT_FAILURE);
00691 }
00692
00693 phymir_state<realnum,chi2_type,LIMPAR,STDLEN> phymir;
00694
00695
00696 (void)remove(STATEFILE_BACKUP);
00697 FILE *tmp = open_data( STATEFILE, "r", AS_LOCAL_ONLY_TRY );
00698 if( tmp != NULL )
00699 {
00700 fclose( tmp );
00701
00702 FILE *dest = open_data( STATEFILE_BACKUP, "wb", AS_LOCAL_ONLY_TRY );
00703 if( dest != NULL )
00704 {
00705 append_file( dest, STATEFILE );
00706 fclose( dest );
00707 }
00708 }
00709
00710 phymir_mode mode = optimize.lgParallel ? ( cpu.lgMPI() ? PHYMIR_MPI : PHYMIR_FORK ) : PHYMIR_SEQ;
00711 long nCPU = optimize.lgParallel ? ( cpu.lgMPI() ? cpu.nCPU() : optimize.useCPU ) : 1;
00712 if( optimize.lgOptCont )
00713 {
00714 phymir.continue_from_state( optimize_func, nvarPhymir, STATEFILE, toler,
00715 optimize.nIterOptim, mode, nCPU );
00716 }
00717 else
00718 {
00719 phymir.init_state_file_name( STATEFILE );
00720 phymir.init_strings( t_version::Inst().chDate, t_version::Inst().chVersion,
00721 cpu.host_name() );
00722 phymir.initial_run( optimize_func, nvarPhymir, xc, del, toler,
00723 optimize.nIterOptim, mode, nCPU );
00724 }
00725
00726 phymir.optimize_with_restart();
00727
00728 if( phymir.lgMaxIterExceeded() )
00729 {
00730 fprintf( ioQQQ, " Optimizer exceeding maximum iterations.\n" );
00731 fprintf( ioQQQ, " This can be reset with the OPTIMIZE ITERATIONS command.\n" );
00732 }
00733
00734
00735 optimize.nOptimiz = phymir.noptim();
00736 for( int i=0; i < nvarPhymir; i++ )
00737 {
00738 xc[i] = phymir.xval(i);
00739 optimize.varmax[i] = min(phymir.xmax(i),optimize.varang[i][1]);
00740 optimize.varmin[i] = max(phymir.xmin(i),optimize.varang[i][0]);
00741 }
00742 *ymin = phymir.yval();
00743 return;
00744 }
00745
00747 inline void wr_block(const void *ptr,
00748 size_t len,
00749 const char *fnam)
00750 {
00751 DEBUG_ENTRY( "wr_block()" );
00752
00753 FILE *fdes = open_data( fnam, "wb", AS_LOCAL_ONLY );
00754 if( fwrite(ptr,len,size_t(1),fdes) != 1 ) {
00755 printf( "error writing on file: %s\n",fnam );
00756 fclose(fdes);
00757 cdEXIT(EXIT_FAILURE);
00758 }
00759 fclose(fdes);
00760 return;
00761 }
00762
00764 inline void rd_block(void *ptr,
00765 size_t len,
00766 const char *fnam)
00767 {
00768 DEBUG_ENTRY( "rd_block()" );
00769
00770 FILE *fdes = open_data( fnam, "rb", AS_LOCAL_ONLY );
00771 if( fread(ptr,len,size_t(1),fdes) != 1 ) {
00772 printf( "error reading on file: %s\n",fnam );
00773 fclose(fdes);
00774 cdEXIT(EXIT_FAILURE);
00775 }
00776 fclose(fdes);
00777 return;
00778 }