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 by Thomas Leurent
35 : : // LAST-MOD: 26-Nov-03 at 15:47:42 by Thomas Leurent
36 : : //
37 : : // DESCRIPTION:
38 : : // ============
39 : : /*! \file MsqHessian.cpp
40 : :
41 : : \author Thomas Leurent
42 : :
43 : : */
44 : :
45 : : #include "MsqHessian.hpp"
46 : : #include "MsqTimer.hpp"
47 : :
48 : : #include <cmath>
49 : : #include <iostream>
50 : :
51 : : namespace MBMesquite
52 : : {
53 : :
54 : 64 : MsqHessian::MsqHessian()
55 : : : mEntries( 0 ), mRowStart( 0 ), mColIndex( 0 ), mSize( 0 ), mPreconditioner( 0 ), precondArraySize( 0 ), mR( 0 ),
56 : 64 : mZ( 0 ), mP( 0 ), mW( 0 ), cgArraySizes( 0 ), maxCGiter( 50 )
57 : : {
58 : 64 : }
59 : :
60 : 63 : MsqHessian::~MsqHessian()
61 : : {
62 : 63 : clear();
63 : 63 : }
64 : :
65 : 70 : void MsqHessian::clear()
66 : : {
67 : 70 : mSize = precondArraySize = cgArraySizes = 0;
68 : :
69 [ + + ][ + + ]: 143987 : delete[] mEntries;
70 : 70 : mEntries = 0;
71 [ + + ]: 70 : delete[] mRowStart;
72 : 70 : mRowStart = 0;
73 [ + + ]: 70 : delete[] mColIndex;
74 : 70 : mColIndex = 0;
75 : :
76 [ + + ][ + + ]: 5121 : delete[] mPreconditioner;
77 : 70 : mPreconditioner = 0;
78 : :
79 [ + + ]: 70 : delete[] mR;
80 : 70 : mR = 0;
81 [ + + ]: 70 : delete[] mZ;
82 : 70 : mZ = 0;
83 [ + + ]: 70 : delete[] mP;
84 : 70 : mP = 0;
85 [ + + ]: 70 : delete[] mW;
86 : 70 : mW = 0;
87 : 70 : }
88 : :
89 : : /*! \brief creates a sparse structure for a Hessian, based on the
90 : : connectivity information contained in the PatchData.
91 : : Only the upper triangular part of the Hessian is stored. */
92 : 62 : void MsqHessian::initialize( PatchData& pd, MsqError& err )
93 : : {
94 : : MSQ_FUNCTION_TIMER( "MsqHession::initialize" );
95 [ - + ][ # # ]: 62 : delete[] mEntries;
96 [ - + ]: 62 : delete[] mRowStart;
97 [ - + ]: 62 : delete[] mColIndex;
98 : :
99 : 62 : size_t num_vertices = pd.num_free_vertices();
100 : 62 : size_t num_elements = pd.num_elements();
101 : : size_t const* vtx_list;
102 : : size_t e, r, rs, re, c, cs, ce, nz, nnz, nve, i, j;
103 [ - + ][ # # ]: 62 : MsqMeshEntity* patchElemArray = pd.get_element_array( err );MSQ_CHKERR( err );
104 : :
105 [ - + ]: 62 : if( num_vertices == 0 )
106 : : {
107 [ # # # # ]: 0 : MSQ_SETERR( err )( "No vertices in PatchData", MsqError::INVALID_ARG );
108 : 0 : return;
109 : : }
110 : :
111 : 62 : mSize = num_vertices;
112 : :
113 : : // Calculate the offsets for a CSC representation of the accumulation
114 : : // pattern.
115 : :
116 [ + - ]: 62 : size_t* col_start = new size_t[num_vertices + 1];
117 : : // mAccumElemStart = new size_t[num_elements+1];
118 : : // mAccumElemStart[0] = 0;
119 : :
120 [ + + ]: 27402 : for( i = 0; i < num_vertices; ++i )
121 : : {
122 : 27340 : col_start[i] = 0;
123 : : }
124 : :
125 [ + + ]: 69644 : for( e = 0; e < num_elements; ++e )
126 : : {
127 : 69582 : nve = patchElemArray[e].node_count();
128 : 69582 : vtx_list = patchElemArray[e].get_vertex_index_array();
129 : 69582 : int nfe = 0;
130 : :
131 [ + + ]: 357024 : for( i = 0; i < nve; ++i )
132 : : {
133 : 287442 : r = vtx_list[i];
134 [ + + ]: 287442 : if( r < num_vertices ) ++nfe;
135 : :
136 [ + + ]: 1043228 : for( j = i; j < nve; ++j )
137 : : {
138 : 755786 : c = vtx_list[j];
139 : :
140 [ + + ]: 755786 : if( r <= c )
141 : : {
142 [ + + ]: 558848 : if( c < num_vertices ) ++col_start[c];
143 : : }
144 : : else
145 : : {
146 [ + + ]: 196938 : if( r < num_vertices ) ++col_start[r];
147 : : }
148 : : }
149 : : }
150 : : // mAccumElemStart[e+1] = mAccumElemStart[e] + (nfe+1)*nfe/2;
151 : : }
152 : :
153 : 62 : nz = 0;
154 [ + + ]: 27402 : for( i = 0; i < num_vertices; ++i )
155 : : {
156 : 27340 : j = col_start[i];
157 : 27340 : col_start[i] = nz;
158 : 27340 : nz += j;
159 : : }
160 : 62 : col_start[i] = nz;
161 : :
162 : : // Finished putting matrix into CSC representation
163 : :
164 [ + - ]: 62 : int* row_instr = new int[5 * nz];
165 [ + - ]: 62 : size_t* row_index = new size_t[nz];
166 : :
167 : 62 : nz = 0;
168 [ + + ]: 69644 : for( e = 0; e < num_elements; ++e )
169 : : {
170 : 69582 : nve = patchElemArray[e].node_count();
171 : 69582 : vtx_list = patchElemArray[e].get_vertex_index_array();
172 : :
173 [ + + ]: 357024 : for( i = 0; i < nve; ++i )
174 : : {
175 : 287442 : r = vtx_list[i];
176 : :
177 [ + + ]: 1043228 : for( j = i; j < nve; ++j )
178 : : {
179 : 755786 : c = vtx_list[j];
180 : :
181 [ + + ]: 755786 : if( r <= c )
182 : : {
183 [ + + ]: 558848 : if( c < num_vertices )
184 : : {
185 : 418928 : row_index[col_start[c]] = r;
186 : 418928 : row_instr[col_start[c]] = nz++;
187 : 418928 : ++col_start[c];
188 : : }
189 : : }
190 : : else
191 : : {
192 [ + + ]: 196938 : if( r < num_vertices )
193 : : {
194 : 126782 : row_index[col_start[r]] = c;
195 : : // can't use -nz, but can negate row_instr[col_start[r]]
196 : 126782 : row_instr[col_start[r]] = nz++;
197 : 126782 : row_instr[col_start[r]] = -row_instr[col_start[r]];
198 : 126782 : ++col_start[r];
199 : : }
200 : : }
201 : : }
202 : : }
203 : : }
204 : :
205 [ + + ]: 27340 : for( i = num_vertices - 1; i > 0; --i )
206 : : {
207 : 27278 : col_start[i + 1] = col_start[i];
208 : : }
209 : 62 : col_start[1] = col_start[0];
210 : 62 : col_start[0] = 0;
211 : :
212 : : // cout << "col_start: ";
213 : : // for (int t=0; t<num_vertices+1; ++t)
214 : : // cout << col_start[t] << " ";
215 : : // cout << endl;
216 : : // cout << "row_index: ";
217 : : // for (int t=0; t<nz; ++t)
218 : : // cout << row_index[t] << " ";
219 : : // cout << endl;
220 : : // cout << "row_instr: ";
221 : : // for (int t=0; t<nz; ++t)
222 : : // cout << row_instr[t] << " ";
223 : : // cout << endl;
224 : :
225 : : // Convert CSC to CSR
226 : : // First calculate the offsets in the row
227 : :
228 [ + - ]: 62 : size_t* row_start = new size_t[num_vertices + 1];
229 : :
230 [ + + ]: 27402 : for( i = 0; i < num_vertices; ++i )
231 : : {
232 : 27340 : row_start[i] = 0;
233 : : }
234 : :
235 [ + + ]: 545772 : for( i = 0; i < nz; ++i )
236 : : {
237 : 545710 : ++row_start[row_index[i]];
238 : : }
239 : :
240 : 62 : nz = 0;
241 [ + + ]: 27402 : for( i = 0; i < num_vertices; ++i )
242 : : {
243 : 27340 : j = row_start[i];
244 : 27340 : row_start[i] = nz;
245 : 27340 : nz += j;
246 : : }
247 : 62 : row_start[i] = nz;
248 : :
249 : : // Now calculate the pattern
250 : :
251 [ + - ]: 62 : size_t* col_index = new size_t[nz];
252 [ + - ]: 62 : int* col_instr = new int[nz];
253 : :
254 [ + + ]: 27402 : for( i = 0; i < num_vertices; ++i )
255 : : {
256 : 27340 : cs = col_start[i];
257 : 27340 : ce = col_start[i + 1];
258 : :
259 [ + + ]: 573050 : while( cs < ce )
260 : : {
261 : 545710 : r = row_index[cs];
262 : :
263 : 545710 : col_index[row_start[r]] = i;
264 : 545710 : col_instr[row_start[r]] = row_instr[cs];
265 : :
266 : 545710 : ++row_start[r];
267 : 545710 : ++cs;
268 : : }
269 : : }
270 : :
271 [ + + ]: 27340 : for( i = num_vertices - 1; i > 0; --i )
272 : : {
273 : 27278 : row_start[i + 1] = row_start[i];
274 : : }
275 : 62 : row_start[1] = row_start[0];
276 : 62 : row_start[0] = 0;
277 : :
278 [ + - ]: 62 : delete[] row_index;
279 : :
280 : : // Now that the matrix is CSR
281 : : // Column indices for each row are sorted
282 : :
283 : : // Compaction -- count the number of nonzeros
284 : 62 : mRowStart = col_start; // don't need to reallocate
285 : : // mAccumulation = row_instr; // don't need to reallocate
286 [ + - ]: 62 : delete[] row_instr;
287 : :
288 [ + + ]: 27464 : for( i = 0; i <= num_vertices; ++i )
289 : : {
290 : 27402 : mRowStart[i] = 0;
291 : : }
292 : :
293 : 62 : nnz = 0;
294 [ + + ]: 27402 : for( i = 0; i < num_vertices; ++i )
295 : : {
296 : 27340 : rs = row_start[i];
297 : 27340 : re = row_start[i + 1];
298 : :
299 : 27340 : c = num_vertices;
300 [ + + ]: 573050 : while( rs < re )
301 : : {
302 [ + + ]: 545710 : if( c != col_index[rs] )
303 : : {
304 : : // This is an unseen nonzero
305 : :
306 : 152094 : c = col_index[rs];
307 : 152094 : ++mRowStart[i];
308 : 152094 : ++nnz;
309 : : }
310 : :
311 : : // if (col_instr[rs] >= 0) {
312 : : // mAccumulation[col_instr[rs]] = nnz - 1;
313 : : //}
314 : : // else {
315 : : // mAccumulation[-col_instr[rs]] = 1 - nnz;
316 : : //}
317 : :
318 : 545710 : ++rs;
319 : : }
320 : : }
321 : :
322 : 62 : nnz = 0;
323 [ + + ]: 27402 : for( i = 0; i < num_vertices; ++i )
324 : : {
325 : 27340 : j = mRowStart[i];
326 : 27340 : mRowStart[i] = nnz;
327 : 27340 : nnz += j;
328 : : }
329 : 62 : mRowStart[i] = nnz;
330 : :
331 [ + - ]: 62 : delete[] col_instr;
332 : :
333 : : // Fill in the compacted hessian matrix
334 : :
335 [ + - ]: 62 : mColIndex = new size_t[nnz];
336 : :
337 [ + + ]: 27402 : for( i = 0; i < num_vertices; ++i )
338 : : {
339 : 27340 : rs = row_start[i];
340 : 27340 : re = row_start[i + 1];
341 : :
342 : 27340 : c = num_vertices;
343 [ + + ]: 573050 : while( rs < re )
344 : : {
345 [ + + ]: 545710 : if( c != col_index[rs] )
346 : : {
347 : : // This is an unseen nonzero
348 : :
349 : 152094 : c = col_index[rs];
350 : 152094 : mColIndex[mRowStart[i]] = c;
351 : 152094 : mRowStart[i]++;
352 : : }
353 : 545710 : ++rs;
354 : : }
355 : : }
356 : :
357 [ + + ]: 27340 : for( i = num_vertices - 1; i > 0; --i )
358 : : {
359 : 27278 : mRowStart[i + 1] = mRowStart[i];
360 : : }
361 : 62 : mRowStart[1] = mRowStart[0];
362 : 62 : mRowStart[0] = 0;
363 : :
364 [ + - ]: 62 : delete[] row_start;
365 [ + - ]: 62 : delete[] col_index;
366 : :
367 [ + - ][ + - ]: 152156 : mEntries = new Matrix3D[nnz]; // On Solaris, no initializer allowed for new of an array
[ + + # # ]
368 [ + + ]: 152156 : for( i = 0; i < nnz; ++i )
369 : 152094 : mEntries[i] = 0.; // so we initialize all entries manually.
370 : :
371 : : // origin_pd = &pd;
372 : :
373 : 62 : return;
374 : : }
375 : :
376 : 0 : void MsqHessian::initialize( const MsqHessian& other )
377 : : {
378 [ # # ]: 0 : if( !other.mSize )
379 : : {
380 [ # # ][ # # ]: 0 : delete[] mEntries;
381 [ # # ]: 0 : delete[] mRowStart;
382 [ # # ]: 0 : delete[] mColIndex;
383 : 0 : mEntries = 0;
384 : 0 : mRowStart = 0;
385 : 0 : mColIndex = 0;
386 : 0 : mSize = 0;
387 [ # # ]: 0 : return;
388 : : }
389 : :
390 [ # # ][ # # ]: 0 : if( mSize != other.mSize || mRowStart[mSize] != other.mRowStart[mSize] )
391 : : {
392 [ # # ][ # # ]: 0 : delete[] mEntries;
393 [ # # ]: 0 : delete[] mRowStart;
394 [ # # ]: 0 : delete[] mColIndex;
395 : :
396 : 0 : mSize = other.mSize;
397 : :
398 [ # # ]: 0 : mRowStart = new size_t[mSize + 1];
399 [ # # ][ # # ]: 0 : mEntries = new Matrix3D[other.mRowStart[mSize]];
[ # # # # ]
400 [ # # ]: 0 : mColIndex = new size_t[other.mRowStart[mSize]];
401 : : }
402 : :
403 : 0 : memcpy( mRowStart, other.mRowStart, sizeof( size_t ) * ( mSize + 1 ) );
404 : 0 : memcpy( mColIndex, other.mColIndex, sizeof( size_t ) * mRowStart[mSize] );
405 : : }
406 : :
407 : 0 : void MsqHessian::add( const MsqHessian& other )
408 : : {
409 [ # # ]: 0 : assert( mSize == other.mSize );
410 [ # # ]: 0 : assert( !memcmp( mRowStart, other.mRowStart, sizeof( size_t ) * ( mSize + 1 ) ) );
411 [ # # ]: 0 : assert( !memcmp( mColIndex, other.mColIndex, sizeof( size_t ) * mRowStart[mSize] ) );
412 [ # # ]: 0 : for( unsigned i = 0; i < mRowStart[mSize]; ++i )
413 : 0 : mEntries[i] += other.mEntries[i];
414 : 0 : }
415 : :
416 : : /*! \param diag is an STL vector of size MsqHessian::size() . */
417 : 0 : void MsqHessian::get_diagonal_blocks( std::vector< Matrix3D >& diag, MsqError& /*err*/ ) const
418 : : {
419 : : // make sure we have enough memory, so that no reallocation is needed later.
420 [ # # ]: 0 : if( diag.size() != size() ) { diag.reserve( size() ); }
421 : :
422 [ # # ]: 0 : for( size_t i = 0; i < size(); ++i )
423 : : {
424 : 0 : diag[i] = mEntries[mRowStart[i]];
425 : : }
426 : 0 : }
427 : :
428 : : /*! compute a preconditioner used in the preconditioned conjugate gradient
429 : : algebraic solver. In fact, this computes \f$ M^{-1} \f$ .
430 : : */
431 : 319 : void MsqHessian::compute_preconditioner( MsqError& /*err*/ )
432 : : {
433 : : // reallocates arrays if size of the Hessian has changed too much.
434 [ - + ][ # # ]: 319 : if( mSize > precondArraySize || mSize < precondArraySize / 10 )
435 : : {
436 [ + + ][ + + ]: 28356 : delete[] mPreconditioner;
437 [ + - ][ + - ]: 34136 : mPreconditioner = new Matrix3D[mSize];
[ + + ][ # # ]
438 : : }
439 : :
440 : : Matrix3D* diag_block;
441 : : double sum, tmp;
442 : : size_t m;
443 : : // For each diagonal block, the (inverted) preconditioner is
444 : : // the inverse of the sum of the diagonal entries.
445 [ + + ]: 34136 : for( m = 0; m < mSize; ++m )
446 : : {
447 : 33817 : diag_block = mEntries + mRowStart[m]; // Gets block at position m,m .
448 : :
449 : : #if !DIAGONAL_PRECONDITIONER
450 : : // calculate LDL^T factorization of the diagonal block
451 : : // L = [1 pre[0][1] pre[0][2]]
452 : : // [0 1 pre[1][2]]
453 : : // [0 0 1 ]
454 : : // inv(D) = [pre[0][0] 0 0 ]
455 : : // [0 pre[1][1] 0 ]
456 : : // [0 0 pre[2][2]]
457 : :
458 : : // If the first method of calculating a preconditioner fails,
459 : : // use the diagonal method.
460 : 33817 : bool use_diag = false;
461 : :
462 [ + + ]: 33817 : if( fabs( ( *diag_block )[0][0] ) < DBL_EPSILON ) { use_diag = true; }
463 : : else
464 : : {
465 : 33416 : mPreconditioner[m][0][0] = 1.0 / ( *diag_block )[0][0];
466 : 33416 : mPreconditioner[m][0][1] = ( *diag_block )[0][1] * mPreconditioner[m][0][0];
467 : 33416 : mPreconditioner[m][0][2] = ( *diag_block )[0][2] * mPreconditioner[m][0][0];
468 : 33416 : sum = ( ( *diag_block )[1][1] - ( *diag_block )[0][1] * mPreconditioner[m][0][1] );
469 [ - + ]: 33416 : if( fabs( sum ) <= DBL_EPSILON )
470 : 0 : use_diag = true;
471 : : else
472 : : {
473 : 33416 : mPreconditioner[m][1][1] = 1.0 / sum;
474 : :
475 : 33416 : tmp = ( *diag_block )[1][2] - ( *diag_block )[0][2] * mPreconditioner[m][0][1];
476 : :
477 : 33416 : mPreconditioner[m][1][2] = mPreconditioner[m][1][1] * tmp;
478 : :
479 : 33416 : sum = ( ( *diag_block )[2][2] - ( *diag_block )[0][2] * mPreconditioner[m][0][2] -
480 : 33416 : mPreconditioner[m][1][2] * tmp );
481 : :
482 [ - + ]: 33416 : if( fabs( sum ) <= DBL_EPSILON )
483 : 0 : use_diag = true;
484 : : else
485 : 33416 : mPreconditioner[m][2][2] = 1.0 / sum;
486 : : }
487 : : }
488 [ + + ]: 33817 : if( use_diag )
489 : : #endif
490 : : {
491 : : // Either this is a fixed vertex or the diagonal block is not
492 : : // invertible. Switch to the diagonal preconditioner in this
493 : : // case.
494 : :
495 : 401 : sum = ( *diag_block )[0][0] + ( *diag_block )[1][1] + ( *diag_block )[2][2];
496 [ - + ]: 401 : if( fabs( sum ) > DBL_EPSILON )
497 : 0 : sum = 1 / sum;
498 : : else
499 : 401 : sum = 0.0;
500 : :
501 : 401 : mPreconditioner[m][0][0] = sum;
502 : 401 : mPreconditioner[m][0][1] = 0.0;
503 : 401 : mPreconditioner[m][0][2] = 0.0;
504 : 401 : mPreconditioner[m][1][1] = sum;
505 : 401 : mPreconditioner[m][1][2] = 0.0;
506 : 401 : mPreconditioner[m][2][2] = sum;
507 : : }
508 : : }
509 [ # # ]: 319 : }
510 : :
511 : : /*! uses the preconditionned conjugate gradient algebraic solver
512 : : to find d in \f$ H * d = -g \f$ .
513 : : \param x : the solution, usually the descent direction d.
514 : : \param b : -b will be the right hand side. Usually b is the gradient.
515 : : */
516 : 319 : void MsqHessian::cg_solver( Vector3D x[], Vector3D b[], MsqError& err )
517 : : {
518 : : MSQ_FUNCTION_TIMER( "MsqHessian::cg_solver" );
519 : :
520 : : // reallocates arrays if size of the Hessian has changed too much.
521 [ + + ][ - + ]: 319 : if( mSize > cgArraySizes || mSize < cgArraySizes / 10 )
522 : : {
523 [ - + ]: 55 : delete[] mR;
524 [ - + ]: 55 : delete[] mZ;
525 [ - + ]: 55 : delete[] mP;
526 [ - + ]: 55 : delete[] mW;
527 [ + - ][ + - ]: 5835 : mR = new Vector3D[mSize];
[ + + ]
528 [ + - ][ + - ]: 5835 : mZ = new Vector3D[mSize];
[ + + ]
529 [ + - ][ + - ]: 5835 : mP = new Vector3D[mSize];
[ + + ]
530 [ + - ][ + - ]: 5835 : mW = new Vector3D[mSize];
[ + + ]
531 : 55 : cgArraySizes = mSize;
532 : : }
533 : :
534 : : size_t i;
535 : : double alpha_, alpha, beta;
536 : 319 : double cg_tol = 1e-2; // 1e-2 will give a reasonably good solution (~1%).
537 : 319 : double norm_g = length( b, mSize );
538 : 319 : double norm_r = norm_g;
539 : : double rzm1; // r^T_{k-1} z_{k-1}
540 : : double rzm2; // r^T_{k-2} z_{k-2}
541 : :
542 [ - + ][ # # ]: 319 : this->compute_preconditioner( err );MSQ_CHKERR( err ); // get M^{-1} for diagonal blocks
543 : :
544 [ + + ]: 34136 : for( i = 0; i < mSize; ++i )
545 [ + - ]: 33817 : x[i] = 0.;
546 [ + + ]: 34136 : for( i = 0; i < mSize; ++i )
547 [ + - ]: 33817 : mR[i] = -b[i]; // r = -b because x_0 = 0 and we solve H*x = -b
548 : 319 : norm_g *= cg_tol;
549 : :
550 : 319 : this->apply_preconditioner( mZ, mR, err ); // solve Mz = r (computes z = M^-1 r)
551 [ + + ]: 34136 : for( i = 0; i < mSize; ++i )
552 : 33817 : mP[i] = mZ[i]; // p_1 = z_0
553 : 319 : rzm1 = inner( mZ, mR, mSize ); // inner product r_{k-1}^T z_{k-1}
554 : :
555 : 319 : size_t cg_iter = 0;
556 [ + + ][ + - ]: 2710 : while( ( norm_r > norm_g ) && ( cg_iter < maxCGiter ) )
557 : : {
558 : 2524 : ++cg_iter;
559 : :
560 : 2524 : axpy( mW, mSize, *this, mP, mSize, 0, 0, err ); // w = A * p_k
561 : :
562 : 2524 : alpha_ = inner( mP, mW, mSize ); // alpha_ = p_k^T A p_k
563 [ + + ]: 2524 : if( alpha_ <= 0.0 )
564 : : {
565 [ + + ]: 133 : if( 1 == cg_iter )
566 : : {
567 [ + + ]: 132 : for( i = 0; i < mSize; ++i )
568 : 129 : x[i] += mP[i]; // x_{k+1} = x_k + p_{k+1}
569 : : }
570 : 133 : break; // Newton goes on with this direction of negative curvature
571 : : }
572 : :
573 : 2391 : alpha = rzm1 / alpha_;
574 : :
575 [ + + ]: 236323 : for( i = 0; i < mSize; ++i )
576 [ + - ]: 233932 : x[i] += alpha * mP[i]; // x_{k+1} = x_k + alpha_{k+1} p_{k+1}
577 [ + + ]: 236323 : for( i = 0; i < mSize; ++i )
578 [ + - ]: 233932 : mR[i] -= alpha * mW[i]; // r_{k+1} = r_k - alpha_{k+1} A p_{k+1}
579 : 2391 : norm_r = length( mR, mSize );
580 : :
581 : 2391 : this->apply_preconditioner( mZ, mR, err ); // solve Mz = r (computes z = M^-1 r)
582 : :
583 : 2391 : rzm2 = rzm1;
584 : 2391 : rzm1 = inner( mZ, mR, mSize ); // inner product r_{k-1}^T z_{k-1}
585 : 2391 : beta = rzm1 / rzm2;
586 [ + + ]: 236323 : for( i = 0; i < mSize; ++i )
587 [ + - ][ + - ]: 233932 : mP[i] = mZ[i] + beta * mP[i]; // p_k = z_{k-1} + Beta_k * p_{k-1}
588 : : }
589 : 319 : }
590 : :
591 : 634 : void MsqHessian::product( Vector3D* v, const Vector3D* x ) const
592 : : {
593 : : // zero output array
594 : 634 : memset( v, 0, size() * sizeof( *v ) );
595 : :
596 : : // for each row
597 : 634 : const Matrix3D* m_iter = mEntries;
598 : 634 : const size_t* c_iter = mColIndex;
599 [ + + ]: 2108900 : for( size_t r = 0; r < size(); ++r )
600 : : {
601 : :
602 : : // diagonal entry
603 : 2108266 : plusEqAx( v[r], *m_iter, x[r] );
604 : 2108266 : ++m_iter;
605 : 2108266 : ++c_iter;
606 : :
607 : : // off-diagonal entires
608 [ + + ]: 11456873 : for( size_t c = mRowStart[r] + 1; c != mRowStart[r + 1]; ++c )
609 : : {
610 : 9348607 : plusEqAx( v[r], *m_iter, x[*c_iter] );
611 : 9348607 : plusEqTransAx( v[*c_iter], *m_iter, x[r] );
612 : 9348607 : ++m_iter;
613 : 9348607 : ++c_iter;
614 : : }
615 : : }
616 : 634 : }
617 : :
618 : : /* ------------------ I/O ----------------- */
619 : :
620 : : //! Prints out the MsqHessian blocks.
621 : 0 : std::ostream& operator<<( std::ostream& s, const MsqHessian& h )
622 : : {
623 : : size_t i, j;
624 : 0 : s << "MsqHessian of size: " << h.mSize << "x" << h.mSize << "\n";
625 [ # # ]: 0 : for( i = 0; i < h.mSize; ++i )
626 : : {
627 : 0 : s << " ROW " << i << " ------------------------\n";
628 [ # # ]: 0 : for( j = h.mRowStart[i]; j < h.mRowStart[i + 1]; ++j )
629 : : {
630 : 0 : s << " column " << h.mColIndex[j] << " ----\n";
631 : 0 : s << h.mEntries[j];
632 : : }
633 : : }
634 : 0 : return s;
635 : : }
636 : :
637 [ + - ][ + - ]: 120 : } // namespace MBMesquite
|
