nmrHFTISolver.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-2014 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 _nmrHFTISolver_h
#define _nmrHFTISolver_h

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

// forward declaration of fortran routine for very old cisstNetlib
#ifndef CISSTNETLIB_VERSION
#  define NEEDS_HFTI_FORWARD_DECLARATION
#else
#  if (CISSTNETLIB_VERSION_MAJOR == 3)
#    define NEEDS_HFTI_FORWARD_DECLARATION
#  endif
#endif


#ifdef NEEDS_HFTI_FORWARD_DECLARATION
CISSTNETLIB_INTEGER hfti_(CISSTNETLIB_DOUBLE * a, CISSTNETLIB_INTEGER * mda, CISSTNETLIB_INTEGER * m, CISSTNETLIB_INTEGER * n,
                          CISSTNETLIB_DOUBLE * b, CISSTNETLIB_INTEGER * mdb, CISSTNETLIB_INTEGER * nb, CISSTNETLIB_DOUBLE * tau,
                          CISSTNETLIB_INTEGER * krank, CISSTNETLIB_DOUBLE * rnorm, CISSTNETLIB_DOUBLE * h__, CISSTNETLIB_DOUBLE * g,
                          CISSTNETLIB_INTEGER * ip);
#  undef NEEDS_HFTI_FORWARD_DECLARATION
#endif

class nmrHFTISolver {

protected:
    CISSTNETLIB_INTEGER M;
    CISSTNETLIB_INTEGER N;
    CISSTNETLIB_INTEGER NB;
    vctDynamicMatrix<CISSTNETLIB_DOUBLE> A;
    vctDynamicMatrix<CISSTNETLIB_DOUBLE> B;
    CISSTNETLIB_DOUBLE tau;
    CISSTNETLIB_INTEGER krank;
    vctDynamicMatrix<CISSTNETLIB_DOUBLE> RNORM;
    vctDynamicMatrix<CISSTNETLIB_DOUBLE> G;
    vctDynamicMatrix<CISSTNETLIB_DOUBLE> H;
    vctDynamicMatrix<CISSTNETLIB_INTEGER> IP;
    vctDynamicMatrix<CISSTNETLIB_DOUBLE> X;

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


    nmrHFTISolver(CISSTNETLIB_INTEGER m, CISSTNETLIB_INTEGER n, CISSTNETLIB_INTEGER nb) {
        Allocate(m, n, nb);
    }


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


    inline void Allocate(CISSTNETLIB_INTEGER m, CISSTNETLIB_INTEGER n, CISSTNETLIB_INTEGER nb) {
        M = m;
        N = n;
        NB = nb;
        vctDynamicMatrix<CISSTNETLIB_DOUBLE> RNORM;
        vctDynamicMatrix<CISSTNETLIB_DOUBLE> G;
        vctDynamicMatrix<CISSTNETLIB_DOUBLE> H;
        vctDynamicMatrix<CISSTNETLIB_INTEGER> IP;
        A.SetSize(M, N, VCT_COL_MAJOR);
        B.SetSize(M, NB, VCT_COL_MAJOR);
        RNORM.SetSize(NB, 1, VCT_COL_MAJOR);
        G.SetSize(N, 1, VCT_COL_MAJOR);
        H.SetSize(N, 1, VCT_COL_MAJOR);
        IP.SetSize(N, 1, VCT_COL_MAJOR);
        X.SetSize(M, NB, 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("nmrHFTISolver Solve: Sizes used for Allocate were different"));
        }

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

        hfti_(A.Pointer(), &M, &M, &N, B.Pointer(), &M, &NB, &tau,
              &krank, RNORM.Pointer(),
              H.Pointer(), G.Pointer(), IP.Pointer());
        //error handling??
    }
};


#endif // _nmrHFTISolver_h