![]() |
NFFT
3.3.1
|
00001 /* 00002 * Copyright (c) 2002, 2016 Jens Keiner, Stefan Kunis, Daniel Potts 00003 * 00004 * This program is free software; you can redistribute it and/or modify it under 00005 * the terms of the GNU General Public License as published by the Free Software 00006 * Foundation; either version 2 of the License, or (at your option) any later 00007 * version. 00008 * 00009 * This program is distributed in the hope that it will be useful, but WITHOUT 00010 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 00011 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 00012 * details. 00013 * 00014 * You should have received a copy of the GNU General Public License along with 00015 * this program; if not, write to the Free Software Foundation, Inc., 51 00016 * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 00017 */ 00018 00019 #include "infft.h" 00020 00021 static R cnrm1(const C *x, const INT n) 00022 { 00023 INT k; 00024 R nrm = K(0.0); 00025 00026 for (k = 0; k < n; k++) 00027 nrm += CABS(x[k]); 00028 00029 return nrm; 00030 } 00031 00032 static R nrm1(const R *x, const INT n) 00033 { 00034 INT k; 00035 R nrm = K(0.0); 00036 00037 for (k = 0; k < n; k++) 00038 nrm += FABS(x[k]); 00039 00040 return nrm; 00041 } 00042 00043 static R cerr2(const C *x, const C *y, const INT n) 00044 { 00045 INT k; 00046 R err = K(0.0); 00047 00048 if (!y) 00049 { 00050 for (k = 0; k < n; k++) 00051 err += CONJ(x[k]) * x[k]; 00052 } 00053 else 00054 { 00055 for (k = 0; k < n; k++) 00056 err += CONJ(x[k]-y[k]) * (x[k]-y[k]); 00057 } 00058 00059 return SQRT(err); 00060 } 00061 00062 static R err2(const R *x, const R *y, const INT n) 00063 { 00064 INT k; 00065 R err = K(0.0); 00066 00067 if (!y) 00068 { 00069 for (k = 0; k < n; k++) 00070 err += x[k]*x[k]; 00071 } 00072 else 00073 { 00074 for (k = 0; k < n; k++) 00075 err += (x[k]-y[k]) * (x[k]-y[k]); 00076 } 00077 00078 return SQRT(err); 00079 } 00080 00081 static R cerri(const C *x, const C *y, const INT n) 00082 { 00083 INT k; 00084 R err = K(0.0), t; 00085 00086 if (!y) 00087 { 00088 for (k = 0; k < n; k++) 00089 { 00090 t = CABS(x[k]); 00091 err = MAX(err, t); 00092 } 00093 } 00094 else 00095 { 00096 for (k = 0; k < n; k++) 00097 { 00098 t = CABS(x[k] - y[k]); 00099 err = MAX(err, t); 00100 } 00101 } 00102 00103 return err; 00104 } 00105 00106 static R erri(const R *x, const R *y, const INT n) 00107 { 00108 INT k; 00109 R err = K(0.0), t; 00110 00111 if (!y) 00112 { 00113 for (k = 0; k < n; k++) 00114 { 00115 t = FABS(x[k]); 00116 err = MAX(err, t); 00117 } 00118 } 00119 else 00120 { 00121 for (k = 0; k < n; k++) 00122 { 00123 t = FABS(x[k] - y[k]); 00124 err = MAX(err, t); 00125 } 00126 } 00127 00128 return err; 00129 } 00130 00133 R Y(error_l_infty_complex)(const C *x, const C *y, const INT n) 00134 { 00135 return (cerri(x, y, n)/cerri(x, 0, n)); 00136 } 00137 00140 R Y(error_l_infty_double)(const R *x, const R *y, const INT n) 00141 { 00142 return (erri(x, y, n)/erri(x, 0, n)); 00143 } 00144 00147 R Y(error_l_infty_1_complex)(const C *x, const C *y, const INT n, 00148 const C *z, const INT m) 00149 { 00150 return (cerri(x, y, n)/cnrm1(z, m)); 00151 } 00152 00155 R Y(error_l_infty_1_double)(const R *x, const R *y, const INT n, const R *z, 00156 const INT m) 00157 { 00158 return (erri(x, y, n)/nrm1(z, m)); 00159 } 00160 00163 R Y(error_l_2_complex)(const C *x, const C *y, const INT n) 00164 { 00165 return (cerr2(x, y, n)/cerr2(x, 0, n)); 00166 } 00167 00170 R Y(error_l_2_double)(const R *x, const R *y, const INT n) 00171 { 00172 return (err2(x, y, n)/err2(x, NULL, n)); 00173 }