nmrLSISolver.h

Go to the documentation of this file.

/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-    */
/* ex: set filetype=cpp softtabstop=4 shiftwidth=4 tabstop=4 cindent expandtab: */

/*

  Author(s):  Ankur Kapoor
  Created on: 2004-10-30

  (C) Copyright 2004-2013 Johns Hopkins University (JHU), All Rights
  Reserved.

--- begin cisst license - do not edit ---

This software is provided "as is" under an open source license, with
no warranty.  The complete license can be found in license.txt and
http://www.cisst.org/cisst/license.txt.

--- end cisst license ---
*/


#ifndef _nmrLSISolver_h
#define _nmrLSISolver_h

#include <cisstVector/vctDynamicMatrix.h>
#include <cisstNumerical/nmrSVDRSSolver.h>
#include <cisstNumerical/nmrLDPSolver.h>
#include <cisstNumerical/nmrNetlib.h>
#include <cisstNumerical/nmrExport.h>


class nmrLSISolver {
    /* Uses lawson and hanson's SVD, which gives U^{T}B and V,
       along with singular values */
protected:
    CISSTNETLIB_INTEGER Ma;
    CISSTNETLIB_INTEGER Na;
    CISSTNETLIB_INTEGER Mg;
    vctDynamicMatrix<CISSTNETLIB_DOUBLE> VSiUtb; /* V*S^{-1}*U^{T}*b */
    vctDynamicMatrix<CISSTNETLIB_DOUBLE> GTilde; /* \tilde{G} = G*V*S^{-1} */
    vctDynamicMatrix<CISSTNETLIB_DOUBLE> HTilde; /* \tilde{H} = h-G*V*S^{-1}*U^{T}*b */
    vctDynamicMatrix<CISSTNETLIB_DOUBLE> Z;      /* work area */
    vctDynamicMatrix<CISSTNETLIB_DOUBLE> X;      /* solution */
    vctDynamicMatrix<CISSTNETLIB_DOUBLE> Si;     /* Dim(n, n) matrix of singular values inverse*/
    vctDynamicMatrix<CISSTNETLIB_DOUBLE> VSi;    /* matrix of V*S^{-1} */
    nmrSVDRSSolver SVDRS;
    nmrLDPSolver LDP;

public:
    nmrLSISolver(void):
        Ma(0),
        Na(0),
        Mg(0)
    {
        Allocate(Ma, Na, Mg);
    }


    nmrLSISolver(CISSTNETLIB_INTEGER ma, CISSTNETLIB_INTEGER na, CISSTNETLIB_INTEGER mg) {
        Allocate(ma, na, mg);
    }


    nmrLSISolver(vctDynamicMatrix<CISSTNETLIB_DOUBLE> &C, vctDynamicMatrix<CISSTNETLIB_DOUBLE> &d,
                 vctDynamicMatrix<CISSTNETLIB_DOUBLE> &A, vctDynamicMatrix<CISSTNETLIB_DOUBLE> &b) {
        Allocate(C, d, A, b);
    }


    inline void Allocate(CISSTNETLIB_INTEGER ma, CISSTNETLIB_INTEGER na, CISSTNETLIB_INTEGER mg) {
        Ma = ma;
        Na = na;
        Mg = mg;
        VSiUtb.SetSize(std::max(Ma, Na), 1, VCT_COL_MAJOR);
        GTilde.SetSize(Mg, Na, VCT_COL_MAJOR);
        HTilde.SetSize(Mg, 1, VCT_COL_MAJOR);
        Z.SetSize(Na, 1, VCT_COL_MAJOR);
        X.SetSize(Na, 1, VCT_COL_MAJOR);
        Si.SetSize(Na, Ma, VCT_COL_MAJOR);
        VSi.SetSize(Na, Na, VCT_COL_MAJOR);
        SVDRS.Allocate(Ma, Na);
        LDP.Allocate(Mg, Na);
    }


    inline void Allocate(vctDynamicMatrix<CISSTNETLIB_DOUBLE> &C, vctDynamicMatrix<CISSTNETLIB_DOUBLE> &d,
                         vctDynamicMatrix<CISSTNETLIB_DOUBLE> &A, vctDynamicMatrix<CISSTNETLIB_DOUBLE> &b) {
        Allocate(C.rows(), C.cols(), A.rows());
    }


    CISST_EXPORT void Solve(vctDynamicMatrix<CISSTNETLIB_DOUBLE> &C, vctDynamicMatrix<CISSTNETLIB_DOUBLE> &d,
                            vctDynamicMatrix<CISSTNETLIB_DOUBLE> &A, vctDynamicMatrix<CISSTNETLIB_DOUBLE> &b) throw (std::runtime_error);

    inline const vctDynamicMatrix<CISSTNETLIB_DOUBLE> &GetX(void) const {
        return X;
    }
};


#endif // _nmrLSISolver_h