matrix.hpp

#ifndef MATRIX_CPP_H
#define MATRIX_CPP_H

#include <cstdlib>
#include <exception>
#include <iostream>
#include <fstream>
#include <list>
#include <map>
#include <vector>
#include <memory>
#include <complex>

#include "exceptions.hpp"
#include "shared_ptr.hpp"

typedef std::complex<double> complex;
typedef std::vector<complex> Vector;

typedef STD_TR1::shared_ptr<Vector> Vector_ref;


template <int>
class selector{
    size_t cc;
public:
    selector():cc(0){};
    explicit selector(size_t i):cc(i){
    };
    
    size_t get() const {return cc;};
    void set(size_t c) {cc = c;};
    
    selector& operator++(){ 
    ++cc; return *this;};
    selector operator++(int){
    selector tmp(cc); ++cc; return tmp;};
    selector& operator--(){ 
    --cc; return *this;};
    selector operator--(int){
    selector tmp(cc); --cc; return tmp;};
    selector operator+(const selector c) const{ 
    return selector(cc + c.cc);};
    selector operator-(const selector c) const{ 
    return selector(cc - c.cc);};
    selector& operator+=(const int c){ 
    cc += c; return *this;};
    selector& operator-=(const int c){ 
    cc -= c; return *this;};
    bool operator==(const selector c) const { 
    return (cc == c.cc);};
    bool operator!=(const selector c) const { 
    return (cc != c.cc);};

    bool operator<(const selector c){
      return cc < c.cc;};
    bool operator>(const selector c){
      return cc > c.cc;};
    bool operator<=(const selector c){
      return cc <= c.cc;};
    bool operator>=(const selector c){
      return cc >= c.cc;};
};

typedef selector<0> row;
typedef selector<1> column;    



class Matrix{
public:
  virtual ~Matrix() throw() {};

  void assert_size(row r, column c);
  void assert_size(row r, column c) const;
  
  void Insert(row r, column c, complex v) { insert(r,c,v);};
  void Add(row r, column c, complex v) { add(r,c,v);};
  complex Get(row r, column c) const { return get(r,c);};    

  row nrows() const { return row(num_rows());};
  size_t num_of_rows() const { return num_rows();};

  column ncols() const { return column(num_cols());};
  size_t num_of_cols() const { return num_cols();};

    void add_diagonal( complex e );
    void insert_diagonal( complex e );
    void multiply_by_number( complex v);
    void print(std::ostream&) const;
    void formatted_print(std::ostream&) const;

    bool check_hermitian() const;
    bool check_realsymmetric() const;
    Matrix & operator +=(Matrix const & rhs);

private:    
    virtual void insert(row r, column c, complex v) = 0;
    virtual void add(row r, column c, complex v) = 0;
    virtual complex get(row r, column c) const = 0;    

    virtual void resize(row r, column c) {
    throw (E_Bad_Size(r.get()*c.get()));
    };

    virtual size_t num_rows() const = 0;
    virtual size_t num_cols() const = 0;
};

std::ostream& operator<<(std::ostream&, const Matrix&);



class sub_Matrix: public Matrix{
public:

    sub_Matrix(Matrix& t, row r, column c): 
    myM(&t), my_c(c), my_r(r) {};

    virtual ~sub_Matrix() throw() {};

private:
    Matrix* myM;
    column my_c;
    row my_r;

    sub_Matrix();

    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);
    };  

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

    size_t num_rows() const { return myM->num_of_rows() - my_r.get();};
    size_t num_cols() const { return myM->num_of_cols() - my_c.get();};
};

class row_iterator{
public:
    row_iterator():r(0),c(0){};
    row_iterator(Matrix* mm, row rr, column cc):r(rr),c(cc),m(mm){};

    row_iterator& operator++(){ ++r; return *this;};
    row_iterator operator++(int){ row_iterator ri(*this); ++r; return ri;};

    row_iterator& operator--(){ --r; return *this;};
    row_iterator operator--(int){ row_iterator ri(*this); --r; return ri;};

    bool operator==(const row_iterator& ri){ return ri.r == r;};
    bool operator!=(const row_iterator& ri){ return ri.r != r;};

    std::complex<double> operator*(){ return m->Get(r,c);}; 
  
private:
    row r;
    column c;
    Matrix* m;
};


class transpose_Matrix: public Matrix{
public:
    transpose_Matrix(Matrix* m):myM(m){};
    virtual ~transpose_Matrix() throw() {};

private:
    Matrix* myM;

    transpose_Matrix();

    void insert(row r, column c, complex v){
    myM->Insert(row(c.get()),column(r.get()),conj(v));
    };
    void add(row r, column c, complex v){
    myM->Add(row(c.get()),column(r.get()),conj(v));
    };
    complex get(row r, column c) const{
    return myM->Get(row(c.get()),column(r.get()));
    };

    size_t num_rows() const {return myM->num_of_cols();};
    size_t num_cols() const {return myM->num_of_rows();};
};


class correct_transpose_Matrix: public Matrix{
public:
    correct_transpose_Matrix(Matrix* m):myM(m){};
    virtual ~correct_transpose_Matrix() throw() {};

private:
    Matrix* myM;

    correct_transpose_Matrix();

    void insert(row r, column c, complex v){
    myM->Insert(row(c.get()),column(r.get()),v);
    };
    void add(row r, column c, complex v){
    myM->Add(row(c.get()),column(r.get()),v);
    };
    complex get(row r, column c) const{
    return myM->Get(row(c.get()),column(r.get()));
    };

    size_t num_rows() const {return myM->num_of_cols();};
    size_t num_cols() const {return myM->num_of_rows();};
};

class scaled_Matrix: public Matrix{
public:
    scaled_Matrix(complex s,Matrix* m):s(s),myM(m){};
    virtual ~scaled_Matrix() throw() {};

private:
    complex s;
    Matrix* myM;

    scaled_Matrix();

    void insert(row r, column c, complex v){
    myM->Insert(row(r.get()),column(c.get()),v*s);
    };
    void add(row r, column c, complex v){
    myM->Add(row(r.get()),column(c.get()),v*s);
    };
    complex get(row r, column c) const{
    return myM->Get(row(r.get()),column(c.get()));
    };

    size_t num_rows() const {return myM->num_of_rows();};
    size_t num_cols() const {return myM->num_of_cols();};
};


class conjugate_transpose_Matrix: public Matrix{
public:
    conjugate_transpose_Matrix(Matrix* m):myM(m){};
    virtual ~conjugate_transpose_Matrix() throw() {};

private:
    Matrix* myM;

    conjugate_transpose_Matrix();

    void insert(row r, column c, complex v){
    myM->Insert(row(c.get()),column(r.get()),conj(v));
    };
    void add(row r, column c, complex v){
    myM->Add(row(c.get()),column(r.get()),conj(v));
    };
    complex get(row r, column c) const{
      return conj(myM->Get(row(c.get()),column(r.get())));
    };

    size_t num_rows() const {return myM->num_of_cols();};
    size_t num_cols() const {return myM->num_of_rows();};
};


class calculation_Matrix: public Matrix{
public:

    typedef STD_TR1::shared_ptr<calculation_Matrix> ref;
    typedef STD_TR1::shared_ptr<const calculation_Matrix> const_ref;

    virtual ~calculation_Matrix() throw() {};

    void Apply(const complex* phi, 
           complex* result) const { apply_two(phi, result);};
    void Apply(const Vector& phi, 
           Vector& result) const;

    void Apply(complex* phi) const { apply_one(phi);};    
    void Apply(Vector& phi) const;

    void Scale(complex v) { scale(v);};

private:

// No Overloading for virtual functions!
    virtual 
    void apply_two(const complex* phi, 
           complex* result) const = 0;
    virtual
    void apply_one(complex* phi) const;    

    virtual
    void scale(complex v) = 0;
};


class orbital_dropper_Matrix: public Matrix{
public:

    orbital_dropper_Matrix(Matrix& t,std::vector<int> orbitals_to_drop):
    myM(&t),orbitals_to_drop(orbitals_to_drop)
    {
        size = myM->num_of_cols()<myM->num_of_rows() ? myM->num_of_cols() : myM->num_of_rows();//If the matrix is not square, the list is as long as the smaller dimension.
        correspondence_list = create_correspondence_list(orbitals_to_drop,size);
    };

    virtual ~orbital_dropper_Matrix() throw() {};


    static std::vector<int> create_correspondence_list(std::vector<int> orbitals_to_drop,int n);

private:
    Matrix* myM;

    orbital_dropper_Matrix();

    void insert(row r, column c, complex v){
        if (correspondence_list[r.get()]==-1 || correspondence_list[c.get()]==-1) return;
        else myM->Insert(row(correspondence_list[r.get()]),column(correspondence_list[c.get()]) ,v);
    };

    void add(row r, column c, complex v){
        if (correspondence_list[r.get()]==-1 || correspondence_list[c.get()]==-1) return;
        else myM->Add(row(correspondence_list[r.get()]),column(correspondence_list[c.get()]) ,v);
    };

    complex get(row r, column c) const{
        if (correspondence_list[r.get()]==-1 || correspondence_list[c.get()]==-1)
        {
            std::cout << "orbital_dropper_Matrix: get() of a nonexisting (dropped) item was attempted (this should not happen!). 0 was returned.";
            return 0;
        }
        else return myM->Get(r,c);
    };

    std::vector<int> orbitals_to_drop;
    std::vector<int> correspondence_list;

    int size;

    size_t num_rows() const { return size;};
    size_t num_cols() const { return size;};


};

#endif