Alsvinn  0.5.3
The fast FVM simulator with UQ support
mat.hpp
Go to the documentation of this file.
1 /* Copyright (c) 2018 ETH Zurich, Kjetil Olsen Lye
2  * This program is free software: you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation, either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * This program is distributed in the hope that it will be useful,
8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10  * GNU General Public License for more details.
11  *
12  * You should have received a copy of the GNU General Public License
13  * along with this program. If not, see <http://www.gnu.org/licenses/>.
14  */
15 
16 #pragma once
17 #include <cassert>
18 #include <iostream>
19 #include <string>
20 namespace alsutils {
21 
22 template<class T, size_t NumberOfRows, size_t NumberOfColumns>
23 class matrix {
24 public:
25 
29 
30  for (size_t column = 0; column < NumberOfColumns; ++column) {
31  for (size_t row = 0; row < NumberOfRows; ++row) {
32  data[column][row] = T(0);
33  }
34  }
35  }
36 
38  __device__ __host__ const T& operator()(size_t row, size_t column) const {
39  assert(row < NumberOfRows);
40  assert(column < NumberOfColumns);
41  return data[column][row];
42  }
43 
45  __device__ __host__ T& operator()(size_t row, size_t column) {
46  assert(row < NumberOfRows);
47  assert(column < NumberOfColumns);
48  return data[column][row];
49  }
50 
51 
54  template<class VectorType>
55  __device__ __host__ VectorType operator*(const VectorType& vector) const {
56  static_assert(NumberOfColumns == NumberOfRows,
57  "Matrix-Vector multiplication only supported for quadratic matrices.");
58  static_assert(VectorType::size() == NumberOfColumns,
59  "Matrix vector multiplication given wrong dimensions");
60 
61  VectorType product;
62 
63  for (size_t column = 0; column < NumberOfColumns; ++column) {
64  for (size_t row = 0; row < NumberOfRows; ++row) {
65  product[row] += (*this)(row, column) * vector[column];
66  }
67  }
68 
69  return product;
70  }
71 
72 
73  __device__ __host__ self_type operator*(const self_type& matrix) const {
74  static_assert(NumberOfColumns == NumberOfRows,
75  "Matrix-Matrix multiplication only supported for quadratic matrices.");
76 
77  self_type product;
78 
79  for (size_t row = 0; row < NumberOfRows; ++row) {
80  for (size_t column = 0; column < NumberOfColumns; ++column) {
81  for (size_t i = 0; i < NumberOfRows; ++i) {
82  product(row, column) += (*this)(row, i) * matrix(i, column);
83  }
84  }
85  }
86 
87  return product;
88  }
89 
91  const {
93 
94  for (size_t column = 0; column < NumberOfColumns; ++column) {
95  for (size_t row = 0; row < NumberOfRows; ++row) {
96  transposedMatrix(column, row) = (*this)(row, column);
97  }
98  }
99 
100  return transposedMatrix;
101  }
102 
104  const {
106 
107  for (size_t column = 0; column < NumberOfColumns; ++column) {
108  T norm = 0;
109 
110  for (size_t row = 0; row < NumberOfRows; ++row) {
111  norm += (*this)(row, column) * (*this)(row, column);
112  }
113 
114  norm = sqrtf(norm);
115 
116  for (size_t row = 0; row < NumberOfRows; ++row) {
117  newMatrix(row, column) = (*this)(row, column) / norm;
118  }
119  }
120 
121  return newMatrix;
122  }
123 
125  static_assert(NumberOfColumns == NumberOfRows,
126  "Matrix-Vector multiplication only supported for quadratic matrices.");
128 
129  for (size_t i = 0; i < NumberOfColumns; ++i) {
130  identityMatrix(i, i) = 1;
131  }
132 
133  return identityMatrix;
134  }
135 
136  __host__ std::string str() const {
137  std::stringstream ss;
138 
139  for (size_t i = 0; i < NumberOfRows; ++i) {
140  for (size_t j = 0; j < NumberOfColumns; ++j) {
141  ss << (*this)(i, j);
142 
143  if (j < NumberOfColumns - 1) {
144  ss << ", ";
145  }
146  }
147  }
148 
149  return ss.str();
150  }
151 
152 
153 private:
154  T data[NumberOfColumns][NumberOfRows];
155 };
156 
157 
158 }
159 
160 template<class T, size_t NumberOfRows, size_t NumberOfColumns>
161 std::ostream& operator<<(
162  std::ostream& os,
164  os << "[" << std::endl;
165 
166  for (size_t i = 0; i < NumberOfRows; ++i) {
167  for (size_t j = 0; j < NumberOfColumns; ++j) {
168  os << mat(i, j);
169 
170  if (j < NumberOfColumns - 1) {
171  os << ", ";
172  }
173  }
174 
175  os << std::endl;
176  }
177 
178  os << "]";
179  return os;
180 }
181 
182 
183 
184 
185 
static __device__ __host__ matrix< T, NumberOfColumns, NumberOfRows > identity()
Definition: mat.hpp:124
__device__ __host__ matrix< T, NumberOfColumns, NumberOfRows > transposed() const
Definition: mat.hpp:90
__device__ __host__ VectorType operator*(const VectorType &vector) const
Definition: mat.hpp:55
#define __host__
Definition: types.hpp:46
std::ostream & operator<<(std::ostream &os, const alsutils::matrix< T, NumberOfRows, NumberOfColumns > &mat)
Definition: mat.hpp:161
matrix< T, NumberOfRows, NumberOfColumns > self_type
Definition: mat.hpp:26
__device__ __host__ const T & operator()(size_t row, size_t column) const
Get the matrix element at the given row and column.
Definition: mat.hpp:38
__device__ __host__ matrix()
Creates a matrix initialized to 0.
Definition: mat.hpp:28
Definition: mat.hpp:23
__device__ __host__ matrix< T, NumberOfColumns, NumberOfRows > normalized() const
Definition: mat.hpp:103
Various utilities for mpi and cuda.
Definition: Factory.hpp:3
__host__ std::string str() const
Definition: mat.hpp:136
#define static_assert(x, y)
Definition: types.hpp:52
#define __device__
Definition: types.hpp:45
__device__ __host__ T & operator()(size_t row, size_t column)
Get the matrix element at the given row and column.
Definition: mat.hpp:45
__device__ __host__ self_type operator*(const self_type &matrix) const
Definition: mat.hpp:73