/home/coin/SVN-release/OS-2.2.0/Bonmin/src/Interfaces/Filter/BonFilterSolver.cpp

Go to the documentation of this file.
00001 // (C) Copyright International Business Machines Corporation, Carnegie Mellon University 2006, 2008
00002 // All Rights Reserved.
00003 // This code is published under the Common Public License.
00004 //
00005 // Authors :
00006 // Pierre Bonami, International Business Machines Corporation
00007 //
00008 // Date : 10/02/2006
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 //#define InitializeAll
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 //Static variables
00040 static Ipopt::TNLP * tnlpSolved = NULL;
00041 static fint nnz_h = -1;
00042 
00043 static fint * hStruct = NULL;
00044 
00045 //Permutation to apply to jacobian in order to get it row ordered
00046 static int * permutationJac = NULL;
00047 static int * permutationHess = NULL;
00048 //static int * cache = NULL;
00049 
00050 
00051 extern "C"
00052 {
00053 //Access to filter common bloc
00054   /* common block for problemname */
00055   extern struct {
00056       fint  char_l;
00057       char pname[10];
00058     }
00059   F77_FUNC(cpname,CPNAME);
00060 
00061   /* common block for Hessian storage set to 0, i.e. NO Hessian */
00062   extern struct {
00063       fint phl, phr, phc;
00064     }
00065   F77_FUNC(hessc,HESSC);
00066 
00067   /* common block for upper bound on filter */
00068   extern struct {
00069       real ubd, tt;
00070     }
00071   F77_FUNC(ubdc,UBDC);
00072 
00073   /* common block for infinity & epslon */
00074   extern struct {
00075       real infty, eps;
00076     }
00077   F77_FUNC_(nlp_eps_inf,NLP_EPS_INF);
00078 
00079   /* common block for prfinting from QP solver */
00080   extern struct {
00081       fint n_bqpd_calls, n_bqpd_prfint;
00082     }
00083   F77_FUNC_(bqpd_count,BQPD_COUNT);
00084 
00085   /* common for scaling: scale_mode = 0 (none), 1 (variables), 2 (vars+cons) */
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   /* evaluation of the Hessian of the Lagrangian */
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   // Convert a sparse matrix from triplet format to row ordered packed matrix
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       //DBG_ASSERT(RowJac[permutation2[i]] >= row);
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   // Convert a sparse matrix from triplet format to row ordered packed matrix
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       //DBG_ASSERT(RowJac[permutation2[i]] >= row);
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     //deleteme
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   // Convert a sparse matrix from triplet format to row ordered packed matrix
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     //deleteme
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     //  roption->AddLowerBoundedIntegerOption("timing", "whether to time evaluations (1 = yes)");
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 /* = false */)
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       // Set printlevel for stdout
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     //rescan bounds which may have changed
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 //Permutation to apply to jacobian in order to get it row ordered
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     // 1) Get some dimensions
00483     // 1.a) First from ampl
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     // 1.b) then from options
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     // Setup storage for Filter
00527     int nplusm = n + m;
00528     //Starting point
00529     x = new real [n];
00530 
00531     //tnlp->get_starting_point(n, 1, x, 0, NULL, NULL, m, 0, NULL);
00532     use_warm_start_in_cache_ = false;
00533     //for(int i = 0 ; i < n ; i++) x[i] = 0;
00534     lam = new real [n+m];
00535 #ifdef InitializeAll
00536     for (int i = 0 ; i < n+m ; i++) lam[i] = 0.;
00537 #endif
00538     //bounds
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     // AW: I don't think we need this, it isn't done for ReOptimize either
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;// - (index_style == Ipopt::TNLP::C_STYLE);
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     // Now setup hessian
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     // work arrays
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     // Setup global variables and static variables
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     // Set up scaling
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     //  rho = 1e6;
00755     //  printf("rho = %e\n", rho);
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())//reset initial point and leave
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 }//end namespace Bonmin

Generated on Thu Aug 5 03:02:55 2010 by  doxygen 1.4.7