00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "infft.h"
00022
00023 static R cnrm1(const C *x, const INT n)
00024 {
00025 INT k;
00026 R nrm = K(0.0);
00027
00028 for (k = 0; k < n; k++)
00029 nrm += CABS(x[k]);
00030
00031 return nrm;
00032 }
00033
00034 static R nrm1(const R *x, const INT n)
00035 {
00036 int k;
00037 R nrm = K(0.0);
00038
00039 for (k = 0; k < n; k++)
00040 nrm += FABS(x[k]);
00041
00042 return nrm;
00043 }
00044
00045 static R cerr2(const C *x, const C *y, const INT n)
00046 {
00047 int k;
00048 R err = K(0.0);
00049
00050 if (!y)
00051 {
00052 for (k = 0; k < n; k++)
00053 err += CONJ(x[k]) * x[k];
00054 }
00055 else
00056 {
00057 for (k = 0; k < n; k++)
00058 err += CONJ(x[k]-y[k]) * (x[k]-y[k]);
00059 }
00060
00061 return SQRT(err);
00062 }
00063
00064 static R err2(const R *x, const R *y, const INT n)
00065 {
00066 int k;
00067 R err = K(0.0);
00068
00069 if (!y)
00070 {
00071 for (k = 0; k < n; k++)
00072 err += x[k]*x[k];
00073 }
00074 else
00075 {
00076 for (k = 0; k < n; k++)
00077 err += (x[k]-y[k]) * (x[k]-y[k]);
00078 }
00079
00080 return SQRT(err);
00081 }
00082
00083 static R cerri(const C *x, const C *y, const INT n)
00084 {
00085 int k;
00086 R err = K(0.0), t;
00087
00088 if (!y)
00089 {
00090 for (k = 0; k < n; k++)
00091 {
00092 t = CABS(x[k]);
00093 err = MAX(err, t);
00094 }
00095 }
00096 else
00097 {
00098 for (k = 0; k < n; k++)
00099 {
00100 t = CABS(x[k] - y[k]);
00101 err = MAX(err, t);
00102 }
00103 }
00104
00105 return err;
00106 }
00107
00108 static R erri(const R *x, const R *y, const INT n)
00109 {
00110 int k;
00111 R err = K(0.0), t;
00112
00113 if (!y)
00114 {
00115 for (k = 0; k < n; k++)
00116 {
00117 t = FABS(x[k]);
00118 err = MAX(err, t);
00119 }
00120 }
00121 else
00122 {
00123 for (k = 0; k < n; k++)
00124 {
00125 t = FABS(x[k] - y[k]);
00126 err = MAX(err, t);
00127 }
00128 }
00129
00130 return err;
00131 }
00132
00135 R X(error_l_infty_complex)(const C *x, const C *y, const INT n)
00136 {
00137 return (cerri(x, y, n)/cerri(x, 0, n));
00138 }
00139
00142 R X(error_l_infty_double)(const R *x, const R *y, const INT n)
00143 {
00144 return (erri(x, y, n)/erri(x, 0, n));
00145 }
00146
00149 R X(error_l_infty_1_complex)(const C *x, const C *y, const INT n,
00150 const C *z, const INT m)
00151 {
00152 return (cerri(x, y, n)/cnrm1(z, m));
00153 }
00154
00157 R X(error_l_infty_1_double)(const R *x, const R *y, const INT n, const R *z,
00158 const INT m)
00159 {
00160 return (erri(x, y, n)/nrm1(z, m));
00161 }
00162
00165 R X(error_l_2_complex)(const C *x, const C *y, const INT n)
00166 {
00167 return (cerr2(x, y, n)/cerr2(x, 0, n));
00168 }
00169
00172 R X(error_l_2_double)(const R *x, const R *y, const INT n)
00173 {
00174 return (err2(x, y, n)/err2(x, NULL, n));
00175 }