Greens-code
A modular quantum transport code
mkl_matrix.hpp
1 #ifndef MKL_MATRIX_HEADER_HPP
2 #define MKL_MATRIX_HEADER_HPP
3 
4 #include "matrix.hpp"
5 #include "shared_ptr.hpp"
6 
8 public:
9  typedef STD_TR1::shared_ptr<mkl_Matrix> ref;
10  typedef STD_TR1::shared_ptr<const mkl_Matrix> const_ref;
11 
12  mkl_Matrix():tf(t_n){};
13  ~mkl_Matrix() throw() {};
14 
15 // use with care, only to interact with legacy C and fortran code
16 // that expects pointer to matrix area
17  const complex* raw_storage() const {return get_storage();};
18  complex* raw_storage() {return get_storage();};
19 
20  const complex* raw_storage(const row r, const column c) const { return get_storage(r,c);};
21  complex* raw_storage(const row r, const column c) { return get_storage(r,c);};;
22 
23  enum trans {t_n,t_c,t_t};
24  void set_mkl_t_flag(trans t){tf = t;};
25 
26  void flip_mkl_t_flag();
27 
28  const char* t_flag() const;
29 
30 // leading dimension of the Matrix, usually (tf == t_n) we store
31 // increasing column indizes (constant row) adjacent in memory, and ld is num_of_rows
32  int ld() const {return ld_();};
33 
34 private:
35  static const char* c_t_n;
36  static const char* c_t_c;
37  static const char* c_t_t;
38 
39  trans tf;
40 
41  size_t num_rows() const { if (tf == t_n) return n_rows(); return n_cols();}
42  size_t num_cols() const { if (tf == t_n) return n_cols(); return n_rows();}
43 
44  void scale(complex v);
45 
46  void apply_two(const complex* phi,
47  complex* result) const;
48 
49  virtual size_t n_rows() const = 0;
50  virtual size_t n_cols() const = 0;
51 
52  virtual const complex* get_storage() const = 0;
53  virtual complex* get_storage() = 0;
54 
55  virtual const complex* get_storage(const row r, const column c) const = 0;
56  virtual complex* get_storage(const row r, const column c) = 0;
57 
58  virtual int ld_() const{
59  if (tf == t_n) return num_of_rows();
60  return num_of_cols();
61  };
62 
63 };
64 
65 class sub_mkl_Matrix: public mkl_Matrix{
66 public:
67  sub_mkl_Matrix(mkl_Matrix& t, row r, size_t rowsize,
68  column c, size_t colsize): myM(&t),
69  start(t.raw_storage(r,c)),
70  rows(rowsize),
71  cols(colsize),
72  my_c(c),
73  my_r(r){};
74 
75  virtual ~sub_mkl_Matrix() throw() {};
76 
77 private:
78  mkl_Matrix* myM;
79 
80  complex* start;
81 
82  size_t rows;
83  size_t cols;
84  column my_c;
85  row my_r;
86 
87  size_t n_rows() const {return rows;};
88  size_t n_cols() const {return cols;};
89 
90  const complex* get_storage() const { return start;};
91  complex* get_storage() { return start;};
92 
93  const complex* get_storage(const row r, const column c) const {
94  return myM->raw_storage(r+my_r,c+my_c);};
95  complex* get_storage(const row r, const column c){
96  return myM->raw_storage(r+my_r,c+my_c);};
97 
98  void insert(row r, column c, complex v){
99  myM->Insert(r + my_r,c + my_c,v);
100  };
101 
102  void add(row r, column c, complex v){
103  myM->Add(r + my_r,c + my_c,v);
104  };
105 
106  void scale(complex v);
107 
108  complex get(row r, column c) const{
109  return myM->Get(r + my_r,c + my_c);
110  };
111 
112  int ld_() const { return myM->ld();};
113 };
114 
115 struct Block_Matrix {
116  int n1;
117  int n2;
118 
119  sub_mkl_Matrix b22;
120  sub_mkl_Matrix b11;
121  sub_mkl_Matrix b12;
122  sub_mkl_Matrix b21;
123 
124  Block_Matrix(mkl_Matrix& m, int n_1, int n_2) : n1(n_1), n2(n_2),
125  b22(m, row(n1), n2, column(n1), n2),
126  b11(m, row(0), n1, column(0), n1),
127  b12(m, row(0), n1, column(n1), n2),
128  b21(m, row(n1), n2, column(0), n1) {
129  };
130 };
131 
133 public:
134 
135  wrap_pointer_Matrix(complex* p, row r, column c)
136  : data(p), nc(c.get()), nr(r.get()){};
137 
138  wrap_pointer_Matrix(complex* p, row r)
139  : data(p), nc(r.get()), nr(r.get()){};
140 
141  virtual ~wrap_pointer_Matrix() throw() {};
142 
143 private:
144  complex* data;
145  int nc;
146  int nr;
147 
148 
149  void insert(row r, column c, complex v){
150  data[index(r,c)] = v;
151  };
152 
153  void add(row r, column c, complex v){
154  data[index(r,c)] += v;
155  };
156 
157  complex get(row r, column c) const {
158  return data[index(r,c)];
159  };
160 
161  const complex* get_storage() const { return data;};
162  complex* get_storage(){ return data;};
163 
164  const complex* get_storage(const row r,
165  const column c) const { return &data[index(r,c)];};
166  complex* get_storage(const row r,
167  const column c){ return &data[index(r,c)];};
168 
169  size_t index(row r, column c) const { return r.get() + c.get()*nr;};
170 
171  size_t n_rows() const {return nr;};
172  size_t n_cols() const {return nc;};
173 };
174 
176 public:
177 
178  explicit hamilton_Matrix(int n) : data(n*n,0), N(n){};
179  ~hamilton_Matrix() throw() {};
180  std::vector<complex> get_diag() const;
181  complex* get_diag(complex*) const;
182 
183 private:
184 
185  std::vector<complex> data;
186  int N;
187 
188  void insert(row r, column c, complex v){
189  data[index(r,c)] = v;
190  };
191 
192  void add(row r, column c, complex v){
193  data[index(r,c)] += v;
194  };
195 
196  complex get(row r, column c) const {
197  return data[index(r,c)];
198  };
199 
200  const complex* get_storage() const { return &data[0];};
201  complex* get_storage(){ return &data[0];};
202 
203  const complex* get_storage(const row r,
204  const column c) const { return &data[index(r,c)];};
205  complex* get_storage(const row r,
206  const column c){ return &data[index(r,c)];};
207 
208 
209  size_t index(row r, column c) const { return r.get() + c.get()*N;};
210 
211  size_t n_rows() const {return N;};
212  size_t n_cols() const {return N;};
213 };
214 
215 
217 public:
218 
219  rectangular_Matrix(const row r, const column c)
220  : nc(c.get()), nr(r.get()), data(nc*nr){};
221 
222  ~rectangular_Matrix() throw() {};
223 
224 
225  void resize(const row r, const column c);
226 
227 private:
228  size_t n_rows() const {return nr;};
229  size_t n_cols() const {return nc;};
230 
231  int nc;
232  int nr;
233  std::vector<complex> data;
234 
235  size_t index(row r, column c) const { return r.get() + c.get()*nr;};
236 
237  void insert(row r, column c, complex v){
238  data[index(r,c)] = v;
239  };
240 
241  void add(row r, column c, complex v){
242  data[index(r,c)] += v;
243  };
244 
245  complex get(row r, column c) const {
246  return data[index(r,c)];
247  };
248 
249  const complex* get_storage() const{ return &data[0];};
250  complex* get_storage() { return &data[0];};
251  const complex* get_storage(const row r,
252  const column c) const{ return &data[index(r,c)];};
253  complex* get_storage(const row r,
254  const column c) { return &data[index(r,c)];};
255 
256 };
257 
258 #endif
void Insert(row r, column c, complex v)
Insert an element at a given position (overwritten if already there)
Definition: matrix.hpp:92
complex Get(row r, column c) const
Obtain an element at a given position.
Definition: matrix.hpp:96
Definition: mkl_matrix.hpp:132
Definition: mkl_matrix.hpp:216
Definition: mkl_matrix.hpp:175
Definition: mkl_matrix.hpp:7
Definition: mkl_matrix.hpp:65
This Base-Class offers Matrix-Vector multiplication and scaling.
Definition: matrix.hpp:322
void Add(row r, column c, complex v)
Add an element at a given position.
Definition: matrix.hpp:94
Definition: mkl_matrix.hpp:115