MOAB: Mesh Oriented datABase  (version 5.2.1)
QuasiNewton.cpp
Go to the documentation of this file.
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 QuasiNewton.cpp
00028  *  \brief Port Todd Munson's quasi-Newton solver to Mesquite
00029  *  \author Jason Kraftcheck (Mesquite Port)
00030  */
00031 
00032 #include "Mesquite.hpp"
00033 #include "QuasiNewton.hpp"
00034 #include "MsqDebug.hpp"
00035 #include "MsqError.hpp"
00036 #include "PatchData.hpp"
00037 
00038 namespace MBMesquite
00039 {
00040 
00041 // Force std::vector to release allocated memory
00042 template < typename T >
00043 static inline void free_vector( std::vector< T >& v )
00044 {
00045     std::vector< T > temp;
00046     temp.swap( v );
00047 }
00048 
00049 std::string QuasiNewton::get_name() const
00050 {
00051     return "QuasiNewton";
00052 }
00053 
00054 PatchSet* QuasiNewton::get_patch_set()
00055 {
00056     return PatchSetUser::get_patch_set();
00057 }
00058 
00059 QuasiNewton::QuasiNewton( ObjectiveFunction* of ) : VertexMover( of ), PatchSetUser( true ), mMemento( 0 ) {}
00060 
00061 QuasiNewton::~QuasiNewton()
00062 {
00063     delete mMemento;
00064     mMemento = 0;
00065 }
00066 
00067 void QuasiNewton::initialize( PatchData& pd, MsqError& err )
00068 {
00069     if( !mMemento )
00070     {
00071         mMemento = pd.create_vertices_memento( err );MSQ_CHKERR( err );
00072     }
00073 }
00074 
00075 void QuasiNewton::initialize_mesh_iteration( PatchData& /*pd*/, MsqError& /*err*/ ) {}
00076 
00077 void QuasiNewton::terminate_mesh_iteration( PatchData& /*pd*/, MsqError& /*err*/ ) {}
00078 
00079 void QuasiNewton::cleanup()
00080 {
00081     // release memento
00082     delete mMemento;
00083     mMemento = 0;
00084 
00085     // release coordinates
00086     for( size_t i = 0; i < ( sizeof( w ) / sizeof( w[0] ) ); ++i )
00087         free_vector( w[i] );
00088     // release gradients
00089     for( size_t i = 0; i < ( sizeof( v ) / sizeof( v[0] ) ); ++i )
00090         free_vector( v[i] );
00091 
00092     // release Hessian memory
00093     free_vector( mHess );
00094 
00095     // release temporary array memory
00096     free_vector( x );
00097     free_vector( d );
00098 }
00099 
00100 // Do v += s * x, where v and x are arrays of length n
00101 static inline void plus_eq_scaled( Vector3D* v, double s, const Vector3D* x, size_t n )
00102 {
00103     Vector3D* end = v + n;
00104     for( ; v != end; ++v, ++x )
00105         *v += s * *x;
00106 }
00107 
00108 void QuasiNewton::solve( Vector3D* z_arr, const Vector3D* v_arr ) const
00109 {
00110     SymMatrix3D pd;
00111 
00112     const double small = DBL_EPSILON;
00113     const size_t nn    = mHess.size();
00114     for( size_t i = 0; i < nn; ++i )
00115     {
00116 
00117         // ensure positive definite: perturb a bit if
00118         // diagonal values are zero.
00119         SymMatrix3D dd = mHess[i];
00120         while( fabs( dd[0] ) < small || fabs( dd[3] ) < small || fabs( dd[5] ) < small )
00121             dd += small;
00122 
00123         // factor
00124         pd[0] = 1.0 / dd[0];
00125         pd[1] = dd[1] * pd[0];
00126         pd[2] = dd[2] * pd[0];
00127 
00128         pd[3] = 1.0 / ( dd[3] - dd[1] * pd[1] );
00129         pd[5] = dd[4] - dd[2] * pd[1];
00130         pd[4] = pd[3] * pd[5];
00131         pd[5] = 1.0 / ( dd[5] - dd[2] * pd[2] - pd[4] * pd[5] );
00132 
00133         if( pd[0] <= 0.0 || pd[3] <= 0.0 || pd[5] <= 0.0 )
00134         {
00135             if( dd[0] + dd[3] + dd[5] <= 0 )
00136             {
00137                 // switch to diagonal
00138                 pd[0] = 1.0 / fabs( dd[0] );
00139                 pd[1] = 0.0;
00140                 pd[2] = 0.0;
00141                 pd[3] = 1.0 / fabs( dd[3] );
00142                 pd[4] = 0.0;
00143                 pd[5] = 1.0 / fabs( dd[5] );
00144             }
00145             else
00146             {
00147                 // diagonal preconditioner
00148                 pd[0] = pd[3] = pd[5] = 1.0 / ( dd[0] + dd[3] + dd[5] );
00149                 pd[1] = pd[2] = pd[4] = 0.0;
00150             }
00151         }
00152 
00153         // solve
00154         const Vector3D& vv = v_arr[i];
00155         Vector3D& z        = z_arr[i];
00156         z[0]               = vv[0];
00157         z[1]               = vv[1] - pd[1] * z[0];
00158         z[2]               = vv[2] - pd[2] * z[0] - pd[4] * z[1];
00159 
00160         z[0] *= pd[0];
00161         z[1] *= pd[3];
00162         z[2] *= pd[5];
00163 
00164         z[1] -= pd[4] * z[2];
00165         z[0] -= pd[1] * z[1] + pd[2] * z[2];
00166     }
00167 }
00168 
00169 void QuasiNewton::optimize_vertex_positions( PatchData& pd, MsqError& err )
00170 {
00171     TerminationCriterion& term = *get_inner_termination_criterion();
00172     OFEvaluator& func          = get_objective_function_evaluator();
00173 
00174     const double sigma   = 1e-4;
00175     const double beta0   = 0.25;
00176     const double beta1   = 0.80;
00177     const double tol1    = 1e-8;
00178     const double epsilon = 1e-10;
00179 
00180     // double norm_r; //, norm_g;
00181     double alpha, beta;
00182     double obj, objn;
00183 
00184     size_t i;
00185 
00186     // Initialize stuff
00187     const size_t nn = pd.num_free_vertices();
00188     double a[QNVEC], b[QNVEC], r[QNVEC];
00189     for( i = 0; i < QNVEC; ++i )
00190         r[i] = 0;
00191     for( i = 0; i <= QNVEC; ++i )
00192     {
00193         v[i].clear();
00194         v[i].resize( nn, Vector3D( 0.0 ) );
00195         w[i].clear();
00196         w[i].resize( nn, Vector3D( 0.0 ) );
00197     }
00198     d.resize( nn );
00199     mHess.resize( nn );  // hMesh(mesh);
00200 
00201     bool valid = func.update( pd, obj, v[QNVEC], mHess, err );MSQ_ERRRTN( err );
00202     if( !valid )
00203     {
00204         MSQ_SETERR( err )( "Initial objective function is not valid", MsqError::INVALID_MESH );
00205         return;
00206     }
00207 
00208     while( !term.terminate() )
00209     {
00210         pd.recreate_vertices_memento( mMemento, err );MSQ_ERRRTN( err );
00211         pd.get_free_vertex_coordinates( w[QNVEC] );
00212 
00213         x = v[QNVEC];
00214         for( i = QNVEC; i--; )
00215         {
00216             a[i] = r[i] * inner( &( w[i][0] ), arrptr( x ), nn );
00217             plus_eq_scaled( arrptr( x ), -a[i], &v[i][0], nn );
00218         }
00219 
00220         solve( arrptr( d ), arrptr( x ) );
00221 
00222         for( i = QNVEC; i--; )
00223         {
00224             b[i] = r[i] * inner( &( v[i][0] ), arrptr( d ), nn );
00225             plus_eq_scaled( arrptr( d ), a[i] - b[i], &( w[i][0] ), nn );
00226         }
00227 
00228         alpha = -inner( &( v[QNVEC][0] ), arrptr( d ), nn ); /* direction is negated */
00229         if( alpha > 0.0 )
00230         {
00231             MSQ_SETERR( err )( "No descent.", MsqError::INVALID_MESH );
00232             return;
00233         }
00234 
00235         alpha *= sigma;
00236         beta = 1.0;
00237 
00238         pd.move_free_vertices_constrained( arrptr( d ), nn, -beta, err );MSQ_ERRRTN( err );
00239         valid = func.evaluate( pd, objn, v[QNVEC], err );
00240         if( err.error_code() == err.BARRIER_VIOLATED )
00241             err.clear();  // barrier violated does not represent an actual error here
00242         MSQ_ERRRTN( err );
00243         if( !valid || ( obj - objn < -alpha * beta - epsilon && length( &( v[QNVEC][0] ), nn ) >= tol1 ) )
00244         {
00245 
00246             if( !valid )  // function not defined at trial point
00247                 beta *= beta0;
00248             else  // unacceptable iterate
00249                 beta *= beta1;
00250 
00251             for( ;; )
00252             {
00253                 if( beta < tol1 )
00254                 {
00255                     pd.set_to_vertices_memento( mMemento, err );MSQ_ERRRTN( err );
00256                     MSQ_SETERR( err )( "Newton step not good", MsqError::INTERNAL_ERROR );
00257                     return;
00258                 }
00259 
00260                 pd.set_free_vertices_constrained( mMemento, arrptr( d ), nn, -beta, err );MSQ_ERRRTN( err );
00261                 valid = func.evaluate( pd, objn, err );
00262                 if( err.error_code() == err.BARRIER_VIOLATED )
00263                     err.clear();  // barrier violated does not represent an actual error here
00264                 MSQ_ERRRTN( err );
00265                 if( !valid )  // function undefined at trial point
00266                     beta *= beta0;
00267                 else if( obj - objn < -alpha * beta - epsilon )  // unacceptlable iterate
00268                     beta *= beta1;
00269                 else
00270                     break;
00271             }
00272         }
00273 
00274         for( i = 0; i < QNVEC - 1; ++i )
00275         {
00276             r[i] = r[i + 1];
00277             w[i].swap( w[i + 1] );
00278             v[i].swap( v[i + 1] );
00279         }
00280         w[QNVEC - 1].swap( w[0] );
00281         v[QNVEC - 1].swap( v[0] );
00282 
00283         func.update( pd, obj, v[QNVEC], mHess, err );MSQ_ERRRTN( err );
00284         // norm_r = length_squared( &(v[QNVEC][0]), nn );
00285         // norm_g = sqrt(norm_r);
00286 
00287         // checks stopping criterion
00288         term.accumulate_patch( pd, err );MSQ_ERRRTN( err );
00289         term.accumulate_inner( pd, objn, &v[QNVEC][0], err );MSQ_ERRRTN( err );
00290     }
00291 }
00292 
00293 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines