nmrSVDSolver.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 _nmrSVDSolver_h
#define _nmrSVDSolver_h


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

class nmrSVDSolver {
    // 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 Lda;
    CISSTNETLIB_INTEGER Ldu;
    CISSTNETLIB_INTEGER Ldvt;
    CISSTNETLIB_INTEGER Lwork;
    char Jobu;
    char Jobvt;
    vctDynamicMatrix<CISSTNETLIB_DOUBLE> S;
    vctDynamicMatrix<CISSTNETLIB_DOUBLE> U;
    vctDynamicMatrix<CISSTNETLIB_DOUBLE> Vt;
    vctDynamicMatrix<CISSTNETLIB_DOUBLE> Work;
    CISSTNETLIB_INTEGER Info;
    bool StorageOrder;

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


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

    
    template <class _matrixOwnerType>
    nmrSVDSolver(const vctDynamicMatrixBase<_matrixOwnerType, 
                                            CISSTNETLIB_DOUBLE> &A) 
    { Allocate(A); }
    
    template <vct::size_type _rows, 
              vct::size_type _cols, 
              bool _storageOrder>
    nmrSVDSolver(const vctFixedSizeMatrix<CISSTNETLIB_DOUBLE, 
                                          _rows, 
                                          _cols, 
                                         _storageOrder>& A)
    { Allocate(A); }


    inline void Allocate( CISSTNETLIB_INTEGER m, 
                          CISSTNETLIB_INTEGER n, 
                          bool storageOrder) {
        const CISSTNETLIB_INTEGER one = 1;
        StorageOrder = storageOrder;
        if (storageOrder == VCT_COL_MAJOR) {
            M = m;
            N = n;
        } else {
            M = n;
            N = m;
        }
        Lda = std::max(one, M);
        Ldu = M;
        Ldvt = N;
        Lwork = std::max(3 * std::min(M, N) + std::max(M, N),
                         5 * std::min(M, N));
        Jobu = 'A';
        Jobvt = 'A';
        S.SetSize(std::min(M, N), 1, StorageOrder);
        U.SetSize(Ldu, Ldu, StorageOrder);
        Vt.SetSize(Ldvt, Ldvt, StorageOrder);
        Work.SetSize(Lwork, 1, StorageOrder);
    }


    template <class _matrixOwnerType>
    inline void 
    Allocate(const vctDynamicMatrixBase<_matrixOwnerType, 
                                        CISSTNETLIB_DOUBLE>& A)
    { Allocate(A.rows(), A.cols(), A.IsRowMajor()); }

    template <vct::size_type _rows, 
              vct::size_type _cols, 
              bool _storageOrder>
    inline void 
    Allocate(const vctFixedSizeMatrix<CISSTNETLIB_DOUBLE, 
                                      _rows, 
                                      _cols, 
                                      _storageOrder> & A)
    { Allocate(_rows, _cols, _storageOrder); }


    template <class _matrixOwnerType>
    inline void 
    Solve(vctDynamicMatrixBase<_matrixOwnerType, 
                               CISSTNETLIB_DOUBLE> &A) 
        throw (std::runtime_error){
        /* 
           check that the size and storage order matches with 
           Allocate() 
        */
        if (A.IsRowMajor() != StorageOrder) {
            cmnThrow(std::runtime_error("nmrSVDSolver 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("nmrSVDSolver Solve: Size used for Allocate was different"));
            }
        } 
        else if (A.IsRowMajor()) {
            if ((M != static_cast<CISSTNETLIB_INTEGER>(A.cols())) || 
                (N != static_cast<CISSTNETLIB_INTEGER>(A.rows()))) {
                cmnThrow(std::runtime_error("nmrSVDSolver Solve: Size used for Allocate was different"));
            }
        }
        
        /* check that the matrices are Fortran like */
        if (! A.IsCompact()) {
            cmnThrow(std::runtime_error("nmrSVDSolver Solve: Requires a compact matrix"));
        }

        dgesvd_(&Jobu, &Jobvt, &M, &N,
                A.Pointer(), &Lda, S.Pointer(),
                U.Pointer(), &Ldu,
                Vt.Pointer(), &Ldvt,
                Work.Pointer(), &Lwork, &Info );
    }
    

    template <vct::size_type _rows, 
              vct::size_type _cols, 
              bool _storageOrder>
    inline void 
    Solve(vctFixedSizeMatrix<CISSTNETLIB_DOUBLE, 
                             _rows, 
                             _cols, 
                             _storageOrder>& A) {
        vctDynamicMatrixRef<CISSTNETLIB_DOUBLE> Aref(A);
        Solve(Aref);
    }

    
    inline const vctDynamicMatrix<CISSTNETLIB_DOUBLE> &GetS(void) const
    { return S; }

    inline const vctDynamicMatrix<CISSTNETLIB_DOUBLE> &GetU(void) const
    { return (StorageOrder == VCT_COL_MAJOR) ? U : Vt; }

    inline const vctDynamicMatrix<CISSTNETLIB_DOUBLE> &GetVt(void)const
    { return (StorageOrder == VCT_COL_MAJOR) ? Vt : U; }
};


#ifdef CISST_COMPILER_IS_MSVC
class CISST_DEPRECATED nmrSVDSolver;
#endif // CISST_COMPILER_IS_MSVC


#endif // _nmrSVDSolver_h