18300593117

Committed 06 Oct 2025 10:47PM UTC coverage: 74.862% (-0.2%) from 75.031%

Build # 18300593117

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #759 from wdhawkins/performance

Sweep performance optimizations

Coverage Stats

294 of 302 new or added lines in 15 files covered. (97.35%)

334 existing lines in 80 files now uncovered.

17788 of 23761 relevant lines covered (74.86%)

61852783.95 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

86.31

/framework/data_types/dense_matrix.h

1	// SPDX-FileCopyrightText: 2024 The OpenSn Authors <https://open-sn.github.io/opensn/>
2	// SPDX-License-Identifier: MIT
3
4	#pragma once
5
6	#include "framework/data_types/ndarray.h"
7	#include "framework/data_types/vector.h"
8
9	namespace opensn
10	{
11
12	template <typename TYPE>
13	class DenseMatrix : public NDArray<TYPE, 2>
14	{
15	public:
16	/// Create an empty dense matrix
17	DenseMatrix() : NDArray<TYPE, 2>() {}
18
19	/**
20	* Create a dense matrix with specified number of rows and columns
21	*
22	* \param rows Number of rows
23	* \param cols Number of columns
24	*/
25	DenseMatrix(size_t rows, size_t cols) : NDArray<TYPE, 2>({rows, cols}) {}	58,189,553✔
26
27	/**
28	* Create a dense matrix with specified number of rows and columns and initialize the element to a
29	* given value
30	*
31	* \param rows Number of rows
32	* \param cols Number of columns
33	* \param intitial_value Value to initialize the matrix elements with
34	*/
35	DenseMatrix(size_t rows, size_t cols, TYPE intitial_value)	9,188,035✔
36	: NDArray<TYPE, 2>({rows, cols}, intitial_value)	9,188,035✔
37	{
38	}	9,188,035✔
39
40	/// Return the number of rows
41	size_t Rows() const { return this->dimension()[0]; }	26,549✔
42
43	/// Return the number of columns
44	size_t Columns() const { return this->dimension()[1]; }	91,761✔
45
46	/// Set the elements of the matrix to a specified value
47	void Set(TYPE val) { this->set(val); }	601,903✔
48
49	/// Set the diagonal of the matrix
50	void SetDiagonal(TYPE val)	36✔
51	{
52	auto dims = this->dimension();	36✔
53	auto d = std::min(dims[0], dims[1]);	36✔
54	for (int i = 0; i < d; ++i)	5,754✔
55	(*this)(i, i) = val;	11,436✔
56	}	36✔
57
58	void SetRow(int row, const Vector<TYPE>& values)
59	{
60	assert(Columns() == values.Rows());
61	for (size_t i = 0; i < Columns(); ++i)
62	(*this)(row, i) = values(i);
63	}
64
65	void Add(const DenseMatrix<TYPE>& other)
66	{
67	assert(Rows() == other.Rows());
68	assert(Columns() == other.Columns());
69	for (auto i = 0; i < Rows(); ++i)
70	for (auto j = 0; j < Columns(); ++j)
71	(*this)(i, j) += other(i, j);
72	}
73
74	void Subtract(const DenseMatrix<TYPE>& other)
75	{
76	assert(Rows() == other.Rows());
77	assert(Columns() == other.Columns());
78	for (auto i = 0; i < Rows(); ++i)
79	for (auto j = 0; j < Columns(); ++j)
80	(*this)(i, j) -= other(i, j);
81	}
82
83	/// Multiply matrix with a vector and return resulting vector
84	Vector<TYPE> Mult(const Vector<TYPE>& x) const
85	{
86	auto rows = Rows();
87	auto cols = x.Rows();
88
89	assert(rows > 0);
90	assert(cols == Columns());
91
92	Vector<TYPE> b(rows, 0.0);
93	for (auto i = 0; i < rows; ++i)
94	for (auto j = 0; j < cols; ++j)
95	b(i) += (this)(i, j) x(j);
96	return b;
97	}
98
99	/// Multiply by a matrix and return the resulting matrix
100	DenseMatrix<TYPE> Mult(const DenseMatrix<TYPE>& other) const
101	{
102	auto rows = Rows();
103	assert(rows != 0 and other.Rows() != 0);
104	auto columns = Columns();
105	auto other_cols = other.Columns();
106	assert(columns != 0 and other_cols != 0 and columns == other.Rows());
107	DenseMatrix<TYPE> res(rows, other_cols, 0.);
108	for (auto i = 0; i < rows; ++i)
109	for (auto j = 0; j < other_cols; ++j)
110	for (auto k = 0; k < columns; ++k)
111	res(i, j) += (this)(i, k) other(k, j);
112	return res;
113	}
114
115	/// Returns the transpose of a matrix
116	DenseMatrix<TYPE> Transposed() const
117	{
118	assert(Rows());
119	assert(Columns());
120	DenseMatrix<TYPE> trans(Columns(), Rows());
121	for (auto i = 0; i < Rows(); ++i)
122	for (auto j = 0; j < Columns(); ++j)
123	trans(j, i) = (*this)(i, j);
124	return trans;
125	}
126
127	/// Transpose this matrix
128	void Transpose()
129	{
130	assert(Rows());
131	assert(Columns());
132	DenseMatrix<TYPE> trans(Columns(), Rows());
133	for (auto i = 0; i < Rows(); ++i)
134	for (auto j = 0; j < Columns(); ++j)
135	trans(j, i) = (*this)(i, j);
136	(*this) = trans;
137	}
138
139	/// Prints the matrix to a string and then returns the string.
140	std::string PrintStr() const
141	{
142	std::stringstream out;
143
144	for (int i = 0; i < Rows(); ++i)
145	{
146	for (int j = 0; j < (Columns() - 1); ++j)
147	out << (*this)(i, j) << " ";
148	out << (*this)(i, Columns() - 1);
149
150	if (i < (Rows() - 1))
151	out << "\n";
152	}
153
154	return out.str();
155	}
156	};
157
158	/// Returns the transpose of a matrix.
159	template <typename TYPE>
160	DenseMatrix<TYPE>
161	Transpose(const DenseMatrix<TYPE>& A)
162	{
163	assert(A.Rows());
164	assert(A.Columns());
165	auto AR = A.Rows();
166	auto AC = 0;
167	if (AR)
168	AC = A.Columns();
169
170	DenseMatrix<TYPE> T(AC, AR);
171	for (auto i = 0; i < AR; ++i)
172	for (auto j = 0; j < AC; ++j)
173	T(j, i) = A(i, j);
174	return T;
175	}
176
177	/// Swaps two rows of a matrix.
178	template <typename TYPE>
179	void
180	SwapRows(DenseMatrix<TYPE>& A, size_t r1, size_t r2)	1✔
181	{
182	auto rows = A.Rows();	1✔
183	auto cols = A.Columns();	1✔
184	assert(rows > 0);	1✔
185	assert(cols > 0);	1✔
186	assert(r1 >= 0 and r1 < rows and r2 >= 0 and r2 < rows);	1✔
187
188	for (auto j = 0; j < cols; ++j)	5✔
189	std::swap(A(r1, j), A(r2, j));	12✔
190	}	1✔
191
192	/// Multiply matrix with a constant and return result.
193	template <typename TYPE>
194	DenseMatrix<TYPE>
195	Mult(const DenseMatrix<TYPE>& A, const TYPE c)
196	{
197	auto R = A.Rows();
198	auto C = A.Columns();
199	DenseMatrix<TYPE> B(R, C, 0.);
200	for (auto i = 0; i < R; ++i)
201	for (auto j = 0; j < C; ++j)
202	B(i, j) = A(i, j) * c;
203	return B;
204	}
205
206	/// Multiply matrix with a vector and return resulting vector
207	template <typename TYPE>
208	Vector<TYPE>
209	Mult(const DenseMatrix<TYPE>& A, const Vector<TYPE>& x)	11,044✔
210	{
211	auto R = A.Rows();	11,044✔
212	auto C = x.Rows();	11,044✔
213
214	assert(R > 0);	11,044✔
215	assert(C == A.Columns());	11,044✔
216
217	Vector<TYPE> b(R, 0.0);	11,044✔
218	for (auto i = 0; i < R; ++i)	643,486✔
219	{
220	for (auto j = 0; j < C; ++j)	101,991,234✔
221	b(i) += A(i, j) * x(j);	405,435,168✔
222	}
223
224	return b;	11,044✔
225	}
226
227	/// Mutliply two matrices and return result.
228	template <typename TYPE>
229	DenseMatrix<TYPE>
230	Mult(const DenseMatrix<TYPE>& A, const DenseMatrix<TYPE>& B)	35✔
231	{
232	auto AR = A.Rows();	35✔
233
234	assert(AR != 0 and B.Rows() != 0);	35✔
235
236	auto AC = A.Columns();	35✔
237	auto BC = B.Columns();	35✔
238
239	assert(AC != 0 and BC != 0 and AC == B.Rows());	35✔
240
241	auto CR = AR;	35✔
242	auto CC = BC;	35✔
243	auto Cs = AC;	35✔
244
245	DenseMatrix<TYPE> C(CR, CC, 0.);	35✔
246	for (auto i = 0; i < CR; ++i)	5,749✔
247	for (auto j = 0; j < CC; ++j)	965,336✔
248	for (auto k = 0; k < Cs; ++k)	162,175,134✔
249	C(i, j) += A(i, k) * B(k, j);	644,862,048✔
250	return C;	35✔
251	}
252
253	/// Adds two matrices and returns the result.
254	template <typename TYPE>
255	DenseMatrix<TYPE>
256	Add(const DenseMatrix<TYPE>& A, const DenseMatrix<TYPE>& B)
257	{
258	auto AR = A.Rows();
259	auto BR = B.Rows();
260
261	assert(AR != 0 and B.Rows() != 0);
262	assert(AR == BR);
263
264	auto AC = A.Columns();
265	auto BC = B.Columns();
266
267	assert(AC != 0 and BC != 0);
268	assert(AC == BC);
269
270	DenseMatrix<TYPE> C(AR, AC, 0.0);
271	for (auto i = 0; i < AR; ++i)
272	for (auto j = 0; j < AC; ++j)
273	C(i, j) = A(i, j) + B(i, j);
274	return C;
275	}
276
277	/// Subtracts matrix A from B and returns the result.
278	template <typename TYPE>
279	DenseMatrix<TYPE>
280	Subtract(const DenseMatrix<TYPE>& A, const DenseMatrix<TYPE>& B)
281	{
282	auto AR = A.Rows();
283	auto BR = B.Rows();
284
285	assert(AR != 0 and BR != 0);
286	assert(AR == BR);
287
288	auto AC = A.Columns();
289	auto BC = B.Columns();
290
291	assert(AC != 0 and BC != 0);
292	assert(AC == BC);
293
294	DenseMatrix<TYPE> C(AR, AC, 0.0);
295	for (auto i = 0; i < AR; ++i)
296	for (auto j = 0; j < AC; ++j)
297	C(i, j) = A(i, j) - B(i, j);
298	return C;
299	}
300
301	/// Scale the matrix with a constant value
302	template <typename TYPE>
303	void
304	Scale(DenseMatrix<TYPE>& mat, TYPE alpha)	1,899✔
305	{
306	for (auto i = 0; i < mat.Rows(); ++i)	9,493✔
307	for (auto j = 0; j < mat.Columns(); ++j)	37,970✔
308	mat(i, j) *= alpha;	60,752✔
309	}	1,899✔
310
311	/// Scale the matrix with a constant value
312	template <typename TYPE>
313	DenseMatrix<TYPE>
314	Scaled(const DenseMatrix<TYPE>& mat, TYPE alpha)
315	{
316	DenseMatrix<TYPE> res(mat.Rows(), mat.Columns());
317	for (auto i = 0; i < mat.Rows(); ++i)
318	for (auto j = 0; j < mat.Columns(); ++j)
319	res(i, j) = mat(i, j) * alpha;
320	return res;
321	}
322
323	/// Returns a copy of A with removed rows `r` and removed columns `c`
324	template <typename TYPE>
325	DenseMatrix<TYPE>
326	SubMatrix(const DenseMatrix<TYPE>& A, const size_t r, const size_t c)	×
327	{
328	auto rows = A.Rows();	×
329	auto cols = A.Columns();	×
330	assert((r >= 0) and (r < rows) and (c >= 0) and (c < cols));	×
331
332	DenseMatrix<TYPE> B(rows - 1, cols - 1);	×
333	for (size_t i = 0, ii = 0; i < rows; ++i)	×
334	{
335	if (i != r)	×
336	{
337	for (size_t j = 0, jj = 0; j < cols; ++j)	×
338	{
339	if (j != c)	×
340	{
341	B(ii, jj) = A(i, j);	×
342	++jj;	×
343	}
344	}
345	++ii;	×
346	}
347	}
348	return B;	×
349	}
350
351	/// Computes the determinant of a matrix.
352	template <typename TYPE>
353	double
354	Determinant(const DenseMatrix<TYPE>& A)	1,899✔
355	{
356	auto rows = A.Rows();	1,899✔
357
358	if (rows == 1)	1,899✔
359	return A(0, 0);	×
360	else if (rows == 2)	1,899✔
361	{
362	return A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0);	×
363	}
364	else if (rows == 3)	1,899✔
365	{
366	return A(0, 0) * A(1, 1) * A(2, 2) + A(0, 1) * A(1, 2) * A(2, 0) + A(0, 2) * A(1, 0) * A(2, 1) -	9✔
367	A(0, 0) * A(1, 2) * A(2, 1) - A(0, 1) * A(1, 0) * A(2, 2) - A(0, 2) * A(1, 1) * A(2, 0);	10✔
368	}
369	// http://www.cvl.iis.u-tokyo.ac.jp/~Aiyazaki/tech/teche23.htAl
370	else if (rows == 4)	1,898✔
371	{
372	return A(0, 0) * A(1, 1) * A(2, 2) * A(3, 3) + A(0, 0) * A(1, 2) * A(2, 3) * A(3, 1) +	15,184✔
373	A(0, 0) * A(1, 3) * A(2, 1) * A(3, 2) + A(0, 1) * A(1, 0) * A(2, 3) * A(3, 2) +	15,184✔
374	A(0, 1) * A(1, 2) * A(2, 0) * A(3, 3) + A(0, 1) * A(1, 3) * A(2, 2) * A(3, 0) +	15,184✔
375	A(0, 2) * A(1, 0) * A(2, 1) * A(3, 3) + A(0, 2) * A(1, 1) * A(2, 3) * A(3, 0) +	15,184✔
376	A(0, 2) * A(1, 3) * A(2, 0) * A(3, 1) + A(0, 3) * A(1, 0) * A(2, 2) * A(3, 1) +	15,184✔
377	A(0, 3) * A(1, 1) * A(2, 0) * A(3, 2) + A(0, 3) * A(1, 2) * A(2, 1) * A(3, 0) -	15,184✔
378	A(0, 0) * A(1, 1) * A(2, 3) * A(3, 2) - A(0, 0) * A(1, 2) * A(2, 1) * A(3, 3) -	15,184✔
379	A(0, 0) * A(1, 3) * A(2, 2) * A(3, 1) - A(0, 1) * A(1, 0) * A(2, 2) * A(3, 3) -	15,184✔
380	A(0, 1) * A(1, 2) * A(2, 3) * A(3, 0) - A(0, 1) * A(1, 3) * A(2, 0) * A(3, 2) -	15,184✔
381	A(0, 2) * A(1, 0) * A(2, 3) * A(3, 1) - A(0, 2) * A(1, 1) * A(2, 0) * A(3, 3) -	15,184✔
382	A(0, 2) * A(1, 3) * A(2, 1) * A(3, 0) - A(0, 3) * A(1, 0) * A(2, 1) * A(3, 2) -	15,184✔
383	A(0, 3) * A(1, 1) * A(2, 2) * A(3, 0) - A(0, 3) * A(1, 2) * A(2, 0) * A(3, 1);	17,082✔
384	}
385	else
386	{
UNCOV 387	double det = 0;	×
388	for (auto n = 0; n < rows; ++n)	×
389	{
390	auto M = SubMatrix(A, 0, n);	×
391	double pm = ((n + 1) % 2) * 2.0 - 1.0;	×
392	det += pm * A(0, n) * Determinant(M);	×
393	}
UNCOV 394	return det;	×
395	}
396	}
397
398	/// Gauss Elimination without pivoting.
399	template <typename TYPE>
400	void
401	GaussElimination(DenseMatrix<TYPE>& A, Vector<TYPE>& b, size_t n)
402	{
403	// Forward elimination
404	for (auto i = 0; i < n - 1; ++i)
405	{
406	auto bi = b(i);
407	auto factor = 1.0 / A(i, i);
408	for (size_t j = i + 1; j < n; ++j)
409	{
410	auto val = A(j, i) * factor;
411	b(j) -= val * bi;
412	for (size_t k = i + 1; k < n; ++k)
413	A(j, k) -= val * A(i, k);
414	}
415	}
416
417	// Back substitution
418	for (int i = n - 1; i >= 0; --i)
419	{
420	auto bi = b(i);
421	for (size_t j = i + 1; j < n; ++j)
422	bi -= A(i, j) * b(j);
423	b(i) = bi / A(i, i);
424	}
425	}
426
427	/// Computes the inverse of a matrix using Gauss-Elimination with pivoting.
428	template <typename TYPE>
429	DenseMatrix<TYPE>
430	InverseGEPivoting(const DenseMatrix<TYPE>& A)	34✔
431	{
432	assert(A.Rows() == A.Columns());	34✔
433
434	const auto rows = A.Rows();	34✔
435	const auto cols = A.Columns();	34✔
436
437	DenseMatrix<TYPE> M(rows, cols);	34✔
438	M.Set(0.);	34✔
439	M.SetDiagonal(1.);	34✔
440
441	auto B = A;	34✔
442
443	for (auto c = 0; c < rows; ++c)	5,746✔
444	{
445	// Find a row with the largest pivot value
446	auto max_row = c; // nzr = non-zero row	5,712✔
447	for (auto r = c; r < rows; ++r)	488,376✔
448	if (std::fabs(B(r, c)) > std::fabs(B(max_row, c)))	965,328✔
449	max_row = r;	×
450
451	if (max_row != c)	5,712✔
452	{
453	SwapRows(B, max_row, c);	×
454	SwapRows(M, max_row, c);	×
455	}
456
457	// Eliminate non-zero values
458	for (auto r = 0; r < rows; ++r)	965,328✔
459	{
460	if (r != c)	959,616✔
461	{
462	auto g = B(r, c) / B(c, c);	1,907,808✔
463	if (B(r, c) != 0)	953,904✔
464	{
465	for (auto k = 0; k < rows; ++k)	×
466	{
467	B(r, k) -= B(c, k) * g;	×
468	M(r, k) -= M(c, k) * g;	×
469	}
470	B(r, c) = 0;	×
471	}
472	}
473	else
474	{
475	auto g = 1 / B(c, c);	5,712✔
476	for (auto k = 0; k < rows; ++k)	965,328✔
477	{
478	B(r, k) *= g;	959,616✔
479	M(r, k) *= g;	1,919,232✔
480	}
481	}
482	}
483	}
484	return M;	34✔
485	}	34✔
486
487	/// Computes the inverse of a matrix.
488	template <typename TYPE>
489	DenseMatrix<TYPE>
490	Inverse(const DenseMatrix<TYPE>& A)	1,932✔
491	{
492	auto rows = A.Rows();	1,932✔
493	DenseMatrix<double> M(rows, A.Rows());	1,932✔
494	double f = 0.0;	1,932✔
495
496	// Only calculate the determinant if matrix size is less than 5 since
497	// the inverse is directly calculated for larger matrices. Otherwise,
498	// the inverse routine spends all of its time sitting in the determinant
499	// function which is unnecessary.
500	if (rows < 5)	1,932✔
501	{
502	f = Determinant(A);	1,898✔
503	assert(f != 0.0);	1,898✔
504	f = 1.0 / f;	1,898✔
505	}
506
507	if (rows == 1)	1,932✔
508	M(0, 0) = f;	×
509	else if (rows == 2)	1,932✔
510	{
511	M(0, 0) = A(1, 1);	×
512	M(0, 1) = -A(0, 1);	×
513	M(1, 0) = -A(1, 0);	×
514	M(1, 1) = A(0, 0);	×
515	Scale(M, f);	×
516	}
517	else if (rows == 3)	1,932✔
518	{
519	M(0, 0) = A(2, 2) * A(1, 1) - A(2, 1) * A(1, 2);	×
520	M(0, 1) = -(A(2, 2) * A(0, 1) - A(2, 1) * A(0, 2));	×
521	M(0, 2) = A(1, 2) * A(0, 1) - A(1, 1) * A(0, 2);	×
522	M(1, 0) = -(A(2, 2) * A(1, 0) - A(2, 0) * A(1, 2));	×
523	M(1, 1) = A(2, 2) * A(0, 0) - A(2, 0) * A(0, 2);	×
524	M(1, 2) = -(A(1, 2) * A(0, 0) - A(1, 0) * A(0, 2));	×
525	M(2, 0) = A(2, 1) * A(1, 0) - A(2, 0) * A(1, 1);	×
526	M(2, 1) = -(A(2, 1) * A(0, 0) - A(2, 0) * A(0, 1));	×
527	M(2, 2) = A(1, 1) * A(0, 0) - A(1, 0) * A(0, 1);	×
528	Scale(M, f);	×
529	}
530	else if (rows == 4)	1,932✔
531	{
532	// http://www.cvl.iis.u-tokyo.ac.jp/~Aiyazaki/tech/teche23.htAl
533	M(0, 0) = A(1, 1) * A(2, 2) * A(3, 3) + A(1, 2) * A(2, 3) * A(3, 1) +	11,388✔
534	A(1, 3) * A(2, 1) * A(3, 2) - A(1, 1) * A(2, 3) * A(3, 2) -	11,388✔
535	A(1, 2) * A(2, 1) * A(3, 3) - A(1, 3) * A(2, 2) * A(3, 1);	13,286✔
536	M(0, 1) = A(0, 1) * A(2, 3) * A(3, 2) + A(0, 2) * A(2, 1) * A(3, 3) +	11,388✔
537	A(0, 3) * A(2, 2) * A(3, 1) - A(0, 1) * A(2, 2) * A(3, 3) -	11,388✔
538	A(0, 2) * A(2, 3) * A(3, 1) - A(0, 3) * A(2, 1) * A(3, 2);	13,286✔
539	M(0, 2) = A(0, 1) * A(1, 2) * A(3, 3) + A(0, 2) * A(1, 3) * A(3, 1) +	11,388✔
540	A(0, 3) * A(1, 1) * A(3, 2) - A(0, 1) * A(1, 3) * A(3, 2) -	11,388✔
541	A(0, 2) * A(1, 1) * A(3, 3) - A(0, 3) * A(1, 2) * A(3, 1);	13,286✔
542	M(0, 3) = A(0, 1) * A(1, 3) * A(2, 2) + A(0, 2) * A(1, 1) * A(2, 3) +	11,388✔
543	A(0, 3) * A(1, 2) * A(2, 1) - A(0, 1) * A(1, 2) * A(2, 3) -	11,388✔
544	A(0, 2) * A(1, 3) * A(2, 1) - A(0, 3) * A(1, 1) * A(2, 2);	13,286✔
545
546	M(1, 0) = A(1, 0) * A(2, 3) * A(3, 2) + A(1, 2) * A(2, 0) * A(3, 3) +	11,388✔
547	A(1, 3) * A(2, 2) * A(3, 0) - A(1, 0) * A(2, 2) * A(3, 3) -	11,388✔
548	A(1, 2) * A(2, 3) * A(3, 0) - A(1, 3) * A(2, 0) * A(3, 2);	13,286✔
549	M(1, 1) = A(0, 0) * A(2, 2) * A(3, 3) + A(0, 2) * A(2, 3) * A(3, 0) +	11,388✔
550	A(0, 3) * A(2, 0) * A(3, 2) - A(0, 0) * A(2, 3) * A(3, 2) -	11,388✔
551	A(0, 2) * A(2, 0) * A(3, 3) - A(0, 3) * A(2, 2) * A(3, 0);	13,286✔
552	M(1, 2) = A(0, 0) * A(1, 3) * A(3, 2) + A(0, 2) * A(1, 0) * A(3, 3) +	11,388✔
553	A(0, 3) * A(1, 2) * A(3, 0) - A(0, 0) * A(1, 2) * A(3, 3) -	11,388✔
554	A(0, 2) * A(1, 3) * A(3, 0) - A(0, 3) * A(1, 0) * A(3, 2);	13,286✔
555	M(1, 3) = A(0, 0) * A(1, 2) * A(2, 3) + A(0, 2) * A(1, 3) * A(2, 0) +	11,388✔
556	A(0, 3) * A(1, 0) * A(2, 2) - A(0, 0) * A(1, 3) * A(2, 2) -	11,388✔
557	A(0, 2) * A(1, 0) * A(2, 3) - A(0, 3) * A(1, 2) * A(2, 0);	13,286✔
558
559	M(2, 0) = A(1, 0) * A(2, 1) * A(3, 3) + A(1, 1) * A(2, 3) * A(3, 0) +	11,388✔
560	A(1, 3) * A(2, 0) * A(3, 1) - A(1, 0) * A(2, 3) * A(3, 1) -	11,388✔
561	A(1, 1) * A(2, 0) * A(3, 3) - A(1, 3) * A(2, 1) * A(3, 0);	13,286✔
562	M(2, 1) = A(0, 0) * A(2, 3) * A(3, 1) + A(0, 1) * A(2, 0) * A(3, 3) +	11,388✔
563	A(0, 3) * A(2, 1) * A(3, 0) - A(0, 0) * A(2, 1) * A(3, 3) -	11,388✔
564	A(0, 1) * A(2, 3) * A(3, 0) - A(0, 3) * A(2, 0) * A(3, 1);	13,286✔
565	M(2, 2) = A(0, 0) * A(1, 1) * A(3, 3) + A(0, 1) * A(1, 3) * A(3, 0) +	11,388✔
566	A(0, 3) * A(1, 0) * A(3, 1) - A(0, 0) * A(1, 3) * A(3, 1) -	11,388✔
567	A(0, 1) * A(1, 0) * A(3, 3) - A(0, 3) * A(1, 1) * A(3, 0);	13,286✔
568	M(2, 3) = A(0, 0) * A(1, 3) * A(2, 1) + A(0, 1) * A(1, 0) * A(2, 3) +	11,388✔
569	A(0, 3) * A(1, 1) * A(2, 0) - A(0, 0) * A(1, 1) * A(2, 3) -	11,388✔
570	A(0, 1) * A(1, 3) * A(2, 0) - A(0, 3) * A(1, 0) * A(2, 1);	13,286✔
571
572	M(3, 0) = A(1, 0) * A(2, 2) * A(3, 1) + A(1, 1) * A(2, 0) * A(3, 2) +	11,388✔
573	A(1, 2) * A(2, 1) * A(3, 0) - A(1, 0) * A(2, 1) * A(3, 2) -	11,388✔
574	A(1, 1) * A(2, 2) * A(3, 0) - A(1, 2) * A(2, 0) * A(3, 1);	13,286✔
575	M(3, 1) = A(0, 0) * A(2, 1) * A(3, 2) + A(0, 1) * A(2, 2) * A(3, 0) +	11,388✔
576	A(0, 2) * A(2, 0) * A(3, 1) - A(0, 0) * A(2, 2) * A(3, 1) -	11,388✔
577	A(0, 1) * A(2, 0) * A(3, 2) - A(0, 2) * A(2, 1) * A(3, 0);	13,286✔
578	M(3, 2) = A(0, 0) * A(1, 2) * A(3, 1) + A(0, 1) * A(1, 0) * A(3, 2) +	11,388✔
579	A(0, 2) * A(1, 1) * A(3, 0) - A(0, 0) * A(1, 1) * A(3, 2) -	11,388✔
580	A(0, 1) * A(1, 2) * A(3, 0) - A(0, 2) * A(1, 0) * A(3, 1);	13,286✔
581	M(3, 3) = A(0, 0) * A(1, 1) * A(2, 2) + A(0, 1) * A(1, 2) * A(2, 0) +	11,388✔
582	A(0, 2) * A(1, 0) * A(2, 1) - A(0, 0) * A(1, 2) * A(2, 1) -	11,388✔
583	A(0, 1) * A(1, 0) * A(2, 2) - A(0, 2) * A(1, 1) * A(2, 0);	13,286✔
584	Scale(M, f);	1,898✔
585	}
586	else
587	M = InverseGEPivoting(A);	34✔
588
589	return M;	1,932✔
590	}	×
591
592	/**
593	* Performs power iteration to obtain the fundamental eigen mode. The eigen-value of the fundamental
594	* mode is return whilst the eigen-vector is return via reference.
595	*/
596	template <typename TYPE>
597	double
598	PowerIteration(const DenseMatrix<TYPE>& A,	34✔
599	Vector<TYPE>& e_vec,
600	int max_it = 2000,
601	double tol = 1.0e-13)
602	{
603	auto n = A.Rows();	34✔
604	int it_counter = 0;	34✔
605	Vector<double> y(n, 1.0);	34✔
606	double lambda0 = 0.0;	34✔
607
608	// Perform initial iteration outside of loop
609	auto Ay = Mult(A, y);	34✔
610	auto lambda = Dot(y, Ay);	34✔
611	y = Scaled(Ay, 1.0 / Vec2Norm(Ay));	34✔
612	if (lambda < 0.0)	34✔
613	Scale(y, -1.0);	×
614
615	// Perform convergence loop
616	bool converged = false;	34✔
617	while (not converged and it_counter < max_it)	3,587✔
618	{
619	// Update old eigenvalue
620	lambda0 = std::fabs(lambda);	3,553✔
621	// Calculate new eigenvalue/eigenvector
622	Ay = Mult(A, y);	3,553✔
623	lambda = Dot(y, Ay);	3,553✔
624	y = Scaled(Ay, 1.0 / Vec2Norm(Ay));	3,553✔
625
626	// Check if converged or not
627	if (std::fabs(std::fabs(lambda) - lambda0) <= tol)	3,553✔
628	converged = true;	34✔
629	// Update counter
630	++it_counter;	3,553✔
631	}
632
633	if (lambda < 0.0)	34✔
634	Scale(y, -1.0);	×
635
636	// Renormalize eigenvector for the last time
637	y = Scaled(Ay, 1.0 / lambda);	34✔
638
639	// Set eigenvector, return the eigenvalue
640	e_vec = std::move(y);	34✔
641
642	return lambda;	34✔
643	}	34✔
644		7,089✔
645	} // namespace opensn	7,089✔

Open-Sn / opensn / 18300593117

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous