MOAB: Mesh Oriented datABase  (version 5.3.0)
MsqHessian.hpp
Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines