nmrLSNonLinJacobianSolver.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-2007 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 _nmrLSNonLinJacobianSolver_h
#define _nmrLSNonLinJacobianSolver_h


#include <cisstVector/vctDynamicMatrix.h>

#include <cnetlib.h>

class nmrLSNonLinJacobianSolver {
    // 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 Ldfjac;
    CISSTNETLIB_DOUBLE Tolerance;
    CISSTNETLIB_INTEGER Info;
    CISSTNETLIB_INTEGER Lwork;
    vctDynamicVector<CISSTNETLIB_INTEGER> IWork;
    vctDynamicVector<CISSTNETLIB_DOUBLE> Work;

public:
    nmrLSNonLinJacobianSolver(void)
    {
        Allocate(0, 0);
    }


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

    
    nmrLSNonLinJacobianSolver(vctDynamicVector<CISSTNETLIB_DOUBLE> &X, vctDynamicVector<CISSTNETLIB_DOUBLE> &F) {
        Allocate(X, F);
    }


    inline void Allocate(CISSTNETLIB_INTEGER m, CISSTNETLIB_INTEGER n) {
        M = m;
        N = n;
    Lwork = M * N + 5 * N + M;
    Ldfjac = M;
    Tolerance = 1e-6;
        Work.SetSize(Lwork);
    IWork.SetSize(N);
    }


    inline void Allocate(vctDynamicVector<CISSTNETLIB_DOUBLE> &X, vctDynamicVector<CISSTNETLIB_DOUBLE> &F) {
        Allocate(X.size(), F.size());
    }


    template <int __instanceLineF, 
              class __elementTypeF, 
              int __instanceLineJ, 
              class __elementTypeJ>
    inline void Solve(nmrCallBackFunctionF<__instanceLineF, __elementTypeF> &callBackF,
                      nmrCallBackFunctionJ<__instanceLineJ, __elementTypeJ> &callBackJ,
                      vctDynamicVector<CISSTNETLIB_DOUBLE> &X,
                      vctDynamicVector<CISSTNETLIB_DOUBLE> &F, 
                      vctDynamicVector<CISSTNETLIB_DOUBLE> &J, 
                      CISSTNETLIB_DOUBLE tolerance) throw (std::runtime_error) {
    //if ((N != (CISSTNETLIB_INTEGER) X.size()) || (M != (CISSTNETLIB_INTEGER) F.size()) || (Ldfjac != (CISSTNETLIB_INTEGER) J.rows()) || (N != (CISSTNETLIB_INTEGER) J.cols())) {
        if ((N != static_cast<CISSTNETLIB_INTEGER>(X.size())) || (M != static_cast<CISSTNETLIB_INTEGER>(F.size())) || (Ldfjac*N != static_cast<CISSTNETLIB_INTEGER>(J.size()))) {
        cmnThrow(std::runtime_error("nmrLSNonLinJacobianSolver Solve: Size used for Allocate was different"));
    }
    Tolerance = tolerance;
        /* call the MINPACK C function */
    lmder2_(callBackF.FunctionFlmdif, callBackJ.FunctionFlmder2, &M, &N,
            X.Pointer(), F.Pointer(), J.Pointer(),
            &Ldfjac, &Tolerance, &Info,
            IWork.Pointer(), Work.Pointer(), &Lwork);
    }
    
};

#endif // _nmrLSNonLinJacobianSolver_h