Branch data Line data Source code
1 : : /* *****************************************************************
2 : : MESQUITE -- The Mesh Quality Improvement Toolkit
3 : :
4 : : Copyright 2004 Sandia Corporation and Argonne National
5 : : Laboratory. Under the terms of Contract DE-AC04-94AL85000
6 : : with Sandia Corporation, the U.S. Government retains certain
7 : : rights in this software.
8 : :
9 : : This library is free software; you can redistribute it and/or
10 : : modify it under the terms of the GNU Lesser General Public
11 : : License as published by the Free Software Foundation; either
12 : : version 2.1 of the License, or (at your option) any later version.
13 : :
14 : : This library is distributed in the hope that it will be useful,
15 : : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 : : Lesser General Public License for more details.
18 : :
19 : : You should have received a copy of the GNU Lesser General Public License
20 : : (lgpl.txt) along with this library; if not, write to the Free Software
21 : : Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22 : :
23 : : [email protected], [email protected], [email protected],
24 : : [email protected], [email protected], [email protected]
25 : :
26 : : ***************************************************************** */
27 : : // -*- Mode : c++; tab-width: 3; c-tab-always-indent: t; indent-tabs-mode: nil; c-basic-offset: 3
28 : : // -*-
29 : : //
30 : : // AUTHOR: Todd Munson <[email protected]>
31 : : // ORG: Argonne National Laboratory
32 : : // E-MAIL: [email protected]
33 : : //
34 : : // ORIG-DATE: 2-Jan-03 at 11:02:19 bu Thomas Leurent
35 : : // LAST-MOD: 5-Oct-04 by Jason Kraftcheck
36 : : //
37 : : // DESCRIPTION:
38 : : // ============
39 : : /*! \file MsqHessian.hpp
40 : :
41 : : The MsqHessian class stores a sparse hessian for a given objective
42 : : function. The objective function must be C2 and such that its hessian
43 : : has non-zero entries only for the duplet of derivatives corresponding
44 : : to nodes of a same element.
45 : :
46 : : \author Thomas Leurent
47 : : */
48 : :
49 : : #ifndef MsqHessian_hpp
50 : : #define MsqHessian_hpp
51 : :
52 : : #include "Mesquite.hpp"
53 : : #include "Matrix3D.hpp"
54 : : #include "PatchData.hpp"
55 : : #include "MsqTimer.hpp"
56 : :
57 : : #include <iosfwd>
58 : :
59 : : namespace MBMesquite
60 : : {
61 : : class ObjectiveFunction;
62 : :
63 : : /*!
64 : : \class MsqHessian
65 : : \brief Vector3D is the object that effeciently stores the objective function
66 : : Hessian each entry is a Matrix3D object (i.e. a vertex Hessian).
67 : : */
68 : : class MESQUITE_EXPORT MsqHessian
69 : : {
70 : : protected: // data accessed directly in tests.
71 : : Matrix3D* mEntries; //!< CSR block entries. size: nb of nonzero blocks, i.e. mRowStart[mSize] .
72 : : size_t* mRowStart; //!< start of each row in mEntries. size: nb of vertices (mSize).
73 : : size_t* mColIndex; //!< CSR block structure: column indexes of the row entries.
74 : :
75 : : size_t mSize; //!< number of rows (or number of columns, this is a square matrix).
76 : :
77 : : Matrix3D* mPreconditioner;
78 : : size_t precondArraySize;
79 : :
80 : : Vector3D* mR; //!< array used in the CG solver
81 : : Vector3D* mZ; //!< array used in the CG solver
82 : : Vector3D* mP; //!< array used in the CG solver
83 : : Vector3D* mW; //!< array used in the CG solver
84 : : size_t cgArraySizes; //!< size of arrays allocated in the CG solver.
85 : : size_t maxCGiter; //!< max nb of iterations of the CG solver.
86 : :
87 : : public:
88 : : MsqHessian();
89 : : ~MsqHessian();
90 : :
91 : : void initialize( PatchData& pd, MsqError& err );
92 : : void initialize( const MsqHessian& other );
93 : :
94 : : inline void zero_out();
95 : :
96 : 2249970 : size_t size() const
97 : : {
98 : 2249970 : return mSize;
99 : : }
100 : :
101 : : //! returns the diagonal blocks, memory must be allocated before call.
102 : : void get_diagonal_blocks( std::vector< Matrix3D >& diag, MsqError& err ) const;
103 : :
104 : : Matrix3D* get_block( size_t i, size_t j );
105 : : const Matrix3D* get_block( size_t i, size_t j ) const;
106 : :
107 : : // inline void accumulate_entries(PatchData &pd, const size_t &elem_index,
108 : : // Matrix3D* const &mat3d_array, MsqError &err);
109 : :
110 : : void compute_preconditioner( MsqError& err );
111 : : void apply_preconditioner( Vector3D zloc[], Vector3D rloc[], MsqError& err );
112 : :
113 : : void cg_solver( Vector3D x[], Vector3D b[], MsqError& err );
114 : : //! Hessian - vector product, summed with a second vector (optional).
115 : : friend void axpy( Vector3D res[], size_t size_r, const MsqHessian& H, const Vector3D x[], size_t size_x,
116 : : const Vector3D y[], size_t size_y, MsqError& err );
117 : :
118 : : //! r = this * x, where r and x are arrays of length size().
119 : : void product( Vector3D* r, const Vector3D* x ) const;
120 : :
121 : : friend std::ostream& operator<<( std::ostream&, const MsqHessian& );
122 : :
123 : : inline void add( size_t row, size_t col, const Matrix3D& m, MsqError& err );
124 : :
125 : : inline void scale( double value );
126 : :
127 : : void add( const MsqHessian& other );
128 : :
129 : : // Release all storage. Object is invalid until next call
130 : : // to initialize(..).
131 : : void clear();
132 : :
133 : : inline double norm() const; //!< Frobenius norm
134 : :
135 : : private:
136 : : MsqHessian& operator=( const MsqHessian& h );
137 : : MsqHessian( const MsqHessian& copy );
138 : : };
139 : :
140 : : inline double MsqHessian::norm() const
141 : : {
142 : : double sum = 0.0;
143 : : for( size_t i = 0; i < mRowStart[mSize]; ++i )
144 : : sum += Frobenius_2( mEntries[i] );
145 : : return sqrt( sum );
146 : : }
147 : :
148 : : /*! Sets all Hessian entries to zero. This is usually used before
149 : : starting to accumulate elements hessian in the objective function
150 : : hessian. */
151 : 366 : inline void MsqHessian::zero_out()
152 : : {
153 [ - + ]: 366 : if( mSize == 0 ) return; // empty hessian.
154 : :
155 : : size_t i;
156 [ + + ]: 945836 : for( i = 0; i < mRowStart[mSize]; ++i )
157 : : {
158 : 945470 : mEntries[i].zero();
159 : : }
160 : : }
161 : :
162 : 270 : inline void MsqHessian::scale( double value )
163 : : {
164 [ + + ]: 767743 : for( size_t i = 0; i < mRowStart[mSize]; ++i )
165 : 767473 : mEntries[i] *= value;
166 : 270 : }
167 : :
168 : : /*! Accumulates entries of an element hessian into an objective function
169 : : hessian. Make sure to use zero_out() before starting the accumulation
170 : : process.
171 : :
172 : : \param pd: PatchData in that contains the element which Hessian
173 : : we are accumulating in the Hessian matrix. This must be the same
174 : : PatchData that was used in MsqHessian::initialize().
175 : : \param elem_index: index of the element in the PatchData.
176 : : \param mat3d_array: This is the upper triangular part of the element Hessian
177 : : for all nodes, including fixed nodes, for which the entries must be null Matrix3Ds.
178 : : \param nb_mat3d. The size of the mat3d_array: (n+1)n/2, where n is
179 : : the number of nodes in the element.
180 : : */
181 : : // inline void MsqHessian::accumulate_entries(PatchData &pd, const size_t &elem_index,
182 : : // Matrix3D* const &mat3d_array, MsqError &err)
183 : : // {
184 : : // if (&pd != origin_pd) {
185 : : // MSQ_SETERR(err)(
186 : : // "Cannot accumulate elements from a different patch. "
187 : : // "Use MsqHessian::initialize first.",
188 : : // MsqError::INVALID_ARG );
189 : : // return;
190 : : // }
191 : : //
192 : : // size_t nve = pd.get_element_array(err)[elem_index].vertex_count();
193 : : // size_t* ve = pd.get_element_array(err)[elem_index].get_vertex_index_array();
194 : : // size_t e = mAccumElemStart[elem_index];
195 : : // size_t i, j, c = 0;
196 : : // for (i = 0; i < nve; ++i)
197 : : // {
198 : : // for (j = i; j < nve; ++j)
199 : : // {
200 : : // if (ve[i] < mSize && ve[j] < mSize)
201 : : // {
202 : : // int k = mAccumulation[e++];
203 : : // if (k >= 0)
204 : : // mEntries[k] += mat3d_array[c];
205 : : // else
206 : : // mEntries[-k].plus_transpose_equal( mat3d_array[c] );
207 : : // }
208 : : // ++c;
209 : : // }
210 : : // }
211 : : // }
212 : :
213 : 3262527 : inline void MsqHessian::add( size_t row, size_t col, const Matrix3D& m, MsqError& err )
214 : : {
215 [ + + ]: 3262527 : if( row <= col )
216 : : {
217 [ + - ]: 8408683 : for( size_t i = mRowStart[row]; i != mRowStart[row + 1]; ++i )
218 [ + + ]: 8408683 : if( mColIndex[i] == col )
219 : : {
220 : 2945207 : mEntries[i] += m;
221 : 2945207 : return;
222 : : }
223 : : }
224 : : else
225 : : {
226 [ + - ]: 1437153 : for( size_t i = mRowStart[col]; i != mRowStart[col + 1]; ++i )
227 [ + + ]: 1437153 : if( mColIndex[i] == row )
228 : : {
229 : 317320 : mEntries[i].plus_transpose_equal( m );
230 : 317320 : return;
231 : : }
232 : : }
233 : :
234 : : MSQ_SETERR( err )
235 [ # # ]: 3262527 : ( MsqError::INVALID_ARG, "Hessian entry (%lu,%lu) does not exist.", (unsigned long)row, (unsigned long)col );
236 : : }
237 : :
238 : : /*!
239 : : \param res: array of Vector3D in which the result is stored.
240 : : \param size_r: size of the res array.
241 : : \param x: vector multiplied by the Hessian.
242 : : \param size_x: size of the x array.
243 : : \param y: vector added to the Hessian vector product. Set to 0 (NULL) if not needed.
244 : : \param size_y: size of the y array. Set to 0 if not needed.
245 : : */
246 : 2524 : inline void axpy( Vector3D res[], size_t size_r, const MsqHessian& H, const Vector3D x[], size_t size_x,
247 : : const Vector3D y[], size_t size_y, MsqError& /*err*/ )
248 : : {
249 [ + - ][ + - ]: 2524 : if( ( size_r != H.mSize ) || ( size_x != H.mSize ) || ( size_y != H.mSize && size_y != 0 ) )
250 : : {
251 : : // throw an error
252 : : }
253 : :
254 [ + - ][ + - ]: 2524 : Vector3D tmpx, tmpm; // for cache opt.
255 : 2524 : size_t* col = H.mColIndex;
256 : 2524 : const size_t nn = H.mSize;
257 : : size_t rl; // row length
258 : : size_t el; // entries index
259 : : size_t lo;
260 : : size_t c; // column index
261 : : size_t i, j;
262 : :
263 [ - + ]: 2524 : if( y != 0 )
264 : : {
265 [ # # ]: 0 : for( i = 0; i < nn; ++i )
266 : : {
267 [ # # ]: 0 : res[i] = y[i];
268 : : }
269 : : }
270 : : else
271 : : { // y == 0
272 [ + + ]: 244727 : for( i = 0; i < nn; ++i )
273 : : {
274 [ + - ]: 242203 : res[i] = 0.;
275 : : }
276 : : }
277 : :
278 : 2524 : el = 0;
279 [ + + ]: 244727 : for( i = 0; i < nn; ++i )
280 : : {
281 : 242203 : rl = H.mRowStart[i + 1] - H.mRowStart[i];
282 : 242203 : lo = *col++;
283 : :
284 : : // Diagonal entry
285 [ + - ]: 242203 : tmpx = x[i];
286 [ + - ]: 242203 : eqAx( tmpm, H.mEntries[el], tmpx );
287 : 242203 : ++el;
288 : :
289 : : // Non-diagonal entries
290 [ + + ]: 1651323 : for( j = 1; j < rl; ++j )
291 : : {
292 : 1409120 : c = *col++;
293 : : // res[i] += H.mEntries[e] * x[c];
294 [ + - ]: 1409120 : plusEqAx( tmpm, H.mEntries[el], x[c] );
295 : : // res[c] += transpose(H.mEntries[e]) * tmpxi;
296 [ + - ]: 1409120 : plusEqTransAx( res[c], H.mEntries[el], tmpx );
297 : 1409120 : ++el;
298 : : }
299 [ + - ]: 242203 : res[lo] += tmpm;
300 : : }
301 : 2524 : }
302 : :
303 : : /*! Computes \f$ z=M^{-1}r \f$ . */
304 : 2710 : inline void MsqHessian::apply_preconditioner( Vector3D zloc[], Vector3D rloc[], MsqError& /*err*/ )
305 : : {
306 : : size_t m;
307 : :
308 [ + + ]: 270459 : for( m = 0; m < mSize; ++m )
309 : : {
310 : : #ifdef DIAGONAL_PRECONDITIONER
311 : : // preconditioner is identity matrix for now.
312 : : zloc[m][0] = mPreconditioner[m][0][0] * rloc[m][0];
313 : : zloc[m][1] = mPreconditioner[m][1][1] * rloc[m][1];
314 : : zloc[m][2] = mPreconditioner[m][2][2] * rloc[m][2];
315 : : #else
316 : : // z = inv(L^T) * r
317 : 267749 : zloc[m][0] = rloc[m][0];
318 : 267749 : zloc[m][1] = rloc[m][1] - mPreconditioner[m][0][1] * zloc[m][0];
319 : 267749 : zloc[m][2] = rloc[m][2] - mPreconditioner[m][0][2] * zloc[m][0] - mPreconditioner[m][1][2] * zloc[m][1];
320 : :
321 : : // z = inv(D) * z
322 : 267749 : zloc[m][0] *= mPreconditioner[m][0][0];
323 : 267749 : zloc[m][1] *= mPreconditioner[m][1][1];
324 : 267749 : zloc[m][2] *= mPreconditioner[m][2][2];
325 : :
326 : : // z = inv(L) * z
327 : 267749 : zloc[m][2] = zloc[m][2];
328 : 267749 : zloc[m][1] = zloc[m][1] - mPreconditioner[m][1][2] * zloc[m][2];
329 : 267749 : zloc[m][0] = zloc[m][0] - mPreconditioner[m][0][1] * zloc[m][1] - mPreconditioner[m][0][2] * zloc[m][2];
330 : : #endif
331 : : }
332 : 2710 : }
333 : :
334 : : /*! Returns a pointer to the Matrix3D block at position i,j if it exist.
335 : : Returns the NULL pointer if position i,j (0-based) is a NULL entry.
336 : : Note that block i,j must be in the upper triangular part of the
337 : : (symetric) hessian. */
338 : 140342 : inline Matrix3D* MsqHessian::get_block( size_t i, size_t j )
339 : : {
340 : : size_t c;
341 : :
342 [ + - ][ + - ]: 140342 : if( i >= mSize || j >= mSize || j < i ) return NULL;
[ - + ]
343 : :
344 [ + - ]: 140342 : for( c = mRowStart[i]; c < mRowStart[i + 1]; ++c )
345 : : {
346 [ + - ]: 140342 : if( mColIndex[c] == j ) return ( mEntries + c );
347 : : }
348 : :
349 : : // if there is no block at position i,j (zero entry).
350 : 0 : return NULL;
351 : : }
352 : : inline const Matrix3D* MsqHessian::get_block( size_t i, size_t j ) const
353 : : {
354 : : return const_cast< MsqHessian* >( this )->get_block( i, j );
355 : : }
356 : :
357 : : /* ------------------ I/O ----------------- */
358 : :
359 : : std::ostream& operator<<( std::ostream& s, const MsqHessian& h );
360 : :
361 : : } // namespace MBMesquite
362 : :
363 : : #endif // MsqHessian_hpp
|
