NFFT 3.5.3alpha
error.c
1/*
2 * Copyright (c) 2002, 2017 Jens Keiner, Stefan Kunis, Daniel Potts
3 *
4 * This program is free software; you can redistribute it and/or modify it under
5 * the terms of the GNU General Public License as published by the Free Software
6 * Foundation; either version 2 of the License, or (at your option) any later
7 * version.
8 *
9 * This program is distributed in the hope that it will be useful, but WITHOUT
10 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
12 * details.
13 *
14 * You should have received a copy of the GNU General Public License along with
15 * this program; if not, write to the Free Software Foundation, Inc., 51
16 * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17 */
18
19#include "infft.h"
20
21static R cnrm1(const C *x, const INT n)
22{
23 INT k;
24 R nrm = K(0.0);
25
26 for (k = 0; k < n; k++)
27 nrm += CABS(x[k]);
28
29 return nrm;
30}
31
32static R nrm1(const R *x, const INT n)
33{
34 INT k;
35 R nrm = K(0.0);
36
37 for (k = 0; k < n; k++)
38 nrm += FABS(x[k]);
39
40 return nrm;
41}
42
43static R cerr2(const C *x, const C *y, const INT n)
44{
45 INT k;
46 R err = K(0.0);
47
48 if (!y)
49 {
50 for (k = 0; k < n; k++)
51 err += CONJ(x[k]) * x[k];
52 }
53 else
54 {
55 for (k = 0; k < n; k++)
56 err += CONJ(x[k]-y[k]) * (x[k]-y[k]);
57 }
58
59 return SQRT(err);
60}
61
62static R err2(const R *x, const R *y, const INT n)
63{
64 INT k;
65 R err = K(0.0);
66
67 if (!y)
68 {
69 for (k = 0; k < n; k++)
70 err += x[k]*x[k];
71 }
72 else
73 {
74 for (k = 0; k < n; k++)
75 err += (x[k]-y[k]) * (x[k]-y[k]);
76 }
77
78 return SQRT(err);
79}
80
81static R cerri(const C *x, const C *y, const INT n)
82{
83 INT k;
84 R err = K(0.0), t;
85
86 if (!y)
87 {
88 for (k = 0; k < n; k++)
89 {
90 t = CABS(x[k]);
91 err = MAX(err, t);
92 }
93 }
94 else
95 {
96 for (k = 0; k < n; k++)
97 {
98 t = CABS(x[k] - y[k]);
99 err = MAX(err, t);
100 }
101 }
102
103 return err;
104}
105
106static R erri(const R *x, const R *y, const INT n)
107{
108 INT k;
109 R err = K(0.0), t;
110
111 if (!y)
112 {
113 for (k = 0; k < n; k++)
114 {
115 t = FABS(x[k]);
116 err = MAX(err, t);
117 }
118 }
119 else
120 {
121 for (k = 0; k < n; k++)
122 {
123 t = FABS(x[k] - y[k]);
124 err = MAX(err, t);
125 }
126 }
127
128 return err;
129}
130
133R Y(error_l_infty_complex)(const C *x, const C *y, const INT n)
134{
135 return (cerri(x, y, n)/cerri(x, 0, n));
136}
137
140R Y(error_l_infty_double)(const R *x, const R *y, const INT n)
141{
142 return (erri(x, y, n)/erri(x, 0, n));
143}
144
147R Y(error_l_infty_1_complex)(const C *x, const C *y, const INT n,
148 const C *z, const INT m)
149{
150 return (cerri(x, y, n)/cnrm1(z, m));
151}
152
155R Y(error_l_infty_1_double)(const R *x, const R *y, const INT n, const R *z,
156 const INT m)
157{
158 return (erri(x, y, n)/nrm1(z, m));
159}
160
163R Y(error_l_2_complex)(const C *x, const C *y, const INT n)
164{
165 return (cerr2(x, y, n)/cerr2(x, 0, n));
166}
167
170R Y(error_l_2_double)(const R *x, const R *y, const INT n)
171{
172 return (err2(x, y, n)/err2(x, NULL, n));
173}
Internal header file for auxiliary definitions and functions.