TBTK
Need a break? Support the development by playing Polarity Puzzles
SparseMatrix.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_SPARSE_MATRIX
24 #define COM_DAFER45_TBTK_SPARSE_MATRIX
25 
26 #include "TBTK/TBTKMacros.h"
27 
28 #include <algorithm>
29 #include <tuple>
30 #include <vector>
31 
32 namespace TBTK{
33 
34 template<typename DataType>
35 class SparseMatrix{
36 public:
38  enum class StorageFormat {CSR, CSC};
39 
41  SparseMatrix();
42 
44  SparseMatrix(StorageFormat);
45 
47  SparseMatrix(const SparseMatrix &sparseMatrix);
48 
50  SparseMatrix(SparseMatrix &&sparseMatrix);
51 
53  SparseMatrix(
54  StorageFormat storageFormat,
55  unsigned int numRows,
56  unsigned int numCols
57  );
58 
60  SparseMatrix& operator=(const SparseMatrix &sparseMatrix);
61 
63  SparseMatrix& operator=(SparseMatrix &&sparseMatrix);
64 
66  ~SparseMatrix();
67 
69  void add(unsigned int row, unsigned int col, const DataType &value);
70 
72  void setStorageFormat(StorageFormat storageFormat);
73 
75  unsigned int getNumRows() const;
76 
78  unsigned int getNumColumns() const;
79 
81  unsigned int getCSRNumMatrixElements() const;
82 
84  unsigned int getCSCNumMatrixElements() const;
85 
87  const unsigned int* getCSRRowPointers() const;
88 
90  const unsigned int* getCSCColumnPointers() const;
91 
93  const unsigned int* getCSRColumns() const;
94 
96  const unsigned int* getCSCRows() const;
97 
99  const DataType* getCSRValues() const;
100 
102  const DataType* getCSCValues() const;
103 
105  void construct();
106 
113  SparseMatrix& operator+=(const SparseMatrix &rhs);
114 
121  SparseMatrix operator+(const SparseMatrix &rhs) const;
122 
129  SparseMatrix& operator-=(const SparseMatrix &rhs);
130 
137  SparseMatrix operator-(const SparseMatrix &rhs) const;
138 
145  SparseMatrix operator*(const SparseMatrix &rhs) const;
146 
153  SparseMatrix& operator*=(const DataType &rhs);
154 
161  SparseMatrix operator*(const DataType &rhs) const;
162 
170  friend SparseMatrix operator*(
171  const DataType &lhs,
172  const SparseMatrix &rhs
173  ){
174  return rhs*lhs;
175  }
176 
181  SparseMatrix hermitianConjugate() const;
182 
186  DataType trace() const;
187 
189  void print() const;
190 private:
193  std::vector<
194  std::tuple<unsigned int, unsigned int, DataType>
195  > dictionaryOfKeys;
196 
198  int numRows, numCols;
199 
202  bool allowDynamicDimensions;
203 
205  StorageFormat storageFormat;
206 
209  unsigned int *csxXPointers;
210 
212  unsigned int *csxY;
213 
215  DataType *csxValues;
216 
219  int csxNumMatrixElements;
220 public:
222  std::vector<
223  std::vector<std::tuple<unsigned int, DataType>>
224  > constructLIL();
225 
227  void sortLIL(
228  std::vector<
229  std::vector<std::tuple<unsigned int, DataType>>
230  > &listOfLists
231  ) const;
232 
235  void mergeLIL(
236  std::vector<
237  std::vector<std::tuple<unsigned int, DataType>>
238  > &listOfLists
239  ) const;
240 
242  void constructCSX();
243 
248  void convertCSXToLIL();
249 
256  static void multiply(
257  const SparseMatrix &lhs,
258  const SparseMatrix &rhs,
259  SparseMatrix &result
260  );
261 };
262 
263 template<typename DataType>
264 inline SparseMatrix<DataType>::SparseMatrix(){
265  csxXPointers = nullptr;
266  csxY = nullptr;
267  csxValues = nullptr;
268 }
269 
270 template<typename DataType>
271 inline SparseMatrix<DataType>::SparseMatrix(StorageFormat storageFormat){
272  this->storageFormat = storageFormat;
273 
274  numRows = -1;
275  numCols = -1;
276  allowDynamicDimensions = true;
277 
278  csxXPointers = nullptr;
279  csxY = nullptr;
280  csxValues = nullptr;
281  csxNumMatrixElements = -1;
282 }
283 
284 template<typename DataType>
285 inline SparseMatrix<DataType>::SparseMatrix(
286  StorageFormat storageFormat,
287  unsigned int numRows,
288  unsigned int numCols
289 ){
290  this->storageFormat = storageFormat;
291 
292  this->numRows = numRows;
293  this->numCols = numCols;
294  allowDynamicDimensions = false;
295 
296  csxXPointers = nullptr;
297  csxY = nullptr;
298  csxValues = nullptr;
299  csxNumMatrixElements = -1;
300 }
301 
302 template<typename DataType>
303 inline SparseMatrix<DataType>::SparseMatrix(
304  const SparseMatrix &sparseMatrix
305 ){
306  storageFormat = sparseMatrix.storageFormat;
307 
308  numRows = sparseMatrix.numRows;
309  numCols = sparseMatrix.numCols;
310  allowDynamicDimensions = sparseMatrix.allowDynamicDimensions;
311 
312  csxNumMatrixElements = sparseMatrix.csxNumMatrixElements;
313  if(csxNumMatrixElements == -1){
314  TBTKAssert(
315  sparseMatrix.csxXPointers == nullptr
316  && sparseMatrix.csxY == nullptr
317  && sparseMatrix.csxValues == nullptr,
318  "SparseMatrix::SparseMatrix()",
319  "Invalid pointers in the original SparseMatrix.",
320  "This should never happen, contact the developer."
321  );
322 
323  csxXPointers = nullptr;
324  csxY = nullptr;
325  csxValues = nullptr;
326  }
327  else{
328  TBTKAssert(
329  sparseMatrix.csxXPointers != nullptr
330  && sparseMatrix.csxY != nullptr
331  && sparseMatrix.csxValues != nullptr,
332  "SparseMatrix::SparseMatrix()",
333  "Invalid pointers in the original SparseMatrix.",
334  "This should never happen, contact the developer."
335  );
336 
337  switch(storageFormat){
338  case StorageFormat::CSR:
339  csxXPointers = new unsigned int[numRows+1];
340  for(int row = 0; row < numRows+1; row++){
341  csxXPointers[row]
342  = sparseMatrix.csxXPointers[row];
343  }
344  break;
345  case StorageFormat::CSC:
346  csxXPointers = new unsigned int[numCols+1];
347  for(int col = 0; col < numCols+1; col++){
348  csxXPointers[col]
349  = sparseMatrix.csxXPointers[col];
350  }
351  break;
352  default:
353  TBTKExit(
354  "SparseMatrix::SparseMatrix()",
355  "Unknow StorageFormat.",
356  "This should never happen, contact the"
357  << " developer."
358  );
359  }
360 
361  csxY = new unsigned int[csxNumMatrixElements];
362  csxValues = new DataType[csxNumMatrixElements];
363  for(int n = 0; n < csxNumMatrixElements; n++){
364  csxY[n] = sparseMatrix.csxY[n];
365  csxValues[n] = sparseMatrix.csxValues[n];
366  }
367  }
368 }
369 
370 template<typename DataType>
371 inline SparseMatrix<DataType>::SparseMatrix(
372  SparseMatrix &&sparseMatrix
373 ){
374  storageFormat = sparseMatrix.storageFormat;
375 
376  numRows = sparseMatrix.numRows;
377  numCols = sparseMatrix.numCols;
378  allowDynamicDimensions = sparseMatrix.allowDynamicDimensions;
379 
380  csxNumMatrixElements = sparseMatrix.csxNumMatrixElements;
381  if(csxNumMatrixElements == -1){
382  TBTKAssert(
383  sparseMatrix.csxXPointers == nullptr
384  && sparseMatrix.csxY == nullptr
385  && sparseMatrix.csxValues == nullptr,
386  "SparseMatrix::SparseMatrix()",
387  "Invalid pointers in the original SparseMatrix.",
388  "This should never happen, contact the developer."
389  );
390 
391  csxXPointers = nullptr;
392  csxY = nullptr;
393  csxValues = nullptr;
394  }
395  else{
396  TBTKAssert(
397  sparseMatrix.csxXPointers != nullptr
398  && sparseMatrix.csxY != nullptr
399  && sparseMatrix.csxValues != nullptr,
400  "SparseMatrix::SparseMatrix()",
401  "Invalid pointers in the original SparseMatrix.",
402  "This should never happen, contact the developer."
403  );
404 
405  csxXPointers = sparseMatrix.csxXPointers;
406  sparseMatrix.csxXPointers = nullptr;
407 
408  csxY = sparseMatrix.csxY;
409  sparseMatrix.csxY = nullptr;
410 
411  csxValues = sparseMatrix.csxValues;
412  sparseMatrix.csxValues = nullptr;
413  }
414 }
415 
416 template<typename DataType>
417 inline SparseMatrix<DataType>::~SparseMatrix(){
418  if(csxXPointers != nullptr)
419  delete [] csxXPointers;
420  if(csxY != nullptr)
421  delete [] csxY;
422  if(csxValues != nullptr)
423  delete [] csxValues;
424 }
425 
426 template<typename DataType>
427 inline SparseMatrix<DataType>& SparseMatrix<DataType>::operator=(
428  const SparseMatrix &rhs
429 ){
430  if(this != &rhs){
431  storageFormat = rhs.storageFormat;
432 
433  numRows = rhs.numRows;
434  numCols = rhs.numCols;
435  allowDynamicDimensions = rhs.allowDynamicDimensions;
436 
437  if(csxXPointers != nullptr)
438  delete [] csxXPointers;
439  if(csxY != nullptr)
440  delete [] csxY;
441  if(csxValues != nullptr)
442  delete [] csxValues;
443 
444  csxNumMatrixElements = rhs.csxNumMatrixElements;
445  if(csxNumMatrixElements == -1){
446  TBTKAssert(
447  rhs.csxXPointers == nullptr
448  && rhs.csxY == nullptr
449  && rhs.csxValues == nullptr,
450  "SparseMatrix::operator=()",
451  "Invalid pointers in the original SparseMatrix.",
452  "This should never happen, contact the developer."
453  );
454 
455  csxXPointers = nullptr;
456  csxY = nullptr;
457  csxValues = nullptr;
458  }
459  else{
460  TBTKAssert(
461  rhs.csxXPointers != nullptr
462  && rhs.csxY != nullptr
463  && rhs.csxValues != nullptr,
464  "SparseMatrix::SparseMatrix()",
465  "Invalid pointers in the original SparseMatrix.",
466  "This should never happen, contact the developer."
467  );
468 
469  switch(storageFormat){
470  case StorageFormat::CSR:
471  csxXPointers = new unsigned int[numRows+1];
472  for(unsigned int row = 0; row < numRows+1; row++){
473  csxXPointers[row]
474  = rhs.csxXPointers[row];
475  }
476  break;
477  case StorageFormat::CSC:
478  csxXPointers = new unsigned int[numCols+1];
479  for(unsigned int col = 0; col < numCols+1; col++){
480  csxXPointers[col]
481  = rhs.csxXPointers[col];
482  }
483  break;
484  default:
485  TBTKExit(
486  "SparseMatrix::SparseMatrix()",
487  "Unknow StorageFormat.",
488  "This should never happen, contact the"
489  << " developer."
490  );
491  }
492 
493  csxY = new unsigned int[csxNumMatrixElements];
494  csxValues = new DataType[csxNumMatrixElements];
495  for(unsigned int n = 0; n < csxNumMatrixElements; n++){
496  csxY[n] = rhs.csxY[n];
497  csxValues[n] = rhs.csxValues[n];
498  }
499  }
500  }
501 
502  return *this;
503 }
504 
505 template<typename DataType>
506 inline SparseMatrix<DataType>& SparseMatrix<DataType>::operator=(
507  SparseMatrix &&rhs
508 ){
509  if(this != &rhs){
510  storageFormat = rhs.storageFormat;
511 
512  numRows = rhs.numRows;
513  numCols = rhs.numCols;
514  allowDynamicDimensions = rhs.allowDynamicDimensions;
515 
516  if(csxXPointers != nullptr)
517  delete [] csxXPointers;
518  if(csxY != nullptr)
519  delete [] csxY;
520  if(csxValues != nullptr)
521  delete [] csxValues;
522 
523  csxNumMatrixElements = rhs.csxNumMatrixElements;
524  if(csxNumMatrixElements == -1){
525  TBTKAssert(
526  rhs.csxXPointers == nullptr
527  && rhs.csxY == nullptr
528  && rhs.csxValues == nullptr,
529  "SparseMatrix::operator=()",
530  "Invalid pointers in the original SparseMatrix.",
531  "This should never happen, contact the developer."
532  );
533 
534  csxXPointers = nullptr;
535  csxY = nullptr;
536  csxValues = nullptr;
537  }
538  else{
539  TBTKAssert(
540  rhs.csxXPointers != nullptr
541  && rhs.csxY != nullptr
542  && rhs.csxValues != nullptr,
543  "SparseMatrix::SparseMatrix()",
544  "Invalid pointers in the original SparseMatrix.",
545  "This should never happen, contact the developer."
546  );
547 
548  csxXPointers = rhs.csxXPointers;
549  rhs.csxXPointers = nullptr;
550 
551  csxY = rhs.csxY;
552  rhs.csxY = nullptr;
553 
554  csxValues = rhs.csxValues;
555  rhs.csxValues = nullptr;
556  }
557  }
558 
559  return *this;
560 }
561 
562 template<typename DataType>
563 inline void SparseMatrix<DataType>::add(
564  unsigned int row,
565  unsigned int col,
566  const DataType &value
567 ){
568  if(!allowDynamicDimensions){
569  TBTKAssert(
570  (int)row < numRows && (int)col < numCols,
571  "SparseMatrix::add()",
572  "Invalid matrix entry. The matrix was constructed"
573  " specifiying that the matrix dimension is '"
574  << numRows << "x" << numCols << " but tried to add"
575  << " matrix element with index '(" << row << ", "
576  << col << ")'",
577  "Ensure that the matrix elements are in range, or use"
578  << " a constructor which does not fix the matrix"
579  << " dimension."
580  );
581  }
582 
583  dictionaryOfKeys.push_back(std::make_tuple(row, col, value));
584 }
585 
586 template<typename DataType>
587 inline void SparseMatrix<DataType>::setStorageFormat(
588  StorageFormat storageFormat
589 ){
590  if(this->storageFormat != storageFormat){
591  convertCSXToLIL();
592  this->storageFormat = storageFormat;
593  constructCSX();
594  }
595 }
596 
597 template<typename DataType>
598 inline unsigned int SparseMatrix<DataType>::getNumRows() const{
599  TBTKAssert(
600  numRows != -1,
601  "SparseMatrix::getNumRows()",
602  "Number of rows not yet determined.",
603  ""
604  );
605 
606  return (unsigned int)numRows;
607 }
608 
609 template<typename DataType>
610 inline unsigned int SparseMatrix<DataType>::getNumColumns() const{
611  TBTKAssert(
612  numRows != -1,
613  "SparseMatrix::getNumRows()",
614  "Number of rows not yet determined.",
615  ""
616  );
617 
618  return (unsigned int)numCols;
619 }
620 
621 template<typename DataType>
622 inline unsigned int SparseMatrix<DataType>::getCSRNumMatrixElements() const{
623  TBTKAssert(
624  storageFormat == StorageFormat::CSR,
625  "SparseMatrix::getCSRNumMatrixElements()",
626  "Tried to access CSR number of matrix elements, but the matrix"
627  " is not on the CSR storage format.",
628  "Use SparseMatrix::setFormat() to change the storage format."
629  );
630 
631  if(csxNumMatrixElements == -1)
632  return 0;
633  else
634  return (unsigned int)csxNumMatrixElements;
635 }
636 
637 template<typename DataType>
638 inline unsigned int SparseMatrix<DataType>::getCSCNumMatrixElements() const{
639  TBTKAssert(
640  storageFormat == StorageFormat::CSC,
641  "SparseMatrix::getCSCNumMatrixElements()",
642  "Tried to access CSC number of matrix elements, but the matrix"
643  " is not on the CSC storage format.",
644  "Use SparseMatrix::setFormat() to change the storage format."
645  );
646 
647  if(csxNumMatrixElements == -1)
648  return 0;
649  else
650  return (unsigned int)csxNumMatrixElements;
651 }
652 
653 template<typename DataType>
654 inline const unsigned int* SparseMatrix<DataType>::getCSRRowPointers() const{
655  TBTKAssert(
656  storageFormat == StorageFormat::CSR,
657  "SparseMatrix::getCSRRowPointers()",
658  "Tried to access CSR row pointers, but the matrix is not on"
659  " the CSR storage format.",
660  "Use SparseMatrix::setFormat() to change the storage format."
661  );
662 
663  TBTKAssert(
664  csxXPointers != nullptr,
665  "SparseMatrix::getCSRRowPointers()",
666  "Tried to access CSR row pointers, but row pointers have not"
667  << " been constructed yet.",
668  ""
669  );
670 
671  return csxXPointers;
672 }
673 
674 template<typename DataType>
675 inline const unsigned int* SparseMatrix<DataType>::getCSCColumnPointers() const{
676  TBTKAssert(
677  storageFormat == StorageFormat::CSC,
678  "SparseMatrix::getCSCColumnPointers()",
679  "Tried to access CSC row pointers, but the matrix is not on"
680  " the CSC storage format.",
681  "Use SparseMatrix::setFormat() to change the storage format."
682  );
683 
684  TBTKAssert(
685  csxXPointers != nullptr,
686  "SparseMatrix::getCSCColumnPointers()",
687  "Tried to access CSC column pointers, but column pointers have"
688  << " not been constructed yet.",
689  ""
690  );
691 
692  return csxXPointers;
693 }
694 
695 template<typename DataType>
696 inline const unsigned int* SparseMatrix<DataType>::getCSRColumns() const{
697  TBTKAssert(
698  storageFormat == StorageFormat::CSR,
699  "SparseMatrix::getCSRColumns()",
700  "Tried to access CSR columns, but the matrix is not on the CSR"
701  << " storage format.",
702  "Use SparseMatrix::setFormat() to change the storage format."
703  );
704 
705  TBTKAssert(
706  csxY != nullptr,
707  "SparseMatrix::getCSRColumns()",
708  "Tried to access CSR columns, but columns have not been"
709  << " constructed yet.",
710  ""
711  );
712 
713  return csxY;
714 }
715 
716 template<typename DataType>
717 inline const unsigned int* SparseMatrix<DataType>::getCSCRows() const{
718  TBTKAssert(
719  storageFormat == StorageFormat::CSC,
720  "SparseMatrix::getCSCRows()",
721  "Tried to access CSC rows, but the matrix is not on the CSC"
722  << " storage format.",
723  "Use SparseMatrix::setFormat() to change the storage format."
724  );
725 
726  TBTKAssert(
727  csxY != nullptr,
728  "SparseMatrix::getCSCRows()",
729  "Tried to access CSC rows, but rows have not been constructed"
730  << " yet.",
731  ""
732  );
733 
734  return csxY;
735 }
736 
737 template<typename DataType>
738 inline const DataType* SparseMatrix<DataType>::getCSRValues() const{
739  TBTKAssert(
740  storageFormat == StorageFormat::CSR,
741  "SparseMatrix::getCSRValues()",
742  "Tried to access CSR values, but the matrix is not on the CSR"
743  << " storage format.",
744  "Use SparseMatrix::setFormat() to change the storage format."
745  );
746 
747  TBTKAssert(
748  csxValues != nullptr,
749  "SparseMatrix::getCSRValues()",
750  "Tried to access CSR values, but values have not been"
751  << " constructed yet.",
752  ""
753  );
754 
755  return csxValues;
756 }
757 
758 template<typename DataType>
759 inline const DataType* SparseMatrix<DataType>::getCSCValues() const{
760  TBTKAssert(
761  storageFormat == StorageFormat::CSC,
762  "SparseMatrix::getCSCValues()",
763  "Tried to access CSC values, but the matrix is not on the CSC"
764  << " storage format.",
765  "Use SparseMatrix::setFormat() to change the storage format."
766  );
767 
768  TBTKAssert(
769  csxValues != nullptr,
770  "SparseMatrix::getCSCValues()",
771  "Tried to access CSC values, but values have not been"
772  << " constructed yet.",
773  ""
774  );
775 
776  return csxValues;
777 }
778 
779 template<typename DataType>
780 inline void SparseMatrix<DataType>::construct(){
781  constructCSX();
782 }
783 
784 template<typename DataType>
785 inline SparseMatrix<DataType>& SparseMatrix<DataType>::operator+=(
786  const SparseMatrix &rhs
787 ){
788  TBTKAssert(
789  csxNumMatrixElements != -1 && rhs.csxNumMatrixElements != -1,
790  "SparseMatrix::operator+=()",
791  "Unable to add matrices since the matrices have not yet been"
792  << " constructed.",
793  "Ensure that SparseMatrix::construct() has been called for"
794  << " both matrices."
795  );
796  TBTKAssert(
797  storageFormat == rhs.storageFormat,
798  "SparseMatrix::operator+=()",
799  "The left and right hand sides must have the same storage"
800  << " format. But the left hand side has storage format '" << (
801  storageFormat == StorageFormat::CSR
802  ? "StorageFormat::CSR"
803  : "StorageFormat::CSC"
804  ) << "' while the right hand side has storage format '" << (
805  rhs.storageFormat == StorageFormat::CSR
806  ? "StorageFormat::CSR"
807  : "StorageFormat::CSC"
808  ) << "'.",
809  ""
810  );
811  TBTKAssert(
812  allowDynamicDimensions == rhs.allowDynamicDimensions,
813  "SparseMatrix::operator+=()",
814  "The left and right hand sides must either both have dynamic"
815  << " or both not have dynamic dimensions. But the left hand"
816  << " side " << (
817  allowDynamicDimensions
818  ? "has dynamic dimensions "
819  : "does not have dynamic dimensions "
820  ) << " whilte the right hand side " << (
821  allowDynamicDimensions
822  ? "has dynamic dimensions."
823  : "does not have dynamic dimensions."
824  ),
825  "Whether the SparseMatrix has dynamic dimensions or not"
826  << " depends on whether the number of rows and columns are"
827  << " passed to the SparseMatrix constructor or not."
828  );
829 
830  if(!allowDynamicDimensions){
831  TBTKAssert(
832  numRows == rhs.numRows && numCols == rhs.numCols,
833  "SparseMatrix::operator+=()",
834  "The left and right hand sides must have the same"
835  << " dimensions, but the left hand side has dimension"
836  << " '" << numRows << "x" << numCols << "' while the"
837  << " right hand side has dimensions '" << rhs.numRows
838  << "x" << numCols << "'.",
839  "If both matrices have dynamic dimensions their"
840  << " dimensions do not need to agree. To create"
841  << " matrices with dynamic dimensions, do not pass row"
842  << " and column numbers to the SparseMatrix"
843  << " constructor."
844  );
845  }
846 
847  convertCSXToLIL();
848 
849  switch(storageFormat){
850  case StorageFormat::CSR:
851  {
852  for(int row = 0; row < rhs.numRows; row++){
853  for(
854  int n = rhs.csxXPointers[row];
855  n < (int)rhs.csxXPointers[row+1];
856  n++
857  ){
858  add(row, rhs.csxY[n], rhs.csxValues[n]);
859  }
860  }
861  break;
862  }
863  case StorageFormat::CSC:
864  for(int col = 0; col < rhs.numCols; col++){
865  for(
866  int n = rhs.csxXPointers[col];
867  n < (int)rhs.csxXPointers[col+1];
868  n++
869  ){
870  add(rhs.csxY[n], col, rhs.csxValues[n]);
871  }
872  }
873  break;
874  default:
875  TBTKExit(
876  "SparseMatrix::operator+=()",
877  "Unknown storage format.",
878  "This should never happen, contact the developer."
879  );
880  }
881 
882  construct();
883 
884  return *this;
885 }
886 
887 template<typename DataType>
888 inline SparseMatrix<DataType> SparseMatrix<DataType>::operator+(
889  const SparseMatrix &rhs
890 ) const{
891  SparseMatrix sparseMatrix = *this;
892 
893  return sparseMatrix += rhs;
894 }
895 
896 template<typename DataType>
897 inline SparseMatrix<DataType>& SparseMatrix<DataType>::operator-=(
898  const SparseMatrix &rhs
899 ){
900  TBTKAssert(
901  csxNumMatrixElements != -1 && rhs.csxNumMatrixElements != -1,
902  "SparseMatrix::operator-=()",
903  "Unable to subtract matrices since the matrices have not yet"
904  << " been constructed.",
905  "Ensure that SparseMatrix::construct() has been called for"
906  << " both matrices."
907  );
908  TBTKAssert(
909  storageFormat == rhs.storageFormat,
910  "SparseMatrix::operator-=()",
911  "The left and right hand sides must have the same storage"
912  << " format. But the left hand side has storage format '" << (
913  storageFormat == StorageFormat::CSR
914  ? "StorageFormat::CSR"
915  : "StorageFormat::CSC"
916  ) << "' while the right hand side has storage format '" << (
917  rhs.storageFormat == StorageFormat::CSR
918  ? "StorageFormat::CSR"
919  : "StorageFormat::CSC"
920  ) << "'.",
921  ""
922  );
923  TBTKAssert(
924  allowDynamicDimensions == rhs.allowDynamicDimensions,
925  "SparseMatrix::operator-=()",
926  "The left and right hand sides must either both have dynamic"
927  << " or both not have dynamic dimensions. But the left hand"
928  << " side " << (
929  allowDynamicDimensions
930  ? "has dynamic dimensions "
931  : "does not have dynamic dimensions "
932  ) << " whilte the right hand side " << (
933  allowDynamicDimensions
934  ? "has dynamic dimensions."
935  : "does not have dynamic dimensions."
936  ),
937  "Whether the SparseMatrix has dynamic dimensions or not"
938  << " depends on whether the number of rows and columns are"
939  << " passed to the SparseMatrix constructor or not."
940  );
941 
942  if(!allowDynamicDimensions){
943  TBTKAssert(
944  numRows == rhs.numRows && numCols == rhs.numCols,
945  "SparseMatrix::operator-=()",
946  "The left and right hand sides must have the same"
947  << " dimensions, but the left hand side has dimension"
948  << " '" << numRows << "x" << numCols << "' while the"
949  << " right hand side has dimensions '" << rhs.numRows
950  << "x" << numCols << "'.",
951  "If both matrices have dynamic dimensions their"
952  << " dimensions do not need to agree. To create"
953  << " matrices with dynamic dimensions, do not pass row"
954  << " and column numbers to the SparseMatrix"
955  << " constructor."
956  );
957  }
958 
959  convertCSXToLIL();
960 
961  switch(storageFormat){
962  case StorageFormat::CSR:
963  {
964  for(int row = 0; row < rhs.numRows; row++){
965  for(
966  int n = rhs.csxXPointers[row];
967  n < (int)rhs.csxXPointers[row+1];
968  n++
969  ){
970  add(row, rhs.csxY[n], -rhs.csxValues[n]);
971  }
972  }
973  break;
974  }
975  case StorageFormat::CSC:
976  for(int col = 0; col < rhs.numCols; col++){
977  for(
978  int n = rhs.csxXPointers[col];
979  n < (int)rhs.csxXPointers[col+1];
980  n++
981  ){
982  add(rhs.csxY[n], col, -rhs.csxValues[n]);
983  }
984  }
985  break;
986  default:
987  TBTKExit(
988  "SparseMatrix::operator-=()",
989  "Unknown storage format.",
990  "This should never happen, contact the developer."
991  );
992  }
993 
994  construct();
995 
996  return *this;
997 }
998 
999 template<typename DataType>
1000 inline SparseMatrix<DataType> SparseMatrix<DataType>::operator-(
1001  const SparseMatrix &rhs
1002 ) const{
1003  SparseMatrix sparseMatrix = *this;
1004 
1005  return sparseMatrix -= rhs;
1006 }
1007 
1008 template<typename DataType>
1009 inline SparseMatrix<DataType> SparseMatrix<DataType>::operator*(
1010  const SparseMatrix &rhs
1011 ) const{
1012  TBTKAssert(
1013  csxNumMatrixElements != -1 && rhs.csxNumMatrixElements != -1,
1014  "SparseMatrix::operator*=()",
1015  "Unable to multiply matrices since the matrices have not yet"
1016  << " been constructed.",
1017  "Ensure that SparseMatrix::construct() has been called for"
1018  << " both matrices."
1019  );
1020  TBTKAssert(
1021  storageFormat == rhs.storageFormat,
1022  "SparseMatrix::operator*=()",
1023  "The left and right hand sides must have the same storage"
1024  << " format. But the left hand side has storage format '" << (
1025  storageFormat == StorageFormat::CSR
1026  ? "StorageFormat::CSR"
1027  : "StorageFormat::CSC"
1028  ) << "' while the right hand side has storage format '" << (
1029  rhs.storageFormat == StorageFormat::CSR
1030  ? "StorageFormat::CSR"
1031  : "StorageFormat::CSC"
1032  ) << "'.",
1033  ""
1034  );
1035  TBTKAssert(
1036  allowDynamicDimensions == rhs.allowDynamicDimensions,
1037  "SparseMatrix::operator*=()",
1038  "The left and right hand sides must either both have dynamic,"
1039  << " or both not have dynamic dimensions. But the left hand"
1040  << " side " << (
1041  allowDynamicDimensions
1042  ? "has dynamic dimensions "
1043  : "does not have dynamic dimensions "
1044  ) << " whilte the right hand side " << (
1045  allowDynamicDimensions
1046  ? "has dynamic dimensions."
1047  : "does not have dynamic dimensions."
1048  ),
1049  "Whether the SparseMatrix has dynamic dimensions or not"
1050  << " depends on whether the number of rows and columns are"
1051  << " passed to the SparseMatrix constructor or not."
1052  );
1053 
1054  if(!allowDynamicDimensions){
1055  TBTKAssert(
1056  numCols == rhs.numRows,
1057  "SparseMatrix::operator*=()",
1058  "The number of columns for the left hand side must be"
1059  << " equal to the number of rows for the right hand"
1060  << " side. But the left hand side has '" << numCols
1061  << "' columns while the right hand side has '"
1062  << rhs.numRows << "'.",
1063  "If both matrices have dynamic dimensions their"
1064  << " dimensions do not need to agree. To create"
1065  << " matrices with dynamic dimensions, do not pass row"
1066  << " and column numbers to the SparseMatrix"
1067  << " constructor."
1068  );
1069  }
1070 
1071  SparseMatrix result;
1072  if(allowDynamicDimensions)
1073  result = SparseMatrix(storageFormat);
1074  else
1075  result = SparseMatrix(storageFormat, numRows, rhs.numCols);
1076 
1077  switch(storageFormat){
1078  case StorageFormat::CSR:
1079  {
1080  SparseMatrix rhsCSC = rhs;
1081  rhsCSC.setStorageFormat(StorageFormat::CSC);
1082  multiply(*this, rhsCSC, result);
1083  break;
1084  }
1085  case StorageFormat::CSC:
1086  {
1087  SparseMatrix lhsCSR = *this;
1088  lhsCSR.setStorageFormat(StorageFormat::CSR);
1089  multiply(lhsCSR, rhs, result);
1090  break;
1091  }
1092  default:
1093  TBTKExit(
1094  "SparseMatrix::operator+=()",
1095  "Unknown storage format.",
1096  "This should never happen, contact the developer."
1097  );
1098  }
1099 
1100  result.construct();
1101 
1102  return result;
1103 }
1104 
1105 template<typename DataType>
1106 inline SparseMatrix<DataType>& SparseMatrix<DataType>::operator*=(
1107  const DataType &rhs
1108 ){
1109  TBTKAssert(
1110  csxNumMatrixElements != -1,
1111  "SparseMatrix::operator*=()",
1112  "Matrix not yet constructed.",
1113  "First call SparseMatrix::construct()."
1114  );
1115 
1116  for(int row = 0; row < numRows; row++){
1117  for(
1118  int n = csxXPointers[row];
1119  n < (int)csxXPointers[row+1];
1120  n++
1121  ){
1122  csxValues[n] *= rhs;
1123  }
1124  }
1125 
1126  return *this;
1127 }
1128 
1129 template<typename DataType>
1130 inline SparseMatrix<DataType> SparseMatrix<DataType>::operator*(
1131  const DataType &rhs
1132 ) const{
1133  SparseMatrix result = *this;
1134 
1135  return result *= rhs;
1136 }
1137 
1138 template<typename DataType>
1139 inline SparseMatrix<DataType> SparseMatrix<DataType>::hermitianConjugate(
1140 ) const{
1141  TBTKAssert(
1142  csxNumMatrixElements != -1,
1143  "SparseMatrix::hermitianConjugate()",
1144  "Matrix not yet constructed.",
1145  "First call SparseMatrix::construct()."
1146  );
1147 
1148  SparseMatrix result;
1149  if(allowDynamicDimensions)
1150  result = SparseMatrix(storageFormat);
1151  else
1152  result = SparseMatrix(storageFormat, numRows, numCols);
1153 
1154  switch(storageFormat){
1155  case StorageFormat::CSR:
1156  {
1157  for(int row = 0; row < numRows; row++){
1158  for(
1159  int n = csxXPointers[row];
1160  n < (int)csxXPointers[row+1];
1161  n++
1162  ){
1163  result.add(csxY[n], row, conj(csxValues[n]));
1164  }
1165  }
1166 
1167  break;
1168  }
1169  case StorageFormat::CSC:
1170  {
1171  for(int col = 0; col < numCols; col++){
1172  for(
1173  int n = csxXPointers[col];
1174  n < (int)csxXPointers[col+1];
1175  n++
1176  ){
1177  result.add(col, csxY[n], conj(csxValues[n]));
1178  }
1179  }
1180 
1181  break;
1182  }
1183  default:
1184  TBTKExit(
1185  "SparseMatrix::hermitianConjugate()",
1186  "Unknown storage format.",
1187  "This should never happen, contact the developer."
1188  );
1189  }
1190 
1191  result.construct();
1192 
1193  return result;
1194 }
1195 
1196 template<typename DataType>
1197 inline DataType SparseMatrix<DataType>::trace() const{
1198  TBTKAssert(
1199  csxNumMatrixElements != -1,
1200  "SparseMatrix::trace()",
1201  "Matrix not yet constructed.",
1202  "First call SparseMatrix::construct()."
1203  );
1204 
1205  switch(storageFormat){
1206  case StorageFormat::CSR:
1207  {
1208  DataType result = 0;
1209  for(int row = 0; row < numRows; row++){
1210  for(
1211  int n = csxXPointers[row];
1212  n < (int)csxXPointers[row+1];
1213  n++
1214  ){
1215  if((unsigned int)row == csxY[n])
1216  result += csxValues[n];
1217  }
1218  }
1219 
1220  return result;
1221  }
1222  case StorageFormat::CSC:
1223  {
1224  DataType result = 0;
1225  for(int col = 0; col < numCols; col++){
1226  for(
1227  int n = csxXPointers[col];
1228  n < (int)csxXPointers[col+1];
1229  n++
1230  ){
1231  if(csxY[n] == (unsigned int)col)
1232  result += csxValues[n];
1233  }
1234  }
1235 
1236  return result;
1237  }
1238  default:
1239  TBTKExit(
1240  "SparseMatrix::trace()",
1241  "Unknown storage format.",
1242  "This should never happen, contact the developer."
1243  );
1244  }
1245 }
1246 
1247 template<typename DataType>
1248 inline void SparseMatrix<DataType>::print() const{
1249  Streams::out << "### Dictionary of Keys (DOK) ###\n";
1250  if(dictionaryOfKeys.size() == 0)
1251  Streams::out << "-\n";
1252  for(unsigned int n = 0; n < dictionaryOfKeys.size(); n++){
1253  unsigned int row = std::get<0>(dictionaryOfKeys[n]);
1254  unsigned int col = std::get<1>(dictionaryOfKeys[n]);
1255  const DataType &value = std::get<2>(dictionaryOfKeys[n]);
1256  Streams::out << "(" << row << ", " << col << ", " << value << ")\n";
1257  }
1258 
1259  Streams::out << "\n";
1260  switch(storageFormat){
1261  case StorageFormat::CSR:
1262  Streams::out << "### Compressed sparse row (CSR) ###\n";
1263  if(csxNumMatrixElements == -1){
1264  Streams::out << "-\n";
1265  }
1266  else{
1267  Streams::out << "Row pointers:\n";
1268  for(int row = 0; row < numRows+1; row++)
1269  Streams::out << csxXPointers[row] << "\t";
1270  Streams::out << "\nColumns:\n";
1271  for(int n = 0; n < csxNumMatrixElements; n++)
1272  Streams::out << csxY[n] << "\t";
1273  Streams::out << "\nValues:\n";
1274  for(int n = 0; n < csxNumMatrixElements; n++)
1275  Streams::out << csxValues[n] << "\t";
1276  Streams::out << "\n";
1277  }
1278  break;
1279  case StorageFormat::CSC:
1280  Streams::out << "### Compressed sparse column (CSC) ###\n";
1281  if(csxNumMatrixElements == -1){
1282  Streams::out << "-\n";
1283  }
1284  else{
1285  Streams::out << "Column pointers:\n";
1286  for(int col = 0; col < numCols+1; col++)
1287  Streams::out << csxXPointers[col] << "\t";
1288  Streams::out << "\nRows:\n";
1289  for(int n = 0; n < csxNumMatrixElements; n++)
1290  Streams::out << csxY[n] << "\t";
1291  Streams::out << "\nValues:\n";
1292  for(int n = 0; n < csxNumMatrixElements; n++)
1293  Streams::out << csxValues[n] << "\t";
1294  Streams::out << "\n";
1295  }
1296  break;
1297  default:
1298  TBTKExit(
1299  "SparseMatrix::print()",
1300  "Unknow StorageFormat.",
1301  "This should never happen, contact the developer."
1302  );
1303  }
1304 }
1305 
1306 template<typename DataType>
1307 inline std::vector<
1308  std::vector<std::tuple<unsigned int, DataType>>
1309 > SparseMatrix<DataType>::constructLIL(){
1310  unsigned int numRows = 0;
1311  unsigned int numCols = 0;
1312  for(unsigned int n = 0; n < dictionaryOfKeys.size(); n++){
1313  unsigned int row = std::get<0>(dictionaryOfKeys[n]);
1314  unsigned int col = std::get<1>(dictionaryOfKeys[n]);
1315  if(row >= numRows)
1316  numRows = row+1;
1317  if(col >= numCols)
1318  numCols = col+1;
1319  }
1320 
1321  std::vector<
1322  std::vector<std::tuple<unsigned int, DataType>>
1323  > listOfLists;
1324 
1325  switch(storageFormat){
1326  case StorageFormat::CSR:
1327  listOfLists.reserve(numRows);
1328  for(unsigned int row = 0; row < numRows; row++){
1329  listOfLists.push_back(
1330  std::vector<std::tuple<unsigned int, DataType>>()
1331  );
1332  }
1333 
1334  for(unsigned int n = 0; n < dictionaryOfKeys.size(); n++){
1335  unsigned int row = std::get<0>(dictionaryOfKeys[n]);
1336  unsigned int col = std::get<1>(dictionaryOfKeys[n]);
1337  const DataType &value = std::get<2>(dictionaryOfKeys[n]);
1338 
1339  listOfLists[row].push_back(std::make_tuple(col, value));
1340  }
1341  break;
1342  case StorageFormat::CSC:
1343  listOfLists.reserve(numCols);
1344  for(unsigned int col = 0; col < numCols; col++){
1345  listOfLists.push_back(
1346  std::vector<std::tuple<unsigned int, DataType>>()
1347  );
1348  }
1349 
1350  for(unsigned int n = 0; n < dictionaryOfKeys.size(); n++){
1351  unsigned int row = std::get<0>(dictionaryOfKeys[n]);
1352  unsigned int col = std::get<1>(dictionaryOfKeys[n]);
1353  const DataType &value = std::get<2>(dictionaryOfKeys[n]);
1354 
1355  listOfLists[col].push_back(std::make_tuple(row, value));
1356  }
1357  break;
1358  default:
1359  TBTKExit(
1360  "SparseMatrix::constructLIL()",
1361  "Unknow StorageFormat.",
1362  "This should never happen, contact the developer."
1363  );
1364  }
1365 
1366  sortLIL(listOfLists);
1367  mergeLIL(listOfLists);
1368 
1369  dictionaryOfKeys.clear();
1370 
1371  return listOfLists;
1372 }
1373 
1374 template<typename DataType>
1375 inline void SparseMatrix<DataType>::sortLIL(
1376  std::vector<
1377  std::vector<std::tuple<unsigned int, DataType>>
1378  > &listOfLists
1379 ) const{
1380  for(unsigned int x = 0; x < listOfLists.size(); x++){
1381  std::sort(
1382  listOfLists[x].begin(),
1383  listOfLists[x].end(),
1384  [](
1385  const std::tuple<unsigned int, DataType> &t1,
1386  const std::tuple<unsigned int, DataType> &t2
1387  ){
1388  return std::get<0>(t1) < std::get<0>(t2);
1389  }
1390  );
1391  }
1392 }
1393 
1394 template<typename DataType>
1395 inline void SparseMatrix<DataType>::mergeLIL(
1396  std::vector<
1397  std::vector<std::tuple<unsigned int, DataType>>
1398  > &listOfLists
1399 ) const{
1400  for(unsigned int x = 0; x < listOfLists.size(); x++){
1401  for(int y = listOfLists[x].size()-1; y > 0; y--){
1402  unsigned int y1 = std::get<0>(listOfLists[x][y]);
1403  unsigned int y2 = std::get<0>(listOfLists[x][y-1]);
1404 
1405  if(y1 == y2){
1406  std::get<1>(listOfLists[x][y-1])
1407  += std::get<1>(listOfLists[x][y]);
1408  listOfLists[x].erase(
1409  listOfLists[x].begin() + y
1410  );
1411  }
1412  }
1413  }
1414 }
1415 
1416 template<typename DataType>
1417 inline void SparseMatrix<DataType>::constructCSX(){
1418  convertCSXToLIL();
1419 
1420 // if(dictionaryOfKeys.size() != 0){
1421  std::vector<
1422  std::vector<std::tuple<unsigned int, DataType>>
1423  > listOfLists = constructLIL();
1424 
1425  if(allowDynamicDimensions){
1426  switch(storageFormat){
1427  case StorageFormat::CSR:
1428  numRows = listOfLists.size();
1429  numCols = 0;
1430  for(
1431  unsigned int row = 0;
1432  row < listOfLists.size();
1433  row++
1434  ){
1435  if(listOfLists[row].size() == 0)
1436  continue;
1437 
1438  unsigned int maxCol = std::get<0>(
1439  listOfLists[row].back()
1440  );
1441  if((int)maxCol+1 > numCols)
1442  numCols = maxCol + 1;
1443  }
1444  break;
1445  case StorageFormat::CSC:
1446  numCols = listOfLists.size();
1447  numRows = 0;
1448  for(
1449  unsigned int col = 0;
1450  col < listOfLists.size();
1451  col++
1452  ){
1453  if(listOfLists[col].size() == 0)
1454  continue;
1455 
1456  unsigned int maxRow = std::get<0>(
1457  listOfLists[col].back()
1458  );
1459  if((int)maxRow+1 > numRows)
1460  numRows = maxRow + 1;
1461  }
1462  break;
1463  default:
1464  TBTKExit(
1465  "SparseMatrix::constructCSX()",
1466  "Unknow StorageFormat.",
1467  "This should never happen, contact the"
1468  << " developer."
1469  );
1470  }
1471  }
1472 
1473  csxNumMatrixElements = 0;
1474  for(unsigned int x = 0; x < listOfLists.size(); x++)
1475  csxNumMatrixElements += listOfLists[x].size();
1476 
1477  switch(storageFormat){
1478  case StorageFormat::CSR:
1479  csxXPointers = new unsigned int[numRows+1];
1480  break;
1481  case StorageFormat::CSC:
1482  csxXPointers = new unsigned int[numCols+1];
1483  break;
1484  default:
1485  TBTKExit(
1486  "SparseMatrix::constructCSX()",
1487  "Unknow StorageFormat.",
1488  "This should never happen, contact the"
1489  << " developer."
1490  );
1491  }
1492  csxY = new unsigned int[csxNumMatrixElements];
1493  csxValues = new DataType[csxNumMatrixElements];
1494 
1495  csxXPointers[0] = 0;
1496  unsigned int currentMatrixElement = 0;
1497  for(unsigned int x = 0; x < listOfLists.size(); x++){
1498  csxXPointers[x+1]
1499  = csxXPointers[x]
1500  + listOfLists[x].size();
1501 
1502  for(
1503  unsigned int y = 0;
1504  y < listOfLists[x].size();
1505  y++
1506  ){
1507  csxY[currentMatrixElement]
1508  = std::get<0>(listOfLists[x][y]);
1509  csxValues[currentMatrixElement]
1510  = std::get<1>(listOfLists[x][y]);
1511 
1512  currentMatrixElement++;
1513  }
1514  }
1515 
1516  switch(storageFormat){
1517  case StorageFormat::CSR:
1518  for(
1519  int row = listOfLists.size();
1520  row < numRows;
1521  row++
1522  ){
1523  csxXPointers[row+1] = csxXPointers[row];
1524  }
1525  break;
1526  case StorageFormat::CSC:
1527  for(
1528  int col = listOfLists.size();
1529  col < numCols;
1530  col++
1531  ){
1532  csxXPointers[col+1] = csxXPointers[col];
1533  }
1534  break;
1535  default:
1536  TBTKExit(
1537  "SparseMatrix::constructCSX()",
1538  "Unknow StorageFormat.",
1539  "This should never happen, contact the"
1540  << " developer."
1541  );
1542  }
1543 
1544  TBTKAssert(
1545  csxNumMatrixElements == (int)currentMatrixElement,
1546  "SparseMatrix::constructCSX()",
1547  "Invalid number of matrix elements.",
1548  "This should never happen, contact the developer."
1549  );
1550 // }
1551 }
1552 
1553 template<typename DataType>
1554 inline void SparseMatrix<DataType>::convertCSXToLIL(){
1555  if(csxNumMatrixElements != -1){
1556  TBTKAssert(
1557  csxXPointers != nullptr
1558  && csxY != nullptr
1559  && csxValues != nullptr,
1560  "SparseMatrix::convertCSXToLIL()",
1561  "'csxNumMatrixElements' is not -1, but a csx-pointer"
1562  << " is a nullptr.",
1563  "This should never happen, contact the developer."
1564  );
1565 
1566  switch(storageFormat){
1567  case StorageFormat::CSR:
1568  {
1569  unsigned int row = 0;
1570  for(int n = 0; n < csxNumMatrixElements; n++){
1571  if((int)csxXPointers[row+1] == n){
1572  for(
1573  int r = row+1;
1574  r < numRows+1;
1575  r++
1576  ){
1577  row++;
1578  if((int)csxXPointers[r+1] > n)
1579  break;
1580  }
1581  }
1582 
1583  add(row, csxY[n], csxValues[n]);
1584  }
1585 
1586  break;
1587  }
1588  case StorageFormat::CSC:
1589  {
1590  unsigned int col = 0;
1591  for(int n = 0; n < csxNumMatrixElements; n++){
1592  if((int)csxXPointers[col+1] == n){
1593  for(
1594  int c = col+1;
1595  c < numCols+1;
1596  c++
1597  ){
1598  col++;
1599  if((int)csxXPointers[c+1] > n)
1600  break;
1601  }
1602  }
1603 
1604  add(csxY[n], col, csxValues[n]);
1605  }
1606 
1607  break;
1608  }
1609  default:
1610  TBTKExit(
1611  "SparseMatrix::convertCSXToLIL()",
1612  "Unknow StorageFormat.",
1613  "This should never happen, contact the"
1614  << " developer."
1615  );
1616  }
1617 
1618  delete [] csxXPointers;
1619  delete [] csxY;
1620  delete [] csxValues;
1621  csxXPointers = nullptr;
1622  csxY = nullptr;
1623  csxValues = nullptr;
1624  csxNumMatrixElements = -1;
1625  }
1626 }
1627 
1628 template<typename DataType>
1629 inline void SparseMatrix<DataType>::multiply(
1630  const SparseMatrix &lhs,
1631  const SparseMatrix &rhs,
1632  SparseMatrix &result
1633 ){
1634  TBTKAssert(
1635  (
1636  lhs.storageFormat == StorageFormat::CSR
1637  && rhs.storageFormat == StorageFormat::CSC
1638  ),
1639  "SparseMatrix::multiply()",
1640  "Storage format combination not supported.",
1641  "This should never happen, contact the developer."
1642  //This algorithm assumes that the left hand side is on CSR
1643  //format and the right hand side is on CSC format. Appropriate
1644  //checks should be performed in the calling function.
1645  );
1646 
1647  for(int lhsRow = 0; lhsRow < lhs.numRows; lhsRow++){
1648  for(int rhsCol = 0; rhsCol < rhs.numCols; rhsCol++){
1649  DataType scalarProduct = 0;
1650  int lhsElement = lhs.csxXPointers[lhsRow];
1651  int rhsElement = rhs.csxXPointers[rhsCol];
1652  int lhsTerminatingElement = lhs.csxXPointers[lhsRow+1];
1653  int rhsTerminatingElement = rhs.csxXPointers[rhsCol+1];
1654  while(true){
1655  if(
1656  lhsElement == lhsTerminatingElement
1657  || rhsElement == rhsTerminatingElement
1658  ){
1659  break;
1660  }
1661 
1662  int difference = lhs.csxY[lhsElement]
1663  - rhs.csxY[rhsElement];
1664 
1665  if(difference < 0){
1666  lhsElement++;
1667  }
1668  else if(difference > 0){
1669  rhsElement++;
1670  }
1671  else{
1672  scalarProduct += lhs.csxValues[
1673  lhsElement
1674  ]*rhs.csxValues[
1675  rhsElement
1676  ];
1677 
1678  lhsElement++;
1679  rhsElement++;
1680  }
1681  }
1682 
1683  if(scalarProduct != DataType(0))
1684  result.add(lhsRow, rhsCol, scalarProduct);
1685  }
1686  }
1687 }
1688 
1689 }; //End of namesapce TBTK
1690 
1691 #endif
TBTK::SparseMatrix::getCSRColumns
const unsigned int * getCSRColumns() const
Definition: SparseMatrix.h:696
TBTK::SparseMatrix::getCSRValues
const DataType * getCSRValues() const
Definition: SparseMatrix.h:738
TBTK::SparseMatrix::operator*
SparseMatrix operator*(const SparseMatrix &rhs) const
Definition: SparseMatrix.h:1009
TBTK::SparseMatrix::operator*
friend SparseMatrix operator*(const DataType &lhs, const SparseMatrix &rhs)
Definition: SparseMatrix.h:170
TBTK::SparseMatrix
Definition: SparseMatrix.h:35
TBTK::SparseMatrix::getCSCColumnPointers
const unsigned int * getCSCColumnPointers() const
Definition: SparseMatrix.h:675
TBTK::SparseMatrix::sortLIL
void sortLIL(std::vector< std::vector< std::tuple< unsigned int, DataType >> > &listOfLists) const
Definition: SparseMatrix.h:1375
TBTK::SparseMatrix::hermitianConjugate
SparseMatrix hermitianConjugate() const
Definition: SparseMatrix.h:1139
TBTK::SparseMatrix::add
void add(unsigned int row, unsigned int col, const DataType &value)
Definition: SparseMatrix.h:563
TBTK::SparseMatrix::getNumRows
unsigned int getNumRows() const
Definition: SparseMatrix.h:598
TBTK::SparseMatrix::constructCSX
void constructCSX()
Definition: SparseMatrix.h:1417
TBTK::SparseMatrix::operator-=
SparseMatrix & operator-=(const SparseMatrix &rhs)
Definition: SparseMatrix.h:897
TBTK::SparseMatrix::operator*=
SparseMatrix & operator*=(const DataType &rhs)
Definition: SparseMatrix.h:1106
TBTK::SparseMatrix::getCSCNumMatrixElements
unsigned int getCSCNumMatrixElements() const
Definition: SparseMatrix.h:638
TBTK::SparseMatrix::operator+
SparseMatrix operator+(const SparseMatrix &rhs) const
Definition: SparseMatrix.h:888
TBTK::SparseMatrix::construct
void construct()
Definition: SparseMatrix.h:780
TBTK::SparseMatrix::~SparseMatrix
~SparseMatrix()
Definition: SparseMatrix.h:417
TBTK::SparseMatrix::constructLIL
std::vector< std::vector< std::tuple< unsigned int, DataType > > > constructLIL()
TBTK::SparseMatrix::getNumColumns
unsigned int getNumColumns() const
Definition: SparseMatrix.h:610
TBTK::SparseMatrix::trace
DataType trace() const
Definition: SparseMatrix.h:1197
TBTK::SparseMatrix::multiply
static void multiply(const SparseMatrix &lhs, const SparseMatrix &rhs, SparseMatrix &result)
Definition: SparseMatrix.h:1629
TBTK::SparseMatrix::operator=
SparseMatrix & operator=(const SparseMatrix &sparseMatrix)
Definition: SparseMatrix.h:427
TBTK::SparseMatrix::convertCSXToLIL
void convertCSXToLIL()
Definition: SparseMatrix.h:1554
TBTK::SparseMatrix::getCSRRowPointers
const unsigned int * getCSRRowPointers() const
Definition: SparseMatrix.h:654
TBTK::SparseMatrix< std::complex< double > >::StorageFormat
StorageFormat
Definition: SparseMatrix.h:38
TBTK::SparseMatrix::SparseMatrix
SparseMatrix()
Definition: SparseMatrix.h:264
TBTK::SparseMatrix::operator-
SparseMatrix operator-(const SparseMatrix &rhs) const
Definition: SparseMatrix.h:1000
TBTK::SparseMatrix::mergeLIL
void mergeLIL(std::vector< std::vector< std::tuple< unsigned int, DataType >> > &listOfLists) const
Definition: SparseMatrix.h:1395
TBTK::SparseMatrix::operator+=
SparseMatrix & operator+=(const SparseMatrix &rhs)
Definition: SparseMatrix.h:785
TBTK::Streams::out
static std::ostream out
Definition: Streams.h:70
TBTK::SparseMatrix::print
void print() const
Definition: SparseMatrix.h:1248
TBTK::SparseMatrix::getCSCValues
const DataType * getCSCValues() const
Definition: SparseMatrix.h:759
TBTKMacros.h
Precompiler macros.
TBTK::SparseMatrix::setStorageFormat
void setStorageFormat(StorageFormat storageFormat)
Definition: SparseMatrix.h:587
TBTK::SparseMatrix::getCSRNumMatrixElements
unsigned int getCSRNumMatrixElements() const
Definition: SparseMatrix.h:622
TBTK::SparseMatrix::getCSCRows
const unsigned int * getCSCRows() const
Definition: SparseMatrix.h:717