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-2018 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) CISST_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