NFFT Logo 3.2.2
error.c
00001 /*
00002  * Copyright (c) 2002, 2012 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 /* $Id: util.c 3483 2010-04-23 19:02:34Z keiner $ */
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 }

Generated on Fri Oct 12 2012 by Doxygen 1.8.0-20120409