mkl_matrix.hpp

#ifndef MKL_MATRIX_HEADER_HPP
#define MKL_MATRIX_HEADER_HPP

#include "matrix.hpp"
#include "shared_ptr.hpp"

class mkl_Matrix: public calculation_Matrix{
public:
    typedef STD_TR1::shared_ptr<mkl_Matrix> ref;
    typedef STD_TR1::shared_ptr<const mkl_Matrix> const_ref;

    mkl_Matrix():tf(t_n){};
    ~mkl_Matrix() throw() {};

// use with care, only to interact with legacy C and fortran code
// that expects pointer to matrix area
    const complex* raw_storage() const {return get_storage();};
    complex* raw_storage() {return get_storage();};

    const complex* raw_storage(const row r, const column c) const  { return get_storage(r,c);};
    complex* raw_storage(const row r, const column c)  { return get_storage(r,c);};;

    enum trans {t_n,t_c,t_t};
    void set_mkl_t_flag(trans t){tf = t;};

    void flip_mkl_t_flag();

    const char* t_flag() const;

// leading dimension of the Matrix, usually (tf == t_n) we store
// increasing column indizes (constant row) adjacent in memory, and ld is num_of_rows
    int ld() const {return ld_();};

private:
    static const char* c_t_n;
    static const char* c_t_c;
    static const char* c_t_t;
    
    trans tf;

    size_t num_rows() const { if (tf == t_n) return n_rows(); return n_cols();}
    size_t num_cols() const { if (tf == t_n) return n_cols(); return n_rows();}

    void scale(complex v);

    void apply_two(const complex* phi, 
           complex* result) const;

    virtual size_t n_rows() const = 0;
    virtual size_t n_cols() const = 0;

    virtual const complex* get_storage() const  = 0;
    virtual complex* get_storage() = 0;

    virtual const complex* get_storage(const row r, const column c) const  = 0;
    virtual complex* get_storage(const row r, const column c) = 0;

    virtual int ld_() const{ 
    if (tf == t_n) return num_of_rows(); 
        return num_of_cols();
    };

};

class sub_mkl_Matrix: public mkl_Matrix{
public:
    sub_mkl_Matrix(mkl_Matrix& t, row r, size_t rowsize, 
           column c, size_t colsize): myM(&t),
                          start(t.raw_storage(r,c)),
                          rows(rowsize),
                          cols(colsize), 
                          my_c(c),
                          my_r(r){};
    
    virtual ~sub_mkl_Matrix() throw() {};
    
private:
    mkl_Matrix* myM;

    complex* start;

    size_t rows;
    size_t cols;
    column my_c;
    row my_r;

    size_t n_rows() const {return rows;};
    size_t n_cols() const {return cols;};

    const complex* get_storage() const  { return start;};
    complex* get_storage() { return start;};

    const complex* get_storage(const row r, const column c) const { 
    return myM->raw_storage(r+my_r,c+my_c);};
    complex* get_storage(const row r, const column c){
    return myM->raw_storage(r+my_r,c+my_c);};

    void insert(row r, column c, complex v){
    myM->Insert(r + my_r,c + my_c,v);
    };  

    void add(row r, column c, complex v){
    myM->Add(r + my_r,c + my_c,v);
    };  

    void scale(complex v);

    complex get(row r, column c) const{
    return myM->Get(r + my_r,c + my_c);
    };

    int ld_() const { return myM->ld();};
};

struct Block_Matrix {
    int n1;
    int n2;

    sub_mkl_Matrix b22;
    sub_mkl_Matrix b11;
    sub_mkl_Matrix b12;
    sub_mkl_Matrix b21;

    Block_Matrix(mkl_Matrix& m, int n_1, int n_2) : n1(n_1), n2(n_2),
    b22(m, row(n1), n2, column(n1), n2),
    b11(m, row(0), n1, column(0), n1),
    b12(m, row(0), n1, column(n1), n2),
    b21(m, row(n1), n2, column(0), n1) {
    };
};

class wrap_pointer_Matrix: public mkl_Matrix{
public:

    wrap_pointer_Matrix(complex* p, row r, column  c) 
    : data(p), nc(c.get()), nr(r.get()){};

    wrap_pointer_Matrix(complex* p, row r) 
    : data(p), nc(r.get()), nr(r.get()){};

    virtual ~wrap_pointer_Matrix() throw() {};

private:
    complex* data;
    int nc;
    int nr;

    
    void insert(row r, column c, complex v){
    data[index(r,c)] = v;
    };  

    void add(row r, column c, complex v){
    data[index(r,c)] += v;
    };  

    complex get(row r, column c) const { 
    return data[index(r,c)];
    };

    const complex* get_storage() const { return data;};
    complex* get_storage(){ return data;};

    const complex* get_storage(const row r, 
                   const column c) const { return &data[index(r,c)];};
    complex* get_storage(const row r, 
             const column c){ return &data[index(r,c)];};

    size_t index(row r, column c) const { return r.get() + c.get()*nr;};

    size_t n_rows() const {return nr;};
    size_t n_cols() const {return nc;};
};

class hamilton_Matrix: public mkl_Matrix{
public:

    explicit hamilton_Matrix(int n) : data(n*n,0), N(n){};
    ~hamilton_Matrix() throw() {};
    std::vector<complex> get_diag() const; 
    complex* get_diag(complex*) const; 

private:

    std::vector<complex> data;
    int N;

    void insert(row r, column c, complex v){
    data[index(r,c)] = v;
    };  

    void add(row r, column c, complex v){
    data[index(r,c)] += v;
    };  

    complex get(row r, column c) const { 
    return data[index(r,c)];
    };

    const complex* get_storage() const { return &data[0];};
    complex* get_storage(){ return &data[0];};

    const complex* get_storage(const row r, 
                   const column c) const { return &data[index(r,c)];};
    complex* get_storage(const row r, 
             const column c){ return &data[index(r,c)];};


    size_t index(row r, column c) const { return r.get() + c.get()*N;};
    
    size_t n_rows() const {return N;};
    size_t n_cols() const {return N;};
};


class rectangular_Matrix: public mkl_Matrix{
public:

    rectangular_Matrix(const row r, const column  c)
    :  nc(c.get()), nr(r.get()), data(nc*nr){};

    ~rectangular_Matrix() throw() {};
    
    
    void resize(const row r, const column c); 
        
private:
    size_t n_rows() const {return nr;};
    size_t n_cols() const {return nc;};

    int nc;
    int nr;
    std::vector<complex> data;

    size_t index(row r, column c) const { return r.get() + c.get()*nr;};

    void insert(row r, column c, complex v){
    data[index(r,c)] = v;
    };  

    void add(row r, column c, complex v){
    data[index(r,c)] += v;
    };  

    complex get(row r, column c) const { 
    return data[index(r,c)];
    };

    const complex* get_storage() const{ return &data[0];};
    complex* get_storage() { return &data[0];};
    const complex* get_storage(const row r, 
                   const column c) const{ return &data[index(r,c)];};
    complex* get_storage(const row r, 
             const column c) { return &data[index(r,c)];};

};

#endif