MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 /* ***************************************************************** 00002 MESQUITE -- The Mesh Quality Improvement Toolkit 00003 00004 Copyright 2004 Sandia Corporation and Argonne National 00005 Laboratory. Under the terms of Contract DE-AC04-94AL85000 00006 with Sandia Corporation, the U.S. Government retains certain 00007 rights in 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 diachin2@llnl.gov, djmelan@sandia.gov, mbrewer@sandia.gov, 00024 pknupp@sandia.gov, tleurent@mcs.anl.gov, tmunson@mcs.anl.gov 00025 00026 ***************************************************************** */ 00027 // -*- Mode : c++; tab-width: 3; c-tab-always-indent: t; indent-tabs-mode: nil; c-basic-offset: 3 00028 // -*- 00029 // 00030 // AUTHOR: Todd Munson <tmunson@mcs.anl.gov> 00031 // ORG: Argonne National Laboratory 00032 // E-MAIL: tmunson@mcs.anl.gov 00033 // 00034 // ORIG-DATE: 2-Jan-03 at 11:02:19 bu Thomas Leurent 00035 // LAST-MOD: 5-Oct-04 by Jason Kraftcheck 00036 // 00037 // DESCRIPTION: 00038 // ============ 00039 /*! \file MsqHessian.hpp 00040 00041 The MsqHessian class stores a sparse hessian for a given objective 00042 function. The objective function must be C2 and such that its hessian 00043 has non-zero entries only for the duplet of derivatives corresponding 00044 to nodes of a same element. 00045 00046 \author Thomas Leurent 00047 */ 00048 00049 #ifndef MsqHessian_hpp 00050 #define MsqHessian_hpp 00051 00052 #include "Mesquite.hpp" 00053 #include "Matrix3D.hpp" 00054 #include "PatchData.hpp" 00055 #include "MsqTimer.hpp" 00056 00057 #include <iosfwd> 00058 00059 namespace MBMesquite 00060 { 00061 class ObjectiveFunction; 00062 00063 /*! 00064 \class MsqHessian 00065 \brief Vector3D is the object that effeciently stores the objective function 00066 Hessian each entry is a Matrix3D object (i.e. a vertex Hessian). 00067 */ 00068 class MESQUITE_EXPORT MsqHessian 00069 { 00070 protected: // data accessed directly in tests. 00071 Matrix3D* mEntries; //!< CSR block entries. size: nb of nonzero blocks, i.e. mRowStart[mSize] . 00072 size_t* mRowStart; //!< start of each row in mEntries. size: nb of vertices (mSize). 00073 size_t* mColIndex; //!< CSR block structure: column indexes of the row entries. 00074 00075 size_t mSize; //!< number of rows (or number of columns, this is a square matrix). 00076 00077 Matrix3D* mPreconditioner; 00078 size_t precondArraySize; 00079 00080 Vector3D* mR; //!< array used in the CG solver 00081 Vector3D* mZ; //!< array used in the CG solver 00082 Vector3D* mP; //!< array used in the CG solver 00083 Vector3D* mW; //!< array used in the CG solver 00084 size_t cgArraySizes; //!< size of arrays allocated in the CG solver. 00085 size_t maxCGiter; //!< max nb of iterations of the CG solver. 00086 00087 public: 00088 MsqHessian(); 00089 ~MsqHessian(); 00090 00091 void initialize( PatchData& pd, MsqError& err ); 00092 void initialize( const MsqHessian& other ); 00093 00094 inline void zero_out(); 00095 00096 size_t size() const 00097 { 00098 return mSize; 00099 } 00100 00101 //! returns the diagonal blocks, memory must be allocated before call. 00102 void get_diagonal_blocks( std::vector< Matrix3D >& diag, MsqError& err ) const; 00103 00104 Matrix3D* get_block( size_t i, size_t j ); 00105 const Matrix3D* get_block( size_t i, size_t j ) const; 00106 00107 // inline void accumulate_entries(PatchData &pd, const size_t &elem_index, 00108 // Matrix3D* const &mat3d_array, MsqError &err); 00109 00110 void compute_preconditioner( MsqError& err ); 00111 void apply_preconditioner( Vector3D zloc[], Vector3D rloc[], MsqError& err ); 00112 00113 void cg_solver( Vector3D x[], Vector3D b[], MsqError& err ); 00114 //! Hessian - vector product, summed with a second vector (optional). 00115 friend void axpy( Vector3D res[], size_t size_r, const MsqHessian& H, const Vector3D x[], size_t size_x, 00116 const Vector3D y[], size_t size_y, MsqError& err ); 00117 00118 //! r = this * x, where r and x are arrays of length size(). 00119 void product( Vector3D* r, const Vector3D* x ) const; 00120 00121 friend std::ostream& operator<<( std::ostream&, const MsqHessian& ); 00122 00123 inline void add( size_t row, size_t col, const Matrix3D& m, MsqError& err ); 00124 00125 inline void scale( double value ); 00126 00127 void add( const MsqHessian& other ); 00128 00129 // Release all storage. Object is invalid until next call 00130 // to initialize(..). 00131 void clear(); 00132 00133 inline double norm() const; //!< Frobenius norm 00134 00135 private: 00136 MsqHessian& operator=( const MsqHessian& h ); 00137 MsqHessian( const MsqHessian& copy ); 00138 }; 00139 00140 inline double MsqHessian::norm() const 00141 { 00142 double sum = 0.0; 00143 for( size_t i = 0; i < mRowStart[mSize]; ++i ) 00144 sum += Frobenius_2( mEntries[i] ); 00145 return sqrt( sum ); 00146 } 00147 00148 /*! Sets all Hessian entries to zero. This is usually used before 00149 starting to accumulate elements hessian in the objective function 00150 hessian. */ 00151 inline void MsqHessian::zero_out() 00152 { 00153 if( mSize == 0 ) return; // empty hessian. 00154 00155 size_t i; 00156 for( i = 0; i < mRowStart[mSize]; ++i ) 00157 { 00158 mEntries[i].zero(); 00159 } 00160 } 00161 00162 inline void MsqHessian::scale( double value ) 00163 { 00164 for( size_t i = 0; i < mRowStart[mSize]; ++i ) 00165 mEntries[i] *= value; 00166 } 00167 00168 /*! Accumulates entries of an element hessian into an objective function 00169 hessian. Make sure to use zero_out() before starting the accumulation 00170 process. 00171 00172 \param pd: PatchData in that contains the element which Hessian 00173 we are accumulating in the Hessian matrix. This must be the same 00174 PatchData that was used in MsqHessian::initialize(). 00175 \param elem_index: index of the element in the PatchData. 00176 \param mat3d_array: This is the upper triangular part of the element Hessian 00177 for all nodes, including fixed nodes, for which the entries must be null Matrix3Ds. 00178 \param nb_mat3d. The size of the mat3d_array: (n+1)n/2, where n is 00179 the number of nodes in the element. 00180 */ 00181 // inline void MsqHessian::accumulate_entries(PatchData &pd, const size_t &elem_index, 00182 // Matrix3D* const &mat3d_array, MsqError &err) 00183 // { 00184 // if (&pd != origin_pd) { 00185 // MSQ_SETERR(err)( 00186 // "Cannot accumulate elements from a different patch. " 00187 // "Use MsqHessian::initialize first.", 00188 // MsqError::INVALID_ARG ); 00189 // return; 00190 // } 00191 // 00192 // size_t nve = pd.get_element_array(err)[elem_index].vertex_count(); 00193 // size_t* ve = pd.get_element_array(err)[elem_index].get_vertex_index_array(); 00194 // size_t e = mAccumElemStart[elem_index]; 00195 // size_t i, j, c = 0; 00196 // for (i = 0; i < nve; ++i) 00197 // { 00198 // for (j = i; j < nve; ++j) 00199 // { 00200 // if (ve[i] < mSize && ve[j] < mSize) 00201 // { 00202 // int k = mAccumulation[e++]; 00203 // if (k >= 0) 00204 // mEntries[k] += mat3d_array[c]; 00205 // else 00206 // mEntries[-k].plus_transpose_equal( mat3d_array[c] ); 00207 // } 00208 // ++c; 00209 // } 00210 // } 00211 // } 00212 00213 inline void MsqHessian::add( size_t row, size_t col, const Matrix3D& m, MsqError& err ) 00214 { 00215 if( row <= col ) 00216 { 00217 for( size_t i = mRowStart[row]; i != mRowStart[row + 1]; ++i ) 00218 if( mColIndex[i] == col ) 00219 { 00220 mEntries[i] += m; 00221 return; 00222 } 00223 } 00224 else 00225 { 00226 for( size_t i = mRowStart[col]; i != mRowStart[col + 1]; ++i ) 00227 if( mColIndex[i] == row ) 00228 { 00229 mEntries[i].plus_transpose_equal( m ); 00230 return; 00231 } 00232 } 00233 00234 MSQ_SETERR( err ) 00235 ( MsqError::INVALID_ARG, "Hessian entry (%lu,%lu) does not exist.", (unsigned long)row, (unsigned long)col ); 00236 } 00237 00238 /*! 00239 \param res: array of Vector3D in which the result is stored. 00240 \param size_r: size of the res array. 00241 \param x: vector multiplied by the Hessian. 00242 \param size_x: size of the x array. 00243 \param y: vector added to the Hessian vector product. Set to 0 (NULL) if not needed. 00244 \param size_y: size of the y array. Set to 0 if not needed. 00245 */ 00246 inline void axpy( Vector3D res[], size_t size_r, const MsqHessian& H, const Vector3D x[], size_t size_x, 00247 const Vector3D y[], size_t size_y, MsqError& /*err*/ ) 00248 { 00249 if( ( size_r != H.mSize ) || ( size_x != H.mSize ) || ( size_y != H.mSize && size_y != 0 ) ) 00250 { 00251 // throw an error 00252 } 00253 00254 Vector3D tmpx, tmpm; // for cache opt. 00255 size_t* col = H.mColIndex; 00256 const size_t nn = H.mSize; 00257 size_t rl; // row length 00258 size_t el; // entries index 00259 size_t lo; 00260 size_t c; // column index 00261 size_t i, j; 00262 00263 if( y != 0 ) 00264 { 00265 for( i = 0; i < nn; ++i ) 00266 { 00267 res[i] = y[i]; 00268 } 00269 } 00270 else 00271 { // y == 0 00272 for( i = 0; i < nn; ++i ) 00273 { 00274 res[i] = 0.; 00275 } 00276 } 00277 00278 el = 0; 00279 for( i = 0; i < nn; ++i ) 00280 { 00281 rl = H.mRowStart[i + 1] - H.mRowStart[i]; 00282 lo = *col++; 00283 00284 // Diagonal entry 00285 tmpx = x[i]; 00286 eqAx( tmpm, H.mEntries[el], tmpx ); 00287 ++el; 00288 00289 // Non-diagonal entries 00290 for( j = 1; j < rl; ++j ) 00291 { 00292 c = *col++; 00293 // res[i] += H.mEntries[e] * x[c]; 00294 plusEqAx( tmpm, H.mEntries[el], x[c] ); 00295 // res[c] += transpose(H.mEntries[e]) * tmpxi; 00296 plusEqTransAx( res[c], H.mEntries[el], tmpx ); 00297 ++el; 00298 } 00299 res[lo] += tmpm; 00300 } 00301 } 00302 00303 /*! Computes \f$ z=M^{-1}r \f$ . */ 00304 inline void MsqHessian::apply_preconditioner( Vector3D zloc[], Vector3D rloc[], MsqError& /*err*/ ) 00305 { 00306 size_t m; 00307 00308 for( m = 0; m < mSize; ++m ) 00309 { 00310 #ifdef DIAGONAL_PRECONDITIONER 00311 // preconditioner is identity matrix for now. 00312 zloc[m][0] = mPreconditioner[m][0][0] * rloc[m][0]; 00313 zloc[m][1] = mPreconditioner[m][1][1] * rloc[m][1]; 00314 zloc[m][2] = mPreconditioner[m][2][2] * rloc[m][2]; 00315 #else 00316 // z = inv(L^T) * r 00317 zloc[m][0] = rloc[m][0]; 00318 zloc[m][1] = rloc[m][1] - mPreconditioner[m][0][1] * zloc[m][0]; 00319 zloc[m][2] = rloc[m][2] - mPreconditioner[m][0][2] * zloc[m][0] - mPreconditioner[m][1][2] * zloc[m][1]; 00320 00321 // z = inv(D) * z 00322 zloc[m][0] *= mPreconditioner[m][0][0]; 00323 zloc[m][1] *= mPreconditioner[m][1][1]; 00324 zloc[m][2] *= mPreconditioner[m][2][2]; 00325 00326 // z = inv(L) * z 00327 zloc[m][2] = zloc[m][2]; 00328 zloc[m][1] = zloc[m][1] - mPreconditioner[m][1][2] * zloc[m][2]; 00329 zloc[m][0] = zloc[m][0] - mPreconditioner[m][0][1] * zloc[m][1] - mPreconditioner[m][0][2] * zloc[m][2]; 00330 #endif 00331 } 00332 } 00333 00334 /*! Returns a pointer to the Matrix3D block at position i,j if it exist. 00335 Returns the NULL pointer if position i,j (0-based) is a NULL entry. 00336 Note that block i,j must be in the upper triangular part of the 00337 (symetric) hessian. */ 00338 inline Matrix3D* MsqHessian::get_block( size_t i, size_t j ) 00339 { 00340 size_t c; 00341 00342 if( i >= mSize || j >= mSize || j < i ) return NULL; 00343 00344 for( c = mRowStart[i]; c < mRowStart[i + 1]; ++c ) 00345 { 00346 if( mColIndex[c] == j ) return ( mEntries + c ); 00347 } 00348 00349 // if there is no block at position i,j (zero entry). 00350 return NULL; 00351 } 00352 inline const Matrix3D* MsqHessian::get_block( size_t i, size_t j ) const 00353 { 00354 return const_cast< MsqHessian* >( this )->get_block( i, j ); 00355 } 00356 00357 /* ------------------ I/O ----------------- */ 00358 00359 std::ostream& operator<<( std::ostream& s, const MsqHessian& h ); 00360 00361 } // namespace MBMesquite 00362 00363 #endif // MsqHessian_hpp