1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
|
// Copyright (c) 2017-2023, University of Tennessee. All rights reserved.
// SPDX-License-Identifier: BSD-3-Clause
// This program is free software: you can redistribute it and/or modify it under
// the terms of the BSD 3-Clause license. See the accompanying LICENSE file.
#ifndef ERROR_HH
#define ERROR_HH
#include "blas.hh"
#include <vector>
// -----------------------------------------------------------------------------
// returns absolute error, || x - xref ||_2
// TODO: generalize to more arbitrary iterators than std::vector.
// TODO: use LAPACK's lassq algorithm to avoid numerical issues in summing squares.
template< typename T1, typename T2 >
blas::real_type< T1, T2 >
abs_error( std::vector<T1>& x, std::vector<T2>& xref )
{
using real_t = blas::real_type< T1, T2 >;
if (x.size() != xref.size()) {
return std::numeric_limits<real_t>::quiet_NaN();
}
real_t tmp;
real_t diff = 0;
for (size_t i = 0; i < x.size(); ++i) {
tmp = std::abs( x[i] - xref[i] );
diff += tmp*tmp;
}
diff = sqrt( diff );
return diff;
}
// -----------------------------------------------------------------------------
// returns relative error, || x - xref ||_2 / || xref ||_2
template< typename T1, typename T2 >
blas::real_type< T1, T2 >
rel_error( std::vector<T1>& x, std::vector<T2>& xref )
{
using real_t = blas::real_type< T1, T2 >;
if (x.size() != xref.size()) {
return std::numeric_limits<real_t>::quiet_NaN();
}
real_t tmp;
real_t diff = 0;
real_t norm = 0;
for (size_t i = 0; i < x.size(); ++i) {
tmp = std::abs( x[i] - xref[i] );
diff += tmp*tmp;
tmp = std::abs( xref[i] );
norm += tmp*tmp;
}
diff = sqrt( diff );
norm = sqrt( norm );
return diff / norm;
}
#endif // #ifndef ERROR_HH
|
