nmrSVDRSSolver.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 _nmrSVDRSSolver_h
#define _nmrSVDRSSolver_h


#include <cisstVector/vctDynamicMatrix.h>
#include <cisstNumerical/nmrNetlib.h>


class nmrSVDRSSolver {

protected:
    CISSTNETLIB_INTEGER M;                    /* No. of rows of A, B, G */
    CISSTNETLIB_INTEGER N;                    /* No. of cols of A, No of rows and cols of V */
    CISSTNETLIB_INTEGER Lda;                  /* >= max(m, n) */
    CISSTNETLIB_INTEGER Ldb;                  /* >= m */
    CISSTNETLIB_INTEGER Nb;                   /* No. of cols of B and G */
    vctDynamicMatrix<CISSTNETLIB_DOUBLE> S;   /* Singular values, Dim(S)>n */
    vctDynamicMatrix<CISSTNETLIB_DOUBLE> Work;/* Work space */

public:

    nmrSVDRSSolver(void): M(0), N(0)
    { Allocate(M, N); }


    nmrSVDRSSolver( CISSTNETLIB_INTEGER m,
                    CISSTNETLIB_INTEGER n,
                    CISSTNETLIB_INTEGER nb = 1) {
        Allocate(m, n, nb);
    }


    nmrSVDRSSolver( vctDynamicMatrix<CISSTNETLIB_DOUBLE> &A,
                    vctDynamicMatrix<CISSTNETLIB_DOUBLE> &B )
    { Allocate(A, B); }


    inline void Allocate( CISSTNETLIB_INTEGER m,
                          CISSTNETLIB_INTEGER n,
                          CISSTNETLIB_INTEGER nb = 1) {
        M = m;
        N = n;
        Lda = std::max(M, N);
        Ldb = M;
        Nb = nb;
        S.SetSize(N, 1, VCT_COL_MAJOR);
        Work.SetSize(2 * N, 1, VCT_COL_MAJOR);
    }


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


    inline void Solve( vctDynamicMatrix<CISSTNETLIB_DOUBLE> &A,
                       vctDynamicMatrix<CISSTNETLIB_DOUBLE> &B )
        throw (std::runtime_error)
    {
        /* check that the size matches with Allocate() */
        if ((M  != static_cast<CISSTNETLIB_INTEGER>(A.rows())) ||
            (N  != static_cast<CISSTNETLIB_INTEGER>(A.cols())) ||
            (Nb != static_cast<CISSTNETLIB_INTEGER>(B.cols()))) {
            cmnThrow(std::runtime_error("nmrSVDRSSolver Solve: Sizes used for Allocate were different"));
        }

        /* check other dimensions */
        if (M != static_cast<CISSTNETLIB_INTEGER>(B.rows())) {
            cmnThrow(std::runtime_error("nmrSVDRSSolver Solve: Sizes of parameters are incompatible"));
        }

        /* check that the matrices are Fortran like */
        if (! (A.IsFortran() && B.IsFortran())) {
            cmnThrow(std::runtime_error("nmrSVDRSSolver Solve: All parameters must be Fortran compatible"));
        }

#if defined(CISSTNETLIB_VERSION_MAJOR)
#if (CISSTNETLIB_VERSION_MAJOR >= 3)
        cisstNetlib_svdrs_(A.Pointer(), &Lda, &M, &N, B.Pointer(),
                           &Ldb, &Nb, S.Pointer(), Work.Pointer());
#endif
#else // no major version
        svdrs_(A.Pointer(), &Lda, &M, &N, B.Pointer(),
               &Ldb, &Nb, S.Pointer(), Work.Pointer());
#endif // CISSTNETLIB_VERSION
    }


    inline const vctDynamicMatrix<CISSTNETLIB_DOUBLE> & GetS(void) const
    { return S; }

};


#endif // _nmrSVDRSSolver_h