supermatrix_MUMPS.hpp

#ifndef SUPERMATRIX_MUMPS_HPP
#define SUPERMATRIX_MUMPS_HPP

#include <memory>
#include "parallel.hpp"
#include "supermatrix.hpp"
#include "my_mkl.hpp"

#ifdef USECBLAS
namespace mkl{
    struct sparse_matrix_t{
    int nnz;
    int N;
    int indexbase;
    const int* irn;
    const int* jcn;
    complex* pointer;
    };

    const int SPARSE_MATRIX_TYPE_GENERAL = 0;
    const int SPARSE_OPERATION_NON_TRANSPOSE = 0;

    struct matrix_descr{
    int type;
    };

    const int SPARSE_INDEX_BASE_ZERO = 0;
    const int SPARSE_INDEX_BASE_ONE = 1;

    void mkl_sparse_z_create_coo(sparse_matrix_t*,int,int,int,int,const int*,const int*,complex*);
    void mkl_sparse_destroy(sparse_matrix_t&);
    void mkl_sparse_z_mv(int,complex,sparse_matrix_t,const matrix_descr&,const complex*,complex,complex*);
}

#endif

struct mklMultiplier_data{
    size_t N;
    size_t nnz;
    Vector values;
    std::vector<int> irn;
    std::vector<int> jcn;
    mkl::sparse_matrix_t pointer;

    mklMultiplier_data(){
    }
    
    ~mklMultiplier_data(){
//  if (pointer)
        mkl_sparse_destroy(pointer);
    }
};

typedef std::vector<std::pair<int,complex> > Multiplier_column_data;

typedef std::vector<Multiplier_column_data> Multiplier_data;

class super_Matrix_opaque_iterator: public super_Matrix_iterator_factory{    
public:
    size_t nnz() const{ return 0;};
};

class super_Matrix_Multiply: public super_Matrix_strategy{
    std::unique_ptr<Multiplier_data> data;

public:
    super_Matrix_Multiply( super_Matrix*,
               std::unique_ptr<Multiplier_data>);
    ~super_Matrix_Multiply(){};

    void Insert(row, column, complex){
    fail("Insert not valid in Multiplier Phase!\n");};
    void Add(row, column, complex){
    fail("Add not valid in Multiplier Phase!\n");};
    complex Get(row, column) const{
    fail("Get not valid in Multiplier Phase!\n"); return complex(0.,0.);};

    void addmyselfto(Matrix&) const{
    fail("Matrix Add not valid in Multiplier Phase!\n");
    };

    void Scale(complex);

    p_iterator_factory provide_iterator(){
    return p_iterator_factory(new super_Matrix_opaque_iterator);
    };
    p_const_iterator_factory provide_const_iterator() const{
    return p_const_iterator_factory(new super_Matrix_opaque_iterator);
    };

    void Apply(const complex* phi, 
                     complex* result) const;
    void Apply(      complex* result) const;
    void Output(std::ostream& os) const;
    void obtain_Diagonal(complex* diag);

    std::unique_ptr<super_Matrix_strategy> deep_copy(
    STD_TR1::shared_ptr<super_Matrix> new_mat
    ) const { fail("Multiplier Matrix cannot be cloned");
    return std::unique_ptr<super_Matrix_strategy>(
        new super_Matrix_SetUp(&*new_mat));
    };

    status my_status() const{return Multiply;};
    
};


class super_Matrix_mklMultiply_iterator: public super_Matrix_iterator_factory{    
    std::shared_ptr<mklMultiplier_data> data;
public:
    super_Matrix_mklMultiply_iterator(std::shared_ptr<mklMultiplier_data> d):
    data(d){}

    size_t nnz() const{ return data->nnz;};
    
    const mklMultiplier_data& get_data(){ return *data;};

    const complex* get_raw_storage() const { return &(data->values[0]);};
    complex* get_modifiable_raw_storage(){ return &(data->values[0]);};

    const int* get_row_data() const { return &(data->irn[0]);};
    const int* get_col_data() const { return &(data->jcn[0]);};

    int* get_modifiable_row_data() { return &(data->irn[0]);};
    int* get_modifiable_col_data() { return &(data->jcn[0]);};
    
};

class super_Matrix_mklMultiply_const_iterator: public super_Matrix_iterator_factory{    
    std::shared_ptr<const mklMultiplier_data> data;
public:
    super_Matrix_mklMultiply_const_iterator(std::shared_ptr<mklMultiplier_data> d):
    data(d){}

    size_t nnz() const{ return data->nnz;};
    
    const mklMultiplier_data& get_data(){ return *data;};

    const complex* get_raw_storage() const { return &(data->values[0]);};

    const int* get_row_data() const { return &(data->irn[0]);};
    const int* get_col_data() const { return &(data->jcn[0]);};
};

class super_Matrix_mklMultiply: public super_Matrix_strategy{
    std::shared_ptr<mklMultiplier_data> data;

    complex scalefactor;
    
public:
    super_Matrix_mklMultiply( super_Matrix*,
                  std::unique_ptr<mklMultiplier_data>);
    ~super_Matrix_mklMultiply();

    void Insert(row, column, complex){
    fail("Insert not valid in Multiplier Phase!\n");};
    void Add(row, column, complex){
    fail("Add not valid in Multiplier Phase!\n");};
    complex Get(row, column) const{
    fail("Get not valid in Multiplier Phase!\n"); return complex(0.,0.);};

    void addmyselfto(Matrix&) const{
    fail("Matrix Add not valid in Multiplier Phase!\n");
    };

    void Scale(complex);

    p_iterator_factory provide_iterator(){
    return p_iterator_factory( new super_Matrix_mklMultiply_iterator(data));
    };
    p_const_iterator_factory provide_const_iterator() const {
    return p_const_iterator_factory( new super_Matrix_mklMultiply_const_iterator(data));
    };

    void Apply(const complex* phi, 
                     complex* result) const;
    void Apply(      complex* result) const;
    void Output(std::ostream& os) const;
    void obtain_Diagonal(complex* diag);

    std::unique_ptr<super_Matrix_strategy> deep_copy(
    STD_TR1::shared_ptr<super_Matrix> new_mat
    ) const { fail("Multiplier Matrix cannot be cloned");
    return std::unique_ptr<super_Matrix_strategy>(
        new super_Matrix_SetUp(&*new_mat));
    };

    status my_status() const{return Multiply;};
    
};

struct MUMPS_matrix {
    size_t N;
    Vector data;
    std::vector<int> irn;
    std::vector<int> jcn;
};

class super_Matrix_Baked: public super_Matrix_strategy{
private:
    std::unique_ptr<MUMPS_matrix> mp;
    bool distributed;
public:
    super_Matrix_Baked( super_Matrix*,
            bool,
            std::unique_ptr<MUMPS_matrix>);
    ~super_Matrix_Baked();

    bool is_distributed() const { return distributed;};

    p_iterator_factory provide_iterator(){
    return p_iterator_factory(new super_Matrix_opaque_iterator);
    };
    p_const_iterator_factory provide_const_iterator() const{
    return p_const_iterator_factory(new super_Matrix_opaque_iterator);
    };

    void addmyselfto(Matrix&) const{
    fail("Matrix Add not valid in Baked Phase!\n");
    };

    void Insert(row, column, complex);
    void Add(row, column, complex);
    complex Get(row, column) const;

    void obtain_Diagonal(complex* diag);

    void Scale(complex);

    void Apply(const complex* phi, 
                     complex* result) const;
    void Apply(      complex* result) const;
    void Output(std::ostream& os) const;

    std::unique_ptr<super_Matrix_strategy> deep_copy(
    STD_TR1::shared_ptr<super_Matrix> new_mat
    ) const;

    status my_status() const{return Baked;};

    std::unique_ptr<MUMPS_matrix> release(); 
};

struct MUMPS_pimpl;

class super_Matrix_Factorized: public super_Matrix_strategy{

    std::unique_ptr<MUMPS_matrix> mp;

    std::unique_ptr<MUMPS_pimpl> id;

    bool distributed;
    
public:
    super_Matrix_Factorized( super_Matrix*,
                 bool,
                 std::unique_ptr<MUMPS_matrix>);
    ~super_Matrix_Factorized();


    p_iterator_factory provide_iterator(){
    return p_iterator_factory(new super_Matrix_opaque_iterator);
    };
    p_const_iterator_factory provide_const_iterator() const{
    return p_const_iterator_factory(new super_Matrix_opaque_iterator);
    };

    void Insert(row, column, complex){
    fail("Insert not valid in Factorized Phase!\n");};
    void Add(row, column, complex){
    fail("Add not valid in Factorized Phase!\n");};
    complex Get(row, column) const{
    fail("Get not valid in Factorized Phase!\n"); return complex(0.,0.);};
    void Scale(complex){
    fail("Scale not valid in Factorized Phase!\n");};

    void addmyselfto(Matrix&) const{
    fail("Matrix Add not valid in Factorized Phase!\n");
    };

    bool is_distributed() const { return distributed;};

    void obtain_Diagonal(complex* diag);

    void SolveDiagonal();

    void Apply(const complex* phi, 
                     complex* result) const;
    void Apply(      complex* result) const;
    void Output(std::ostream& os) const {
    fail("Output not valid in Factorized Phase!\n");};

    status my_status() const{return Factorized_dist;};

    std::unique_ptr<super_Matrix_strategy> deep_copy(
    STD_TR1::shared_ptr<super_Matrix> new_mat
    ) const { fail("Factorized Matrix cannot be cloned");
    return std::unique_ptr<super_Matrix_strategy>(
        new super_Matrix_SetUp(&*new_mat));};
};

class super_Matrix_SetUp;
class matrix_factory{
public:
    static std::unique_ptr<MUMPS_matrix > get_fact(MPI_Comm, const super_Matrix_SetUp&, bool);
    static std::unique_ptr<Multiplier_data > get_mult(MPI_Comm, const super_Matrix_SetUp&);
    static std::unique_ptr<mklMultiplier_data > get_mkl_mult(MPI_Comm, const super_Matrix_SetUp&);
};

#endif