Peano 4
Loading...
Searching...
No Matches
DynamicMatrix.cpp
Go to the documentation of this file.
1#include "DynamicMatrix.h"
4
5#include <iostream>
6#include <cstring>
7
8
12
14 if (innerProduct) {
15 assertionMsg(false, "not implemented yet");
16 }
17 else {
18 _rows = lhs._rows*rhs._rows;
19 _cols = lhs._cols*rhs._cols;
20
21 _m = tarch::allocateMemory<double>(_rows*_cols,tarch::MemoryLocation::Heap);
22
23 for (int lhsRow=0; lhsRow<lhs._rows; lhsRow++)
24 for (int rhsRow=0; rhsRow<rhs._rows; rhsRow++)
25 for (int lhsCol=0; lhsCol<lhs._cols; lhsCol++)
26 for (int rhsCol=0; rhsCol<rhs._cols; rhsCol++) {
27 _m[ serialise(lhsRow*rhs._rows+rhsRow, lhsCol*rhs._cols+rhsCol) ]
28 = lhs._m[ lhs.serialise(lhsRow,lhsCol) ] * rhs._m[ rhs.serialise(rhsRow,rhsCol) ];
29 }
30 }
31}
32
33tarch::la::DynamicMatrix::DynamicMatrix(int rows, int cols, std::initializer_list< std::initializer_list<double> > values):
34 _cols(cols),
35 _rows(rows),
36 _m(nullptr) {
37 _m = tarch::allocateMemory<double>(rows*cols,tarch::MemoryLocation::Heap);
38
39 int index = 0;
40 for (typename std::initializer_list< std::initializer_list<double> >::const_iterator p = values.begin(); p!=values.end(); p++)
41 for (typename std::initializer_list<double>::const_iterator pp = p->begin(); pp!=p->end(); pp++) {
42 _m[index] = *pp;
43 index++;
44 }
45}
46
48 _cols(rhs._cols),
49 _rows(rhs._rows),
50 _m(nullptr) {
51 _m = tarch::allocateMemory<double>(_rows*_cols,tarch::MemoryLocation::Heap);
52 memcpy(_m, rhs._m, _cols*_rows*sizeof(double));
53}
54
56 _cols(rhs._cols),
57 _rows(rhs._rows),
58 _m(rhs._m) {
59 rhs._m = nullptr;
60 rhs._cols = 0;
61 rhs._rows = 0;
62}
63
64tarch::la::DynamicMatrix& tarch::la::DynamicMatrix::operator=( std::initializer_list< std::initializer_list<double> > values ) {
65 int index = 0;
66 for (typename std::initializer_list< std::initializer_list<double> >::const_iterator p = values.begin(); p!=values.end(); p++)
67 for (typename std::initializer_list<double>::const_iterator pp = p->begin(); pp!=p->end(); pp++) {
68 assertion3(index<_cols*_rows, _cols, _rows, index);
69 _m[index] = *pp;
70 index++;
71 }
72 return *this;
73}
74
76 _cols(cols),
77 _rows(rows),
78 _m(nullptr) {
79 _m = tarch::allocateMemory<double>(rows*cols,tarch::MemoryLocation::Heap);
80 std::fill_n(_m,rows*cols,0.0);
81}
82
83/*
84tarch::la::DynamicMatrix::DynamicMatrix(const DynamicMatrix&& value):
85 _cols(value._cols),
86 _rows(value._rows),
87 _m(value._m) {
88}
89*/
90
96
98 bool result = (_rows==matrix._rows and _cols==matrix._cols);
99
100 for (int col=0; col<_cols; col++)
101 for (int row=0; row<_rows; row++) {
102 result &= tarch::la::equals( _m[ serialise(row,col) ], matrix._m[ serialise(row,col) ] );
103 }
104
105 return result;
106}
107
109 assertion4(row>=0, row, col, _rows, _cols);
110 assertion4(col>=0, row, col, _rows, _cols);
111 assertion4(row<_rows, row, col, _rows, _cols);
112 assertion4(col<_cols, row, col, _rows, _cols);
113 return _m[ serialise(row,col) ];
114}
115
117 assertion(row>=0);
118 assertion(col>=0);
119 assertion(row<_rows);
120 assertion(col<_cols);
121 return _m[ serialise(row,col) ];
122}
123
125 tarch::la::DynamicMatrix result(rows,rows);
126 for (int i=0; i<rows; i++) {
127 result(i,i) = 1.0;
128 }
129 return result;
130}
131
133 for (int col=0; col<_cols; col++)
134 for (int row=0; row<_rows; row++) {
135 _m[ serialise(row,col) ] *= value;
136 }
137}
138
139void tarch::la::DynamicMatrix::multiply( double* result, double* x ) {
140 std::fill_n(result,_rows,0.0);
141
142 for (int col=0; col<_cols; col++)
143 for (int row=0; row<_rows; row++) {
144 result[row] += _m[ serialise(row,col) ] * x[row];
145 }
146}
147
148std::string tarch::la::DynamicMatrix::toString(bool addLineBreaks) const {
149 std::ostringstream msg;
150 msg << "(rows=" << _rows << ",cols=" << _cols << ",{";
151 if (addLineBreaks) msg << std::endl;
152 for (int row=0; row<_rows; row++) {
153 if (addLineBreaks)
154 msg << std::endl;
155 else {
156 if (row!=0) msg << ",";
157 msg << "{";
158 }
159
160 for (int col=0; col<_cols; col++) {
161 if (addLineBreaks)
162 msg << " ";
163 else if (col!=0)
164 msg << ",";
165 msg << _m[serialise(row,col)];
166 }
167
168 if (not addLineBreaks)
169 msg << "}";
170 }
171 if (addLineBreaks) msg << std::endl;
172 msg << "})";
173 return msg.str();
174}
175
176std::string tarch::la::DynamicMatrix::vectorToString( double* values, int entries, bool addLineBreaks ) {
177 std::ostringstream msg;
178 msg << "(entries=" << entries << ",{";
179 for (int i=0; i<entries; i++) {
180 if (addLineBreaks)
181 msg << std::endl;
182 else {
183 if (i!=0) msg << ",";
184 msg << values[i];
185 }
186 }
187 msg << "})";
188 return msg.str();
189}
190
192 double* __restrict__ result,
193 const double* __restrict__ x,
194 int batchCount
195) {
196 batchedMultiplyAoS(result, x, batchCount, _rows, 0);
197}
198
200 double* __restrict__ result,
201 const double* __restrict__ x,
202 int batchSize,
203 int resultSize,
204 int firstRow
205) {
206 assertion3( batchSize>0, batchSize, resultSize, firstRow );
207 assertion3( resultSize>0, batchSize, resultSize, firstRow );
208 assertion3( firstRow>=0, batchSize, resultSize, firstRow );
209
210 std::fill_n(result,resultSize*batchSize,0.0);
211
212 for (int col=0; col<_cols; col++) {
213 for (int rowInResult=0; rowInResult<resultSize; rowInResult++)
214 #ifdef SharedOMP
215 #pragma omp simd
216 #endif
217 for (int i=0; i<batchSize; i++)
218 {
219 assertion3(rowInResult + firstRow<_rows, batchSize, resultSize, firstRow);
220 result[rowInResult*batchSize+i] += _m[ serialise(rowInResult + firstRow,col) ] * x[col*batchSize+i];
221 }
222 }
223}
224
225void tarch::la::DynamicMatrix::insertEmptyColumns( int number, int where, int repeatEveryKColumns ) {
226 assertion4( number>0, number, where, repeatEveryKColumns, toString() );
227 assertion4( where>=0, number, where, repeatEveryKColumns, toString() );
228 assertion4( where<=_cols, number, where, repeatEveryKColumns, toString() );
229 assertion4( repeatEveryKColumns>=0, number, where, repeatEveryKColumns, toString() );
230
231 double* oldData = _m;
232 int oldRows = _rows;
233 int oldCols = _cols;
234
235 const int numberOfInsertions = repeatEveryKColumns==0 ? number : (_cols/repeatEveryKColumns) * number;
236 _cols+=numberOfInsertions;
237
238 _m = tarch::allocateMemory<double>(_rows*_cols,tarch::MemoryLocation::Heap);
239 std::fill_n(_m,_rows*_cols,0.0);
240
241 int imageColumn = 0;
242 for (int col=0; col<oldCols; col++) {
243 if (col==where) {
244 imageColumn += number;
245 if (repeatEveryKColumns>=1) {
246 where += repeatEveryKColumns;
247 }
248 }
249 for (int row=0; row<oldRows; row++) {
250 _m[ serialise(row,imageColumn) ] = oldData[ serialise(row,col,oldRows,oldCols) ];
251 }
252 imageColumn++;
253 }
254
256}
257
258void tarch::la::DynamicMatrix::insertEmptyRows( int number, int where, int repeatEveryKRows ) {
259 assertion4( number>0, number, where, repeatEveryKRows, toString() );
260 assertion4( where>=0, number, where, repeatEveryKRows, toString() );
261 assertion4( where<_cols, number, where, repeatEveryKRows, toString() );
262 assertion4( repeatEveryKRows>=0, number, where, repeatEveryKRows, toString() );
263
264 double* oldData = _m;
265 int oldRows = _rows;
266 int oldCols = _cols;
267
268 const int numberOfInsertions = repeatEveryKRows==0 ? number : (_rows/repeatEveryKRows) * number;
269 _rows+=numberOfInsertions;
270
271 _m = tarch::allocateMemory<double>(_rows*_cols,tarch::MemoryLocation::Heap);
272 std::fill_n(_m,_rows*_cols,0.0);
273
274 int imageRow = 0;
275 for (int row=0; row<oldRows; row++) {
276 if (row==where) {
277 imageRow += number;
278 if (repeatEveryKRows>=1) {
279 where += repeatEveryKRows;
280 }
281 }
282 for (int col=0; col<oldCols; col++) {
283 _m[ serialise(imageRow,col) ] = oldData[ serialise(row,col,oldRows,oldCols) ];
284 }
285 imageRow++;
286 }
287
289}
290
292 double* oldData = _m;
293 int oldRows = _rows;
294 int oldCols = _cols;
295 _cols--;
296
297
298 _m = tarch::allocateMemory<double>(_rows*_cols,tarch::MemoryLocation::Heap);
299 std::fill_n(_m,_rows*_cols,0.0);
300
301 for (int col=0; col<_cols; col++)
302 for (int row=0; row<_rows; row++) {
303 int preImageColumn = col<number ? col : col+1;
304 _m[ serialise(row,col) ] = oldData[ serialise(row,preImageColumn,oldRows,oldCols) ];
305 }
306
308}
309
311 double* oldData = _m;
312 _m = tarch::allocateMemory<double>(_rows*_cols,tarch::MemoryLocation::Heap);
313 std::fill_n(_m,_rows*_cols,0.0);
314
315 for (int col=0; col<_cols; col++)
316 for (int row=0; row<_rows; row++) {
317 int imageRow = row + shift;
318 if (wrap) imageRow %= _rows;
319 if (imageRow>=0 and imageRow<_rows) {
320 _m[ serialise(imageRow,col) ] = oldData[ serialise(row,col) ];
321 }
322 }
323
325}
326
328 double* oldData = _m;
329 _m = tarch::allocateMemory<double>(_rows*_cols,tarch::MemoryLocation::Heap);
330 std::fill_n(_m,_rows*_cols,0.0);
331
332 for (int col=0; col<_cols; col++)
333 for (int row=0; row<_rows; row++) {
334 int imageCol = col+ shift;
335 if (wrap) imageCol %= _cols;
336 if (imageCol>=0 and imageCol<_cols) {
337 _m[ serialise(row,imageCol) ] = oldData[ serialise(row,col) ];
338 }
339 }
340
342}
343
345 return _rows;
346}
347
349 return _cols;
350}
351
352void tarch::la::DynamicMatrix::replicateRows( int blockSize, int numberOfReplications, int shiftAfterEveryReplication, bool extendColumnsToAccommodateShifts ) {
353 assertion3( blockSize>=1, blockSize, numberOfReplications, shiftAfterEveryReplication );
354 assertion3( _rows%blockSize==0, blockSize, numberOfReplications, shiftAfterEveryReplication );
355 assertion3( numberOfReplications>=2, blockSize, numberOfReplications, shiftAfterEveryReplication );
356
357 double* oldData = _m;
358 int oldRows = _rows;
359 int oldCols = _cols;
360
361 assertion2(numberOfReplications>=1,numberOfReplications,shiftAfterEveryReplication);
362 assertion2(shiftAfterEveryReplication>=0,numberOfReplications,shiftAfterEveryReplication);
363
364 const int numberOfBlocks = _rows / blockSize;
365 _rows = _rows*numberOfReplications;
366 if (extendColumnsToAccommodateShifts) {
367 _cols = _cols+shiftAfterEveryReplication*(numberOfReplications-1);
368 }
369
370 _m = tarch::allocateMemory<double>(_rows*_cols,tarch::MemoryLocation::Heap);
371
372 std::fill_n(_m,_rows*_cols,0.0);
373
374 for (int block=0; block<numberOfBlocks; block++) {
375 for (int replication=0; replication<numberOfReplications; replication++) {
376 for (int blockRow=0; blockRow<blockSize; blockRow++) {
377 const int destRow = (block*numberOfReplications+replication)*blockSize+blockRow;
378 const int srcRow = block*blockSize+blockRow;
379 for (int col=0; col<oldCols; col++) {
380 const int destCol = col+replication*shiftAfterEveryReplication;
381 if (destCol<_cols) {
382 _m[ serialise( destRow,destCol) ] = oldData[ serialise(srcRow,col,oldRows,oldCols) ];
383 }
384 }
385 }
386 }
387 }
388
390}
391
392void tarch::la::DynamicMatrix::replicateCols( int blockSize, int numberOfReplications, int shiftAfterEveryReplication, bool extendColumnsToAccommodateShifts ) {
393 assertion3( blockSize>=1, blockSize, numberOfReplications, shiftAfterEveryReplication );
394 assertion3( _rows%blockSize==0, blockSize, numberOfReplications, shiftAfterEveryReplication );
395 assertion3( numberOfReplications>=2, blockSize, numberOfReplications, shiftAfterEveryReplication );
396
397 double* oldData = _m;
398 int oldRows = _rows;
399 int oldCols = _cols;
400
401 assertion2(numberOfReplications>=1,numberOfReplications,shiftAfterEveryReplication);
402 assertion2(shiftAfterEveryReplication>=0,numberOfReplications,shiftAfterEveryReplication);
403
404 const int numberOfBlocks = _cols / blockSize;
405 _cols = _cols*numberOfReplications;
406 if (extendColumnsToAccommodateShifts) {
407 _rows = _rows+shiftAfterEveryReplication*(numberOfReplications-1);
408 }
409
410 _m = tarch::allocateMemory<double>(_rows*_cols,tarch::MemoryLocation::Heap);
411
412 std::fill_n(_m,_rows*_cols,0.0);
413
414 for (int block=0; block<numberOfBlocks; block++) {
415 for (int replication=0; replication<numberOfReplications; replication++) {
416 for (int blockCol=0; blockCol<blockSize; blockCol++) {
417 const int destCol = (block*numberOfReplications+replication)*blockSize+blockCol;
418 const int srcCol = block*blockSize+blockCol;
419 for (int row=0; row<oldRows; row++) {
420 const int destRow = row+replication*shiftAfterEveryReplication;
421 if (destRow<_rows) {
422 _m[ serialise( destRow,destCol) ] = oldData[ serialise(row,srcCol,oldRows,oldCols) ];
423 }
424 }
425 }
426 }
427 }
428
430}
431
433 tarch::la::DynamicMatrix result(A.rows(),B.cols());
434
435 assertionEquals( A.cols(), B.rows() );
436
437 for (int row=0; row<A.rows(); row++)
438 for (int col=0; col<B.cols(); col++)
439 for (int i=0; i<A.cols(); i++) {
440 result(row,col) += A(row,i) * B(i,col);
441 }
442
443 return result;
444}
#define assertion2(expr, param0, param1)
#define assertion4(expr, param0, param1, param2, param3)
#define assertion3(expr, param0, param1, param2)
#define assertionEquals(lhs, rhs)
#define assertionMsg(expr, message)
#define assertion(expr)
And from this we can write down f$ nabla phi_i nabla phi_i dx but since we are constructing matrix let s investigate the f$ our matrix elements will nabla phi_i dx f By this will be a sparse as these basis functions are chosen to not overlap with each other almost everywhere In other they have only local support We can read off the right hand side values
And from this we can write down f$ nabla phi_i nabla phi_i dx but since we are constructing matrix let s investigate the f$ our matrix elements will nabla phi_i dx f By this will be a sparse matrix
tarch::la::DynamicMatrix kroneckerProduct(const tarch::la::DynamicMatrix &lhs, const tarch::la::DynamicMatrix &rhs)
Wrapper around static routine, so I don't have to use full-qualified name.
My standard matrix is a matrix where the size is fixed at compile time.
double & operator()(int row, int col)
~DynamicMatrix()
Free the array on the heap.
void multiply(double *result, double *x)
void batchedMultiplyAoS(double *__restrict__ result, const double *__restrict__ x, int batchCount, int resultSize, int firstRow)
This operation assumes that x holds a whole batch of vectors in AoS format.
void insertEmptyRows(int number, int where, int repeatEveryKColumns=0)
static std::string vectorToString(double *values, int entries, bool addLineBreaks=false)
I often need this in combination with the toString() operation above.
void removeColumn(int number)
void shiftColumnsRight(int shift, bool wrap=false)
DynamicMatrix & operator=(const DynamicMatrix &)=delete
bool operator==(double values[][Cols]) const
void shiftRowsDown(int shift, bool wrap=false)
Shift the rows to the right.
void insertEmptyColumns(int number, int where, int repeatEveryKColumns=0)
Insert zero columns.
void replicateRows(int blockSize, int numberOfReplications, int shiftAfterEveryReplication, bool extendColumnsToAccommodateShifts)
Split the matrix into blocks of rows of size blockSize.
DynamicMatrix(int rows, int cols)
Create empty matrix.
void replicateCols(int blockSize, int numberOfReplications, int shiftAfterEveryReplication, bool extendColumnsToAccommodateShifts)
std::string toString(bool addLineBreaks=false) const
static tarch::la::DynamicMatrix id(int rows)
Create (square) identify matrix with rows rows and column.
static int serialise(int row, int col, int Rows, int Cols)
void scale(double value)
Scale all entries.
Vector< Rows, Scalar > col(const Matrix< Rows, Cols, Scalar > &matrix, int whichColumn)
Extract row from matrix.
Matrix< Rows, Cols, Scalar > operator*(const Matrix< Rows, X, Scalar > &lhs, const Matrix< X, Cols, Scalar > &rhs)
bool equals(const Matrix< Rows, Cols, Scalar > &lhs, const Matrix< Rows, Cols, Scalar > &rhs, const Scalar &tolerance=NUMERICAL_ZERO_DIFFERENCE)
Compares to matrices on equality by means of a numerical accuracy.
Vector< Cols, Scalar > row(const Matrix< Rows, Cols, Scalar > &matrix, int whichRow)
Extract row from matrix.
std::string toString(MemoryLocation value)
void freeMemory(void *data, MemoryLocation location)
@ Heap
Create data on the heap of the local device.