MOAB: Mesh Oriented datABase  (version 5.4.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[],
00116                       size_t size_r,
00117                       const MsqHessian& H,
00118                       const Vector3D x[],
00119                       size_t size_x,
00120                       const Vector3D y[],
00121                       size_t size_y,
00122                       MsqError& err );
00123 
00124     //! r = this * x, where r and x are arrays of length size().
00125     void product( Vector3D* r, const Vector3D* x ) const;
00126 
00127     friend std::ostream& operator<<( std::ostream&, const MsqHessian& );
00128 
00129     inline void add( size_t row, size_t col, const Matrix3D& m, MsqError& err );
00130 
00131     inline void scale( double value );
00132 
00133     void add( const MsqHessian& other );
00134 
00135     // Release all storage.  Object is invalid until next call
00136     // to initialize(..).
00137     void clear();
00138 
00139     inline double norm() const;  //!< Frobenius norm
00140 
00141   private:
00142     MsqHessian& operator=( const MsqHessian& h );
00143     MsqHessian( const MsqHessian& copy );
00144 };
00145 
00146 inline double MsqHessian::norm() const
00147 {
00148     double sum = 0.0;
00149     for( size_t i = 0; i < mRowStart[mSize]; ++i )
00150         sum += Frobenius_2( mEntries[i] );
00151     return sqrt( sum );
00152 }
00153 
00154 /*! Sets all Hessian entries to zero. This is usually used before
00155   starting to accumulate elements hessian in the objective function
00156   hessian. */
00157 inline void MsqHessian::zero_out()
00158 {
00159     if( mSize == 0 ) return;  // empty hessian.
00160 
00161     size_t i;
00162     for( i = 0; i < mRowStart[mSize]; ++i )
00163     {
00164         mEntries[i].zero();
00165     }
00166 }
00167 
00168 inline void MsqHessian::scale( double value )
00169 {
00170     for( size_t i = 0; i < mRowStart[mSize]; ++i )
00171         mEntries[i] *= value;
00172 }
00173 
00174 /*! Accumulates entries of an element hessian into an objective function
00175   hessian. Make sure to use zero_out() before starting the accumulation
00176   process.
00177 
00178   \param pd: PatchData in that contains the element which Hessian
00179   we are accumulating in the Hessian matrix. This must be the same
00180   PatchData that was used in MsqHessian::initialize().
00181   \param elem_index: index of the element in the PatchData.
00182   \param mat3d_array: This is the upper triangular part of the element Hessian
00183   for all nodes, including fixed nodes, for which the entries must be null Matrix3Ds.
00184   \param nb_mat3d. The size of the mat3d_array: (n+1)n/2, where n is
00185   the number of nodes in the element.
00186 */
00187 //  inline void MsqHessian::accumulate_entries(PatchData &pd, const size_t &elem_index,
00188 //                                             Matrix3D* const &mat3d_array, MsqError &err)
00189 //    {
00190 //      if (&pd != origin_pd) {
00191 //        MSQ_SETERR(err)(
00192 //                    "Cannot accumulate elements from a different patch. "
00193 //                    "Use MsqHessian::initialize first.",
00194 //                    MsqError::INVALID_ARG );
00195 //        return;
00196 //      }
00197 //
00198 //      size_t nve = pd.get_element_array(err)[elem_index].vertex_count();
00199 //      size_t* ve = pd.get_element_array(err)[elem_index].get_vertex_index_array();
00200 //      size_t e = mAccumElemStart[elem_index];
00201 //      size_t i, j, c = 0;
00202 //      for (i = 0; i < nve; ++i)
00203 //      {
00204 //        for (j = i; j < nve; ++j)
00205 //        {
00206 //          if (ve[i] < mSize && ve[j] < mSize)
00207 //          {
00208 //            int k = mAccumulation[e++];
00209 //            if (k >= 0)
00210 //              mEntries[k] += mat3d_array[c];
00211 //            else
00212 //              mEntries[-k].plus_transpose_equal( mat3d_array[c] );
00213 //          }
00214 //          ++c;
00215 //        }
00216 //      }
00217 //    }
00218 
00219 inline void MsqHessian::add( size_t row, size_t col, const Matrix3D& m, MsqError& err )
00220 {
00221     if( row <= col )
00222     {
00223         for( size_t i = mRowStart[row]; i != mRowStart[row + 1]; ++i )
00224             if( mColIndex[i] == col )
00225             {
00226                 mEntries[i] += m;
00227                 return;
00228             }
00229     }
00230     else
00231     {
00232         for( size_t i = mRowStart[col]; i != mRowStart[col + 1]; ++i )
00233             if( mColIndex[i] == row )
00234             {
00235                 mEntries[i].plus_transpose_equal( m );
00236                 return;
00237             }
00238     }
00239 
00240     MSQ_SETERR( err )
00241     ( MsqError::INVALID_ARG, "Hessian entry (%lu,%lu) does not exist.", (unsigned long)row, (unsigned long)col );
00242 }
00243 
00244 /*!
00245   \param res: array of Vector3D in which the result is stored.
00246   \param size_r: size of the res array.
00247   \param x: vector multiplied by the Hessian.
00248   \param size_x: size of the x array.
00249   \param y: vector added to the Hessian vector product. Set to 0 (NULL) if not needed.
00250   \param size_y: size of the y array. Set to 0 if not needed.
00251 */
00252 inline void axpy( Vector3D res[],
00253                   size_t size_r,
00254                   const MsqHessian& H,
00255                   const Vector3D x[],
00256                   size_t size_x,
00257                   const Vector3D y[],
00258                   size_t size_y,
00259                   MsqError& /*err*/ )
00260 {
00261     if( ( size_r != H.mSize ) || ( size_x != H.mSize ) || ( size_y != H.mSize && size_y != 0 ) )
00262     {
00263         // throw an error
00264     }
00265 
00266     Vector3D tmpx, tmpm;  // for cache opt.
00267     size_t* col     = H.mColIndex;
00268     const size_t nn = H.mSize;
00269     size_t rl;  // row length
00270     size_t el;  // entries index
00271     size_t lo;
00272     size_t c;  // column index
00273     size_t i, j;
00274 
00275     if( y != 0 )
00276     {
00277         for( i = 0; i < nn; ++i )
00278         {
00279             res[i] = y[i];
00280         }
00281     }
00282     else
00283     {  // y == 0
00284         for( i = 0; i < nn; ++i )
00285         {
00286             res[i] = 0.;
00287         }
00288     }
00289 
00290     el = 0;
00291     for( i = 0; i < nn; ++i )
00292     {
00293         rl = H.mRowStart[i + 1] - H.mRowStart[i];
00294         lo = *col++;
00295 
00296         // Diagonal entry
00297         tmpx = x[i];
00298         eqAx( tmpm, H.mEntries[el], tmpx );
00299         ++el;
00300 
00301         // Non-diagonal entries
00302         for( j = 1; j < rl; ++j )
00303         {
00304             c = *col++;
00305             //          res[i] += H.mEntries[e] * x[c];
00306             plusEqAx( tmpm, H.mEntries[el], x[c] );
00307             //          res[c] += transpose(H.mEntries[e]) * tmpxi;
00308             plusEqTransAx( res[c], H.mEntries[el], tmpx );
00309             ++el;
00310         }
00311         res[lo] += tmpm;
00312     }
00313 }
00314 
00315 /*! Computes \f$ z=M^{-1}r \f$ . */
00316 inline void MsqHessian::apply_preconditioner( Vector3D zloc[], Vector3D rloc[], MsqError& /*err*/ )
00317 {
00318     size_t m;
00319 
00320     for( m = 0; m < mSize; ++m )
00321     {
00322 #ifdef DIAGONAL_PRECONDITIONER
00323         // preconditioner is identity matrix for now.
00324         zloc[m][0] = mPreconditioner[m][0][0] * rloc[m][0];
00325         zloc[m][1] = mPreconditioner[m][1][1] * rloc[m][1];
00326         zloc[m][2] = mPreconditioner[m][2][2] * rloc[m][2];
00327 #else
00328         // z = inv(L^T) * r
00329         zloc[m][0] = rloc[m][0];
00330         zloc[m][1] = rloc[m][1] - mPreconditioner[m][0][1] * zloc[m][0];
00331         zloc[m][2] = rloc[m][2] - mPreconditioner[m][0][2] * zloc[m][0] - mPreconditioner[m][1][2] * zloc[m][1];
00332 
00333         // z = inv(D) * z
00334         zloc[m][0] *= mPreconditioner[m][0][0];
00335         zloc[m][1] *= mPreconditioner[m][1][1];
00336         zloc[m][2] *= mPreconditioner[m][2][2];
00337 
00338         // z = inv(L) * z
00339         zloc[m][2] = zloc[m][2];
00340         zloc[m][1] = zloc[m][1] - mPreconditioner[m][1][2] * zloc[m][2];
00341         zloc[m][0] = zloc[m][0] - mPreconditioner[m][0][1] * zloc[m][1] - mPreconditioner[m][0][2] * zloc[m][2];
00342 #endif
00343     }
00344 }
00345 
00346 /*! Returns a pointer to the Matrix3D block at position i,j if it exist.
00347   Returns the NULL pointer if position i,j (0-based) is a NULL entry.
00348   Note that block i,j must be in the upper triangular part of the
00349   (symetric) hessian. */
00350 inline Matrix3D* MsqHessian::get_block( size_t i, size_t j )
00351 {
00352     size_t c;
00353 
00354     if( i >= mSize || j >= mSize || j < i ) return NULL;
00355 
00356     for( c = mRowStart[i]; c < mRowStart[i + 1]; ++c )
00357     {
00358         if( mColIndex[c] == j ) return ( mEntries + c );
00359     }
00360 
00361     // if there is no block at position i,j (zero entry).
00362     return NULL;
00363 }
00364 inline const Matrix3D* MsqHessian::get_block( size_t i, size_t j ) const
00365 {
00366     return const_cast< MsqHessian* >( this )->get_block( i, j );
00367 }
00368 
00369 /* ------------------ I/O ----------------- */
00370 
00371 std::ostream& operator<<( std::ostream& s, const MsqHessian& h );
00372 
00373 }  // namespace MBMesquite
00374 
00375 #endif  // MsqHessian_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines