MOAB: Mesh Oriented datABase  (version 5.2.1)
MBMesquite::MsqHessian Class Reference

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.
Matrix3Dget_block (size_t i, size_t j)
const Matrix3Dget_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

Matrix3DmEntries
 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).
Matrix3DmPreconditioner
size_t precondArraySize
Vector3DmR
 array used in the CG solver
Vector3DmZ
 array used in the CG solver
Vector3DmP
 array used in the CG solver
Vector3DmW
 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

MsqHessianoperator= (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

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

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 213 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 304 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 516 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}
    }
}

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

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

Definition at line 431 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 338 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 352 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:
diagis 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]];
    }
}

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().

Referenced by ObjectiveFunctionTests::compare_hessian_diagonal(), ObjectiveFunctionTests::compare_hessian_gradient(), ObjectiveFunctionTests::compare_numerical_hessian(), ObjectiveFunctionTests::evaluate_internal(), MBMesquite::CompositeOFAdd::evaluate_with_Hessian(), MBMesquite::ObjectiveFunction::evaluate_with_Hessian_diagonal(), CompositeOFTest::get_hessians(), MBMesquite::TrustRegion::optimize_vertex_positions(), MBMesquite::FeasibleNewton::optimize_vertex_positions(), ObjectiveFunctionTest::test_compute_ana_hessian_tet(), ObjectiveFunctionTest::test_compute_ana_hessian_tet_scaled(), ObjectiveFunctionTests::test_handles_invalid_qm(), ObjectiveFunctionTests::test_handles_qm_error(), PMeanPTemplateTest::test_Hessian(), StdDevTemplateTest::test_hessian_fails(), StdDevTemplateTest::test_hessian_fails_sqr(), CompositeOFTest::test_multiply_hessian(), and ObjectiveFunctionTests::test_negate_flag().

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

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

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

Definition at line 151 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 246 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 621 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;
}

Member Data Documentation

size of arrays allocated in the CG solver.

Definition at line 84 of file MsqHessian.hpp.

Referenced by cg_solver(), and clear().

max nb of iterations of the CG solver.

Definition at line 85 of file MsqHessian.hpp.

Referenced by cg_solver().

size_t* MBMesquite::MsqHessian::mColIndex [protected]

CSR block structure: column indexes of the row entries.

Definition at line 73 of file MsqHessian.hpp.

Referenced by add(), MBMesquite::axpy(), clear(), get_block(), initialize(), MBMesquite::operator<<(), and product().

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

Definition at line 71 of file MsqHessian.hpp.

Referenced by add(), MBMesquite::axpy(), clear(), compute_preconditioner(), get_block(), get_diagonal_blocks(), initialize(), norm(), MBMesquite::operator<<(), product(), scale(), and zero_out().

array used in the CG solver

Definition at line 82 of file MsqHessian.hpp.

Referenced by cg_solver(), and clear().

array used in the CG solver

Definition at line 80 of file MsqHessian.hpp.

Referenced by cg_solver(), and clear().

size_t* MBMesquite::MsqHessian::mRowStart [protected]

start of each row in mEntries. size: nb of vertices (mSize).

Definition at line 72 of file MsqHessian.hpp.

Referenced by add(), MBMesquite::axpy(), clear(), compute_preconditioner(), get_block(), get_diagonal_blocks(), initialize(), norm(), MBMesquite::operator<<(), product(), scale(), and zero_out().

size_t MBMesquite::MsqHessian::mSize [protected]

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

Definition at line 75 of file MsqHessian.hpp.

Referenced by add(), apply_preconditioner(), MBMesquite::axpy(), cg_solver(), clear(), compute_preconditioner(), get_block(), initialize(), norm(), MBMesquite::operator<<(), scale(), and zero_out().

array used in the CG solver

Definition at line 83 of file MsqHessian.hpp.

Referenced by cg_solver(), and clear().

array used in the CG solver

Definition at line 81 of file MsqHessian.hpp.

Referenced by cg_solver(), and clear().

Definition at line 78 of file MsqHessian.hpp.

Referenced by clear(), and compute_preconditioner().

List of all members.


The documentation for this class was generated from the following files:
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines