nmrNNLSSolver.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 _nmrNNLSSolver_h
#define _nmrNNLSSolver_h

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

class nmrNNLSSolver {

protected:
    CISSTNETLIB_INTEGER M;
    CISSTNETLIB_INTEGER N;
    CISSTNETLIB_INTEGER Mda;
    CISSTNETLIB_INTEGER Mode;
    CISSTNETLIB_DOUBLE RNorm;
    vctDynamicMatrix<CISSTNETLIB_DOUBLE> X;
    vctDynamicMatrix<CISSTNETLIB_DOUBLE> W;
    vctDynamicMatrix<CISSTNETLIB_DOUBLE> Zz;
    vctDynamicMatrix<CISSTNETLIB_INTEGER> Index;

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


    nmrNNLSSolver(CISSTNETLIB_INTEGER m, CISSTNETLIB_INTEGER n) {
        Allocate(m, n);
    }


    nmrNNLSSolver(vctDynamicMatrix<CISSTNETLIB_DOUBLE> &C, vctDynamicMatrix<CISSTNETLIB_DOUBLE> &d) {
        Allocate(C, d);
    }


    inline void Allocate(CISSTNETLIB_INTEGER m, CISSTNETLIB_INTEGER n) {
        M = m;
        N = n;
        Mda = M;
        X.SetSize(N, 1, VCT_COL_MAJOR);
        W.SetSize(N,1, VCT_COL_MAJOR);
        Zz.SetSize(M,1, VCT_COL_MAJOR);
        Index.SetSize(N, 1, VCT_COL_MAJOR);
    }


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


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

        /* check other dimensions */
        if (C.rows() != d.rows()) {
            cmnThrow(std::runtime_error("nmrNNLSSolver Solve: Sizes of parameters are incompatible"));
        }

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

#if defined(CISSTNETLIB_VERSION_MAJOR)
#if (CISSTNETLIB_VERSION_MAJOR >= 3)
        cisstNetlib_nnls_(C.Pointer(), &Mda, &M, &N, d.Pointer(), X.Pointer(), &RNorm,
              W.Pointer(), Zz.Pointer(), Index.Pointer(), &Mode);
#endif
#else // no major version
        nnls_(C.Pointer(), &Mda, &M, &N, d.Pointer(), X.Pointer(), &RNorm,
              W.Pointer(), Zz.Pointer(), Index.Pointer(), &Mode);
#endif // CISSTNETLIB_VERSION

        //error handling??
    }


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


    /* Get W.  This method must be used after Solve(). */
    inline const vctDynamicMatrix<CISSTNETLIB_DOUBLE> &GetDual(void) const {
        return W;
    }


    /* Get RNorm.  This method must be used after Solve(). */
    inline CISSTNETLIB_DOUBLE GetRNorm(void) const {
        return RNorm;
    }
};


#endif // _nmrNNLSSolver_h