root/library/bdm/math/chmat.cpp

#include "chmat.h"

//using std::endl;

namespace bdm {

void chmat::add ( const chmat &A2, double w ) {
    bdm_assert_debug ( dim == A2.dim, "Matrices of unequal dimension" );
    mat pre = concat_vertical ( Ch, sqrt ( w ) * A2.Ch );
    mat post = zeros ( pre.rows(), pre.cols() );
    if ( !qr ( pre, post ) ) {
        bdm_warning ( "Unstable QR in chmat add" );
    }
    Ch = post ( 0, dim - 1, 0, dim - 1 );
};

void chmat::opupdt ( const vec &v, double w ) {
//TODO see cholupdt in lhotse
    mat Z;
    mat R;
    mat V ( 1, v.length() );
    V.set_row ( 0, v*sqrt ( w ) );
    Z = concat_vertical ( Ch, V );
    qr ( Z, R );
    Ch = R ( 0, Ch.rows() - 1, 0, Ch.cols() - 1 );
}

mat chmat::to_mat() const {
    mat F = Ch.T() * Ch;
    return F;
}

void chmat::mult_sym ( const mat &C ) {
    bdm_assert_debug ( C.rows() == dim, "Wrong dimension of U" );
    if ( !qr ( Ch*C.T(), Ch ) ) {
        bdm_warning ( "QR unstable in chmat mult_sym" );
    }
}

void chmat::mult_sym ( const mat &C , chmat &U ) const {
    bdm_assert_debug ( C.rows() == U.dim, "Wrong dimension of U" );
    if ( !qr ( Ch*C.T(), U.Ch ) ) {
        bdm_warning ( "QR unstable in chmat mult_sym" );
    }
}

void chmat::mult_sym_t ( const mat &C ) {
    bdm_assert_debug ( C.cols() == dim, "Wrong dimension of U" );
    if ( !qr ( Ch*C, Ch ) ) {
        bdm_warning ( "QR unstable in chmat mult_sym" );
    }
}

void chmat::mult_sym_t ( const mat &C, chmat &U ) const {
    bdm_assert_debug ( C.cols() == U.dim, "Wrong dimension of U" );
    if ( !qr ( Ch*C, U.Ch ) ) {
        bdm_warning ( "QR unstable in chmat mult_sym" );
    }
}

double chmat::logdet() const {
    double ldet = 0.0;
    int i;
    //sum of logs of (possibly negative!) diagonal entries
    for ( i = 0; i < Ch.rows(); i++ ) {
        ldet += log ( std::fabs ( Ch ( i, i ) ) );
    }
    return 2*ldet; //compensate for Ch being sqrt()
}

//TODO can be done more efficiently using BLAS, see triangular matrices
vec chmat::sqrt_mult ( const vec &v ) const {
    vec pom;
    pom = Ch.T() * v;
    return pom;
}

double chmat::qform ( const vec &v ) const {
    vec pom;
    pom = Ch * v;
    return pom*pom;
}

double chmat::invqform ( const vec &v ) const {
    vec pom ( v.length() );
    forward_substitution ( Ch.T(), v, pom );
    return pom*pom;
}

void chmat::clear() {
    Ch.clear();
}

}