Vector3D is the object that effeciently stores the objective function Hessian each entry is a Matrix3D object (i.e. a vertex Hessian). More...

#include <MsqHessian.hpp>

Collaboration diagram for MsqHessian:

[legend]

Public Member Functions
	MsqHessian ()

	~MsqHessian ()

void	initialize (PatchData &pd, MsqError &err)
	creates a sparse structure for a Hessian, based on the connectivity information contained in the PatchData. Only the upper triangular part of the Hessian is stored. More...

void	zero_out ()

size_t	size ()

void	get_diagonal_blocks (msq_std::vector< Matrix3D > &diag, MsqError &err)
	returns the diagonal blocks, memory must be allocated before call. More...

Matrix3D *	get_block (size_t i, size_t j)

void	accumulate_entries (PatchData &pd, const size_t &elem_index, Matrix3D *const &mat3d_array, MsqError &err)

void	compute_preconditioner (MsqError &err)

void	apply_preconditioner (Vector3D zloc[], Vector3D rloc[], MsqError &err)

void	cg_solver (Vector3D x[], Vector3D b[], MsqError &err)

	MsqHessian ()

	~MsqHessian ()

void	initialize (PatchData &pd, MsqError &err)

void	zero_out ()

size_t	size ()

void	get_diagonal_blocks (msq_std::vector< Matrix3D > &diag, MsqError &err)
	returns the diagonal blocks, memory must be allocated before call. More...

Matrix3D *	get_block (size_t i, size_t j)

void	accumulate_entries (PatchData &pd, const size_t &elem_index, Matrix3D *const &mat3d_array, MsqError &err)

void	compute_preconditioner (MsqError &err)

void	apply_preconditioner (Vector3D zloc[], Vector3D rloc[], MsqError &err)

void	cg_solver (Vector3D x[], Vector3D b[], MsqError &err)

Protected Attributes
PatchData *	origin_pd

MsqMeshEntity *	patchElemArray
	stored once during initialization for More...

Matrix3D *	mEntries
	CSR block entries. size: nb of nonzero blocks, i.e. mRowStart[mSize] . More...

size_t *	mRowStart
	start of each row in mEntries. size: nb of vertices (mSize). More...

size_t *	mColIndex
	CSR block structure: column indexes of the row entries. More...

int *	mAccumulation
	accumulation pattern instructions More...

size_t *	mAccumElemStart
	Starting index in mAccumulation for element i, i=1,... More...

size_t	mSize
	number of rows (or number of columns, this is a square matrix). More...

Matrix3D *	mPreconditioner

size_t	precondArraySize

Vector3D *	mR
	array used in the CG solver More...

Vector3D *	mZ
	array used in the CG solver More...

Vector3D *	mP
	array used in the CG solver More...

Vector3D *	mW
	array used in the CG solver More...

size_t	cgArraySizes
	size of arrays allocated in the CG solver. More...

size_t	maxCGiter
	max nb of iterations of the CG solver. More...

Friends
class	ObjectiveFunction

void	axpy (Vector3D res[], size_t size_r, const MsqHessian &H, const Vector3D x[], size_t size_x, const Vector3D y[], size_t size_y, MsqError &err)
	Hessian - vector product, summed with a second vector (optional). More...

msq_stdio::ostream &	operator<< (msq_stdio::ostream &, const MsqHessian &)
	Prints out the MsqHessian blocks. More...

void	axpy (Vector3D res[], size_t size_r, const MsqHessian &H, const Vector3D x[], size_t size_x, const Vector3D y[], size_t size_y, MsqError &err)
	Hessian - vector product, summed with a second vector (optional). More...

msq_stdio::ostream &	operator<< (msq_stdio::ostream &, const MsqHessian &)
	Prints out the MsqHessian blocks. More...

Detailed Description

Vector3D is the object that effeciently stores the objective function Hessian each entry is a Matrix3D object (i.e. a vertex Hessian).

Definition at line 75 of file includeLinks/MsqHessian.hpp.

Constructor & Destructor Documentation

MsqHessian ( )

Definition at line 62 of file Misc/MsqHessian.cpp.

                        :
   origin_pd(0), mEntries(0), mRowStart(0), mColIndex(0), 
   mAccumulation(0), mAccumElemStart(0), mSize(0), 
   mPreconditioner(0), precondArraySize(0),
   mR(0), mZ(0), mP(0), mW(0), cgArraySizes(0), maxCGiter(50)
 { }

~MsqHessian ( )

Definition at line 70 of file Misc/MsqHessian.cpp.

References MsqHessian::mAccumElemStart, MsqHessian::mAccumulation, MsqHessian::mColIndex, MsqHessian::mEntries, MsqHessian::mP, MsqHessian::mPreconditioner, MsqHessian::mR, MsqHessian::mRowStart, MsqHessian::mW, and MsqHessian::mZ.

 {
   delete[] mEntries;    
   delete[] mRowStart;   
   delete[] mColIndex; 
 
   delete[] mAccumulation;
   delete[] mAccumElemStart;
 
   delete[] mPreconditioner;
 
   delete[] mR;
   delete[] mZ;
   delete[] mP;
   delete[] mW;
 }

MsqHessian ( )

~MsqHessian ( )

Member Function Documentation

void accumulate_entries	(	PatchData &	pd,
		const size_t &	elem_index,
		Matrix3D *const &	mat3d_array,
		MsqError &	err
	)

inline

Accumulates entries of an element hessian into an objective function hessian. Make sure to use zero_out() before starting the accumulation process.

Parameters

pd,:	PatchData in that contains the element which Hessian we are accumulating in the Hessian matrix. This must be the same PatchData that was used in MsqHessian::initialize().
elem_index,:	index of the element in the PatchData.
mat3d_array,:	This is the upper triangular part of the element Hessian for all nodes, including fixed nodes, for which the entries must be null Matrix3Ds.
nb_mat3d.	The size of the mat3d_array: (n+1)n/2, where n is the number of nodes in the element.

Definition at line 153 of file includeLinks/MsqHessian.hpp.

References PatchData::get_element_array(), i, MsqError::INVALID_ARG, j, MsqHessian::mAccumElemStart, MsqHessian::mAccumulation, MsqHessian::mEntries, MSQ_SETERR, MsqHessian::origin_pd, Matrix3D::plus_transpose_equal(), and MsqMeshEntity::vertex_count().

Referenced by LPtoPTemplate::compute_analytical_hessian().

     {
       if (&pd != origin_pd) {
         MSQ_SETERR(err)( 
                     "Cannot accumulate elements from a different patch. "
                     "Use MsqHessian::initialize first.",
                     MsqError::INVALID_ARG ); 
         return;
       }
 
       size_t nve = pd.get_element_array(err)[elem_index].vertex_count(); 
       const size_t nb_mat3d = (nve+1)*nve/2;
     
       size_t e = mAccumElemStart[elem_index];
       size_t i;
       int j;
       for (i = 0; i < nb_mat3d; ++i) {
         j = mAccumulation[e++];
         if (j >= 0)
           mEntries[j] += mat3d_array[i];
         else
           mEntries[-j].plus_transpose_equal(mat3d_array[i]);
       }
     }

Here is the call graph for this function:

Here is the caller graph for this function:

void accumulate_entries	(	PatchData &	pd,
		const size_t &	elem_index,
		Matrix3D *const &	mat3d_array,
		MsqError &	err
	)

inline

void apply_preconditioner	(	Vector3D	zloc[],
		Vector3D	rloc[],
		MsqError &	err
	)

inline

Computes $z=M^{-1}r$ .

Definition at line 240 of file includeLinks/MsqHessian.hpp.

References MsqHessian::mPreconditioner, and MsqHessian::mSize.

Referenced by MsqHessian::cg_solver().

     {
       size_t m;
 
       for (m=0; m<mSize; ++m) {
 #ifdef DIAGONAL_PRECONDITIONER
         // preconditioner is identity matrix for now.
         zloc[m][0] = mPreconditioner[m][0][0] * rloc[m][0]; 
         zloc[m][1] = mPreconditioner[m][1][1] * rloc[m][1]; 
         zloc[m][2] = mPreconditioner[m][2][2] * rloc[m][2]; 
 #else
         // z = inv(L^T) * r
         zloc[m][0] = rloc[m][0];
         zloc[m][1] = rloc[m][1] - mPreconditioner[m][0][1] * zloc[m][0];
         zloc[m][2] = rloc[m][2] - mPreconditioner[m][0][2] * zloc[m][0] - mPreconditioner[m][1][2] * zloc[m][1];
 
         // z = inv(D) * z
         zloc[m][0] *= mPreconditioner[m][0][0];
         zloc[m][1] *= mPreconditioner[m][1][1];
         zloc[m][2] *= mPreconditioner[m][2][2];
 
         // z = inv(L) * z
         zloc[m][2] = zloc[m][2];
         zloc[m][1] = zloc[m][1] - mPreconditioner[m][1][2] * zloc[m][2];
         zloc[m][0] = zloc[m][0] - mPreconditioner[m][0][1] * zloc[m][1] - mPreconditioner[m][0][2] * zloc[m][2];
 #endif
       }
     }

Here is the caller graph for this function:

void apply_preconditioner	(	Vector3D	zloc[],
		Vector3D	rloc[],
		MsqError &	err
	)

void cg_solver	(	Vector3D	x[],
		Vector3D	b[],
		MsqError &	err
	)

uses the preconditionned conjugate gradient algebraic solver to find d in $ H * d = -g $ .

Parameters

x	: the solution, usually the descent direction d.
b	: -b will be the right hand side. Usually b is the gradient.

Definition at line 440 of file Misc/MsqHessian.cpp.

References MsqHessian::apply_preconditioner(), MsqHessian::axpy, MsqHessian::cgArraySizes, MsqHessian::compute_preconditioner(), i, Mesquite::inner(), Mesquite::length(), MsqHessian::maxCGiter, MsqHessian::mP, MsqHessian::mR, MsqHessian::mSize, MSQ_CHKERR, MSQ_FUNCTION_TIMER, MsqHessian::mW, and MsqHessian::mZ.

Referenced by FeasibleNewton::optimize_vertex_positions().

 {
   MSQ_FUNCTION_TIMER( "MsqHessian::cg_solver" );
   
   // reallocates arrays if size of the Hessian has changed too much.
   if (mSize > cgArraySizes || mSize < cgArraySizes/10 ) {
     delete[] mR;
     delete[] mZ;
     delete[] mP;
     delete[] mW;
     mR = new Vector3D[mSize];
     mZ = new Vector3D[mSize];
     mP = new Vector3D[mSize];
     mW = new Vector3D[mSize];
     cgArraySizes = mSize;
   }
 
   size_t i;
   double alpha_, alpha, beta; 
   double cg_tol = 1e-2; // 1e-2 will give a reasonably good solution (~1%). 
   double norm_g = length(b, mSize);
   double norm_r = norm_g;
   double rzm1; // r^T_{k-1} z_{k-1}
   double rzm2; // r^T_{k-2} z_{k-2}
 
   this->compute_preconditioner(err); MSQ_CHKERR(err); // get M^{-1} for diagonal blocks
 
   for (i=0; i<mSize; ++i)  x[i] = 0. ;  
   for (i=0; i<mSize; ++i)  mR[i] = -b[i] ;  // r = -b because x_0 = 0 and we solve H*x = -b
   norm_g *= cg_tol;
 
   this->apply_preconditioner(mZ, mR, err); // solve Mz = r (computes z = M^-1 r)
   for (i=0; i<mSize; ++i)  mP[i] = mZ[i] ; // p_1 = z_0  
   rzm1 = inner(mZ,mR,mSize); // inner product r_{k-1}^T z_{k-1} 
     
   size_t cg_iter = 0;
   while ((norm_r > norm_g) && (cg_iter < maxCGiter)) {
     ++cg_iter;
       
     axpy(mW, mSize, *this, mP, mSize, 0,0,err); // w = A * p_k
       
     alpha_ = inner(mP,mW,mSize); // alpha_ = p_k^T A p_k
     if (alpha_ <= 0.0) {
       if (1 == cg_iter) {
         for (i=0; i<mSize; ++i)  x[i] += mP[i]; // x_{k+1} = x_k + p_{k+1} 
       }
       break; // Newton goes on with this direction of negative curvature 
     }
 
     alpha = rzm1 / alpha_;
       
     for (i=0; i<mSize; ++i)  x[i] += alpha*mP[i]; // x_{k+1} = x_k + alpha_{k+1} p_{k+1} 
     for (i=0; i<mSize; ++i)  mR[i] -= alpha*mW[i]; // r_{k+1} = r_k - alpha_{k+1} A p_{k+1} 
     norm_r = length(mR, mSize);
  
     this->apply_preconditioner(mZ, mR, err); // solve Mz = r (computes z = M^-1 r)
       
     rzm2 = rzm1;
     rzm1 = inner(mZ,mR,mSize); // inner product r_{k-1}^T z_{k-1} 
     beta = rzm1 / rzm2;
     for (i=0; i<mSize; ++i)  mP[i] = mZ[i] + beta*mP[i]; // p_k = z_{k-1} + Beta_k * p_{k-1}
   }
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void cg_solver	(	Vector3D	x[],
		Vector3D	b[],
		MsqError &	err
	)

void compute_preconditioner ( MsqError & err )

compute a preconditioner used in the preconditioned conjugate gradient algebraic solver. In fact, this computes $M^{-1}$ .

Definition at line 358 of file Misc/MsqHessian.cpp.

References MsqHessian::mEntries, MsqHessian::mPreconditioner, MsqHessian::mRowStart, MsqHessian::mSize, and MsqHessian::precondArraySize.

Referenced by MsqHessian::cg_solver().

 {
   // reallocates arrays if size of the Hessian has changed too much.
   if (mSize > precondArraySize || mSize < precondArraySize/10 ) {
     delete[] mPreconditioner;
     mPreconditioner = new Matrix3D[mSize];
   }
 
   Matrix3D* diag_block;
   double sum, tmp;
   size_t m;
   // For each diagonal block, the (inverted) preconditioner is
   // the inverse of the sum of the diagonal entries.
   for (m=0; m<mSize; ++m) {
     diag_block = mEntries + mRowStart[m]; // Gets block at position m,m .
 
 #if DIAGONAL_PRECONDITIONER
     // find sum, and computes inverse, or 0 if sum = 0 .
     sum = (*diag_block)[0][0] + (*diag_block)[1][1] + (*diag_block)[2][2];
     double inv_sum;
     if (sum != 0.) 
       inv_sum = 1 / sum;
     else
       inv_sum = 0.;
     
     mPreconditioner[m][0][0] = inv_sum;
     mPreconditioner[m][1][1] = inv_sum;
     mPreconditioner[m][2][2] = inv_sum;
 #else
     // calculate LDL^T factorization of the diagonal block
     //  L = [1 pre[0][1] pre[0][2]]
     //      [0 1         pre[1][2]]
     //      [0 0         1        ]
     //  inv(D) = [pre[0][0] 0         0        ]
     //           [0         pre[1][1] 0        ]
     //           [0         0         pre[2][2]]
 
     if ((*diag_block)[0][0] == 0.0) {
       // Either this is a fixed vertex or the diagonal block is not
       // invertible.  Switch to the diagonal preconditioner in this
       // case.
 
       sum = (*diag_block)[0][0] + (*diag_block)[1][1] + (*diag_block)[2][2];
       if (sum != 0.0) 
         sum = 1 / sum;
       
       mPreconditioner[m][0][0] = sum;
       mPreconditioner[m][0][1] = 0.0;
       mPreconditioner[m][0][2] = 0.0;
       mPreconditioner[m][1][1] = sum;
       mPreconditioner[m][1][2] = 0.0;
       mPreconditioner[m][2][2] = sum;
     }
     else {
       mPreconditioner[m][0][0] = 1.0 / (*diag_block)[0][0];
       mPreconditioner[m][0][1] = (*diag_block)[0][1] * mPreconditioner[m][0][0];
       mPreconditioner[m][0][2] = (*diag_block)[0][2] * mPreconditioner[m][0][0];
 
       mPreconditioner[m][1][1] = 
         1.0 / ((*diag_block)[1][1] - 
                (*diag_block)[0][1] * mPreconditioner[m][0][1]);
      
       tmp = (*diag_block)[1][2] - 
             (*diag_block)[0][2] * mPreconditioner[m][0][1];
 
       mPreconditioner[m][1][2] = mPreconditioner[m][1][1] * tmp;
 
       mPreconditioner[m][2][2] = 
         1.0 / ((*diag_block)[2][2] - 
                (*diag_block)[0][2]*mPreconditioner[m][0][2] - 
                mPreconditioner[m][1][2]*tmp);
     }
 #endif
   }
 }

Here is the caller graph for this function:

Matrix3D * get_block	(	size_t	i,
		size_t	j
	)

inline

Returns a pointer to the Matrix3D block at position i,j if it exist. Returns the NULL pointer if position i,j (0-based) is a NULL entry. Note that block i,j must be in the upper triangular part of the (symetric) hessian.

Definition at line 276 of file includeLinks/MsqHessian.hpp.

References MsqHessian::mColIndex, MsqHessian::mEntries, MsqHessian::mRowStart, and MsqHessian::mSize.

     {
       size_t c;
       
       if (i >= mSize || j >= mSize || j < i)
         return NULL;
 
       for (c=mRowStart[i]; c<mRowStart[i+1]; ++c) {
         if (mColIndex[c] == j)
           return ( mEntries + c );
       }
     
       // if there is no block at position i,j (zero entry).
       return NULL;
     }

Matrix3D* get_block	(	size_t	i,
		size_t	j
	)

void get_diagonal_blocks	(	msq_std::vector< Matrix3D > &	diag,
		MsqError &	err
	)

returns the diagonal blocks, memory must be allocated before call.

Parameters

diag	is an STL vector of size MsqHessian::size() .

Definition at line 341 of file Misc/MsqHessian.cpp.

References i, MsqHessian::mEntries, MsqHessian::mRowStart, and MsqHessian::size().

 {
   // make sure we have enough memory, so that no reallocation is needed later.
   if (diag.size() != size()) {
     diag.reserve(size());
   }
 
   for (size_t i=0; i<size(); ++i) {
     diag[i] = mEntries[mRowStart[i]];
   }
 }

Here is the call graph for this function:

void get_diagonal_blocks	(	msq_std::vector< Matrix3D > &	diag,
		MsqError &	err
	)

returns the diagonal blocks, memory must be allocated before call.

void initialize	(	PatchData &	pd,
		MsqError &	err
	)

creates a sparse structure for a Hessian, based on the connectivity information contained in the PatchData. Only the upper triangular part of the Hessian is stored.

Definition at line 91 of file Misc/MsqHessian.cpp.

References PatchData::get_element_array(), MsqMeshEntity::get_vertex_index_array(), i, MsqError::INVALID_ARG, j, MsqHessian::mAccumElemStart, MsqHessian::mAccumulation, MsqHessian::mColIndex, MsqHessian::mEntries, MsqHessian::mRowStart, MsqHessian::mSize, MSQ_CHKERR, MSQ_FUNCTION_TIMER, MSQ_SETERR, PatchData::num_elements(), PatchData::num_vertices(), MsqHessian::origin_pd, MsqHessian::patchElemArray, rs(), and MsqMeshEntity::vertex_count().

Referenced by FeasibleNewton::optimize_vertex_positions().

 {
   MSQ_FUNCTION_TIMER( "MsqHession::initialize" );
   delete[] mEntries;
   delete[] mRowStart;
   delete[] mColIndex;
   delete[] mAccumulation;
   delete[] mAccumElemStart;
   
   size_t num_vertices = pd.num_vertices();
   size_t num_elements = pd.num_elements();
   size_t const * vtx_list;
   size_t e, r, rs, re, c, cs, ce, nz, nnz, nve, i, j;
   patchElemArray = pd.get_element_array(err); MSQ_CHKERR(err);
 
   if (num_vertices == 0) {
     MSQ_SETERR( err )( "No vertices in PatchData", MsqError::INVALID_ARG);
     return;
   }
 
   mSize = num_vertices;
 
   // Calculate the offsets for a CSC representation of the accumulation
   // pattern.
 
   size_t* col_start = new size_t[num_vertices + 1];
   mAccumElemStart = new size_t[num_elements+1];
   mAccumElemStart[0] = 0;
   
   for (i = 0; i < num_vertices; ++i) {
     col_start[i] = 0;
   }
 
   for (e = 0; e < num_elements; ++e) {
     nve = patchElemArray[e].vertex_count();
     vtx_list = patchElemArray[e].get_vertex_index_array();
     mAccumElemStart[e+1] = mAccumElemStart[e] + (nve+1)*nve/2;
     
     for (i = 0; i < nve; ++i) {
       r = vtx_list[i];
       
       for (j = i; j < nve; ++j) {
         c = vtx_list[j];
 
         if (r <= c) {
           col_start[c]++;
         }
         else {
           col_start[r]++;
         }
       }
     }
   }
 
   nz = 0;
   for (i = 0; i < num_vertices; ++i) {
     j = col_start[i];
     col_start[i] = nz;
     nz += j;
   }
   col_start[i] = nz;
 
   // Finished putting matrix into CSC representation
 
   int* row_instr = new int[5*nz];
   size_t* row_index = new size_t[nz];
 
   nz = 0;
   for (e = 0; e < num_elements; ++e) {
     nve = patchElemArray[e].vertex_count();
     vtx_list = patchElemArray[e].get_vertex_index_array();
 
     for (i = 0; i < nve; ++i) {
       r = vtx_list[i];
 
       for (j = i; j < nve; ++j) {
         c = vtx_list[j];
 
         if (r <= c) {
           row_index[col_start[c]] = r;
           row_instr[col_start[c]] = nz;
           ++col_start[c];
         }
         else {
           row_index[col_start[r]] = c;
           row_instr[col_start[r]] = -nz;
           ++col_start[r];
         }
         
         ++nz;
       }
     }
   }
 
   for (i = num_vertices-1; i > 0; --i) {
     col_start[i+1] = col_start[i];
   }
   col_start[1] = col_start[0];
   col_start[0] = 0;
 
   //   cout << "col_start: ";
   //   for (int t=0; t<num_vertices+1; ++t)
   //     cout << col_start[t] << " ";
   //   cout << endl;
   //   cout << "row_index: ";
   //   for (int t=0; t<nz; ++t)
   //     cout << row_index[t] << " ";
   //   cout << endl;
   //   cout << "row_instr: ";
   //   for (int t=0; t<nz; ++t)
   //     cout << row_instr[t] << " ";
   //   cout << endl;
   
   
   // Convert CSC to CSR
   // First calculate the offsets in the row
 
   size_t* row_start = new size_t[num_vertices + 1];
     
   for (i = 0; i < num_vertices; ++i) {
     row_start[i] = 0;
   }
 
   for (i = 0; i < nz; ++i) {
     ++row_start[row_index[i]];
   }
 
   nz = 0;
   for (i = 0; i < num_vertices; ++i) {
     j = row_start[i];
     row_start[i] = nz;
     nz += j;
   }
   row_start[i] = nz;
     
   // Now calculate the pattern
 
   size_t* col_index = new size_t[nz];
   int* col_instr = new int[nz];
 
   for (i = 0; i < num_vertices; ++i) {
     cs = col_start[i];
     ce = col_start[i+1];
 
     while(cs < ce) {
       r = row_index[cs];
 
       col_index[row_start[r]] = i;
       col_instr[row_start[r]] = row_instr[cs];
 
       ++row_start[r];
       ++cs;
     }
   }
 
   for (i = num_vertices-1; i > 0; --i) {
     row_start[i+1] = row_start[i];
   }
   row_start[1] = row_start[0];
   row_start[0] = 0;
 
   delete[] row_index;
 
   // Now that the matrix is CSR
   // Column indices for each row are sorted
 
   // Compaction -- count the number of nonzeros
   mRowStart = col_start;   // don't need to reallocate
   mAccumulation = row_instr;   // don't need to reallocate
 
   for (i = 0; i <= num_vertices; ++i) {
     mRowStart[i] = 0;
   }
 
   nnz = 0;
   for (i = 0; i < num_vertices; ++i) {
     rs = row_start[i];
     re = row_start[i+1];
 
     c = num_vertices;
     while(rs < re) {
       if (c != col_index[rs]) {
         // This is an unseen nonzero
 
         c = col_index[rs];
         ++mRowStart[i];
         ++nnz;
       }
 
       if (col_instr[rs] >= 0) {
         mAccumulation[col_instr[rs]] = nnz - 1;
       }
       else {
         mAccumulation[-col_instr[rs]] = 1 - nnz;
       }
       
       ++rs;
     }
   }
 
   nnz = 0;
   for (i = 0; i < num_vertices; ++i) {
     j = mRowStart[i];
     mRowStart[i] = nnz;
     nnz += j;
   }
   mRowStart[i] = nnz;
 
   delete [] col_instr;
 
   // Fill in the compacted hessian matrix
 
   mColIndex = new size_t[nnz];
 
   for (i = 0; i < num_vertices; ++i) {
     rs = row_start[i];
     re = row_start[i+1];
 
     c = num_vertices;
     while(rs < re) {
       if (c != col_index[rs]) {
         // This is an unseen nonzero
 
         c = col_index[rs];
         mColIndex[mRowStart[i]] = c;
         mRowStart[i]++;
       }
       ++rs;
     }
   }
 
   for (i = num_vertices-1; i > 0; --i) {
     mRowStart[i+1] = mRowStart[i];
   }
   mRowStart[1] = mRowStart[0];
   mRowStart[0] = 0;
   
   delete [] row_start;
   delete [] col_index;
 
   mEntries = new Matrix3D[nnz]; // On Solaris, no initializer allowed for new of an array 
   for (i=0;i<nnz;++i) mEntries[i] = 0.; // so we initialize all entries manually. 
 
   origin_pd = &pd;
 
   return;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void initialize	(	PatchData &	pd,
		MsqError &	err
	)

size_t size ( )

inline

Definition at line 107 of file src/Misc/MsqHessian.hpp.

References MsqHessian::mSize.

107 {return mSize;}

Mesquite::MsqHessian::mSize

size_t mSize

number of rows (or number of columns, this is a square matrix).

Definition: includeLinks/MsqHessian.hpp:89

size_t size ( )

inline

Definition at line 107 of file includeLinks/MsqHessian.hpp.

References MsqHessian::mSize.

Referenced by MsqHessian::get_diagonal_blocks().

107 {return mSize;}

Mesquite::MsqHessian::mSize

size_t mSize

number of rows (or number of columns, this is a square matrix).

Definition: includeLinks/MsqHessian.hpp:89

Here is the caller graph for this function:

void zero_out ( )

inline

void zero_out ( )

inline

Sets all Hessian entries to zero. This is usually used before starting to accumulate elements hessian in the objective function hessian.

Definition at line 129 of file includeLinks/MsqHessian.hpp.

References i, MsqHessian::mEntries, MsqHessian::mRowStart, MsqHessian::mSize, and Matrix3D::zero().

Referenced by LPtoPTemplate::compute_analytical_hessian().

     {
       if (mSize==0) return; // empty hessian.
     
       size_t i;
       for (i=0; i<mRowStart[mSize]; ++i) {
         mEntries[i].zero();
       }
     }

Here is the call graph for this function:

Here is the caller graph for this function:

Friends And Related Function Documentation

void axpy	(	Vector3D	res[],
		size_t	size_r,
		const MsqHessian &	H,
		const Vector3D	x[],
		size_t	size_x,
		const Vector3D	y[],
		size_t	size_y,
		MsqError &	err
	)

friend

Hessian - vector product, summed with a second vector (optional).

Parameters

res,:	array of Vector3D in which the result is stored.
size_r,:	size of the res array.
x,:	vector multiplied by the Hessian.
size_x,:	size of the x array.
y,:	vector added to the Hessian vector product. Set to 0 (NULL) if not needed.
size_y,:	size of the y array. Set to 0 if not needed.

Definition at line 187 of file includeLinks/MsqHessian.hpp.

Referenced by MsqHessian::cg_solver().

     {
       if ((size_r != H.mSize) || (size_x != H.mSize) ||
           (size_y != H.mSize && size_y != 0)) {
         // throw an error
       }
 
       Vector3D tmpx, tmpm; // for cache opt.
       size_t* col = H.mColIndex;
       const size_t nn = H.mSize;
       size_t rl; // row length
       size_t el; // entries index
       size_t lo;
       size_t c;  // column index
       size_t i, j;
      
       if (y != 0) {
         for (i = 0; i < nn; ++i) {
           res[i] = y[i];
         }
       } 
       else {          // y == 0
         for (i = 0; i < nn; ++i) {
           res[i] = 0.;
         }
       }
       
       el = 0;
       for (i = 0; i < nn; ++i) {
         rl = H.mRowStart[i+1] - H.mRowStart[i];
         lo = *col++;
 
         // Diagonal entry
         tmpx = x[i];
         eqAx(tmpm, H.mEntries[el], tmpx);
         ++el;
         
         //Non-diagonal entries
         for (j = 1; j < rl; ++j) {
           c = *col++;
           //          res[i] += H.mEntries[e] * x[c];
           plusEqAx(tmpm, H.mEntries[el], x[c]);
           //          res[c] += transpose(H.mEntries[e]) * tmpxi;
           plusEqTransAx(res[c], H.mEntries[el], tmpx);
           ++el;
         }
         res[lo] += tmpm;
       }
     }

void axpy	(	Vector3D	res[],
		size_t	size_r,
		const MsqHessian &	H,
		const Vector3D	x[],
		size_t	size_x,
		const Vector3D	y[],
		size_t	size_y,
		MsqError &	err
	)

friend

Hessian - vector product, summed with a second vector (optional).

Parameters

res,:	array of Vector3D in which the result is stored.
size_r,:	size of the res array.
x,:	vector multiplied by the Hessian.
size_x,:	size of the x array.
y,:	vector added to the Hessian vector product. Set to 0 (NULL) if not needed.
size_y,:	size of the y array. Set to 0 if not needed.

Definition at line 187 of file includeLinks/MsqHessian.hpp.

     {
       if ((size_r != H.mSize) || (size_x != H.mSize) ||
           (size_y != H.mSize && size_y != 0)) {
         // throw an error
       }
 
       Vector3D tmpx, tmpm; // for cache opt.
       size_t* col = H.mColIndex;
       const size_t nn = H.mSize;
       size_t rl; // row length
       size_t el; // entries index
       size_t lo;
       size_t c;  // column index
       size_t i, j;
      
       if (y != 0) {
         for (i = 0; i < nn; ++i) {
           res[i] = y[i];
         }
       } 
       else {          // y == 0
         for (i = 0; i < nn; ++i) {
           res[i] = 0.;
         }
       }
       
       el = 0;
       for (i = 0; i < nn; ++i) {
         rl = H.mRowStart[i+1] - H.mRowStart[i];
         lo = *col++;
 
         // Diagonal entry
         tmpx = x[i];
         eqAx(tmpm, H.mEntries[el], tmpx);
         ++el;
         
         //Non-diagonal entries
         for (j = 1; j < rl; ++j) {
           c = *col++;
           //          res[i] += H.mEntries[e] * x[c];
           plusEqAx(tmpm, H.mEntries[el], x[c]);
           //          res[c] += transpose(H.mEntries[e]) * tmpxi;
           plusEqTransAx(res[c], H.mEntries[el], tmpx);
           ++el;
         }
         res[lo] += tmpm;
       }
     }