TBTK
Matrix.h
Go to the documentation of this file.
1 /* Copyright 2017 Kristofer Bj√∂rnson
2  *
3  * Licensed under the Apache License, Version 2.0 (the "License");
4  * you may not use this file except in compliance with the License.
5  * You may obtain a copy of the License at
6  *
7  * http://www.apache.org/licenses/LICENSE-2.0
8  *
9  * Unless required by applicable law or agreed to in writing, software
10  * distributed under the License is distributed on an "AS IS" BASIS,
11  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12  * See the License for the specific language governing permissions and
13  * limitations under the License.
14  */
15 
23 #ifndef COM_DAFER45_TBTK_MATRIX
24 #define COM_DAFER45_TBTK_MATRIX
25 
26 #include "TBTK/TBTKMacros.h"
27 
28 #include <complex>
29 
30 namespace TBTK{
31 
32 template<typename DataType, unsigned int ROWS = 0, unsigned int COLS = 0>
33 class Matrix{
34 public:
36  Matrix();
37 
39  Matrix(const Matrix<DataType, ROWS, COLS> &matrix);
40 
43 
45  ~Matrix();
46 
50  );
51 
55  );
56 
58  const DataType& at(unsigned int row, unsigned int col) const;
59 
61  DataType& at(unsigned int row, unsigned int col);
62 
64  unsigned int getNumRows() const;
65 
67  unsigned int getNumCols() const;
68 private:
70  DataType data[ROWS*COLS];
71 };
72 
73 template<typename DataType>
74 class Matrix<DataType, 0, 0>{
75 public:
77  Matrix(unsigned int rows, unsigned int cols);
78 
80  Matrix(const Matrix<DataType, 0, 0> &matrix);
81 
84 
86  ~Matrix();
87 
90 
93 
95  const DataType& at(unsigned int row, unsigned int col) const;
96 
98  DataType& at(unsigned int row, unsigned int col);
99 
101  unsigned int getNumRows() const;
102 
104  unsigned int getNumCols() const;
105 
108  const Matrix<DataType, 0, 0> &rhs
109  ) const;
110 private:
112  DataType *data;
113 
115  unsigned int rows;
116 
118  unsigned int cols;
119 };
120 
121 template<>
122 class Matrix<std::complex<double>, 0, 0>{
123 public:
125  Matrix(unsigned int rows, unsigned int cols);
126 
128  Matrix(const Matrix<std::complex<double>, 0, 0> &matrix);
129 
131  Matrix(Matrix<std::complex<double>, 0, 0> &&matrix);
132 
134  ~Matrix();
135 
138  const Matrix<std::complex<double>, 0, 0> &rhs
139  );
140 
143  Matrix<std::complex<double>, 0, 0> &&rhs
144  );
145 
147  const std::complex<double>& at(
148  unsigned int row,
149  unsigned int col
150  ) const;
151 
153  std::complex<double>& at(unsigned int row, unsigned int col);
154 
156  unsigned int getNumRows() const;
157 
159  unsigned int getNumCols() const;
160 
163  const Matrix<std::complex<double>, 0, 0> &rhs
164  ) const;
165 
167  void invert();
168 
170  std::complex<double> determinant();
171 private:
173  std::complex<double> *data;
174 
176  unsigned int rows;
177 
179  unsigned int cols;
180 };
181 
182 template<typename DataType, unsigned int ROWS, unsigned int COLS>
184 }
185 
186 template<typename DataType, unsigned int ROWS, unsigned int COLS>
188  for(unsigned int n = 0; n < ROWS*COLS; n++)
189  data[n] = matrix.data[n];
190 }
191 
192 template<typename DataType, unsigned int ROWS, unsigned int COLS>
194  for(unsigned int n = 0; n < ROWS*COLS; n++)
195  data[n] = matrix.data[n];
196 }
197 
198 template<typename DataType>
199 Matrix<DataType, 0, 0>::Matrix(unsigned int rows, unsigned int cols){
200  this->rows = rows;
201  this->cols = cols;
202  data = new DataType[rows*cols];
203 }
204 
205 template<typename DataType>
207  rows = matrix.rows;
208  cols = matrix.cols;
209  data = new DataType[rows*cols];
210  for(unsigned int n = 0; n < rows*cols; n++)
211  data[n] = matrix.data[n];
212 }
213 
214 template<typename DataType>
216  rows = matrix.rows;
217  cols = matrix.cols;
218  data = matrix.data;
219  matrix.data = nullptr;
220 }
221 
222 inline Matrix<std::complex<double>, 0, 0>::Matrix(unsigned int rows, unsigned int cols){
223  this->rows = rows;
224  this->cols = cols;
225  data = new std::complex<double>[rows*cols];
226 }
227 
229  const Matrix<std::complex<double>, 0, 0> &matrix
230 ){
231  rows = matrix.rows;
232  cols = matrix.cols;
233  data = new std::complex<double>[rows*cols];
234  for(unsigned int n = 0; n < rows*cols; n++)
235  data[n] = matrix.data[n];
236 }
237 
239  Matrix<std::complex<double>, 0, 0> &&matrix
240 ){
241  rows = matrix.rows;
242  cols = matrix.cols;
243  data = matrix.data;
244  matrix.data = nullptr;
245 }
246 
247 template<typename DataType, unsigned int ROWS, unsigned int COLS>
249 }
250 
251 template<typename DataType>
253  if(data != nullptr)
254  delete [] data;
255 }
256 
258  if(data != nullptr)
259  delete [] data;
260 }
261 
262 template<typename DataType, unsigned int ROWS, unsigned int COLS>
265 ){
266  if(this != &rhs){
267  for(unsigned int n = 0; n < ROWS*COLS; n++)
268  data[n] = rhs.data[n];
269  }
270 
271  return *this;
272 }
273 
274 template<typename DataType, unsigned int ROWS, unsigned int COLS>
277 ){
278  if(this != &rhs){
279  for(unsigned int n = 0; n < ROWS*COLS; n++)
280  data[n] = rhs.data[n];
281  }
282 
283  return *this;
284 }
285 
286 template<typename DataType>
288  const Matrix<DataType, 0, 0> &rhs
289 ){
290  if(this != &rhs){
291  rows = rhs.rows;
292  cols = rhs.cols;
293 
294  if(data != nullptr)
295  delete [] data;
296 
297  data = new DataType[rows*cols];
298  for(unsigned int n = 0; n < rows*cols; n++)
299  data[n] = rhs.data[n];
300  }
301 
302  return *this;
303 }
304 
305 template<typename DataType>
308 ){
309  if(this != &rhs){
310  rows = rhs.rows;
311  cols = rhs.cols;
312 
313  if(data != nullptr)
314  delete [] data;
315 
316  data = rhs.data;
317  rhs.data = nullptr;
318  }
319 
320  return *this;
321 }
322 
324  const Matrix<std::complex<double>, 0, 0> &rhs
325 ){
326  if(this != &rhs){
327  rows = rhs.rows;
328  cols = rhs.cols;
329 
330  if(data != nullptr)
331  delete [] data;
332 
333  data = new std::complex<double>[rows*cols];
334  for(unsigned int n = 0; n < rows*cols; n++)
335  data[n] = rhs.data[n];
336  }
337 
338  return *this;
339 }
340 
341 inline Matrix<std::complex<double>, 0, 0>& Matrix<std::complex<double>, 0, 0>::operator=(
342  Matrix<std::complex<double>, 0, 0> &&rhs
343 ){
344  if(this != &rhs){
345  rows = rhs.rows;
346  cols = rhs.cols;
347 
348  if(data != nullptr)
349  delete [] data;
350 
351  data = rhs.data;
352  rhs.data = nullptr;
353  }
354 
355  return *this;
356 }
357 
358 template<typename DataType, unsigned int ROWS, unsigned int COLS>
360  unsigned int row,
361  unsigned int col
362 ) const{
363  return data[row + ROWS*col];
364 }
365 
366 template<typename DataType>
368  unsigned int row,
369  unsigned int col
370 ) const{
371  return data[row + rows*col];
372 }
373 
374 inline const std::complex<double>& Matrix<std::complex<double>, 0, 0>::at(
375  unsigned int row,
376  unsigned int col
377 ) const{
378  return data[row + rows*col];
379 }
380 
381 template<typename DataType, unsigned int ROWS, unsigned int COLS>
383  unsigned int row,
384  unsigned int col
385 ){
386  return data[row + ROWS*col];
387 }
388 
389 template<typename DataType>
391  unsigned int row,
392  unsigned int col
393 ){
394  return data[row + rows*col];
395 }
396 
397 inline std::complex<double>& Matrix<std::complex<double>, 0, 0>::at(
398  unsigned int row,
399  unsigned int col
400 ){
401  return data[row + rows*col];
402 }
403 
404 template<typename DataType, unsigned int ROWS, unsigned int COLS>
406  return ROWS;
407 }
408 
409 template<typename DataType>
411  return rows;
412 }
413 
414 inline unsigned int Matrix<std::complex<double>, 0, 0>::getNumRows() const{
415  return rows;
416 }
417 
418 template<typename DataType, unsigned int ROWS, unsigned int COLS>
420  return COLS;
421 }
422 
423 template<typename DataType>
425  return cols;
426 }
427 
428 inline unsigned int Matrix<std::complex<double>, 0, 0>::getNumCols() const{
429  return cols;
430 }
431 
432 template<typename DataType>
434  const Matrix<DataType, 0, 0> &rhs
435 ) const{
436  TBTKAssert(
437  cols == rhs.rows,
438  "Matrix::operator*()",
439  "Incompatible matrix dimensions.",
440  "The matrix dimensions are " << rows << "x" << cols << " and "
441  << rhs.rows << "x" << rhs.cols << "\n"
442  );
443 
444  Matrix<DataType> result(rows, rhs.cols);
445  for(unsigned int row = 0; row < rows; row++){
446  for(unsigned int col = 0; col < rhs.cols; col++){
447  result.at(row, col) = 0.;
448 
449  for(unsigned int n = 0; n < cols; n++){
450  result.at(row, col)
451  += at(row, n)*rhs.at(n, col);
452  }
453  }
454  }
455 
456  return result;
457 }
458 
459 inline const Matrix<std::complex<double>, 0, 0> Matrix<std::complex<double>, 0, 0>::operator*(
460  const Matrix<std::complex<double>, 0, 0> &rhs
461 ) const{
462  TBTKAssert(
463  cols == rhs.rows,
464  "Matrix::operator*()",
465  "Incompatible matrix dimensions.",
466  "The matrix dimensions are " << rows << "x" << cols << " and "
467  << rhs.rows << "x" << rhs.cols << "\n"
468  );
469 
470  Matrix<std::complex<double>> result(rows, rhs.cols);
471  for(unsigned int row = 0; row < rows; row++){
472  for(unsigned int col = 0; col < rhs.cols; col++){
473  result.at(row, col) = 0.;
474 
475  for(unsigned int n = 0; n < cols; n++){
476  result.at(row, col)
477  += at(row, n)*rhs.at(n, col);
478  }
479  }
480  }
481 
482  return result;
483 }
484 
485 extern "C"{
486  void zgetrf_(
487  int *M,
488  int *N,
489  std::complex<double> *A,
490  int *lda,
491  int *ipiv,
492  int *info
493  );
494  void zgetri_(
495  int *N,
496  std::complex<double> *A,
497  int *lda,
498  int *ipiv,
499  std::complex<double> *work,
500  int *lwork,
501  int *info
502  );
503 };
504 
505 inline void Matrix<std::complex<double>, 0, 0>::invert(){
506  TBTKAssert(
507  rows == cols,
508  "Matrix::invert()",
509  "Invalid matrix dimension. Only square matrices can be"
510  << " inverted, but the matrix has size " << rows << "x" << cols
511  << "\n",
512  ""
513  );
514 
515  int *ipiv = new int[std::min(rows, cols)];
516  int lwork = rows*cols;
517  std::complex<double> *work = new std::complex<double>[lwork];
518  int info;
519 
520  zgetrf_((int*)&rows, (int*)&cols, data, (int*)&rows, ipiv, &info);
521 
522  if(info < 0){
523  TBTKExit(
524  "Matrix::invert()",
525  "Argument '" << -info << "' to zgetrf_() is invlid.",
526  "This should never happen, contact the developer."
527  );
528  }
529  else if(info > 0){
530  TBTKExit(
531  "Matrix::invert()",
532  "Unable to invert matrix since it is signular.",
533  ""
534  );
535  }
536 
537  zgetri_((int*)&rows, data, (int*)&rows, ipiv, work, &lwork, &info);
538  TBTKAssert(
539  info == 0,
540  "Matrix::invert()",
541  "Inversion failed with error code 'INFO = " << info << "'.",
542  "See the documentation for the lapack function zgetri_() for"
543  << " further information."
544  );
545 
546  delete [] ipiv;
547  delete [] work;
548 }
549 
550 inline std::complex<double> Matrix<std::complex<double>, 0, 0>::determinant(){
551  TBTKAssert(
552  rows == cols,
553  "Matrix::determinant()",
554  "Invalid matrix dimension. The determinant can only be"
555  << " calculated for square matrices, but the matrix has size "
556  << rows << "x" << cols << "\n",
557  ""
558  );
559 
560  std::complex<double> *copy = new std::complex<double>[rows*cols];
561  for(unsigned int n = 0; n < rows*cols; n++)
562  copy[n] = data[n];
563 
564  int *ipiv = new int[std::min(rows, cols)];
565 // int lwork = rows*cols;
566 // std::complex<double> *work = new std::complex<double>[lwork];
567  int info;
568 
569  zgetrf_((int*)&rows, (int*)&cols, copy, (int*)&rows, ipiv, &info);
570 
571  if(info < 0){
572  TBTKExit(
573  "Matrix::determinant()",
574  "Argument '" << -info << "' to zgetrf_() is invlid.",
575  "This should never happen, contact the developer."
576  );
577  }
578  else if(info > 0){
579  TBTKExit(
580  "Matrix::determinant()",
581  "Unable to invert matrix since it is signular.",
582  ""
583  );
584  }
585 
586  std::complex<double> det = 1.;
587  for(unsigned int n = 0; n < rows; n++){
588  Streams::out << ipiv[n] << "\n";
589  if(ipiv[n]-1 == (int)n)
590  det *= copy[rows*n + n];
591  else
592  det *= -copy[rows*n + n];
593  }
594 
595  delete [] copy;
596  delete [] ipiv;
597 
598  return det;
599 }
600 
601 }; //End namespace TBTK
602 
603 #endif
FockStateRuleSet operator*(const LadderOperator< BitRegister > &ladderOperator, const FockStateRuleSet &fockStateRuleSet)
Definition: FockStateRuleSet.h:103
Precompiler macros.
unsigned int getNumRows() const
Definition: Matrix.h:405
const DataType & at(unsigned int row, unsigned int col) const
Definition: Matrix.h:367
Definition: Matrix.h:74
static std::ostream out
Definition: Streams.h:36
const DataType & at(unsigned int row, unsigned int col) const
Definition: Matrix.h:359
unsigned int getNumCols() const
Definition: Matrix.h:419
Definition: Matrix.h:33
Matrix< DataType, ROWS, COLS > & operator=(const Matrix< DataType, ROWS, COLS > &rhs)
Definition: Matrix.h:263
Definition: ModelFactory.h:35
Matrix()
Definition: Matrix.h:183
~Matrix()
Definition: Matrix.h:248