00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include "BonminConfig.h"
00011
00012 #include "BonFilterSolver.hpp"
00013 #include "BonFilterWarmStart.hpp"
00014
00015 #include <fstream>
00016
00017 #include "CoinTime.hpp"
00018 #include<algorithm>
00019 typedef Bonmin::FilterSolver::fint fint;
00020 typedef Bonmin::FilterSolver::real real;
00021
00022
00023
00024
00025 typedef long ftnlen;
00026 extern "C"
00027 {
00028 void F77_FUNC(filtersqp,FILTERSQP)(
00029 fint *n, fint *m, fint *kmax, fint *maxa,
00030 fint *maxf, fint *mlp, fint *mxwk, fint *mxiwk,
00031 fint *iprint, fint *nout, fint *ifail, real *rho,
00032 real *x, real *c, real *f, real *fmin, real *bl,
00033 real *bu, real *s, real *a, fint *la, real *ws,
00034 fint *lws, real *lam, char *cstype, real *user,
00035 fint *iuser, fint *maxiter, fint *istat,
00036 real *rstat, ftnlen cstype_len);
00037 }
00038
00039
00040 static Ipopt::TNLP * tnlpSolved = NULL;
00041 static fint nnz_h = -1;
00042
00043 static fint * hStruct = NULL;
00044
00045
00046 static int * permutationJac = NULL;
00047 static int * permutationHess = NULL;
00048
00049
00050
00051 extern "C"
00052 {
00053
00054
00055 extern struct {
00056 fint char_l;
00057 char pname[10];
00058 }
00059 F77_FUNC(cpname,CPNAME);
00060
00061
00062 extern struct {
00063 fint phl, phr, phc;
00064 }
00065 F77_FUNC(hessc,HESSC);
00066
00067
00068 extern struct {
00069 real ubd, tt;
00070 }
00071 F77_FUNC(ubdc,UBDC);
00072
00073
00074 extern struct {
00075 real infty, eps;
00076 }
00077 F77_FUNC_(nlp_eps_inf,NLP_EPS_INF);
00078
00079
00080 extern struct {
00081 fint n_bqpd_calls, n_bqpd_prfint;
00082 }
00083 F77_FUNC_(bqpd_count,BQPD_COUNT);
00084
00085
00086 extern struct {
00087 fint scale_mode, phe;
00088 }
00089 F77_FUNC(scalec,SCALEC);
00090 }
00091
00092 extern "C"
00093 {
00094
00096 void F77_FUNC(objfun,OBJFUN)(real *x, fint *n, real * f, real *user, fint * iuser, fint * errflag) {
00097 (*errflag) = !tnlpSolved->eval_f(*n, x, 1, *f);
00098 }
00099
00101 void
00102 F77_FUNC(confun,CONFUN)(real * x, fint * n , fint *m, real *c, real *a, fint * la, real * user, fint * iuser,
00103 fint * errflag) {
00104 (*errflag) = !tnlpSolved->eval_g(*n, x, 1, *m, c);
00105 }
00106
00107 void
00108 F77_FUNC(gradient,GRADIENT)(fint *n, fint *m, fint * mxa, real * x, real *a, fint * la,
00109 fint * maxa, real * user, fint * iuser, fint * errflag) {
00110 (*errflag) = !tnlpSolved->eval_grad_f(*n, x, 1, a);
00112 int nnz = la[0] - *n - 1;
00113 double * values = new double [nnz];
00114 (*errflag) = !tnlpSolved->eval_jac_g(*n, x, 1, *m, nnz, NULL, NULL, values) || (*errflag);
00115 a+= *n;
00116 for (int i = 0 ; i < nnz ; i++) {
00117 int indice = permutationJac[i];
00118 if (indice > nnz) {
00119 #ifndef NDEBUG
00120 std::cout<<"Error in gradient computation, i: "<<i
00121 <<" in row order "<<permutationJac[i]<<std::endl;
00122 #endif
00123 }
00124 *a++ = values[indice];
00125 }
00126 delete [] values;
00127 }
00128
00129
00130 void
00131 F77_FUNC(hessian,HESSIAN)(real *x, fint *n, fint *m, fint *phase, real *lam,
00132 real *ws, fint *lws, real *user, fint *iuser,
00133 fint *l_hess, fint *li_hess, fint *errflag) {
00134 Number obj_factor = (*phase == 1)? 0. : 1.;
00135 fint end = nnz_h + (*n) + 2;
00136
00137 for (int i = 0 ; i < end ; i++) {
00138 lws[i] = hStruct[i];
00139 }
00140 *l_hess = nnz_h;
00141 *li_hess = nnz_h + *n + 3;
00142 Number * mlam = NULL;
00143 if (*m > 0) {
00144 mlam = new Number[*m];
00145 }
00146 for (int i = 0; i<*m; i++) {
00147 mlam[i] = -lam[*n+i];
00148 }
00149 Number * values = new Number [nnz_h];
00150 (*errflag) = !tnlpSolved->eval_h(*n, x, 1, obj_factor, *m, mlam ,1, hStruct[0] - 1, NULL, NULL, values);
00151 delete [] mlam;
00152 for (int i = 0 ; i < nnz_h ; i++) ws[i] = values[permutationHess[i]];
00153 delete [] values;
00154 }
00155
00156 }
00157
00158 namespace Bonmin
00159 {
00160
00161 struct Transposer
00162 {
00163 const Index* rowIndices;
00164 const Index* colIndices;
00165 bool operator()(int i, int j)
00166 {
00167 return rowIndices[i]<rowIndices[j] ||
00168 (rowIndices[i]==rowIndices[j] && colIndices[i] < colIndices[j]);
00169 }
00170 };
00171
00172
00173 void FilterSolver::TMat2RowPMat(bool symmetric, fint n, fint m, int nnz,
00174 const Index* iRow,
00175 const Index* iCol, int * permutation2,
00176 fint * lws, int nnz_offset, int n_offset,
00177 Ipopt::TNLP::IndexStyleEnum index_style)
00178 {
00179 for (int i = 0 ; i < nnz ; i++)
00180 permutation2[i] = i;
00181
00182 Transposer lt;
00183 if (symmetric) {
00184 Index* tmpRow = new Index[nnz];
00185 Index* tmpCol = new Index[nnz];
00186 for (int i=0; i<nnz; i++) {
00187 const Index& irow = iRow[i];
00188 const Index& jcol = iCol[i];
00189 if (irow > jcol) {
00190 tmpRow[i] = irow;
00191 tmpCol[i] = jcol;
00192 }
00193 else {
00194 tmpRow[i] = jcol;
00195 tmpCol[i] = irow;
00196 }
00197 }
00198 lt.rowIndices = tmpRow;
00199 lt.colIndices = tmpCol;
00200 }
00201 else {
00202 lt.rowIndices = iRow;
00203 lt.colIndices = iCol;
00204 }
00205
00206 std::sort(permutation2, &permutation2[nnz], lt);
00207
00208 const int idx_offset = (index_style == Ipopt::TNLP::C_STYLE);
00209 fint row = 1-idx_offset;
00210 lws[0] = nnz + nnz_offset + 1;
00211 fint * inds = lws + nnz_offset + 1;
00212 fint * start = inds + nnz + n_offset;
00213 *start++ = 1 + nnz_offset;
00214 for (fint i = 0 ; i < nnz ; i++) {
00215 inds[i] = lt.colIndices[permutation2[i]] + idx_offset;
00216
00217 if (lt.rowIndices[permutation2[i]] > row) {
00218 for (;row < lt.rowIndices[permutation2[i]] ; row++)
00219 *start++ = i + nnz_offset + 1;
00220 }
00221 }
00222 for (;row <= m-idx_offset ; row++)
00223 *start++ = nnz + nnz_offset +1;
00224
00225 #if 0
00226 for (int i = 0; i<nnz_offset+1; i++)
00227 printf("lws[%3d] = %3d\n", i, lws[i]);
00228 for (int i = nnz_offset+1; i<nnz_offset+nnz+1; i++)
00229 printf("lws[%3d] = %3d [%3d,%3d]\n", i, lws[i], lt.rowIndices[permutation2[i-nnz_offset-1]], lt.colIndices[permutation2[i-nnz_offset-1]]);
00230 for (int i = nnz_offset+nnz+1; i<lws[0]+m+2; i++)
00231 printf("lws[%3d] = %3d\n", i, lws[i]);
00232 #endif
00233
00234 if (symmetric) {
00235 delete [] lt.rowIndices;
00236 delete [] lt.colIndices;
00237 }
00238
00239
00240 }
00241
00242 #if 0
00243
00244 void TMat2RowPMat_(fint n, fint m, int nnz, const Index* iRow,
00245 const Index* iCol, int * permutation2,
00246 fint * lws, int offset)
00247 {
00248 for (int i = 0 ; i < nnz ; i++)
00249 permutation2[i] = i;
00250
00251 Transposer lt;
00252 lt.rowIndices = iRow;
00253 lt.colIndices = iCol;
00254
00255 std::sort(permutation2, &permutation2[nnz], lt);
00256
00257 fint row = 1;
00258 lws[0] = nnz + offset + 1;
00259 fint * inds = lws + offset + 1;
00260 fint * start = inds + nnz + 1;
00261
00262 for (fint i = 0 ; i < nnz ; i++) {
00263 inds[i] = iCol[permutation2[i]];
00264
00265 if (iRow[permutation2[i]] >= row) {
00266 for (;row <= iRow[permutation2[i]] ; row++)
00267 *start++ = i + offset + 1;
00268 }
00269 }
00270 for (;row <= m+1 ; row++)
00271 *start++ = nnz + offset +1;
00272
00273
00274 for (int i = 0; i<offset+1; i++)
00275 printf("alws[%3d] = %3d\n", i, lws[i]);
00276 for (int i = offset+1; i<offset+nnz+1; i++)
00277 printf("alws[%3d] = %3d [%3d,%3d]\n", i, lws[i], lt.rowIndices[permutation2[i-offset-1]], lt.colIndices[permutation2[i-offset-1]]);
00278 for (int i = offset+nnz+1; i<lws[0]+m+2; i++)
00279 printf("alws[%3d] = %3d\n", i, lws[i]);
00280 }
00281
00282
00283
00284 void TMat2ColPMat_(fint n, fint m, int nnz, const Index* iRow,
00285 const Index* iCol,
00286 int* permutationHess2, fint * lws, int offset)
00287 {
00288 for (int i = 0 ; i < nnz ; i++)
00289 permutationHess2[i] = i;
00290
00291 fint col = 1;
00292 lws[0] = nnz + offset + 1;
00293 fint * inds = lws + 1;
00294 fint * start = inds + nnz + offset;
00295
00296 Transposer lt;
00297 lt.rowIndices = iCol;
00298 lt.colIndices = iRow;
00299
00300 std::sort(permutationHess2, permutationHess2 + nnz, lt);
00301
00302 for (fint i = 0 ; i < nnz ; i++) {
00303 inds[offset + i] = iRow[permutationHess2[i]];
00304 if (iCol[permutationHess2[i]] >= col) {
00305 for (;col <= iCol[permutationHess2[i]] ; col++)
00306 *start++ = i + offset + 1;
00307 }
00308 }
00309 for (;col <= n+1 ; col++)
00310 *start++ = nnz + offset +1;
00311
00312
00313 for (int i = 0; i<offset+1; i++)
00314 printf("blws[%3d] = %3d\n", i, lws[i]);
00315 for (int i = offset+1; i<offset+nnz+1; i++)
00316 printf("blws[%3d] = %3d [%3d,%3d]\n", i, lws[i], lt.rowIndices[permutationHess2[i-offset-1]], lt.colIndices[permutationHess2[i-offset-1]]);
00317 for (int i = offset+nnz+1; i<lws[0]+m+2; i++)
00318 printf("blws[%3d] = %3d\n", i, lws[i]);
00319 }
00320 #endif
00321
00322
00323 std::string FilterSolver::solverName_ = "filter SQP";
00324
00325 void
00326 FilterSolver::registerOptions(Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions)
00327 {
00328 roptions->SetRegisteringCategory("FilterSQP options", RegisteredOptions::FilterCategory);
00329 roptions->AddLowerBoundedNumberOption("eps", "Tolerance for SQP solver",
00330 0., 1, 1e-08, "");
00331 roptions->AddLowerBoundedNumberOption("infty","A large number (1E20)",0.,1, 1e20, "");
00332 roptions->AddBoundedIntegerOption("iprint", "Print level (0=silent, 3=verbose)", 0,6,0);
00333 roptions->AddLowerBoundedIntegerOption("kmax", "Dimension of null-space",
00334 -1, -1, "");
00335 roptions->AddLowerBoundedIntegerOption("maxf","Maximum filter length",0,100);
00336 roptions->AddLowerBoundedIntegerOption("maxiter", "Maximum number of iterations",0,1000);
00337 roptions->AddLowerBoundedIntegerOption("mlp","Maximum level for degeneracy (bqpd)",0, 1000);
00338 roptions->AddLowerBoundedIntegerOption("mxlws", "FINTEGER workspace increment", 0, 500000);
00339 roptions->AddLowerBoundedIntegerOption("mxws", "REAL workspace increment",
00340 0,2000000);
00341 roptions->AddLowerBoundedNumberOption("rho_init", "Initial trust region size",0,1,10.);
00342
00343 roptions->AddLowerBoundedNumberOption("tt", "Parameter for upper bound on filter",0,1, 125e-2);
00344 roptions->AddLowerBoundedNumberOption("ubd", "Parameter for upper bound on filter", 0 , 1,1e2);
00345
00346 }
00347
00348
00349 FilterSolver::FilterSolver(bool createEmpty )
00350 :
00351 TNLPSolver(),
00352 warmF_(NULL),
00353 cached_(NULL)
00354 {}
00355
00356 FilterSolver::FilterSolver(Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions,
00357 Ipopt::SmartPtr<Ipopt::OptionsList> options,
00358 Ipopt::SmartPtr<Ipopt::Journalist> journalist,
00359 const std::string & prefix):
00360 TNLPSolver(roptions, options, journalist, prefix),
00361 warmF_(NULL),
00362 cached_(NULL)
00363 {}
00364
00365 FilterSolver::FilterSolver(Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions,
00366 Ipopt::SmartPtr<Ipopt::OptionsList> options,
00367 Ipopt::SmartPtr<Ipopt::Journalist> journalist):
00368 TNLPSolver(roptions, options, journalist, "bonmin."),
00369 warmF_(NULL),
00370 cached_(NULL)
00371 {}
00372
00373
00374 FilterSolver::FilterSolver(const FilterSolver & other):
00375 TNLPSolver(other),
00376 warmF_(NULL),
00377 cached_(NULL)
00378 {
00379 warmF_ = (other.warmF_.IsValid()) ? dynamic_cast<FilterWarmStart *>(other.warmF_->clone()):NULL;
00380 }
00381
00382 Ipopt::SmartPtr <TNLPSolver>
00383 FilterSolver::clone()
00384 {
00385 Ipopt::SmartPtr<FilterSolver> retval = new FilterSolver(*this);
00386 return GetRawPtr(retval);
00387 }
00388
00389 FilterSolver::~FilterSolver()
00390 {}
00391
00392 bool
00393 FilterSolver::Initialize(std::string optFile)
00394 {
00395 std::ifstream is;
00396 if (optFile != "") {
00397 try {
00398 is.open(optFile.c_str());
00399 }
00400 catch (std::bad_alloc) {
00401 journalist_->Printf(Ipopt::J_SUMMARY, Ipopt::J_MAIN, "\nEXIT: Not enough memory.\n");
00402 return false;
00403 }
00404 catch (...) {
00405 Ipopt::IpoptException E("Unknown Exception caught in ipopt", "Unknown File", -1);
00406 E.ReportException(*journalist_);
00407 return false;
00408 }
00409 }
00410 bool retval = Initialize(is);
00411 if (is) {
00412 is.close();
00413 }
00414 if(!options_->GetIntegerValue("print_level",default_log_level_,""))
00415 default_log_level_ = 1;
00416 return retval;
00417 }
00418
00419 bool
00420 FilterSolver::Initialize(std::istream &is)
00421 {
00422
00423 Index ivalue;
00424 options_->GetIntegerValue("print_level", ivalue, "");
00425 EJournalLevel print_level = (EJournalLevel)ivalue;
00426 SmartPtr<Journal> stdout_jrnl = journalist_->GetJournal("console");
00427 if (IsValid(stdout_jrnl)) {
00428
00429 stdout_jrnl->SetAllPrintLevels(print_level);
00430 stdout_jrnl->SetPrintLevel(J_DBG, J_NONE);
00431 }
00432
00433 if (is.good()) {
00434 options_->ReadFromStream(*journalist_, is);
00435 }
00436 return true;
00437 }
00438
00440 TNLPSolver::ReturnStatus
00441 FilterSolver::OptimizeTNLP(const Ipopt::SmartPtr<Ipopt::TNLP> & tnlp)
00442 {
00443 if (cached_.IsNull() || !cached_->use_warm_start_in_cache_) {
00444 cached_ = new cachedInfo(tnlp, options_);
00445 }
00446 cached_->load_ws(warmF_);
00447 return callOptimizer();
00448 }
00449
00451 TNLPSolver::ReturnStatus
00452 FilterSolver::ReOptimizeTNLP(const Ipopt::SmartPtr<Ipopt::TNLP> & tnlp)
00453 {
00454 assert(tnlp == cached_->tnlp_);
00455 cached_->load_ws(warmF_);
00456
00457 assert(cached_->bounds);
00458 int n = cached_->n;
00459 int m = cached_->m;
00460 tnlp->get_bounds_info(n, cached_->bounds, &cached_->bounds[n+m],
00461 m, &cached_->bounds[n], &cached_->bounds[2*n + m]);
00462
00463 tnlpSolved = static_cast<Ipopt::TNLP *>(Ipopt::GetRawPtr(tnlp));
00464 nnz_h = cached_->nnz_h_;
00465
00466 hStruct = cached_->hStruct_;
00467
00468
00469 permutationJac = cached_->permutationJac_;
00470 permutationHess = cached_->permutationHess_;
00471
00472
00473 return callOptimizer();
00474 }
00475
00476
00477
00478 void
00479 FilterSolver::cachedInfo::initialize(const Ipopt::SmartPtr<Ipopt::TNLP> & tnlp,
00480 Ipopt::SmartPtr<Ipopt::OptionsList>& options)
00481 {
00482
00483
00484 int nnz_jac_g;
00485
00486 Ipopt::TNLP::IndexStyleEnum index_style;
00487 Ipopt::Index nv, nc, nnz_j, nnz_hess;
00488 tnlp->get_nlp_info( nv, nc,
00489 nnz_j, (Ipopt::Index&) nnz_hess,
00490 index_style);
00491 n = nv;
00492 m = nc;
00493 nnz_jac_g = nnz_j;
00494 nnz_h_ = nnz_hess;
00495
00496 nnz_h = nnz_h_;
00497
00498
00499
00500 Ipopt::Index kmax_ipt;
00501 options->GetIntegerValue("kmax", kmax_ipt, "filter.");
00502 if (kmax_ipt == -1) {
00503 kmax = n;
00504 }
00505 else {
00506 kmax = kmax_ipt;
00507 kmax = std::min(kmax,n);
00508 }
00509 Ipopt::Index mlp_ipt;
00510 options->GetIntegerValue("mlp", mlp_ipt,"filter.");
00511 mlp = mlp_ipt;
00512
00513 Ipopt::Index maxf_ipt;
00514 options->GetIntegerValue("maxf", maxf_ipt,"filter.");
00515 maxf = maxf_ipt;
00516
00517 fint mxwk0;
00518 Ipopt::Index mxwk0_ipt;
00519 options->GetIntegerValue("mxws", mxwk0_ipt, "filter.");
00520 mxwk0 = mxwk0_ipt;
00521
00522 fint mxiwk0;
00523 Ipopt::Index mxiwk0_ipt;
00524 options->GetIntegerValue("mxlws", mxiwk0_ipt, "filter.");
00525 mxiwk0 = mxiwk0_ipt;
00526
00527 int nplusm = n + m;
00528
00529 x = new real [n];
00530
00531
00532 use_warm_start_in_cache_ = false;
00533
00534 lam = new real [n+m];
00535 #ifdef InitializeAll
00536 for (int i = 0 ; i < n+m ; i++) lam[i] = 0.;
00537 #endif
00538
00539 bounds = new real [2*n + 2*m];
00540
00541 tnlp->get_bounds_info(n, bounds, bounds + nplusm, m, bounds + n, bounds + n + nplusm);
00542
00543 #if 0
00544 double infty = F77_FUNC_(nlp_eps_inf,NLP_EPS_INF).infty;
00545
00546 for (int i = 0 ; i < nplusm ; i++) {
00547 if (bounds[i] < -infty) bounds[i] = - infty;
00548 }
00549
00550 real * ubounds = bounds + nplusm;
00551 for (int i = 0 ; i < nplusm ; i++) {
00552 if (ubounds[i] > infty) ubounds[i] = infty;
00553 }
00554 #endif
00555
00556 maxa = n + nnz_jac_g;
00557 fint maxia = n + nnz_jac_g + m + 3;
00558 a = new real[maxa];
00559 la = new fint [maxia];
00560
00561 int * RowJac = new int [nnz_jac_g];
00562 int * ColJac = new int [nnz_jac_g];
00563
00564 la[nnz_jac_g + n + 1] = 1;
00565
00566 for (fint i = 1; i <= n ; i++)
00567 la[i] = i;
00568 tnlp->eval_jac_g( nv, NULL, 0, nc , nnz_j, RowJac, ColJac, NULL);
00569
00570 permutationJac = permutationJac_ = new int [nnz_jac_g];
00571 TMat2RowPMat(false, n, m, nnz_jac_g, RowJac, ColJac, permutationJac,
00572 la, n, 1, index_style);
00573
00574 delete [] RowJac;
00575 delete [] ColJac;
00576
00577
00578 permutationHess = permutationHess_ = new int[nnz_h];
00579 hStruct_ = new fint[nnz_h + n + 3];
00580 int * cache = new int[2*nnz_h + 1];
00581 F77_FUNC(hessc,HESSC).phl = 1;
00582 tnlp->eval_h((Ipopt::Index&) n, NULL, 0, 1., (Ipopt::Index&) m, NULL, 0, (Ipopt::Index&) nnz_h, cache + nnz_h, cache , NULL);
00583
00584 TMat2RowPMat(true, n, n, nnz_h, cache, cache + nnz_h, permutationHess,
00585 hStruct_, 0, 0, index_style);
00586
00587 delete [] cache;
00588
00589 fint lh1 = nnz_h + 8 + 2 * n + m;
00590 maxWk = 21*n + 8*m + mlp + 8*maxf + lh1 + kmax*(kmax+9)/2 + mxwk0;
00591 maxiWk = 13*n + 4*m + mlp + lh1 + kmax + 113 + mxiwk0;
00592
00593 ws = new real[maxWk];
00594 lws = new fint[maxiWk];
00595 #ifdef InitializeAll
00596 for (int i = 0 ; i < maxWk ; i++) ws[i] = 0;
00597 for (int i = 0 ; i < maxiWk ; i++) lws[i] = 0;
00598 #endif
00599
00600
00601 hStruct = hStruct_;
00602 tnlpSolved = static_cast<Ipopt::TNLP *>(Ipopt::GetRawPtr(tnlp));
00603
00604 options->GetNumericValue("ubd",F77_FUNC(ubdc,UBDC).ubd, "filter.");
00605 options->GetNumericValue("tt", F77_FUNC(ubdc,UBDC).tt, "filter.");
00606 options->GetNumericValue("eps", F77_FUNC_(nlp_eps_inf,NLP_EPS_INF).eps, "filter.");
00607 options->GetNumericValue("infty", F77_FUNC_(nlp_eps_inf,NLP_EPS_INF).infty, "filter.");
00608 rho = 10.;
00609 maxiter = 1000;
00610 options->GetIntegerValue("maxiter", (Ipopt::Index &) maxiter, "filter.");
00611 options->GetNumericValue("rho_init",rho,"filter.");
00612
00613
00614
00615 F77_FUNC(scalec,SCALEC).scale_mode = 0;
00616 s = new real [n+m];
00617
00618 istat = new fint[14];
00619 rstat = new real[7];
00620 #ifdef InitializeAll
00621 for (int i=0; i<14; i++) {
00622 istat[0] = 43;
00623 }
00624 for (int i=0; i<7; i++) {
00625 rstat[0] = 42.;
00626 }
00627 #endif
00628
00629 fmin = -1e100;
00630 Ipopt::Index bufy;
00631 options->GetIntegerValue("iprint",bufy, "filter.");
00632 iprint = bufy;
00633 nout = 6;
00634 cstype = new char[m];
00635 Ipopt::TNLP::LinearityType * const_types =
00636 new Ipopt::TNLP::LinearityType[m];
00637 tnlp->get_constraints_linearity(m, const_types);
00638 for (int i = 0 ; i < m ; i++) {
00639 if (const_types[i] == Ipopt::TNLP::LINEAR) {
00640 cstype[i] = 'L';
00641 }
00642 else
00643 cstype[i] = 'N';
00644 }
00645 delete [] const_types;
00646 c = new double[m];
00647 tnlp_ = Ipopt::GetRawPtr(tnlp);
00648 }
00649
00651 TNLPSolver::ReturnStatus
00652 FilterSolver::callOptimizer()
00653 {
00654 cached_->optimize();
00655
00656 TNLPSolver::ReturnStatus optimizationStatus = TNLPSolver::exception;
00657 Ipopt::SolverReturn status = Ipopt::INTERNAL_ERROR;
00658 fint ifail = cached_->ifail;
00659 switch (ifail) {
00660 case 0:
00661 optimizationStatus = TNLPSolver::solvedOptimal;
00662 status = Ipopt::SUCCESS;
00663 break;
00664 case 1:
00665 optimizationStatus = TNLPSolver::unbounded;
00666 status = Ipopt::DIVERGING_ITERATES;
00667 case 2:
00668 case 3:
00669 case 4:
00670 optimizationStatus = TNLPSolver::provenInfeasible;
00671 status = Ipopt::LOCAL_INFEASIBILITY;
00672 break;
00673 case 5:
00674 case 6:
00675 case 8:
00676 optimizationStatus = TNLPSolver::iterationLimit;
00677 status = Ipopt::MAXITER_EXCEEDED;
00678 break;
00679 case 7:
00680 optimizationStatus = TNLPSolver::externalException;
00681 status = Ipopt::INTERNAL_ERROR;
00682 break;
00683 case 9:
00684 case 10:
00685 optimizationStatus = TNLPSolver::exception;
00686 status = Ipopt::INTERNAL_ERROR;
00687 break;
00688 }
00689
00690 Number* mlam = NULL;
00691 if (cached_->m>0) {
00692 mlam = new Number[cached_->m];
00693 }
00694 for (int i = 0; i<cached_->m; i++) {
00695 mlam[i] = -cached_->lam[cached_->n+i];
00696 }
00697 Number* z_L = new Number[cached_->n];
00698 Number* z_U = new Number[cached_->n];
00699 const int os = cached_->n+cached_->m;
00700 for (int i=0; i<cached_->n; i++) {
00701 if (cached_->x[i] == cached_->bounds[i]) {
00702 z_L[i] = Max(0.,cached_->lam[i]);
00703 }
00704 else {
00705 z_L[i] = 0.;
00706 }
00707 if (cached_->x[i] == cached_->bounds[os+i]) {
00708 z_U[i] = Max(0.,-cached_->lam[i]);
00709 }
00710 else {
00711 z_U[i] = 0.;
00712 }
00713 }
00714 cached_->tnlp_->finalize_solution(status, cached_->n,
00715 cached_->x, z_L, z_U,
00716 cached_->m, cached_->c, mlam,
00717 cached_->f, NULL, NULL);
00718 delete [] mlam;
00719 delete [] z_L;
00720 delete [] z_U;
00721 return optimizationStatus;
00722 }
00724 void
00725 FilterSolver::cachedInfo::load_ws(Coin::SmartPtr<FilterWarmStart> warmF){
00726 if(!warmF.IsValid()) return;
00727 const fint xsize = warmF->primalSize();
00728 const real* xarray = warmF->primal();
00729 for (int i = 0; i<xsize; i++) {
00730 x[i] = xarray[i];
00731 }
00732 CoinCopyN(warmF->dual(), warmF->dualSize(), lam);
00733 CoinCopyN(warmF->lwsArray(), warmF->lwsSize(), lws);
00734 for (int i = 0 ; i < 14 ; i ++) {
00735 istat[i] = warmF->istat()[i];
00736 }
00737 use_warm_start_in_cache_ = true;
00738 }
00740 void
00741 FilterSolver::cachedInfo::optimize()
00742 {
00743 if (use_warm_start_in_cache_) {
00744 ifail = -1;
00745 use_warm_start_in_cache_ = false;
00746 }
00747 else {
00748 tnlp_->get_starting_point(n, 1, x, 0, NULL, NULL, m, 0, NULL);
00749 ifail = 0;
00750 }
00751 cpuTime_ = - CoinCpuTime();
00752 fint cstype_len = 1;
00753 rho = 10;
00754
00755
00756 #if 0
00757 printf("========= 3333333333333333 =============\n");
00758 for (int i=0; i<n; i++) {
00759 printf("xL[%3d] = %15.8e xU[%3d] = %15.8e\n", i, bounds[i], i, bounds[m+n+i]);
00760 }
00761 for (int i=0; i<m; i++) {
00762 printf("gL[%3d] = %15.8e gU[%3d] = %15.8e\n", i, bounds[n+i], i, bounds[m+2*n+i]);
00763 }
00764 #endif
00765 #if 0
00766 for (int i=0; i<n; i++) {
00767 printf("fxstart[%2d] = %23.16e\n", i, x[i]);
00768 }
00769 #endif
00770 F77_FUNC(filtersqp,FILTERSQP)(&n, &m, &kmax, & maxa, &maxf, &mlp, &maxWk,
00771 &maxiWk, &iprint, &nout, &ifail, &rho, x,
00772 c, &f, &fmin, bounds,
00773 bounds + n + m,
00774 s, a, la,
00775 ws, lws, lam, cstype,
00776 NULL, NULL,
00777 &maxiter, istat, rstat,
00778 cstype_len);
00779 #if 0
00780 for (int i=0; i<n; i++) {
00781 printf("fxsol[%2d] = %23.16e\n", i, x[i]);
00782 }
00783 #endif
00784 #if 0
00785 printf("final f = %e\n", f);
00786 printf("ifail = %d\n", ifail);
00787 #endif
00788 cpuTime_ += CoinCpuTime();
00789 }
00790
00791 std::string
00792 FilterSolver::UnsolvedFilterError::errorNames_[1] =
00793 {"Internal error in Filter SQP."};
00794
00795 std::string
00796 FilterSolver::UnsolvedFilterError::solverName_ =
00797 "filterSqp";
00798
00799 const std::string&
00800 FilterSolver::UnsolvedFilterError::errorName() const
00801 {
00802 return errorNames_[0];
00803 }
00804
00805 const std::string&
00806 FilterSolver::UnsolvedFilterError::solverName() const
00807 {
00808 return solverName_;
00809 }
00810
00811 bool
00812 FilterSolver::setWarmStart(const CoinWarmStart * warm,
00813 Ipopt::SmartPtr<TMINLP2TNLP> tnlp)
00814 {
00815 if (warm == NULL || cached_.IsNull()) {
00816 cached_ = new cachedInfo(GetRawPtr(tnlp), options_);
00817 }
00818 if(warm == NULL) return 1;
00819 const FilterWarmStart * warmF = dynamic_cast<const FilterWarmStart *> (warm);
00820 assert(warmF);
00821 if (warmF->empty())
00822 {
00823 warmF_ = NULL;
00824 disableWarmStart();
00825 return 1;
00826 }
00827 enableWarmStart();
00828 warmF_ = dynamic_cast<FilterWarmStart *> (warmF->clone());
00829 return true;
00830 }
00831
00832 CoinWarmStart *
00833 FilterSolver::getWarmStart(Ipopt::SmartPtr<TMINLP2TNLP> tnlp) const
00834 {
00835 return new FilterWarmStart(cached_->n, cached_->x,
00836 cached_->n+cached_->m, cached_->lam,
00837 cached_->maxiWk, cached_->lws, cached_->istat);
00838 }
00839
00840 CoinWarmStart *
00841 FilterSolver::getEmptyWarmStart() const
00842 {
00843 return new FilterWarmStart;
00844 }
00845
00846
00848 bool
00849 FilterSolver::warmStartIsValid(const CoinWarmStart * ws) const{
00850 const FilterWarmStart* fws = dynamic_cast<const FilterWarmStart*>(ws);
00851 if (fws && ! fws->empty()) {
00852 return true;
00853 }
00854 return false;
00855 }
00856
00857
00858 }