Ipopt
trunk
|
00001 // Copyright (C) 2004, 2008 International Business Machines and others. 00002 // All Rights Reserved. 00003 // This code is published under the Eclipse Public License. 00004 // 00005 // $Id$ 00006 // 00007 // Authors: Carl Laird, Andreas Waechter IBM 2004-08-13 00008 00009 #ifndef __IPVECTOR_HPP__ 00010 #define __IPVECTOR_HPP__ 00011 00012 #include "IpTypes.hpp" 00013 #include "IpTaggedObject.hpp" 00014 #include "IpCachedResults.hpp" 00015 #include "IpSmartPtr.hpp" 00016 #include "IpJournalist.hpp" 00017 #include "IpException.hpp" 00018 00019 #include <vector> 00020 00021 namespace Ipopt 00022 { 00025 DECLARE_STD_EXCEPTION(UNIMPLEMENTED_LINALG_METHOD_CALLED); 00026 00027 /* forward declarations */ 00028 class VectorSpace; 00029 00047 class Vector : public TaggedObject 00048 { 00049 public: 00055 inline 00056 Vector(const VectorSpace* owner_space); 00057 00059 inline 00060 virtual ~Vector(); 00062 00064 inline 00065 Vector* MakeNew() const; 00066 00068 inline 00069 Vector* MakeNewCopy() const; 00070 00077 inline 00078 void Copy(const Vector& x); 00079 00081 void Scal(Number alpha); 00082 00084 inline 00085 void Axpy(Number alpha, const Vector &x); 00086 00088 inline 00089 Number Dot(const Vector &x) const; 00090 00092 inline 00093 Number Nrm2() const; 00094 00096 inline 00097 Number Asum() const; 00098 00100 inline 00101 Number Amax() const; 00103 00110 inline 00111 void Set(Number alpha); 00112 00114 inline 00115 void ElementWiseDivide(const Vector& x); 00116 00118 inline 00119 void ElementWiseMultiply(const Vector& x); 00120 00122 inline 00123 void ElementWiseMax(const Vector& x); 00124 00126 inline 00127 void ElementWiseMin(const Vector& x); 00128 00130 inline 00131 void ElementWiseReciprocal(); 00132 00134 inline 00135 void ElementWiseAbs(); 00136 00138 inline 00139 void ElementWiseSqrt(); 00140 00144 inline 00145 void ElementWiseSgn(); 00146 00148 inline 00149 void AddScalar(Number scalar); 00150 00152 inline 00153 Number Max() const; 00154 00156 inline 00157 Number Min() const; 00158 00160 inline 00161 Number Sum() const; 00162 00164 inline 00165 Number SumLogs() const; 00167 00175 inline 00176 void AddOneVector(Number a, const Vector& v1, Number c); 00177 00180 inline void AddTwoVectors(Number a, const Vector& v1, 00181 Number b, const Vector& v2, Number c); 00185 inline 00186 Number FracToBound(const Vector& delta, Number tau) const; 00188 inline 00189 void AddVectorQuotient(Number a, const Vector& z, const Vector& s, 00190 Number c); 00192 00195 inline 00196 bool HasValidNumbers() const; 00197 00201 inline 00202 Index Dim() const; 00203 00205 inline 00206 SmartPtr<const VectorSpace> OwnerSpace() const; 00208 00215 void Print(SmartPtr<const Journalist> jnlst, 00216 EJournalLevel level, 00217 EJournalCategory category, 00218 const std::string& name, 00219 Index indent=0, 00220 const std::string& prefix="") const; 00221 void Print(const Journalist& jnlst, 00222 EJournalLevel level, 00223 EJournalCategory category, 00224 const std::string& name, 00225 Index indent=0, 00226 const std::string& prefix="") const; 00228 00229 protected: 00235 virtual void CopyImpl(const Vector& x)=0; 00236 00238 virtual void ScalImpl(Number alpha)=0; 00239 00241 virtual void AxpyImpl(Number alpha, const Vector &x)=0; 00242 00244 virtual Number DotImpl(const Vector &x) const =0; 00245 00247 virtual Number Nrm2Impl() const =0; 00248 00250 virtual Number AsumImpl() const =0; 00251 00253 virtual Number AmaxImpl() const =0; 00254 00256 virtual void SetImpl(Number alpha)=0; 00257 00259 virtual void ElementWiseDivideImpl(const Vector& x)=0; 00260 00262 virtual void ElementWiseMultiplyImpl(const Vector& x)=0; 00263 00265 virtual void ElementWiseMaxImpl(const Vector& x)=0; 00266 00268 virtual void ElementWiseMinImpl(const Vector& x)=0; 00269 00271 virtual void ElementWiseReciprocalImpl()=0; 00272 00274 virtual void ElementWiseAbsImpl()=0; 00275 00277 virtual void ElementWiseSqrtImpl()=0; 00278 00280 virtual void ElementWiseSgnImpl()=0; 00281 00283 virtual void AddScalarImpl(Number scalar)=0; 00284 00286 virtual Number MaxImpl() const=0; 00287 00289 virtual Number MinImpl() const=0; 00290 00292 virtual Number SumImpl() const=0; 00293 00295 virtual Number SumLogsImpl() const=0; 00296 00299 virtual void AddTwoVectorsImpl(Number a, const Vector& v1, 00300 Number b, const Vector& v2, Number c); 00301 00303 virtual Number FracToBoundImpl(const Vector& delta, Number tau) const; 00304 00306 virtual void AddVectorQuotientImpl(Number a, const Vector& z, 00307 const Vector& s, Number c); 00308 00312 virtual bool HasValidNumbersImpl() const; 00313 00315 virtual void PrintImpl(const Journalist& jnlst, 00316 EJournalLevel level, 00317 EJournalCategory category, 00318 const std::string& name, 00319 Index indent, 00320 const std::string& prefix) const =0; 00322 00323 private: 00333 Vector(); 00334 00336 Vector(const Vector&); 00337 00339 Vector& operator=(const Vector&); 00341 00343 const SmartPtr<const VectorSpace> owner_space_; 00344 00348 mutable CachedResults<Number> dot_cache_; 00349 00350 mutable TaggedObject::Tag nrm2_cache_tag_; 00351 mutable Number cached_nrm2_; 00352 00353 mutable TaggedObject::Tag asum_cache_tag_; 00354 mutable Number cached_asum_; 00355 00356 mutable TaggedObject::Tag amax_cache_tag_; 00357 mutable Number cached_amax_; 00358 00359 mutable TaggedObject::Tag max_cache_tag_; 00360 mutable Number cached_max_; 00361 00362 mutable TaggedObject::Tag min_cache_tag_; 00363 mutable Number cached_min_; 00364 00365 mutable TaggedObject::Tag sum_cache_tag_; 00366 mutable Number cached_sum_; 00367 00368 mutable TaggedObject::Tag sumlogs_cache_tag_; 00369 mutable Number cached_sumlogs_; 00370 00371 mutable TaggedObject::Tag valid_cache_tag_; 00372 mutable bool cached_valid_; 00373 00374 // AW: I removed this cache since it gets in the way for the 00375 // quality function search 00376 // /** Cache for FracToBound */ 00377 // mutable CachedResults<Number> frac_to_bound_cache_; 00379 00380 }; 00381 00390 class VectorSpace : public ReferencedObject 00391 { 00392 public: 00398 VectorSpace(Index dim); 00399 00401 virtual ~VectorSpace() 00402 {} 00404 00408 virtual Vector* MakeNew() const=0; 00409 00411 Index Dim() const 00412 { 00413 return dim_; 00414 } 00415 00416 private: 00426 VectorSpace(); 00427 00429 VectorSpace(const VectorSpace&); 00430 00432 VectorSpace& operator=(const VectorSpace&); 00434 00436 const Index dim_; 00437 }; 00438 00439 /* inline methods */ 00440 inline 00441 Vector::~Vector() 00442 {} 00443 00444 inline 00445 Vector::Vector(const VectorSpace* owner_space) 00446 : 00447 TaggedObject(), 00448 owner_space_(owner_space), 00449 dot_cache_(10), 00450 nrm2_cache_tag_(0), 00451 asum_cache_tag_(0), 00452 amax_cache_tag_(0), 00453 max_cache_tag_(0), 00454 min_cache_tag_(0), 00455 sum_cache_tag_(0), 00456 sumlogs_cache_tag_(0), 00457 cached_valid_(0) 00458 { 00459 DBG_ASSERT(IsValid(owner_space_)); 00460 } 00461 00462 inline 00463 Vector* Vector::MakeNew() const 00464 { 00465 return owner_space_->MakeNew(); 00466 } 00467 00468 inline 00469 Vector* Vector::MakeNewCopy() const 00470 { 00471 // ToDo: We can probably copy also the cached values for Norms etc here 00472 Vector* copy = MakeNew(); 00473 copy->Copy(*this); 00474 return copy; 00475 } 00476 00477 inline 00478 void Vector::Copy(const Vector& x) 00479 { 00480 CopyImpl(x); 00481 ObjectChanged(); 00482 // Also copy any cached scalar values from the original vector 00483 // ToDo: Check if that is too much overhead 00484 TaggedObject::Tag x_tag = x.GetTag(); 00485 if (x_tag == x.nrm2_cache_tag_) { 00486 nrm2_cache_tag_ = GetTag(); 00487 cached_nrm2_ = x.cached_nrm2_; 00488 } 00489 if (x_tag == x.asum_cache_tag_) { 00490 asum_cache_tag_ = GetTag(); 00491 cached_asum_ = x.cached_asum_; 00492 } 00493 if (x_tag == x.amax_cache_tag_) { 00494 amax_cache_tag_ = GetTag(); 00495 cached_amax_ = x.cached_amax_; 00496 } 00497 if (x_tag == x.max_cache_tag_) { 00498 max_cache_tag_ = GetTag(); 00499 cached_max_ = x.cached_max_; 00500 } 00501 if (x_tag == x.min_cache_tag_) { 00502 min_cache_tag_ = GetTag(); 00503 cached_min_ = x.cached_min_; 00504 } 00505 if (x_tag == x.sum_cache_tag_) { 00506 sum_cache_tag_ = GetTag(); 00507 cached_sum_ = x.cached_sum_; 00508 } 00509 if (x_tag == x.sumlogs_cache_tag_) { 00510 sumlogs_cache_tag_ = GetTag(); 00511 cached_sumlogs_ = x.cached_sumlogs_; 00512 } 00513 } 00514 00515 inline 00516 void Vector::Axpy(Number alpha, const Vector &x) 00517 { 00518 AxpyImpl(alpha, x); 00519 ObjectChanged(); 00520 } 00521 00522 inline 00523 Number Vector::Dot(const Vector &x) const 00524 { 00525 // The current implementation of the caching doesn't allow to have 00526 // a dependency of something with itself. Therefore, we use the 00527 // Nrm2 method if the dot product is to be taken with the vector 00528 // itself. Might be more efficient anyway. 00529 if (this==&x) { 00530 Number nrm2 = Nrm2(); 00531 return nrm2*nrm2; 00532 } 00533 Number retValue; 00534 if (!dot_cache_.GetCachedResult2Dep(retValue, this, &x)) { 00535 retValue = DotImpl(x); 00536 dot_cache_.AddCachedResult2Dep(retValue, this, &x); 00537 } 00538 return retValue; 00539 } 00540 00541 inline 00542 Number Vector::Nrm2() const 00543 { 00544 if (nrm2_cache_tag_ != GetTag()) { 00545 cached_nrm2_ = Nrm2Impl(); 00546 nrm2_cache_tag_ = GetTag(); 00547 } 00548 return cached_nrm2_; 00549 } 00550 00551 inline 00552 Number Vector::Asum() const 00553 { 00554 if (asum_cache_tag_ != GetTag()) { 00555 cached_asum_ = AsumImpl(); 00556 asum_cache_tag_ = GetTag(); 00557 } 00558 return cached_asum_; 00559 } 00560 00561 inline 00562 Number Vector::Amax() const 00563 { 00564 if (amax_cache_tag_ != GetTag()) { 00565 cached_amax_ = AmaxImpl(); 00566 amax_cache_tag_ = GetTag(); 00567 } 00568 return cached_amax_; 00569 } 00570 00571 inline 00572 Number Vector::Sum() const 00573 { 00574 if (sum_cache_tag_ != GetTag()) { 00575 cached_sum_ = SumImpl(); 00576 sum_cache_tag_ = GetTag(); 00577 } 00578 return cached_sum_; 00579 } 00580 00581 inline 00582 Number Vector::SumLogs() const 00583 { 00584 if (sumlogs_cache_tag_ != GetTag()) { 00585 cached_sumlogs_ = SumLogsImpl(); 00586 sumlogs_cache_tag_ = GetTag(); 00587 } 00588 return cached_sumlogs_; 00589 } 00590 00591 inline 00592 void Vector::ElementWiseSgn() 00593 { 00594 ElementWiseSgnImpl(); 00595 ObjectChanged(); 00596 } 00597 00598 inline 00599 void Vector::Set(Number alpha) 00600 { 00601 // Could initialize caches here 00602 SetImpl(alpha); 00603 ObjectChanged(); 00604 } 00605 00606 inline 00607 void Vector::ElementWiseDivide(const Vector& x) 00608 { 00609 ElementWiseDivideImpl(x); 00610 ObjectChanged(); 00611 } 00612 00613 inline 00614 void Vector::ElementWiseMultiply(const Vector& x) 00615 { 00616 ElementWiseMultiplyImpl(x); 00617 ObjectChanged(); 00618 } 00619 00620 inline 00621 void Vector::ElementWiseReciprocal() 00622 { 00623 ElementWiseReciprocalImpl(); 00624 ObjectChanged(); 00625 } 00626 00627 inline 00628 void Vector::ElementWiseMax(const Vector& x) 00629 { 00630 // Could initialize some caches here 00631 ElementWiseMaxImpl(x); 00632 ObjectChanged(); 00633 } 00634 00635 inline 00636 void Vector::ElementWiseMin(const Vector& x) 00637 { 00638 // Could initialize some caches here 00639 ElementWiseMinImpl(x); 00640 ObjectChanged(); 00641 } 00642 00643 inline 00644 void Vector::ElementWiseAbs() 00645 { 00646 // Could initialize some caches here 00647 ElementWiseAbsImpl(); 00648 ObjectChanged(); 00649 } 00650 00651 inline 00652 void Vector::ElementWiseSqrt() 00653 { 00654 ElementWiseSqrtImpl(); 00655 ObjectChanged(); 00656 } 00657 00658 inline 00659 void Vector::AddScalar(Number scalar) 00660 { 00661 // Could initialize some caches here 00662 AddScalarImpl(scalar); 00663 ObjectChanged(); 00664 } 00665 00666 inline 00667 Number Vector::Max() const 00668 { 00669 if (max_cache_tag_ != GetTag()) { 00670 cached_max_ = MaxImpl(); 00671 max_cache_tag_ = GetTag(); 00672 } 00673 return cached_max_; 00674 } 00675 00676 inline 00677 Number Vector::Min() const 00678 { 00679 if (min_cache_tag_ != GetTag()) { 00680 cached_min_ = MinImpl(); 00681 min_cache_tag_ = GetTag(); 00682 } 00683 return cached_min_; 00684 } 00685 00686 inline 00687 void Vector::AddOneVector(Number a, const Vector& v1, Number c) 00688 { 00689 AddTwoVectors(a, v1, 0., v1, c); 00690 } 00691 00692 inline 00693 void Vector::AddTwoVectors(Number a, const Vector& v1, 00694 Number b, const Vector& v2, Number c) 00695 { 00696 AddTwoVectorsImpl(a, v1, b, v2, c); 00697 ObjectChanged(); 00698 } 00699 00700 inline 00701 Number Vector::FracToBound(const Vector& delta, Number tau) const 00702 { 00703 /* AW: I avoid the caching here, since it leads to overhead in the 00704 quality function search. Caches for this are in 00705 CalculatedQuantities. 00706 Number retValue; 00707 std::vector<const TaggedObject*> tdeps(1); 00708 tdeps[0] = δ 00709 std::vector<Number> sdeps(1); 00710 sdeps[0] = tau; 00711 if (!frac_to_bound_cache_.GetCachedResult(retValue, tdeps, sdeps)) { 00712 retValue = FracToBoundImpl(delta, tau); 00713 frac_to_bound_cache_.AddCachedResult(retValue, tdeps, sdeps); 00714 } 00715 return retValue; 00716 */ 00717 return FracToBoundImpl(delta, tau); 00718 } 00719 00720 inline 00721 void Vector::AddVectorQuotient(Number a, const Vector& z, 00722 const Vector& s, Number c) 00723 { 00724 AddVectorQuotientImpl(a, z, s, c); 00725 ObjectChanged(); 00726 } 00727 00728 inline 00729 bool Vector::HasValidNumbers() const 00730 { 00731 if (valid_cache_tag_ != GetTag()) { 00732 cached_valid_ = HasValidNumbersImpl(); 00733 valid_cache_tag_ = GetTag(); 00734 } 00735 return cached_valid_; 00736 } 00737 00738 inline 00739 Index Vector::Dim() const 00740 { 00741 return owner_space_->Dim(); 00742 } 00743 00744 inline 00745 SmartPtr<const VectorSpace> Vector::OwnerSpace() const 00746 { 00747 return owner_space_; 00748 } 00749 00750 inline 00751 VectorSpace::VectorSpace(Index dim) 00752 : 00753 dim_(dim) 00754 {} 00755 00756 } // namespace Ipopt 00757 00758 // Macro definitions for debugging vectors 00759 #if COIN_IPOPT_VERBOSITY == 0 00760 # define DBG_PRINT_VECTOR(__verbose_level, __vec_name, __vec) 00761 #else 00762 # define DBG_PRINT_VECTOR(__verbose_level, __vec_name, __vec) \ 00763 if (dbg_jrnl.Verbosity() >= (__verbose_level)) { \ 00764 if (dbg_jrnl.Jnlst()!=NULL) { \ 00765 (__vec).Print(dbg_jrnl.Jnlst(), \ 00766 J_ERROR, J_DBG, \ 00767 __vec_name, \ 00768 dbg_jrnl.IndentationLevel()*2, \ 00769 "# "); \ 00770 } \ 00771 } 00772 #endif //if COIN_IPOPT_VERBOSITY == 0 00773 00774 #endif