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>

Inheritance diagram for MBMesquite::MsqHessian:

Collaboration diagram for MBMesquite::MsqHessian:

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.
void	initialize (const MsqHessian &other)
void	zero_out ()
size_t	size () const
void	get_diagonal_blocks (std::vector< Matrix3D > &diag, MsqError &err) const
	returns the diagonal blocks, memory must be allocated before call.
Matrix3D *	get_block (size_t i, size_t j)
const Matrix3D *	get_block (size_t i, size_t j) const
void	compute_preconditioner (MsqError &err)
void	apply_preconditioner (Vector3D zloc[], Vector3D rloc[], MsqError &err)
void	cg_solver (Vector3D x[], Vector3D b[], MsqError &err)
void	product (Vector3D r, const Vector3D x) const
	r = this * x, where r and x are arrays of length size().
void	add (size_t row, size_t col, const Matrix3D &m, MsqError &err)
void	scale (double value)
void	add (const MsqHessian &other)
void	clear ()
double	norm () const
	Frobenius norm.
Protected Attributes
Matrix3D *	mEntries
	CSR block entries. size: nb of nonzero blocks, i.e. mRowStart[mSize] .
size_t *	mRowStart
	start of each row in mEntries. size: nb of vertices (mSize).
size_t *	mColIndex
	CSR block structure: column indexes of the row entries.
size_t	mSize
	number of rows (or number of columns, this is a square matrix).
Matrix3D *	mPreconditioner
size_t	precondArraySize
Vector3D *	mR
	array used in the CG solver
Vector3D *	mZ
	array used in the CG solver
Vector3D *	mP
	array used in the CG solver
Vector3D *	mW
	array used in the CG solver
size_t	cgArraySizes
	size of arrays allocated in the CG solver.
size_t	maxCGiter
	max nb of iterations of the CG solver.
Private Member Functions
MsqHessian &	operator= (const MsqHessian &h)
	MsqHessian (const MsqHessian &copy)
Friends
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).
std::ostream &	operator<< (std::ostream &, const MsqHessian &)
	Prints out the MsqHessian blocks.

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 68 of file MsqHessian.hpp.

Constructor & Destructor Documentation

MBMesquite::MsqHessian::MsqHessian ( )

Definition at line 54 of file MsqHessian.cpp.

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

MBMesquite::MsqHessian::~MsqHessian ( )

Definition at line 60 of file MsqHessian.cpp.

References clear().

{
    clear();
}

MBMesquite::MsqHessian::MsqHessian ( const MsqHessian & copy ) [private]

Member Function Documentation

void MBMesquite::MsqHessian::add	(	size_t	row,
		size_t	col,
		const Matrix3D &	m,
		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 219 of file MsqHessian.hpp.

References MBMesquite::MsqError::INVALID_ARG, mColIndex, mEntries, mRowStart, MSQ_SETERR, and MBMesquite::Matrix3D::plus_transpose_equal().

Referenced by MBMesquite::PatchPowerMeanP::evaluate_with_Hessian(), MBMesquite::CompositeOFAdd::evaluate_with_Hessian(), MBMesquite::LPtoPTemplate::evaluate_with_Hessian(), and MBMesquite::PMeanPTemplate::evaluate_with_Hessian().

{
    if( row <= col )
    {
        for( size_t i = mRowStart[row]; i != mRowStart[row + 1]; ++i )
            if( mColIndex[i] == col )
            {
                mEntries[i] += m;
                return;
            }
    }
    else
    {
        for( size_t i = mRowStart[col]; i != mRowStart[col + 1]; ++i )
            if( mColIndex[i] == row )
            {
                mEntries[i].plus_transpose_equal( m );
                return;
            }
    }

    MSQ_SETERR( err )
    ( MsqError::INVALID_ARG, "Hessian entry (%lu,%lu) does not exist.", (unsigned long)row, (unsigned long)col );
}

void MBMesquite::MsqHessian::add ( const MsqHessian & other )

Definition at line 407 of file MsqHessian.cpp.

References mColIndex, mEntries, mRowStart, and mSize.

{
    assert( mSize == other.mSize );
    assert( !memcmp( mRowStart, other.mRowStart, sizeof( size_t ) * ( mSize + 1 ) ) );
    assert( !memcmp( mColIndex, other.mColIndex, sizeof( size_t ) * mRowStart[mSize] ) );
    for( unsigned i = 0; i < mRowStart[mSize]; ++i )
        mEntries[i] += other.mEntries[i];
}

void MBMesquite::MsqHessian::apply_preconditioner	(	Vector3D	zloc[],
		Vector3D	rloc[],
		MsqError &	err
	)		`[inline]`

Computes \( z=M^{-1}r \) .

Definition at line 316 of file MsqHessian.hpp.

References mPreconditioner, and mSize.

Referenced by 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
    }
}

void MBMesquite::MsqHessian::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 522 of file MsqHessian.cpp.

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

Referenced by MBMesquite::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}
    }
}

void MBMesquite::MsqHessian::clear ( )

Definition at line 65 of file MsqHessian.cpp.

References cgArraySizes, mColIndex, mEntries, mP, mPreconditioner, mR, mRowStart, mSize, mW, mZ, and precondArraySize.

Referenced by MBMesquite::TrustRegion::cleanup(), and ~MsqHessian().

{
    mSize = precondArraySize = cgArraySizes = 0;

    delete[] mEntries;
    mEntries = 0;
    delete[] mRowStart;
    mRowStart = 0;
    delete[] mColIndex;
    mColIndex = 0;

    delete[] mPreconditioner;
    mPreconditioner = 0;

    delete[] mR;
    mR = 0;
    delete[] mZ;
    mZ = 0;
    delete[] mP;
    mP = 0;
    delete[] mW;
    mW = 0;
}

void MBMesquite::MsqHessian::compute_preconditioner ( MsqError & err )

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

Definition at line 434 of file MsqHessian.cpp.

References mEntries, mPreconditioner, mRowStart, mSize, precondArraySize, and moab::sum().

Referenced by cg_solver(), and MsqHessianTest::test_cholesky_preconditioner().

{
    // 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
        // 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 the first method of calculating a preconditioner fails,
        // use the diagonal method.
        bool use_diag = false;

        if( fabs( ( *diag_block )[0][0] ) < DBL_EPSILON )
        {
            use_diag = true;
        }
        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];
            sum                      = ( ( *diag_block )[1][1] - ( *diag_block )[0][1] * mPreconditioner[m][0][1] );
            if( fabs( sum ) <= DBL_EPSILON )
                use_diag = true;
            else
            {
                mPreconditioner[m][1][1] = 1.0 / sum;

                tmp = ( *diag_block )[1][2] - ( *diag_block )[0][2] * mPreconditioner[m][0][1];

                mPreconditioner[m][1][2] = mPreconditioner[m][1][1] * tmp;

                sum = ( ( *diag_block )[2][2] - ( *diag_block )[0][2] * mPreconditioner[m][0][2] -
                        mPreconditioner[m][1][2] * tmp );

                if( fabs( sum ) <= DBL_EPSILON )
                    use_diag = true;
                else
                    mPreconditioner[m][2][2] = 1.0 / sum;
            }
        }
        if( use_diag )
#endif
        {
            // 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( fabs( sum ) > DBL_EPSILON )
                sum = 1 / sum;
            else
                sum = 0.0;

            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;
        }
    }
}

Matrix3D * MBMesquite::MsqHessian::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 350 of file MsqHessian.hpp.

References mColIndex, mEntries, mRowStart, and mSize.

Referenced by ObjectiveFunctionTests::compare_hessian_diagonal(), ObjectiveFunctionTests::compare_numerical_hessian(), MBMesquite::TrustRegion::compute_preconditioner(), MBMesquite::ObjectiveFunction::evaluate_with_Hessian_diagonal(), get_block(), CompositeOFTest::test_add_hessian(), ObjectiveFunctionTest::test_compute_ana_hessian_tet(), ObjectiveFunctionTest::test_compute_ana_hessian_tet_scaled(), PMeanPTemplateTest::test_Hessian(), ObjectiveFunctionTests::test_negate_flag(), CompositeOFTest::test_scalar_add_hessian(), and CompositeOFTest::test_scalar_multiply_hessian().

{
    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;
}

const Matrix3D * MBMesquite::MsqHessian::get_block	(	size_t	i,
		size_t	j
	)		const `[inline]`

Definition at line 364 of file MsqHessian.hpp.

References get_block().

{
    return const_cast< MsqHessian* >( this )->get_block( i, j );
}

void MBMesquite::MsqHessian::get_diagonal_blocks	(	std::vector< Matrix3D > &	diag,
		MsqError &	err
	)		const

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

Parameters:

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

Definition at line 417 of file MsqHessian.cpp.

References mEntries, mRowStart, and 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]];
    }
}

void MBMesquite::MsqHessian::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 92 of file MsqHessian.cpp.

References MBMesquite::PatchData::get_element_array(), MBMesquite::MsqMeshEntity::get_vertex_index_array(), MBMesquite::MsqError::INVALID_ARG, mColIndex, mEntries, mRowStart, mSize, MSQ_CHKERR, MSQ_FUNCTION_TIMER, MSQ_SETERR, MBMesquite::MsqMeshEntity::node_count(), MBMesquite::PatchData::num_elements(), and MBMesquite::PatchData::num_free_vertices().

{
    MSQ_FUNCTION_TIMER( "MsqHession::initialize" );
    delete[] mEntries;
    delete[] mRowStart;
    delete[] mColIndex;

    size_t num_vertices = pd.num_free_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;
    MsqMeshEntity* 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].node_count();
        vtx_list = patchElemArray[e].get_vertex_index_array();
        int nfe  = 0;

        for( i = 0; i < nve; ++i )
        {
            r = vtx_list[i];
            if( r < num_vertices ) ++nfe;

            for( j = i; j < nve; ++j )
            {
                c = vtx_list[j];

                if( r <= c )
                {
                    if( c < num_vertices ) ++col_start[c];
                }
                else
                {
                    if( r < num_vertices ) ++col_start[r];
                }
            }
        }
        // mAccumElemStart[e+1] = mAccumElemStart[e] + (nfe+1)*nfe/2;
    }

    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].node_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 )
                {
                    if( c < num_vertices )
                    {
                        row_index[col_start[c]] = r;
                        row_instr[col_start[c]] = nz++;
                        ++col_start[c];
                    }
                }
                else
                {
                    if( r < num_vertices )
                    {
                        row_index[col_start[r]] = c;
                        // can't use -nz, but can negate row_instr[col_start[r]]
                        row_instr[col_start[r]] = nz++;
                        row_instr[col_start[r]] = -row_instr[col_start[r]];
                        ++col_start[r];
                    }
                }
            }
        }
    }

    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
    delete[] row_instr;

    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;
}

void MBMesquite::MsqHessian::initialize ( const MsqHessian & other )

Definition at line 376 of file MsqHessian.cpp.

References mColIndex, mEntries, mRowStart, and mSize.

{
    if( !other.mSize )
    {
        delete[] mEntries;
        delete[] mRowStart;
        delete[] mColIndex;
        mEntries  = 0;
        mRowStart = 0;
        mColIndex = 0;
        mSize     = 0;
        return;
    }

    if( mSize != other.mSize || mRowStart[mSize] != other.mRowStart[mSize] )
    {
        delete[] mEntries;
        delete[] mRowStart;
        delete[] mColIndex;

        mSize = other.mSize;

        mRowStart = new size_t[mSize + 1];
        mEntries  = new Matrix3D[other.mRowStart[mSize]];
        mColIndex = new size_t[other.mRowStart[mSize]];
    }

    memcpy( mRowStart, other.mRowStart, sizeof( size_t ) * ( mSize + 1 ) );
    memcpy( mColIndex, other.mColIndex, sizeof( size_t ) * mRowStart[mSize] );
}

double MBMesquite::MsqHessian::norm ( ) const [inline]

Frobenius norm.

Definition at line 146 of file MsqHessian.hpp.

References MBMesquite::Frobenius_2(), mEntries, mRowStart, mSize, and moab::sum().

Referenced by MBMesquite::FeasibleNewton::optimize_vertex_positions().

{
    double sum = 0.0;
    for( size_t i = 0; i < mRowStart[mSize]; ++i )
        sum += Frobenius_2( mEntries[i] );
    return sqrt( sum );
}

MsqHessian& MBMesquite::MsqHessian::operator= ( const MsqHessian & h ) [private]

void MBMesquite::MsqHessian::product	(	Vector3D *	r,
		const Vector3D *	x
	)		const

r = this * x, where r and x are arrays of length size().

Definition at line 597 of file MsqHessian.cpp.

References mColIndex, mEntries, mRowStart, MBMesquite::plusEqAx(), MBMesquite::plusEqTransAx(), and size().

Referenced by MBMesquite::TrustRegion::optimize_vertex_positions().

{
    // zero output array
    memset( v, 0, size() * sizeof( *v ) );

    // for each row
    const Matrix3D* m_iter = mEntries;
    const size_t* c_iter   = mColIndex;
    for( size_t r = 0; r < size(); ++r )
    {

        // diagonal entry
        plusEqAx( v[r], *m_iter, x[r] );
        ++m_iter;
        ++c_iter;

        // off-diagonal entires
        for( size_t c = mRowStart[r] + 1; c != mRowStart[r + 1]; ++c )
        {
            plusEqAx( v[r], *m_iter, x[*c_iter] );
            plusEqTransAx( v[*c_iter], *m_iter, x[r] );
            ++m_iter;
            ++c_iter;
        }
    }
}

void MBMesquite::MsqHessian::scale ( double value ) [inline]

Definition at line 168 of file MsqHessian.hpp.

References mEntries, mRowStart, and mSize.

Referenced by MBMesquite::PatchPowerMeanP::evaluate_with_Hessian(), MBMesquite::CompositeOFScalarMultiply::evaluate_with_Hessian(), MBMesquite::LPtoPTemplate::evaluate_with_Hessian(), and MBMesquite::PMeanPTemplate::evaluate_with_Hessian().

{
    for( size_t i = 0; i < mRowStart[mSize]; ++i )
        mEntries[i] *= value;
}

size_t MBMesquite::MsqHessian::size ( ) const [inline]

Definition at line 96 of file MsqHessian.hpp.

Referenced by ObjectiveFunctionTests::compare_hessian_diagonal(), MBMesquite::TrustRegion::compute_preconditioner(), MBMesquite::ObjectiveFunction::evaluate_with_Hessian_diagonal(), get_diagonal_blocks(), product(), CompositeOFTest::test_add_hessian(), PMeanPTemplateTest::test_Hessian(), ObjectiveFunctionTests::test_negate_flag(), CompositeOFTest::test_scalar_add_hessian(), and CompositeOFTest::test_scalar_multiply_hessian().

    {
        return mSize;
    }

void MBMesquite::MsqHessian::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 157 of file MsqHessian.hpp.

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

Referenced by MBMesquite::PatchPowerMeanP::evaluate_with_Hessian(), MBMesquite::LPtoPTemplate::evaluate_with_Hessian(), MBMesquite::PMeanPTemplate::evaluate_with_Hessian(), and MsqHessianTest::test_cholesky_preconditioner().

{
    if( mSize == 0 ) return;  // empty hessian.

    size_t i;
    for( i = 0; i < mRowStart[mSize]; ++i )
    {
        mEntries[i].zero();
    }
}

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 252 of file MsqHessian.hpp.

Referenced by 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;
    }
}

std::ostream& operator<<	(	std::ostream &	,
		const MsqHessian &
	)		`[friend]`

Prints out the MsqHessian blocks.

Definition at line 627 of file MsqHessian.cpp.

{
    size_t i, j;
    s << "MsqHessian of size: " << h.mSize << "x" << h.mSize << "\n";
    for( i = 0; i < h.mSize; ++i )
    {
        s << " ROW " << i << " ------------------------\n";
        for( j = h.mRowStart[i]; j < h.mRowStart[i + 1]; ++j )
        {
            s << "   column " << h.mColIndex[j] << " ----\n";
            s << h.mEntries[j];
        }
    }
    return s;
}