mathtools.h

/***************************************************************************
                          mathtools.h  -  description
                             -------------------
    author               : Armando Navarro-Vazquez
    email                : armando.deus@gmail.com
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/

#ifndef MAGNESMATHTOOLS_H
#define MAGNESMATHTOOLS_H

#include <cstring>
#include <assert.h>
#include <string>
#include <iostream>
#include "exception.h"

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

#include "arrays_defines.h"

#include "toolsexport.h"

namespace magnes
{

template <class X,int N> class D1StaticArray
{
private:
    X m_data[N];
public:
    D1StaticArray()
    {
        assert(N>0);
    }
    D1StaticArray(const X& x)
    {
        for(size_t i=0;i<N;++i)
        {
            m_data[i]=x;
        }
    }
    const X& operator()(size_t i) const
    {
        assert(i<N);
        return m_data[i];
    }

    X& operator()(size_t i)
    {
        return m_data[i];
    }
    operator X*()
    {
        return &m_data[0];
    }
    operator const X*() const
    {
        return &m_data[0];
    }

    size_t size() const
    {
        return N;
    }

};
template <class X> class D1Array
{

private:
    X* m_data;
public:
    D1Array() : m_size ( 0 )
    {
        m_data=NULL;
    }
    explicit D1Array ( size_t i ) : m_size ( i )
    {
        m_data= new X[i];
    }
    D1Array ( size_t i,const X& x );
    D1Array ( const D1Array& x );
    void Initialize ( size_t i )
    {
        if ( m_data )
            delete [] m_data;
        m_data = new X[i];
        m_size=i;
    }
    void Initialize ( size_t i, const X& x );
    ~D1Array()
    {
        if ( m_data )
            delete [] m_data;
    }
    void Clear()
    {
        if ( m_data )
            delete [] m_data;
        m_data=NULL;
        m_size=0;
    }
    X& operator () ( size_t i );
    const X& operator() ( size_t i ) const;
    operator X*()
    {
        return m_data;
    }
    operator const X*() const
    {
        return m_data;
    }
    D1Array& operator= ( const D1Array& x );

    void operator*=(const X& x);
    void operator/=(const X& x);
    size_t size() const
    {
        return m_size;
    }

private:
    size_t m_size;

};


template<class X> D1Array<X>::D1Array ( size_t i, const X& x ) : m_size ( i )
{
    m_data = new X[i];
    for ( size_t k=0;k<i;k++ )
        m_data[k]=x;
}
template <class X> D1Array<X>::D1Array ( const D1Array<X>& x )
{

    m_size=x.size();
    if ( x.m_data !=NULL )
    {
        m_data= new X[m_size];

        memcpy ( m_data,x,m_size*sizeof ( X ) );
    }
    else m_data = NULL;

}

template <class X> D1Array<X>& D1Array<X>::operator= ( const D1Array<X>& x )
{
    if ( m_data )
        delete [] m_data;
    m_size=x.size();
    if ( x.m_data != NULL )
    {
        m_data= new X[m_size];
        memcpy ( m_data,x,m_size*sizeof ( X ) );
    }
    else m_data = NULL;
    return *this;

}
template<class X> void D1Array<X>::Initialize ( size_t i, const X& x )
{
    if ( m_data )
        delete [] m_data;
    m_size=i;
    m_data = new X[i];
    for ( size_t k=0;k<i;k++ )
        m_data[k]=x;
}

template <class X> X& D1Array<X>::operator () ( size_t i )
{
    assert ( i < m_size );
    return m_data[i];
}

template <class X> const X& D1Array<X>::operator () ( size_t i ) const
{
    assert ( i < m_size );
    return m_data[i];
}

template <class X> void D1Array<X>::operator*=(const X& x)
{
    for(size_t i=0;i<this->size();++i)
    {
        (*this)[i]*=x;
    }
}

template <class X> void D1Array<X>::operator/=(const X& x)
{
    for(size_t i=0;i<this->size();++i)
    {
        (*this)[i]/=x;
    }
}


template<class X> class D2Array
{
private:
    X* m_data;
public:
    D2Array()
    {
        m_rows=m_columns=0;
        m_data=NULL;
    }
    explicit D2Array ( size_t i, size_t j ) : m_rows ( i ), m_columns ( j )
    {
        m_data = new X[i*j];
    }
    D2Array ( size_t i,size_t j, const X& x );
    D2Array ( const D2Array& x );
    D2Array& operator= ( const D2Array& x );
    void operator+=(const D2Array& x);
    void operator-=(const D2Array& x);
    void operator*=(const X& x);
    void operator /=(const X& x);
    void Initialize ( size_t i, size_t j )
    {
        assert( i > 0 && j > 0 );
        m_rows=i;
        m_columns=j;
        if ( m_data )
            delete [] m_data;
        m_data= new X[i*j];
    }
    bool Empty() const
    {
        return ( m_data == NULL );
    }
    void Initialize ( size_t i, size_t j,const X& x );
    ~D2Array()
    {
        if ( m_data )
            delete [] m_data;
    }
    void Clear()
    {
        if ( m_data )
            delete [] m_data;
        m_data=NULL;
        m_rows=m_columns=0;
    }
    X& operator() ( size_t i, size_t j );
    const X& operator() ( size_t i, size_t j ) const;
    operator X*()
    {
        return m_data;
    }
    operator const X*() const
    {
        return m_data;
    }
    //arithmetic operators
    size_t NRows() const
    {
        return m_rows;
    }
    size_t NColumns() const
    {
        return m_columns;
    }
    void SwapRows ( size_t i, size_t j );
    void SwapColumns ( size_t i, size_t j );
    D2Array<X> RemoveRow ( size_t row );
    D2Array<X> RemoveColumn ( size_t column );
    D1Array<X> GetRow ( size_t row ) const;
    D1Array<X> GetColumn ( size_t column ) const;
    D2Array<X> SymmetricPart() const;
    D2Array<X> AntiSymmetricPart() const;
    void Transpose();
    X Trace() const;

private:
    size_t m_rows;
    size_t m_columns;
};



template<class X> D2Array<X>::D2Array ( size_t i, size_t j, const X& x ) : m_rows ( i ), m_columns ( j )
{
    size_t total=i*j;
    m_data = new X[total];
    for ( size_t k=0;k<total;++k )
        m_data[k]=x;
}

template<class X> void D2Array<X>::Initialize ( size_t i, size_t j, const X& x )
{
    if ( m_data )
        delete [] m_data;
    m_rows=i;
    m_columns=j;
    size_t total=i*j;
    m_data = new X[total];
    for ( size_t k=0;k<total;++k )
        m_data[k]=x;

}

template <class X> D2Array<X>::D2Array ( const D2Array<X>& x )
{
    m_rows=x.NRows();
    m_columns=x.NColumns();
    if ( x.m_data != NULL )
    {
        m_data= new X[m_rows*m_columns];
        memcpy ( m_data,x,m_rows*m_columns*sizeof ( X ) );
    }
    else m_data = NULL;
}

template <class X> D2Array<X>& D2Array<X>::operator= ( const D2Array& x )
{
    if ( m_data )
        delete [] m_data;
    m_rows=x.NRows();
    m_columns=x.NColumns();
    if ( x.m_data == NULL )
    {
        m_data=NULL;
    }
    else
    {
        m_data= new X[m_rows*m_columns];
        memcpy ( m_data,x,m_rows*m_columns*sizeof ( X ) );
    }
    return *this;
}

template <class X> void D2Array<X>::operator += ( const D2Array& x)
{
    assert ( this->m_rows == x.m_rows && this->m_columns == x.m_columns );
    for(size_t i=0;i<m_rows;++i)
        for(size_t j=0;j<m_columns;++j)
            (*this)(i,j)+=x(i,j);
}

template <class X> void D2Array<X>::operator -= ( const D2Array& x)
{
    assert ( this->m_rows == x.m_rows && this->m_columns == x.m_columns );
    for(size_t i=0;i<m_rows;++i)
        for(size_t j=0;j<m_columns;++j)
            (*this)(i,j)-=x(i,j);
}

template <class X> void D2Array<X>::operator *=(const X& x)
{
    for(size_t i=0;i<this->NRows();++i)
        for(size_t j=0;j<this->NColumns();++j)
        {
            (*this)(i,j)*=x;
        }
}

template <class X> void D2Array<X>::operator /=(const X& x)
{
    for(size_t i=0;i<this->NRows();++i)
        for(size_t j=0;j<this->NColumns();++j)
        {
            (*this)(i,j)/=x;
        }
}

template <class X> X& D2Array<X>::operator() ( size_t i, size_t j )
{
    return m_data[j*m_rows+i];
}

template <class X>
const X& D2Array<X>::operator() ( size_t i, size_t j ) const
{
    return m_data[j*m_rows+i];
}

template <class X> void D2Array<X>::SwapRows ( size_t i, size_t j )
{

    for ( size_t k=0;k<m_columns;++k )
    {
        X tmp=m_data[k*m_rows+j];
        m_data[k*m_rows+j]=m_data[k*m_rows+i];
        m_data[k*m_rows+i]=tmp;
    }
}

template <class X> void D2Array<X>::SwapColumns ( size_t i, size_t j )
{
    for ( size_t k=0;k<m_rows;++k )
    {
        X tmp= ( *this ) ( k,i );
        ( *this ) ( k,i ) = ( *this ) ( k,j );
        ( *this ) ( k,j ) =tmp;
    }
}


template <class X> D2Array<X> D2Array<X>::RemoveRow ( size_t row )
{
    D2Array<X> nmat ( m_rows-1,m_columns );
    size_t k=0;
    for ( size_t i=0;i<m_rows;i++ )
    {
        if ( i != row )
        {
            for ( size_t j=0;j<m_columns;j++ )
                nmat ( k,j ) = ( *this ) ( i,j );
            k++;
        }

    }

    return nmat;
}


template <class X> D2Array<X> D2Array<X>::RemoveColumn ( size_t column )
{
    D2Array<X> nmat ( m_rows,m_columns-1 );

    for ( size_t i=0;i<m_rows;i++ )
    {
        size_t k=0;
        for ( size_t j=0;j<m_columns;j++ )
        {
            if ( j!=column )
                nmat ( i,k++ ) = ( *this ) ( i,j );
        }
    }
    return nmat;
}

template <class X> D1Array<X> D2Array<X>::GetRow ( size_t row ) const
{
    D1Array<X> fila ( m_columns );
    for ( size_t i=0;i<m_columns;i++ )
        fila ( i ) = ( *this ) ( row,i );

    return fila;
}

template <class X> D1Array<X> D2Array<X>::GetColumn ( size_t column ) const
{
    D1Array<X> columna ( m_rows );
    for ( size_t i=0;i<m_rows;i++ )
        columna ( i ) = ( *this ) ( i,column );

    return columna;
}

template <class X> D2Array<X> D2Array<X>::SymmetricPart() const
{
    D2Array<X> t(*this);
    for(size_t i=0;i<this->NRows();++i)
        for(size_t j=i+1;j<this->NColumns();++j)
        {
            X mean=(t(i,j)+t(j,i))/2;
            t(i,j)=t(j,i)=mean;
        }
    return t;
}

template <class X> D2Array<X> D2Array<X>::AntiSymmetricPart() const
{
    D2Array<X> t(*this);
    for(size_t i=0;i<this->NRows();++i)
        for(size_t j=i;j<this->NColumns();++j)
        {
            if ( i == j )
                t(i,i)=0;
            else
            {
                X mean=(t(i,j)-t(j,i))/2;
                t(i,j)=mean;
                t(j,i)=-mean;
            }
        }
    return t;
}


template <class X> void D2Array<X>::Transpose()
{
    D2Array<X>& a=(*this);
    assert ( a.NRows() == a.NColumns() );
    size_t n=a.NRows();
    for(size_t i=0;i<n;++i)
    {
        for(size_t j=i+1;j<n;++j)
        {
            X tmp=a(i,j);
            a(i,j)=a(j,i);
            a(j,i)=tmp;
        }
    }
}

template <class X> X D2Array<X>::Trace() const
{
    assert ( this->NRows() == this->NColumns() );
    X trace=0;
    for(size_t i=0;i<this->NRows();++i)
    {
        trace+=(*this)(i,i);
    }
    return trace;
}

inline double operator * (const D1Array<double>& a, const D1Array<double>& b)
{
    assert ( a.size() == b.size() );
    double c=0;
    size_t n=a.size();
    for(size_t i=0;i<n;++i)
        c+=a[i]*b[i];
    return c;
}

inline D2Array<double> operator * ( const D2Array<double>& a,const D2Array<double>& b )
{
    assert ( a.NColumns() == b.NRows() );
    D2Array<double> c ( a.NRows(),b.NColumns() );

#ifdef WITH_MKL 
    int stride;
    if ( a.NRows() > a.NColumns() ) stride= a.NRows();
    else stride=a.NColumns();
    cblas_dgemm ( CblasColMajor,CblasNoTrans,CblasNoTrans,a.NRows(),b.NColumns(),a.NRows(),1.0,a,stride,b,b.NColumns(),0.0,c,stride );
#endif

#ifdef WITH_VECLIB
    long int stride;
    if ( a.NRows() > a.NColumns() ) stride= a.NRows();
    else stride=a.NColumns();
    cblas_dgemm ( CblasColMajor,CblasNoTrans,CblasNoTrans,a.NRows(),b.NColumns(),a.NRows(),1.0,a,stride,b,b.NColumns(),0.0,c,stride );
#endif
#ifdef WITH_LAPACK

    char transa='N';
    char transb='N';
    int m=a.NRows();
    int n=b.NColumns();
    int k=a.NColumns();
    double alpha=1.0;
    int lda=a.NRows();
    int ldb=k;
    int ldc=m;
    double beta=0;

    dgemm(&transa,&transb,&m,&n,&k,&alpha,const_cast<double*> ( static_cast<const double*> ( a )),&lda,const_cast<double*> ( static_cast<const double*> ( b )),&ldb,&beta,c,&ldc);
#endif

    return c;

}

inline D1Array<double> operator * ( const D2Array<double>& a,const D1Array<double>& b )
{
    assert ( a.NRows() == b.size() );
    D1Array<double> c ( a.NRows() );


#ifdef WITH_MKL 
    cblas_dgemv ( CblasColMajor,CblasNoTrans,a.NRows(),a.NColumns(),1.0,a,a.NRows(),b,1,0.0,c,1 );
#endif

#ifdef WITH_VECLIB
    cblas_dgemv ( CblasColMajor,CblasNoTrans,a.NRows(),a.NColumns(),1.0,a,a.NRows(),b,1,0.0,c,1 );
#endif
#ifdef WITH_LAPACK
    char trans='N';
    int nrows=a.NRows();
    int ncolumns=a.NColumns();
    double alpha=1.0;
    int lda=a.NRows();
    int incx=1;
    double beta=0.0;
    int incy=1;
    dgemv(&trans,&nrows,&ncolumns,&alpha,const_cast<double*> ( static_cast<const double*> ( a ) ),&lda,const_cast<double*> ( static_cast<const double*> ( b ) ),&incx,&beta,c,&incy);
#endif

    return c;
}
#if defined WITH_MKL || defined WITH_VECLIB || defined WITH_LAPACK

inline int Eigen ( D2Array<double>& a, D1Array<double>& e )
{
    char jobz='V';
    char uplo='U';
#ifdef WITH_VECLIB
    __CLPK_integer ldwork=a.NColumns() *a.NColumns();
    D1Array<double> work ( ldwork );
    __CLPK_integer info;
    __CLPK_integer n=a.NColumns();
    __CLPK_integer lda=a.NColumns();
    dsyev_ ( &jobz,&uplo,&n,a,&lda,e,work,&ldwork,&info );
#else
    int ldwork=a.NColumns() *a.NColumns();
    D1Array<double> work ( ldwork );
    int info;
    int n=a.NColumns();
    int lda=a.NColumns();
    dsyev ( &jobz,&uplo,&n,a,&lda,e,work,&ldwork,&info );
#endif
    return info;
}

inline void SVDFit(const D1Array<double>& res, const D2Array<double>& modelmatrix, D1Array<double>& solution)
{

#if defined WITH_MKL || defined WITH_LAPACK
    D2Array<double> A=modelmatrix;
    D1Array<double> b=res;
    int m=A.NRows();
    int n=A.NColumns();
    int nrhs=1;
    int lda=A.NRows();
    int ldb=res.size();
    D1Array<double> work ( 7*m );
    int ldwork=7*m;
    int info;
    D1Array<double> s ( A.NColumns() );
    int rank;
    double rcond=-1; //use machine precision
    dgelss ( &m,&n,&nrhs,A,&lda,b,&ldb,s,&rcond,&rank,work,&ldwork,&info );

    if ( info < 0 )
        throw magnes::Exception("Illegal entry to dgelss routine in SVDFit");
    if ( info > 0 )
        throw magnes::Exception("dgelss routine has not converged in SVDFit");

    for(size_t i=0;i<A.NColumns();++i)
        solution(i)=b(i);

    return;

#endif

#ifdef WITH_VECLIB
    D2Array<double> A=modelmatrix;
    D1Array<double> b=res;
    __CLPK_integer m=A.NRows();
    __CLPK_integer n=A.NColumns();
    __CLPK_integer nrhs=1;
    __CLPK_integer lda=A.NRows();
    __CLPK_integer ldb=b.size();
    D1Array<double> work ( 7*m );
    __CLPK_integer ldwork=7*m;
    __CLPK_integer info;
    D1Array<double> s ( A.NColumns() );
    __CLPK_integer rank;
    double rcond=-1; //use machine precision
    dgelss_ ( &m,&n,&nrhs,A,&lda,b,&ldb,s,&rcond,&rank,work,&ldwork,&info );
    if ( info < 0 )
        throw magnes::Exception("Illegal entry to dgelss routine in SVDFit");
    if ( info > 0 )
        throw magnes::Exception("dgelss routine has not converged in SVDFit");

    for(size_t i=0;i<A.NColumns();++i)
        solution(i)=b(i);

    return;
#endif
    throw magnes::Exception("SVD fit not implemented");
}

inline int Invert ( D2Array<double>& a )
{

#ifdef WITH_VECLIB
    __CLPK_integer m=a.NRows();
    __CLPK_integer n=a.NColumns();
    __CLPK_integer lda=m;
    __CLPK_integer info;
    __CLPK_integer* ipiv = new  __CLPK_integer[m];
    dgetrf_ ( &m,&n,a,&lda,ipiv,&info );
#else
    int m=a.NRows();
    int n=a.NColumns();
    int lda=m;
    int info;
    int* ipiv = new int[m];
    dgetrf ( &m,&n,a,&lda,ipiv,&info );
#endif
    //  std::cout << "info from dgetrf " << info << std::endl;
    D1Array<double> work ( n );
#ifdef WITH_VECLIB
    __CLPK_integer lwork=n;
    dgetri_ ( &n,a,&lda,ipiv,work,&lwork,&info );
#else
    int lwork=n;
    dgetri ( &n,a,&lda,ipiv,work,&lwork,&info );
#endif
    delete [] ipiv;
    //std::cout << "info from dgetri " << info << std::endl;
    return info;

}


#endif

template <class X> class  Scalar
{
public:
    Scalar() :  m_value(NULL) {}
    ~Scalar() { delete m_value; }
    X Value() const { return *m_value; }

    explicit Scalar(X v)
    {
        m_value =new X(v);
    }
    Scalar(const Scalar& v)
    {

        if (  v )
        {
            m_value = new X(v.Value());

        }
        else m_value=NULL;
        //   if (!m_value && !e)  //Do nothing


    }
    Scalar& operator = (const Scalar& e)
    {
        if ( &e != this)
        {
            if ( m_value && e )
            {
                *m_value=e.Value();

            }
            if ( m_value && !e)
            {
                delete m_value;
                m_value=NULL;
            }
            if ( !m_value && e)
            {
                m_value = new X(e.Value());

            }
            //   if (!m_value && !e)  //Do nothing

        }
        return *this;
    }


    Scalar& operator = (X v)
    {
        {
            if ( m_value ) *m_value=v;
            else
                m_value = new X(v);
        }
        return *this;
    }
private:
    typedef void (Scalar::*bool_type) () const;
    void helperfunction() const {}
private:
    //static unity m_unity;
    X* m_value;
public:
    operator bool_type () const
    {
        if ( m_value != 0 ) return &Scalar::helperfunction;
        else return 0;
    }
};

template <class X> class Tensor
{
public:
    Tensor() {}
    const D2Array<double>& LabFrame() const { return m_lftensor; }
    D2Array<double>& LabFrame() { return m_lftensor; }
    D2Array<double>& PrincipalFrame() const { return m_pftensor; }
    D2Array<double>& PrincipalFrame() { return m_pftensor; }
    D1Array<double>& EigenValues() { return m_eigenvalues; }
    D1Array<double>& EigenValues() const { return m_eigenvalues; }
    void Clear()
    {
        m_lftensor.Clear();
        m_pftensor.Clear();
        m_eigenvalues.Clear();
    }
    bool Empty() const { return m_lftensor.Empty() && m_pftensor.Empty(); }
private:
    D2Array<X> m_lftensor;
    D2Array<X> m_pftensor;
    D1Array<X> m_eigenvalues;

    void SetPrincipalFrame();
    //void SetLabFrame();

};
template <class X> void Tensor<X>::SetPrincipalFrame()
{
    Eigen(this->PrincipalFrame(),this->EigenValues());
    this->PrincipalFrame().Transpose();
}

inline std::string MathInfo()
{
    std::string library;
#ifdef WITH_LAPACK
    library="lapack";
    return library;
#endif

#ifdef WITH_VECLIB
    library="veclib macosx";
    return library;
#endif

#ifdef WITH_MKL
    library="Intel MKL";
    return;
#endif
    library="undefined";
    return library;
}

}

#endif //MAGNESMATHTOOLS_H