MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
SmallMatrix_impl.h
Go to the documentation of this file.
1 /*
2  This file forms part of hpGEM. This package has been developed over a number of years by various people at the University of Twente and a full list of contributors can be found at
3  http://hpgem.org/about-the-code/team
4 
5  This code is distributed using BSD 3-Clause License. A copy of which can found below.
6 
7 
8  Copyright (c) 2014, University of Twente
9  All rights reserved.
10 
11  Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
12 
13  1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
14 
15  2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
16 
17  3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
18 
19  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
20  */
21 
22 //Copyright (c) 2013-2020, The MercuryDPM Developers Team. All rights reserved.
23 //For the list of developers, see <http://www.MercuryDPM.org/Team>.
24 //
25 //Redistribution and use in source and binary forms, with or without
26 //modification, are permitted provided that the following conditions are met:
27 // * Redistributions of source code must retain the above copyright
28 // notice, this list of conditions and the following disclaimer.
29 // * Redistributions in binary form must reproduce the above copyright
30 // notice, this list of conditions and the following disclaimer in the
31 // documentation and/or other materials provided with the distribution.
32 // * Neither the name MercuryDPM nor the
33 // names of its contributors may be used to endorse or promote products
34 // derived from this software without specific prior written permission.
35 //
36 //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
37 //ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
38 //WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
39 //DISCLAIMED. IN NO EVENT SHALL THE MERCURYDPM DEVELOPERS TEAM BE LIABLE FOR ANY
40 //DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
41 //(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
42 //LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
43 //ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
44 //(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
45 //SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
46 
47 //Note: This code is copied and adapted from hpGEM (see license above), version 22th of January 2016. It has been
48 //integrated into MercuryDPM at 16th of March 2017.
49 
50 
51 #include "SmallMatrix.h"
52 
53 
54 extern "C"
55 {
56 
58 void dgemv_(const char* trans, int* m, int* n, double* alpha, double* A, int* LDA, double* x, int* incx, double* beta,
59  double* y, int* incy);
60 
62 int
63 dgemm_(const char* transA, const char* transB, int* M, int* N, int* k, double* alpha, double* A, int* LDA, double* B,
64  int* LDB, double* beta, double* C, int* LDC);
65 
67 int daxpy_(int* N, double* DA, double* DX, int* INCX, double* DY, int* INCY);
68 
70 void dgetrf_(int* M, int* N, double* A, int* lda, int* IPIV, int* INFO);
71 
73 void dgetri_(int* N, double* A, int* lda, int* IPIV, double* WORK, int* lwork, int* INFO);
74 
76 void dgesv_(int* N, int* NRHS, double* A, int* lda, int* IPIV, double* B, int* LDB, int* INFO);
77 }
78 
79 
80 template<unsigned int numberOfRows, unsigned int numberOfColumns>
82 {
83  if (numberOfRows == 0)
84  {
85  logger(WARN, "Trying to multiply a vector with a matrix without any rows.");
87  }
88  if (numberOfColumns == 0)
89  {
90  logger(WARN, "Trying to multiply a vector with a matrix without any columns.");
92  }
93  int nr = numberOfRows;
94  int nc = numberOfColumns;
95 
96  int i_one = 1;
97  double d_one = 1.0;
98  double d_zero = 0.0;
99 
101 
102  logger(DEBUG, "Matrix size: % x % \n Vector size: %", nr, nc, right.size());
103 
104  dgemv_("N", &nr, &nc, &d_one, this->data(), &nr, right.data(), &i_one, &d_zero, result.data(), &i_one);
105  return result;
106 }
107 
108 template<unsigned int numberOfRows, unsigned int numberOfColumns>
111 {
112  if (numberOfRows == 0)
113  {
114  logger(WARN, "Trying to multiply a vector with a matrix without any rows.");
115  return SmallVector<numberOfRows>();
116  }
117  if (numberOfColumns == 0)
118  {
119  logger(WARN, "Trying to multiply a vector with a matrix without any columns.");
120  return SmallVector<numberOfRows>();
121  }
122  int nr = numberOfRows;
123  int nc = numberOfColumns;
124 
125  int i_one = 1;
126  double d_one = 1.0;
127  double d_zero = 0.0;
128 
130 
131  logger(DEBUG, "Matrix size: % x % \n Vector size: %", nr, nc, right.size());
132 
133  dgemv_("N", &nr, &nc, &d_one, (const_cast<double*>(this->data())), &nr, right.data(), &i_one, &d_zero,
134  result.data(), &i_one);
135  return result;
136 }
137 
138 template<unsigned int numberOfRows, unsigned int numberOfColumns>
139 template<unsigned int K>
142 {
143  int i = numberOfRows;
144  int j = numberOfColumns;
145  int k = K;
146 
147  if (numberOfColumns == 0)
148  {
149  logger(WARN, "Trying to multiply a matrix with a matrix without any columns.");
151  }
152  //The result of the matrix is left.numberOfRows, right.numberOfColumns()
154 
155  double d_one = 1.0;
156  double d_zero = 0.0;
157 
158  //Let the actual multiplication be done by Fortran
159  dgemm_("N", "N", &i, &k, &j, &d_one, this->data(), &i, const_cast<double*>(other.data()), &j, &d_zero, C.data(),
160  &i);
161 
162  return C;
163 }
164 
165 template<unsigned int numberOfRows, unsigned int numberOfColumns>
166 template<unsigned int K>
169 {
170  int i = numberOfRows;
171  int j = numberOfColumns;
172  int k = K;
173 
174  if (numberOfColumns == 0)
175  {
176  logger(WARN, "Trying to multiply a matrix with a matrix without any columns.");
178  }
179  //The result of the matrix is left.Nrows, right.NCols()
181 
182  double d_one = 1.0;
183  double d_zero = 0.0;
184 
185  //Let the actual multiplication be done by Fortran
186  dgemm_("N", "N", &i, &k, &j, &d_one, const_cast<double*>(this->data()), &i, const_cast<double*>(other.data()), &j,
187  &d_zero, C.data(), &i);
188 
189  return C;
190 }
191 
192 template<unsigned int numberOfRows, unsigned int numberOfColumns>
195 {
196  //blas does not support in-place multiply
197  return (*this) = (*this) * other;
198 }
199 
200 template<unsigned int numberOfRows, unsigned int numberOfColumns>
202 {
203  //copied from MiddleSizeMatrix to prevent constructing a temporary MiddleSizeMatrix
204  logger.assert(numberOfColumns == numberOfRows - 1,
205  "Matrix has wrong dimensions to construct the wedge stuff vector");
207 
208  switch (numberOfRows)
209  {
210  case 2:
211  result[0] = -(*this)(1, 0);
212  result[1] = +(*this)(0, 0);
213  break;
214  case 3:
215  result[0] = (*this)(1, 0) * (*this)(2, 1) - (*this)(2, 0) * (*this)(1, 1);
216  result[1] = (*this)(0, 1) * (*this)(2, 0) - (*this)(0, 0) * (*this)(2, 1); // includes minus sign already!
217  result[2] = (*this)(0, 0) * (*this)(1, 1) - (*this)(1, 0) * (*this)(0, 1);
218  break;
219  case 4:
220  result[0] = (*this)(1, 0) * (-(*this)(2, 1) * (*this)(3, 2) + (*this)(3, 1) * (*this)(2, 2)) +
221  (*this)(2, 0) * ((*this)(1, 1) * (*this)(3, 2) - (*this)(3, 1) * (*this)(1, 2)) +
222  (*this)(3, 0) * (-(*this)(1, 1) * (*this)(2, 2) + (*this)(2, 1) * (*this)(1, 2));
223 
224  result[1] = (*this)(0, 0) * ((*this)(2, 1) * (*this)(3, 2) - (*this)(3, 1) * (*this)(2, 2)) +
225  (*this)(2, 0) * (-(*this)(0, 1) * (*this)(3, 2) + (*this)(3, 1) * (*this)(0, 2)) +
226  (*this)(3, 0) * ((*this)(0, 1) * (*this)(2, 2) - (*this)(2, 1) * (*this)(0, 2));
227  result[2] = (*this)(0, 0) * (-(*this)(1, 1) * (*this)(3, 2) + (*this)(3, 1) * (*this)(1, 2)) +
228  (*this)(1, 0) * ((*this)(0, 1) * (*this)(3, 2) - (*this)(3, 1) * (*this)(0, 2)) +
229  (*this)(3, 0) * (-(*this)(0, 1) * (*this)(1, 2) + (*this)(1, 1) * (*this)(0, 2));
230  result[3] = (*this)(0, 0) * ((*this)(1, 1) * (*this)(2, 2) - (*this)(2, 1) * (*this)(1, 2)) +
231  (*this)(1, 0) * (-(*this)(0, 1) * (*this)(2, 2) + (*this)(2, 1) * (*this)(0, 2)) +
232  (*this)(2, 0) * ((*this)(0, 1) * (*this)(1, 2) - (*this)(1, 1) * (*this)(0, 2));
233  break;
234  default:
235  logger(ERROR, "Wedge product not implemented for this dimension");
236  } //end switch
237 
238  return (result);
239 }
240 
241 template<unsigned int numberOfRows, unsigned int numberOfColumns>
243 {
244  int nr = numberOfRows;
245  int nc = numberOfColumns;
246  int nPivot = std::min(numberOfRows, numberOfColumns);
247  int iPivot[nPivot];
248 
249  SmallMatrix result(*this);
250 
251  int info;
252 
253  dgetrf_(&nr, &nc, result.data(), &nr, iPivot, &info);
254 
255  return result;
256 }
257 
258 //class template specialization for this one function is a waste of code duplication
259 //just let the compiler figure out which case it needs
260 template<unsigned int numberOfRows, unsigned int numberOfColumns>
262 {
263  logger.assert(numberOfRows == numberOfColumns, "Matrix should be square to have a determinant!");
264 
265  switch (numberOfRows)
266  {
267  case 0:
268  return 1;
269  case 1:
270  return (*this)(0, 0);
271  case 2:
272  return (*this)(0, 0) * (*this)(1, 1) - (*this)(0, 1) * (*this)(1, 0);
273 
274  case 3:
275  return (*this)(0, 0) * ((*this)(1, 1) * (*this)(2, 2) - (*this)(1, 2) * (*this)(2, 1)) -
276  (*this)(0, 1) * ((*this)(1, 0) * (*this)(2, 2) - (*this)(2, 0) * (*this)(1, 2)) +
277  (*this)(0, 2) * ((*this)(1, 0) * (*this)(2, 1) - (*this)(2, 0) * (*this)(1, 1));
278 
279  case 4:
280  return ((*this)(3, 0) * (*this)(2, 1) * (*this)(0, 3) - (*this)(2, 0) * (*this)(3, 1) * (*this)(0, 3)) *
281  (*this)(1, 2) +
282  (-(*this)(3, 0) * (*this)(0, 3) * (*this)(2, 2) + (*this)(2, 0) * (*this)(0, 3) * (*this)(3, 2)) *
283  (*this)(1, 1) +
284  ((*this)(3, 1) * (*this)(0, 3) * (*this)(2, 2) - (*this)(2, 1) * (*this)(0, 3) * (*this)(3, 2)) *
285  (*this)(1, 0) +
286  (-(*this)(3, 0) * (*this)(2, 1) * (*this)(1, 3) + (*this)(2, 0) * (*this)(3, 1) * (*this)(1, 3) +
287  (-(*this)(2, 0) * (*this)(3, 3) + (*this)(3, 0) * (*this)(2, 3)) * (*this)(1, 1) +
288  ((*this)(2, 1) * (*this)(3, 3) - (*this)(3, 1) * (*this)(2, 3)) * (*this)(1, 0)) * (*this)(0, 2) +
289  ((*this)(3, 0) * (*this)(1, 3) * (*this)(2, 2) - (*this)(2, 0) * (*this)(1, 3) * (*this)(3, 2) +
290  ((*this)(2, 0) * (*this)(3, 3) - (*this)(3, 0) * (*this)(2, 3)) * (*this)(1, 2) +
291  (-(*this)(2, 2) * (*this)(3, 3) + (*this)(2, 3) * (*this)(3, 2)) * (*this)(1, 0)) * (*this)(0, 1) +
292  (-(*this)(3, 1) * (*this)(1, 3) * (*this)(2, 2) + (*this)(2, 1) * (*this)(1, 3) * (*this)(3, 2) +
293  ((*this)(3, 1) * (*this)(2, 3) - (*this)(2, 1) * (*this)(3, 3)) * (*this)(1, 2) +
294  (*this)(1, 1) * ((*this)(2, 2) * (*this)(3, 3) - (*this)(2, 3) * (*this)(3, 2))) * (*this)(0, 0);
295  // ... says Maple; this can possibly be done more efficiently,
296  // maybe even with LU (with pivoting, though...)
297  default:
298  logger(ERROR, "Computing the Determinant for size % is not implemented", numberOfRows);
299  break;
300  }
301  return 0;
302 }
303 
304 template<unsigned int numberOfRows, unsigned int numberOfColumns>
306 {
307  logger.assert(numberOfRows == numberOfColumns, "Cannot invert a non-square matrix");
309 
310  int nr = numberOfRows;
311  int nc = numberOfColumns;
312 
313  int nPivot = numberOfRows;
314  int iPivot[nPivot];
315 
316  int info = 0;
317 
318  dgetrf_(&nr, &nc, result.data(), &nr, iPivot, &info);
319 
320  int lwork = numberOfRows * numberOfColumns;
321  SmallMatrix work;
322  dgetri_(&nc, result.data(), &nc, iPivot, work.data(), &lwork, &info);
323 
324  return result;
325 }
326 
327 template<unsigned int numberOfRows, unsigned int numberOfColumns>
328 template<unsigned int numberOfRightHandSideColumns>
330 {
331  logger.assert(numberOfRows == numberOfColumns, "can only solve for square matrixes");
332 
333  int n = numberOfRows;
334  int nrhs = numberOfRightHandSideColumns;
335  int info;
336 
337  int IPIV[numberOfRows];
339  dgesv_(&n, &nrhs, matThis.data(), &n, IPIV, B.data(), &n, &info);
340 }
341 
342 template<unsigned int numberOfRows, unsigned int numberOfColumns>
344 {
345  logger.assert(numberOfRows == numberOfColumns, "can only solve for square matrixes");
346 
347  int n = numberOfRows;
348  int nrhs = 1;
349  int info;
350 
351  int IPIV[numberOfRows];
352  SmallMatrix matThis = *this;
353  dgesv_(&n, &nrhs, matThis.data(), &n, IPIV, b.data(), &n, &info);
354 }
355 
356 template<unsigned int numberOfRows, unsigned int numberOfColumns>
358 {
359  if (numberOfColumns == 0)
360  {
361  logger(WARN, "Trying to multiply a vector with a matrix without any columns.");
363  }
364  int nr = numberOfRows;
365  int nc = numberOfColumns;
366 
367  int i_one = 1;
368  double d_one = 1.0;
369  double d_zero = 0.0;
370 
372 
373  logger(DEBUG, "Matrix size: % x % \n Vector size: %", nr, nc, vec.size());
374 
375  dgemv_("T", &nr, &nc, &d_one, mat.data(), &nr, vec.data(), &i_one, &d_zero, result.data(), &i_one);
376  return result;
377 }
378 
Mdouble * data()
Definition: SmallMatrix.h:351
LL< Log::DEBUG > DEBUG
Debug information.
Definition: Logger.cc:58
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
double Mdouble
Definition: GeneralDefine.h:34
LL< Log::INFO > INFO
Info log level.
Definition: Logger.cc:55
Mdouble beta(Mdouble z, Mdouble w)
This is the beta function, returns the approximation based on cmath's implementation of ln(gamma) ...
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
void dgesv_(int *N, int *NRHS, double *A, int *lda, int *IPIV, double *B, int *LDB, int *INFO)
This is used for solve Ax=B for x. Again this is from LAPACK.
Implementation of a 3D vector (by Vitaliy).
Definition: SmallVector.h:61
LL< Log::ERROR > ERROR
Error log level.
Definition: Logger.cc:53
Mdouble determinant() const
const Mdouble * data() const
Definition: SmallVector.h:238
LL< Log::WARN > WARN
Warning log level.
Definition: Logger.cc:54
void solve(SmallMatrix< numberOfRows, numberOfRightHandSideColumns > &B) const
solves Ax=B where A is the current matrix and B is passed in. The result is returned in B...
SmallMatrix inverse() const
return the inverse in the vector result. The size of result matches the matrix.
SmallVector< numberOfRows > computeWedgeStuffVector() const
computeWedgeStuffVector.
SmallMatrix LUfactorisation() const
Return the LUfactorisation of the matrix.
unsigned int size() const
Definition: SmallVector.h:233
void dgemv_(const char *trans, int *m, int *n, double *alpha, double *A, int *LDA, double *x, int *incx, double *beta, double *y, int *incy)
This does matrix times vector and is from blas level 2.
int daxpy_(int *N, double *DA, double *DX, int *INCX, double *DY, int *INCY)
This is the gerneral scalar times vector + vector from blas, hence from blas level 1...
SmallMatrix & operator*=(const Mdouble &scalar)
Does matrix A_ij=scalar*A_ij.
Definition: SmallMatrix.h:220
void dgetri_(int *N, double *A, int *lda, int *IPIV, double *WORK, int *lwork, int *INFO)
This is the inverse calulation also from LAPACK. Calculates inverse if you pass it the LU factorisati...
Data type for small dense matrix.
Definition: SmallMatrix.h:67
SmallVector< numberOfColumns > operator*(SmallVector< numberOfRows > &vec, SmallMatrix< numberOfRows, numberOfColumns > &mat)
Multiplies a matrix with a vector.
void dgetrf_(int *M, int *N, double *A, int *lda, int *IPIV, int *INFO)
This is LU factorisation of the matrix A. This has been taken from LAPACK.
int dgemm_(const char *transA, const char *transB, int *M, int *N, int *k, double *alpha, double *A, int *LDA, double *B, int *LDB, double *beta, double *C, int *LDC)
This is the gernal matrix multiplication from blas level 3.
SmallVector< numberOfRows > operator*(SmallVector< numberOfColumns > &right)
Defines Matrix A times vector B and return vector C i.e. C_,j= A_ij B_,j.