C++ LAPACK

#ifndef _FMATRIX_H
#define _FMATRIX_H

#include<iostream>

template <class T>
class Fmatrix
{
public:
    //default
    Fmatrix(): m(0), n(0), flag(0), mat(new T[0]){ this->settype(); }
    //blank
    Fmatrix(size_t _m): m(_m), n(1), flag(0), mat(new T[_m]){ this->settype(); }
    Fmatrix(size_t _m, size_t _n): m(_m), n(_n), flag(0), mat(new T[_m*_n])
    { this->settype(); }
    //copy
    Fmatrix(const Fmatrix<T>& B) :
    m(B.numrows()), n(B.numcols()), mat(new T[B.numelms()])
    {
        for(size_t i=0; i<B.numelms(); i++){ this->mat<i> = B<i>; }
        this->type=B.gettype();
    }

    ~Fmatrix(){ if(mat) delete[] mat; }

    inline T& operator() (size_t i, size_t j)
    { return this->mat[j*(this->m)+i]; }
    inline T  operator() (size_t i, size_t j) const
    { return this->mat[j*(this->m)+i]; }

    inline T& operator[] (size_t i)       { return mat<i>; }
    inline T  operator[] (size_t i) const { return mat<i>; }

    Fmatrix operator = (const Fmatrix<T>& B)
    {
        if(&B!=this){
            delete[] mat;
            this->m=B.numrows();
            this->n=B.numcols();
            this->flag=B.error();
            mat=new T[(this->m)*(this->n)];
            for(size_t i=0; i<B.numelms(); i++){ this->mat<i> = B<i>;}
        }
        return *this;
    }
    //these (+-*/) won't work for complex
    Fmatrix operator += (const Fmatrix<T>& B)
    {  for(size_t i=0; i< this->numelms(); i++){ this->mat<i>+=B<i>; }  }
    Fmatrix operator -= (const Fmatrix<T>& B)
    {  for(size_t i=0; i< this->numelms(); i++){ this->mat<i>-=B<i>; }  }
    Fmatrix operator *= (const T alpha)
    {  for(size_t i=0; i< this->numelms(); i++){ this->mat<i>*=alpha; }  }
    Fmatrix operator /= (const T beta)
    {  for(size_t i=0; i< this->numelms(); i++){ this->mat<i>/=beta; }  }

    inline T* begin(){ return mat; }

    inline const size_t numrows() const { return this->m; }
    inline const size_t numcols() const { return this->n; }
    inline const size_t numelms() const { return (this->m)*(this->n); }
    inline const bool   error()   const { return this->flag; }
    inline const char   gettype() const { return this->type; }
    
    void print(std::ostream &os=std::cout, char delim='\t')
    {
        for(size_t i=0;i<this->m;i++){for(size_t j=0;j<this->n;j++){
            os << (*this)(i,j) << delim;
        } os << std::endl; } os << std::endl;
    }
    
    void read(std::istream &is=std::cin)
    {
        for(size_t i=0;i<this->m;i++){for(size_t j=0;j<this->n;j++){
            is >> (*this)(i,j);
        } }
    }

private:
    size_t m,n;
    T *mat;
    char type;
    bool flag;

    inline const char settype()
    {
        char tstr[20], otyp[]="czu";
        strcpy(tstr,typeid(T).name());
        if(tstr[0]=='f'||tstr[0]=='d'||tstr[0]=='e') type=tstr[0];
        else if( strcmp(tstr,"7complex")==0 )        type=otyp[0];
        else if( strcmp(tstr,"13doublecomplex")==0 ) type=otyp[1];
        else                                         type=otyp[2];
    }
};

#endif /* _FMATRIX_H */

#ifndef _LINALG_H
#define _LINALG_H

#include <cstdlib>
#include <cstdio>
#include <cmath>
#include <iostream>
#include <fstream>

extern "C" {
    #include "f2c.h"
    #include "blaswrap.h"
    #include "clapack.h"
}

#include "Fmatrix.h"

//these scalar ops (+-*/) won't work for complex
template <class T>
const Fmatrix<T> operator + (const Fmatrix<T>& A, const Fmatrix<T>& B)
{  return Fmatrix<T>(A)+=B;  }
template <class T>
const Fmatrix<T> operator - (const Fmatrix<T>& A, const Fmatrix<T>& B)
{  return Fmatrix<T>(A)+=B;  }
template <class T>
const Fmatrix<T> operator * (const T alpha, const Fmatrix<T>& A)
{  return Fmatrix<T>(A)*=alpha;  }
template <class T>
const Fmatrix<T> operator / (const Fmatrix<T>& A, const T beta)
{  return Fmatrix<T>(A)/=beta;  }
template <class T>
const Fvector<T> operator * ( Fmatrix<T>& A, Fvector<T>& x)
{
    Fvector<T> y(A.numrows());
    gmv(A,x,y);
    return y;
}
template <class T>
const Fmatrix<T> operator * ( Fmatrix<T>& A, Fmatrix<T>& B)
{
    Fmatrix<T> C(A.numrows(),B.numcols());
    gmm(A,B,C);
    return C;
}

template <class T>
int gmv(Fmatrix<T>& A, Fvector<T>& x, Fvector<T>& y, char TRANSA='N',
            T alpha=1, T beta=0, integer INCX=1, integer INCY=1)
{
    integer M, N, LDA;

    M=(TRANSA=='N') ? A.numrows():A.numcols();
    N=(TRANSA=='N') ? A.numcols():A.numrows();
    LDA=M;

    //this still needs a switch on data type
    dgemv_(&TRANSA, &M, &N,
                        (doublereal*)&alpha,
                        (doublereal*)A.begin(), &LDA,
                        (doublereal*)x.begin(), &INCX,
                        (doublereal*)&beta,
                        (doublereal*)y.begin(), &INCY);


}


template <class T>
int gmm(Fmatrix<T>& A, Fmatrix<T>& B, Fmatrix<T>& C,
            char TRANSA='N', char TRANSB='N', T alpha=1, T beta=0)
{
    integer M, N, K, LDA, LDB, LDC;

    M=(TRANSA=='N') ? A.numrows():A.numcols();
    N=(TRANSB=='N') ? B.numcols():B.numrows();
    K=(TRANSA=='N') ? A.numcols():A.numrows();
    LDA=(TRANSA=='N') ? M:K;
    LDB=(TRANSB=='N') ? K:N;
    LDC=C.numrows();

    switch(A.gettype()){//cast explicitly otherwise compiler is cranky
        case 'f':
            sgemm_(&TRANSA, &TRANSB, &M, &N, &K,
                        (real*)&alpha,
                        (real*)A.begin(), &LDA,
                        (real*)B.begin(), &LDB,
                        (real*)&beta,
                        (real*)C.begin(), &LDC);
        break;
        case 'd':
            dgemm_(&TRANSA, &TRANSB, &M, &N, &K,
                        (doublereal*)&alpha,
                        (doublereal*)A.begin(), &LDA,
                        (doublereal*)B.begin(), &LDB,
                        (doublereal*)&beta,
                        (doublereal*)C.begin(), &LDC);
        break;
        case 'c':
            cgemm_(&TRANSA, &TRANSB, &M, &N, &K,
                        (complex*)&alpha,
                        (complex*)A.begin(), &LDA,
                        (complex*)B.begin(), &LDB,
                        (complex*)&beta,
                        (complex*)C.begin(), &LDC);
        break;
        case 'z':
            zgemm_(&TRANSA, &TRANSB, &M, &N, &K,
                        (doublecomplex*)&alpha,
                        (doublecomplex*)A.begin(), &LDA,
                        (doublecomplex*)B.begin(), &LDB,
                        (doublecomplex*)&beta,
                        (doublecomplex*)C.begin(), &LDC);
        break;
        default:
            std::cout << "\nERROR: Something wrong in GMM switch.\n\n";
    }
    return 0;
}

template <class T>
int pca(Fmatrix<T>& A, Fmatrix<T>& U, Fmatrix<T>& S, Fmatrix<T>& PC)
{
   std::ofstream PCoutx("pcx"), PCouty("pcy"), PCoutz("pcz"), PCoutw("pcw");
        PCoutx.precision(16); PCouty.precision(16);
        PCoutz.precision(16); PCoutw.precision(16);
        Fmatrix<T> Y(A.numcols(),A.numrows());
        T s,*mean;

        mean = new T[A.numrows()];
        for(int i=0; i<A.numrows(); i++) mean<i>=0;

        //find means
        for(int i=0; i<A.numrows(); i++){ for(int j=0; j<A.numcols(); j++){
            mean<i>+=A(i,j);
        } }
        for(int i=0; i<A.numrows(); i++) mean<i>/=A.numcols();

        //subtract means from data
        for(int i=0; i<A.numrows(); i++){ for(int j=0; j<A.numcols(); j++){
            A(i,j)-=mean<i>;
        } }

        //form matrix for SVD
        s=sqrt(A.numcols()-1);
        for(int i=0; i<A.numrows(); i++){ for(int j=0; j<A.numcols(); j++){
            Y(j,i)=A(i,j)/s;
        } }

        return svd(Y,U,S,PC);
}

template <class T>
int svd(Fmatrix<T>& A, Fmatrix<T>& U, Fmatrix<T>& S, Fmatrix<T>& V,
            char JOBU='A', char JOBVT='A')
{
    float  *wkf,*rwkc;
    double *wkd,*rwkz;
    complex *wkc;
    doublecomplex *wkz;
    integer M, N, mn, MN, LDA, LDU, LDVT, LWORK, INFO;

    M=A.numrows(); N=A.numcols(); LDA=M; LDU=M; LDVT=N;
    mn=min(M,N);
    MN=max(M,N);
    LWORK=2*max(3*mn+MN,5*mn);

    switch(A.gettype()){//cast explicitly otherwise compiler is cranky
        case 'f':
            wkf = new real[LWORK];
            sgesvd_(&JOBU, &JOBVT, &M, &N,
                        (real*)A.begin(), &LDA,
                        (real*)S.begin(),
                        (real*)U.begin(), &LDU,
                        (real*)V.begin(), &LDVT,
                        wkf, &LWORK, &INFO);
            delete[] wkf;
        break;
        case 'd':
            wkd = new doublereal[LWORK];
            dgesvd_(&JOBU, &JOBVT, &M, &N,
                        (doublereal*)A.begin(), &LDA,
                        (doublereal*)S.begin(),
                        (doublereal*)U.begin(), &LDU,
                        (doublereal*)V.begin(), &LDVT,
                        wkd, &LWORK, &INFO);
            delete[] wkd;
        break;
        case 'c':
            wkc = new complex[LWORK];
            rwkc = new real[5*mn];
            cgesvd_(&JOBU, &JOBVT, &M, &N,
                        (complex*)A.begin(), &LDA,
                        (real*)S.begin(),
                        (complex*)U.begin(), &LDU,
                        (complex*)V.begin(), &LDVT,
                        wkc, &LWORK, rwkc, &INFO);
            delete[] wkc; delete[] rwkc;
        break;
        case 'z':
            wkz = new doublecomplex[LWORK];
            rwkz = new doublereal[5*mn];
            zgesvd_(&JOBU, &JOBVT, &M, &N,
                        (doublecomplex*)A.begin(), &LDA,
                        (doublereal*)S.begin(),
                        (doublecomplex*)U.begin(), &LDU,
                        (doublecomplex*)V.begin(), &LDVT,
                        wkz, &LWORK, rwkz, &INFO);
            delete[] wkz; delete[] rwkz;
        break;
        default:
            std::cout << "\nERROR: Something wrong in SVD switch.\n\n";
    }

    return INFO;
}

#endif /* _LINALG_H */

#include <cstdlib>
#include <cstdio>
#include <cmath>
#include <iostream>

#include "Flinalg.h"

using namespace std;

int main(int argc, char** argv)
{
    cout.precision(2);
    size_t M,N;

    cout << "input numrows: ";
    cin >> M;
    cout << "input numcols: ";
    cin >> N;
    Fmatrix<double> A(M,N); //instantiate MxN matrix
    Fmatrix<double> v(N),u(M); //instantiate Mx1 matrix i.e. vector
    cout << "A is " <<  A.numrows() << 'x' << A.numcols() << ", numelms: " << A.numelms() << endl;
    cout << "v is " <<  v.numrows() << 'x' << v.numcols() << ", numelms: " << v.numelms() << endl;
    cout << "u is " <<  u.numrows() << 'x' << u.numcols() << ", numelms: " << u.numelms() << endl << endl;

    cout << "Input data for matrix A by rows (" << A.numelms() << " numbers):" << endl;
    A.read();
    cout << endl << "A =" << endl;
    A.print();
    
    //Matrix-vector multiply by operator overload
    cout << endl << "Setting v to ones vector..." << endl;
    for(int i=0;i<v.numelms();i++)
        v<i>=1;
    cout << "A*v =" << endl;
    u=A*v;
    u.print();
    
    //generic matrix-matrix multiply
    Fmatrix<double> C(A.numrows(),A.numrows());
    gmm(A,A,C,'N','T');
    cout << "A*A' =" << endl; 
    C.print();
    
    cout << "computing svd..." << endl;
    Fmatrix<double> U(M,M),S(min(M,N)),V(N,N); //min is defined in clapack.h
    cout << "INFO: " << svd(A,U,S,V) << endl; //do svd, output error code
    cout << "singular values:" << endl;
    for (size_t i=0; i < min(M,N); i++)
    { cout << "S[ " << i+1 << " ]\t= " << S<i> << endl; } cout << endl;
    cout << "singular vectors:" << endl;     //rows of V so print transpose
    for(int i=0; i<V.numrows(); i++){ cout << "v" << i+1 <<'\t'; } cout << endl;
    for(int i=0; i<V.numrows(); i++){ for(size_t j=0;j<V.numcols();j++){
        cout << V(j,i) << '\t';
    } cout << endl; }

    return 0;
}