MOAB: Mesh Oriented datABase
(version 5.3.1)
|
00001 /* ***************************************************************** 00002 MESQUITE -- The Mesh Quality Improvement Toolkit 00003 00004 Copyright 2007 Sandia National Laboratories. Developed at the 00005 University of Wisconsin--Madison under SNL contract number 00006 624796. The U.S. Government and the University of Wisconsin 00007 retain certain rights to this software. 00008 00009 This library is free software; you can redistribute it and/or 00010 modify it under the terms of the GNU Lesser General Public 00011 License as published by the Free Software Foundation; either 00012 version 2.1 of the License, or (at your option) any later version. 00013 00014 This library is distributed in the hope that it will be useful, 00015 but WITHOUT ANY WARRANTY; without even the implied warranty of 00016 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00017 Lesser General Public License for more details. 00018 00019 You should have received a copy of the GNU Lesser General Public License 00020 (lgpl.txt) along with this library; if not, write to the Free Software 00021 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00022 00023 (2008) kraftche@cae.wisc.edu 00024 00025 ***************************************************************** */ 00026 00027 /** \file TrustRegion.cpp 00028 * \brief Port Todd Munson's trust region solver to Mesquite 00029 * \author Jason Kraftcheck (Mesquite Port) 00030 */ 00031 00032 #include "Mesquite.hpp" 00033 #include "TrustRegion.hpp" 00034 #include "MsqDebug.hpp" 00035 #include "MsqError.hpp" 00036 #include "PatchData.hpp" 00037 00038 #define USE_FN_PC1 // Use 1st preconditioner from Todd's code 00039 // (alternate is whatever is in MsqHessian already) 00040 #undef DO_STEEP_DESC // Jason's apparently broken hack to fall back to 00041 // steepest descent search direction 00042 00043 namespace MBMesquite 00044 { 00045 00046 // Force std::vector to release allocated memory 00047 template < typename T > 00048 static inline void free_vector( std::vector< T >& v ) 00049 { 00050 std::vector< T > temp; 00051 temp.swap( v ); 00052 } 00053 00054 std::string TrustRegion::get_name() const 00055 { 00056 return "TrustRegion"; 00057 } 00058 00059 PatchSet* TrustRegion::get_patch_set() 00060 { 00061 return PatchSetUser::get_patch_set(); 00062 } 00063 00064 TrustRegion::TrustRegion( ObjectiveFunction* of ) : VertexMover( of ), PatchSetUser( true ), mMemento( 0 ) {} 00065 00066 TrustRegion::~TrustRegion() 00067 { 00068 delete mMemento; 00069 mMemento = 0; 00070 } 00071 00072 void TrustRegion::initialize( PatchData& pd, MsqError& err ) 00073 { 00074 mMemento = pd.create_vertices_memento( err );MSQ_CHKERR( err ); 00075 } 00076 00077 void TrustRegion::initialize_mesh_iteration( PatchData& /*pd*/, MsqError& /*err*/ ) {} 00078 00079 void TrustRegion::terminate_mesh_iteration( PatchData& /*pd*/, MsqError& /*err*/ ) {} 00080 00081 void TrustRegion::cleanup() 00082 { 00083 // release Memento 00084 delete mMemento; 00085 mMemento = 0; 00086 // release temporary array memory 00087 mHess.clear(); 00088 free_vector( mGrad ); 00089 free_vector( wVect ); 00090 free_vector( zVect ); 00091 free_vector( dVect ); 00092 free_vector( pVect ); 00093 free_vector( rVect ); 00094 free_vector( preCond ); 00095 } 00096 00097 static inline void negate( Vector3D* out, const Vector3D* in, size_t nn ) 00098 { 00099 for( size_t i = 0; i < nn; ++i ) 00100 out[i] = -in[i]; 00101 } 00102 00103 // Do v += s * x, where v and x are arrays of length n 00104 static inline void plus_eq_scaled( Vector3D* v, double s, const Vector3D* x, size_t n ) 00105 { 00106 Vector3D* end = v + n; 00107 for( ; v != end; ++v, ++x ) 00108 *v += s * *x; 00109 } 00110 00111 // Do v = s*v - x, where v and x are arrays of length n 00112 static inline void times_eq_minus( Vector3D* v, double s, const Vector3D* x, size_t n ) 00113 { 00114 Vector3D* end = v + n; 00115 for( ; v != end; ++v, ++x ) 00116 { 00117 *v *= s; 00118 *v -= *x; 00119 } 00120 } 00121 00122 void TrustRegion::compute_preconditioner( MsqError& 00123 #ifndef USE_FN_PC1 00124 err 00125 #endif 00126 ) 00127 { 00128 #ifndef USE_FN_PC1 00129 mHessian.calculate_preconditioner( err ); 00130 #else 00131 double dia; 00132 preCond.resize( mHess.size() ); 00133 for( size_t i = 0; i < mHess.size(); ++i ) 00134 { 00135 const Matrix3D& m = *mHess.get_block( i, i ); 00136 dia = m[0][0] + m[1][1] + m[2][2]; 00137 preCond[i] = dia < DBL_EPSILON ? 1.0 : 1.0 / dia; 00138 } 00139 #endif 00140 } 00141 00142 void TrustRegion::apply_preconditioner( Vector3D* z, Vector3D* r, 00143 MsqError& 00144 #ifndef USE_FN_PC1 00145 err 00146 #endif 00147 ) 00148 { 00149 #ifndef USE_FN_PC1 00150 mHessian.apply_preconditioner( z, r, err ); 00151 #else 00152 for( size_t i = 0; i < preCond.size(); ++i ) 00153 z[i] = preCond[i] * r[i]; 00154 #endif 00155 } 00156 00157 void TrustRegion::optimize_vertex_positions( PatchData& pd, MsqError& err ) 00158 { 00159 TerminationCriterion& term = *get_inner_termination_criterion(); 00160 OFEvaluator& func = get_objective_function_evaluator(); 00161 00162 const double cg_tol = 1e-2; 00163 const double eta_1 = 0.01; 00164 const double eta_2 = 0.90; 00165 const double tr_incr = 10; 00166 const double tr_decr_def = 0.25; 00167 const double tr_decr_undef = 0.25; 00168 const double tr_num_tol = 1e-6; 00169 const int max_cg_iter = 10000; 00170 00171 double radius = 1000; /* delta*delta */ 00172 00173 const int nn = pd.num_free_vertices(); 00174 wVect.resize( nn ); 00175 Vector3D* w = arrptr( wVect ); 00176 zVect.resize( nn ); 00177 Vector3D* z = arrptr( zVect ); 00178 dVect.resize( nn ); 00179 Vector3D* d = arrptr( dVect ); 00180 pVect.resize( nn ); 00181 Vector3D* p = arrptr( pVect ); 00182 rVect.resize( nn ); 00183 Vector3D* r = arrptr( rVect ); 00184 00185 double norm_r, norm_g; 00186 double alpha, beta, kappa; 00187 double rz, rzm1; 00188 double dMp, norm_d, norm_dp1, norm_p; 00189 double obj, objn; 00190 00191 int cg_iter; 00192 bool valid; 00193 00194 mHess.initialize( pd, err ); // hMesh(mesh); 00195 valid = func.update( pd, obj, mGrad, mHess, err );MSQ_ERRRTN( err ); 00196 if( !valid ) 00197 { 00198 MSQ_SETERR( err )( "Initial objective function is not valid", MsqError::INVALID_MESH ); 00199 return; 00200 } 00201 compute_preconditioner( err );MSQ_ERRRTN( err ); 00202 pd.recreate_vertices_memento( mMemento, err );MSQ_ERRRTN( err ); 00203 00204 while( !term.terminate() && ( radius > 1e-20 ) ) 00205 { 00206 00207 norm_r = length_squared( arrptr( mGrad ), nn ); 00208 norm_g = sqrt( norm_r ); 00209 00210 memset( d, 0, 3 * sizeof( double ) * nn ); 00211 memcpy( r, arrptr( mGrad ), 00212 nn * sizeof( Vector3D ) ); // memcpy(r, mesh->g, 3*sizeof(double)*nn); 00213 norm_g *= cg_tol; 00214 00215 apply_preconditioner( z, r, err );MSQ_ERRRTN( err ); // prec->apply(z, r, prec, mesh); 00216 negate( p, z, nn ); 00217 rz = inner( r, z, nn ); 00218 00219 dMp = 0; 00220 norm_p = rz; 00221 norm_d = 0; 00222 00223 cg_iter = 0; 00224 while( ( sqrt( norm_r ) > norm_g ) && 00225 #ifdef DO_STEEP_DESC 00226 ( norm_d > tr_num_tol ) && 00227 #endif 00228 ( cg_iter < max_cg_iter ) ) 00229 { 00230 ++cg_iter; 00231 00232 memset( w, 0, 3 * sizeof( double ) * nn ); 00233 // matmul(w, mHess, p); //matmul(w, mesh, p); 00234 mHess.product( w, p ); 00235 00236 kappa = inner( p, w, nn ); 00237 if( kappa <= 0.0 ) 00238 { 00239 alpha = ( sqrt( dMp * dMp + norm_p * ( radius - norm_d ) ) - dMp ) / norm_p; 00240 plus_eq_scaled( d, alpha, p, nn ); 00241 break; 00242 } 00243 00244 alpha = rz / kappa; 00245 00246 norm_dp1 = norm_d + 2.0 * alpha * dMp + alpha * alpha * norm_p; 00247 if( norm_dp1 >= radius ) 00248 { 00249 alpha = ( sqrt( dMp * dMp + norm_p * ( radius - norm_d ) ) - dMp ) / norm_p; 00250 plus_eq_scaled( d, alpha, p, nn ); 00251 break; 00252 } 00253 00254 plus_eq_scaled( d, alpha, p, nn ); 00255 plus_eq_scaled( r, alpha, w, nn ); 00256 norm_r = length_squared( r, nn ); 00257 00258 apply_preconditioner( z, r, err );MSQ_ERRRTN( err ); // prec->apply(z, r, prec, mesh); 00259 00260 rzm1 = rz; 00261 rz = inner( r, z, nn ); 00262 beta = rz / rzm1; 00263 times_eq_minus( p, beta, z, nn ); 00264 00265 dMp = beta * ( dMp + alpha * norm_p ); 00266 norm_p = rz + beta * beta * norm_p; 00267 norm_d = norm_dp1; 00268 } 00269 00270 #ifdef DO_STEEP_DESC 00271 if( norm_d <= tr_num_tol ) 00272 { 00273 norm_g = length( arrptr( mGrad ), nn ); 00274 double ll = 1.0; 00275 if( norm_g < tr_num_tol ) break; 00276 if( norm_g > radius ) ll = radius / nurm_g; 00277 for( int i = 0; i < nn; ++i ) 00278 d[i] = ll * mGrad[i]; 00279 } 00280 #endif 00281 00282 alpha = inner( arrptr( mGrad ), d, nn ); // inner(mesh->g, d, nn); 00283 00284 memset( p, 0, 3 * sizeof( double ) * nn ); 00285 // matmul(p, mHess, d); //matmul(p, mesh, d); 00286 mHess.product( p, d ); 00287 beta = 0.5 * inner( p, d, nn ); 00288 kappa = alpha + beta; 00289 00290 /* Put the new point into the locations */ 00291 pd.move_free_vertices_constrained( d, nn, 1.0, err );MSQ_ERRRTN( err ); 00292 00293 valid = func.evaluate( pd, objn, err ); 00294 if( err.error_code() == err.BARRIER_VIOLATED ) 00295 err.clear(); // barrier violated does not represent an actual error here 00296 MSQ_ERRRTN( err ); 00297 00298 if( !valid ) 00299 { 00300 /* Function not defined at trial point */ 00301 radius *= tr_decr_undef; 00302 pd.set_to_vertices_memento( mMemento, err );MSQ_ERRRTN( err ); 00303 continue; 00304 } 00305 00306 if( ( fabs( kappa ) <= tr_num_tol ) && ( fabs( objn - obj ) <= tr_num_tol ) ) { kappa = 1; } 00307 else 00308 { 00309 kappa = ( objn - obj ) / kappa; 00310 } 00311 00312 if( kappa < eta_1 ) 00313 { 00314 /* Iterate is unacceptable */ 00315 radius *= tr_decr_def; 00316 pd.set_to_vertices_memento( mMemento, err );MSQ_ERRRTN( err ); 00317 continue; 00318 } 00319 00320 /* Iterate is acceptable */ 00321 if( kappa >= eta_2 ) 00322 { 00323 /* Iterate is a very good step, increase radius */ 00324 radius *= tr_incr; 00325 if( radius > 1e20 ) { radius = 1e20; } 00326 } 00327 00328 func.update( pd, obj, mGrad, mHess, err ); 00329 compute_preconditioner( err );MSQ_ERRRTN( err ); 00330 pd.recreate_vertices_memento( mMemento, err );MSQ_ERRRTN( err ); 00331 00332 // checks stopping criterion 00333 term.accumulate_patch( pd, err );MSQ_ERRRTN( err ); 00334 term.accumulate_inner( pd, objn, arrptr( mGrad ), err );MSQ_ERRRTN( err ); 00335 } 00336 } 00337 00338 } // namespace MBMesquite