nmrLSSolver.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-26

  (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 _nmrLSSolver_h
#define _nmrLSSolver_h


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

class nmrLSSolver {
    // we have this class so that we reserve memory only one would
    // help if svd of a same size matrix (or a matrix) that doesnt
    // change much is desired.

protected:
    CISSTNETLIB_INTEGER M;
    CISSTNETLIB_INTEGER N;
    CISSTNETLIB_INTEGER NRHS;
    CISSTNETLIB_INTEGER Lda;
    CISSTNETLIB_INTEGER Ldb;
    CISSTNETLIB_INTEGER Lwork;
    char Trans;
    vctDynamicMatrix<double> Work;
    CISSTNETLIB_INTEGER Info;
    bool StorageOrder;

public:
    nmrLSSolver(void):
        M(0),
        N(0),
        NRHS(0),
        StorageOrder(VCT_COL_MAJOR)
    {
        Allocate(M, N, NRHS, StorageOrder);
    }


    nmrLSSolver(CISSTNETLIB_INTEGER m, CISSTNETLIB_INTEGER n, CISSTNETLIB_INTEGER nrhs, bool storageOrder) {
        Allocate(m, n, nrhs, storageOrder);
    }


    nmrLSSolver(vctDynamicMatrix<double> &A, vctDynamicMatrix<double> &B) {
        Allocate(A, B);
    }


    inline void Allocate(CISSTNETLIB_INTEGER m, CISSTNETLIB_INTEGER n, CISSTNETLIB_INTEGER nrhs, bool storageOrder) {
        const CISSTNETLIB_INTEGER one = 1;
        StorageOrder = storageOrder;
        M = m;
        N = n;
        NRHS = nrhs;
        Lda = std::max(one, M);
        Ldb = std::max(one, std::max(M, N));
        CISSTNETLIB_INTEGER MN = std::min(M, N);
        Lwork = std::max (one, MN + std::max (MN, NRHS));
        Trans = 'N';
        Work.SetSize(Lwork, 1, StorageOrder);
    }


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


    template <class _matrixOwnerType>
    inline void Solve(vctDynamicMatrixBase<_matrixOwnerType, double> &A, vctDynamicMatrixBase<_matrixOwnerType, double> &B) throw (std::runtime_error) {
        /* check that the size and storage order matches with Allocate() */
        if (A.IsRowMajor() != StorageOrder) {
            cmnThrow(std::runtime_error("nmrLSSolver Solve: Storage order used for Allocate was different"));
        }
        if (B.IsRowMajor() != StorageOrder) {
            cmnThrow(std::runtime_error("nmrLSSolver Solve: Storage order used for Allocate was different"));
        }

        /* check sizes based on storage order, there is a more compact
           expression for this test but I find this easier to read and
           debug (Anton) */
        if (A.IsColMajor()) {
            if ((M != static_cast<CISSTNETLIB_INTEGER>(A.rows())) || (N != static_cast<CISSTNETLIB_INTEGER>(A.cols()))) {
                cmnThrow(std::runtime_error("nmrLSSolver Solve: Size used for Allocate was different"));
            }
        }
        if (B.IsColMajor()) {
            if ((M != static_cast<CISSTNETLIB_INTEGER>(B.rows())) || (NRHS != static_cast<CISSTNETLIB_INTEGER>(B.cols()))) {
                cmnThrow(std::runtime_error("nmrLSSolver Solve: Size used for Allocate was different"));
            }
        }

        /* check that the matrices are Fortran like */
        if (! A.IsCompact()) {
            cmnThrow(std::runtime_error("nmrLSSolver Solve: Requires a compact matrix"));
        }
        if (! B.IsCompact()) {
            cmnThrow(std::runtime_error("nmrLSSolver Solve: Requires a compact matrix"));
        }

        /* call the LAPACK C function */
#if defined(CISSTNETLIB_VERSION_MAJOR)
#if (CISSTNETLIB_VERSION_MAJOR >= 3)
        cisstNetlib_dgels_(&Trans, &M, &N, &NRHS,
                           A.Pointer(), &Lda,
                           B.Pointer(), &Ldb,
                           Work.Pointer(), &Lwork, &Info);
#endif
#else // no major version
        dgels_(&Trans, &M, &N, &NRHS,
               A.Pointer(), &Lda,
               B.Pointer(), &Ldb,
               Work.Pointer(), &Lwork, &Info);
#endif // CISSTNETLIB_VERSION
    }

};


#endif // _nmrLSSolver_h