/dbfs1/fh22/linux/src/FMMLIB/tools.h

Go to the documentation of this file.

#ifndef TOOLS_H
#define TOOLS_H



#include <iostream>
#include <fstream>
#include <math.h>
#include <complex>
#include <sstream>
#include <string>
#include <vector>
#include <limits>
#include<ctime>


#include <tbci/matrix.h>
#include <tbci/vector.h>


using std::cout;
using std::complex;
using std::string;

const complex<double> I=complex<double>(0,1.);

template < typename T >
inline std::string to_str(T t)
{
        std::ostringstream ss;
        ss<<t;
        return ss.str();
}



template < typename T >
inline string to_str(T *V, int dim1)
{
        string s="";
        for (int i=0;i<dim1;i++)
        {
                s+=to_str(V[i])+" ";
        }
        return s;
}

template < typename T >
inline string to_str(T **M, int dim1, int dim2=-1)
{
        string s="";
        if (dim2==-1) dim2=dim1;
        for (int i=0;i<dim1;i++)
        {
                for (int j=0;j<dim2;j++) s+=to_str(real(M[i][j]))+" ";
                s+="\n";
        }
        return s;
}


template < typename T >
T str2number(std::string s)
{
   std::stringstream ss(s.c_str());
   T result;
   if (!(ss >> result))
      throw std::exception();
   return result;
} 

inline std::vector<string> strsplit(string str, string delim=" ")
{
        std::vector<string> results;
        int cutAt;
        while( (cutAt = str.find_first_of(delim)) != signed(str.npos) )
        {
                if(cutAt > 0)
                {
                        results.push_back(str.substr(0,cutAt));
                }
                str = str.substr(cutAt+1);
        }
        if((str.length() > 0))
        {
                results.push_back(str);
        }
        return results;
}

inline std::vector<string> readFile(string filename)
{
        std::vector<string> data(0);
        std::ifstream file;
        string line;
        //bool validMatrix=false;

        file.open(filename.c_str());
        if (file.is_open())
        {
                while (! file.eof() )
                {
                        getline(file,line);
                        if (line.length()>1)
                        {
                                if (line.at(0)!='#')
                                {
                                        data.push_back(line);
                                }
                        }
                }
                file.close();
                return data;
        }
        else 
        {
                cout << "Unable to open file"; 
                return data;
        }
}


inline int delta(int i,int j)
{
        return (i==j?1:0);
}


inline string timestamp(void){
        string timestamp;
        time_t ltime;
        tm *t;
        time(&ltime);
        t = gmtime(&ltime);
        std::stringstream h,min,s,d,mon,y;
        
        h<<std::setw(2)<<std::setfill('0') <<  t->tm_hour; 
        min<<std::setw(2)<<std::setfill('0')<<t->tm_min; 
        s<<std::setw(2)<<std::setfill('0')<<t->tm_sec; 
        d<<std::setw(2)<<std::setfill('0')<<t->tm_mday; 
        mon<<std::setw(2)<<std::setfill('0')<<t->tm_mon; 
        y<<(t->tm_year+1900);
        timestamp = y.str()+""+mon.str()+""+d.str()+"-"+h.str()+""+min.str()+""+s.str();
        return timestamp;
}


// Methods for TBCI Matrices and vectors


template <typename T>
inline TBCI::Matrix< T > readMatrix(string filename)
{
        TBCI::Matrix< T > M;
        std::vector<string> stringData(readFile(filename)), line;
        std::vector< std::vector<string> > data(0);
        bool validMatrix=true;

        for (int i=0;i<signed(stringData.size());i++)
        {
                data.push_back(strsplit(stringData[i]));
                if (i>0) if (data[i].size()!=data[i-1].size()) validMatrix=false;
        }
        if (validMatrix)
        {
                M.resize(data.size(), data[0].size());
                for (int i=0;i<signed(data.size());i++)
                {
                        for (int j=0;j<signed(data[i].size());j++) 
                        {
                                M(i,j)=str2number< T >(data[i][j]);
                        }
                }
                return M;
        }
        else 
        {
                cout<<"File does not contain a valid matrix\n";
                return M;
        }
}

template < typename T >
TBCI::Vector< T > readVector(string filename)
{
        std::vector<string> v(strsplit(readFile(filename)[0]));
        TBCI::Vector< T > result(0.,v.size());

        for (int i=0;i<signed(v.size());i++) 
        {
                result[i]=str2number<T>(v[i]);
        }
        return result;
}


template < typename T >
T min(TBCI::Matrix< T > M)
{
        T min=std::numeric_limits< T >::max();
        for (int i=0;i<M.rows();i++)
        {
                for (int j=0;j<M.columns();j++)
                {
                        if (M(i,j)<min) min=M(i,j);
                }
        }
        return min;
}
template < typename T >


T max(TBCI::Matrix< T > M)
{
        T max=std::numeric_limits< T >::min();
        for (int i=0;i<M.rows();i++)
        {
                for (int j=0;j<M.columns();j++)
                {
                        if (M(i,j)<max) max=M(i,j);
                }
        }
        return max;
}


template < typename T >
inline TBCI::Matrix< T > invert(TBCI::Matrix< T > M)
{
        //TBCI::Matrix< T > invM(M.rows(),M.columns());
        //TBCI::gaussj(invM, M);
        //return invM;
        return TBCI::Matrix<T>(inv(TBCI::F_Matrix<T>(M)));
}



template < typename T >
inline TBCI::Matrix< T > multFrontBack(TBCI::Matrix< T > M, TBCI::Vector< T > v)
{
        for (unsigned int row=0;row<M.rows();row++)
        {
                for (unsigned int column=0;column<M.columns();column++)
                {
                        M(row,column)*=v(row)*v(column);
                }
        }
        return M;
}

template < typename T >
inline TBCI::Matrix< T > multFront(TBCI::Matrix< T > M, TBCI::Vector< T > v)
{
        for (unsigned int row=0;row<M.rows();row++)
        {
                for (unsigned int column=0;column<M.columns();column++)
                {
                        M(row,column)*=v(column);
                }
        }
        return M;
}


template < typename T >
TBCI::Matrix<T> real(TBCI::Matrix< complex<T> > M)
{
        TBCI::Matrix< T > N(0.,M.rows(),M.columns());

        for (int i=0;i<signed(N.rows());i++) 
        {
                for (int j=0;j<signed(N.columns());j++) 
                {
                        N(i,j)=real(M(i,j));
                }
        }
        return N;
}
template < typename T >
TBCI::Matrix<T> imag(TBCI::Matrix< complex<T> > M)
{
        TBCI::Matrix< T > N(0.,M.rows(),M.columns());

        for (int i=0;i<signed(N.rows());i++) 
        {
                for (int j=0;j<signed(N.columns());j++) 
                {
                        N(i,j)=imag(M(i,j));
                }
        }
        return N;
}
template < typename T >
TBCI::Matrix<T> arg(TBCI::Matrix< complex<T> > M)
{
        TBCI::Matrix< T > N(0.,M.rows(),M.columns());

        for (int i=0;i<signed(N.rows());i++) 
        {
                for (int j=0;j<signed(N.columns());j++) 
                {
                        N(i,j)=arg(M(i,j));
                }
        }
        return N;
}
template < typename T >
TBCI::Matrix<T> abs(TBCI::Matrix< complex<T> > M)
{
        TBCI::Matrix< T > N(0.,M.rows(),M.columns());

        for (int i=0;i<signed(N.rows());i++) 
        {
                for (int j=0;j<signed(N.columns());j++) 
                {
                        N(i,j)=abs(M(i,j));
                }
        }
        return N;
}
template < typename T >
TBCI::Matrix<T> norm(TBCI::Matrix< complex<T> > M)
{
        TBCI::Matrix< T > N(0.,M.rows(),M.columns());

        for (int i=0;i<signed(N.rows());i++) 
        {
                for (int j=0;j<signed(N.columns());j++) 
                {
                        N(i,j)=norm(M(i,j));
                }
        }
        return N;
}





/************************************************************************************
**************** ALLOCATE / DELETE A 2D ARRAY ****************************************
************************************************************************************/

template < typename T >
T **Allocate2DArray(int nRows, int nCols)
{
        T **ptr2d;
        ptr2d = new T*[nRows];
        for (int i=0;i<nRows;i++)
        {
                ptr2d[i]=new T[nCols];
        }
        return ptr2d;
}


template < typename T >
void Free2DArray(T** Array, int nRows)
{
        for (int i=0;i<nRows;i++)
        {
                delete [] Array[i];
        }
        delete [] Array;
}

/************************************************************************************
**************** ALLOCATE / DELETE A 3D ARRAY ****************************************
************************************************************************************/

 template < typename T >
 T ***Allocate3DArray(int n1, int n2, int n3)
 {
        T ***ptr = new T**[n1];
        for (int i=0;i<n1;i++)
        {
                *(ptr+i) = new T*[n2];
                for (int j=0;j<n2;j++)
                {
                        *(*(ptr+i)+j)=new T[n3];
                }
        }
        return ptr;
 }
 
template < typename T >
void Free3DArray(T*** Array)
{
        delete [] **Array;
        delete [] *Array;
        delete [] Array;
}

/*************************************************************************
*********** CALCULATE INVERSE OF A COMPLEX MATRIX *****
*************************************************************************/
#ifndef CMPLX
#define CMPLX
struct cmplx{double re; double im;};
#endif

extern "C" void zgetri_(int *, cmplx *, int *, int *, cmplx *, int *, int *);
extern "C" void zgetrf_(int *, int*, cmplx *, int *, int *, int *);

template < typename T >
void inv_c(T **lhsMat, T **inv, int dim)
{
        int *IPIV;
        cmplx *lhsMat_cplx, *WORK;
        int INFO;
        
        lhsMat_cplx = new cmplx[dim*dim];//(cmplx *) malloc(sizeof(cmplx)*dim*dim);
        WORK = new cmplx[dim];//(cmplx *)malloc(sizeof(cmplx)*dim);
        IPIV = new int[dim];//(int *)malloc(sizeof(int)*dim);

        for (int i=0;i<dim;i++){
                for (int j=0;j<dim;j++){
                        lhsMat_cplx[i+j*dim].re=real(lhsMat[i][j]);
                        lhsMat_cplx[i+j*dim].im=imag(lhsMat[i][j]);
                }
        }
        // LU factorization (not banded)
        zgetrf_(&dim, &dim,lhsMat_cplx,&dim,IPIV,&INFO);
        zgetri_(&dim, lhsMat_cplx, &dim, IPIV, WORK, &dim, &INFO);      
        for (int i=0;i<dim;i++)
        {
                for (int j=0;j<dim;j++) inv[i][j]= T(lhsMat_cplx[i+j*dim].re,lhsMat_cplx[i+j*dim].im);
        }
        delete [] lhsMat_cplx;
        delete [] WORK;
        delete [] IPIV;
        
}
#endif

---