1 //
2 //  Copyright (C) 2004-2006 Rational Discovery LLC
3 //
4 //   @@ All Rights Reserved @@
5 //  This file is part of the RDKit.
6 //  The contents are covered by the terms of the BSD license
7 //  which is included in the file license.txt, found at the root
8 //  of the RDKit source tree.
9 //
10 #include <RDGeneral/export.h>
11 #ifndef __RD_SQUARE_MATRIX_H__
12 #define __RD_SQUARE_MATRIX_H__
13 
14 #include "Matrix.h"
15 
16 namespace RDNumeric {
17 template <typename TYPE>
18 class SquareMatrix : public Matrix<TYPE> {
19  public:
20   //! brief Square matrix of size N
SquareMatrix()21   SquareMatrix(){};
22 
SquareMatrix(unsigned int N)23   explicit SquareMatrix(unsigned int N) : Matrix<TYPE>(N, N){};
24 
SquareMatrix(unsigned int N,TYPE val)25   SquareMatrix(unsigned int N, TYPE val) : Matrix<TYPE>(N, N, val){};
26 
SquareMatrix(unsigned int N,typename Matrix<TYPE>::DATA_SPTR data)27   SquareMatrix(unsigned int N, typename Matrix<TYPE>::DATA_SPTR data)
28       : Matrix<TYPE>(N, N, data){};
29 
30   // inline unsigned int size() const {
31   //  return d_nRows;
32   //};
33 
34   virtual SquareMatrix<TYPE> &operator*=(TYPE scale) {
35     Matrix<TYPE>::operator*=(scale);
36     return *this;
37   }
38 
39   //! In place matrix multiplication
40   virtual SquareMatrix<TYPE> &operator*=(const SquareMatrix<TYPE> &B) {
41     CHECK_INVARIANT(this->d_nCols == B.numRows(),
42                     "Size mismatch during multiplication");
43 
44     const TYPE *bData = B.getData();
45     TYPE *newData = new TYPE[this->d_dataSize];
46     unsigned int i, j, k;
47     unsigned int idA, idAt, idC, idCt, idB;
48     TYPE *data = this->d_data.get();
49     for (i = 0; i < this->d_nRows; i++) {
50       idA = i * this->d_nRows;
51       idC = idA;
52       for (j = 0; j < this->d_nCols; j++) {
53         idCt = idC + j;
54         newData[idCt] = (TYPE)(0.0);
55         for (k = 0; k < this->d_nCols; k++) {
56           idAt = idA + k;
57           idB = k * this->d_nRows + j;
58           newData[idCt] += (data[idAt] * bData[idB]);
59         }
60       }
61     }
62     boost::shared_array<TYPE> tsptr(newData);
63     this->d_data.swap(tsptr);
64     return (*this);
65   }
66 
67   //! In place matrix transpose
transposeInplace()68   virtual SquareMatrix<TYPE> &transposeInplace() {
69     unsigned int i, j;
70     unsigned int id1, id1t, id2;
71     TYPE temp;
72     TYPE *data = this->d_data.get();
73     for (i = 1; i < this->d_nRows; i++) {
74       id1 = i * this->d_nCols;
75       for (j = 0; j < i; j++) {
76         id1t = id1 + j;
77         id2 = j * this->d_nCols + i;
78         temp = data[id1t];
79         data[id1t] = data[id2];
80         data[id2] = temp;
81       }
82     }
83     return (*this);
84   }
85 };
86 typedef SquareMatrix<double> DoubleSquareMatrix;
87 }  // namespace RDNumeric
88 
89 #endif
90