Mercurial > forge

annotate extra/NaN/src/sumskipnan_mex.cpp @ 12691:6d6285a2a633 octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
use macro ISNAN() instead of C++'s isnan() - because it supports all floating point formats not just double
author schloegl
date Sat, 12 Sep 2015 14:16:39 +0000
parents f26b1170ea90
children 13815b367946
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
7992
db5092052107 OpenMP support included
schloegl
parents: 7889
diff changeset
1
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
2 //-------------------------------------------------------------------
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
3 // C-MEX implementation of SUMSKIPNAN - this function is part of the NaN-toolbox.
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
4 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
5 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
6 // This program is free software; you can redistribute it and/or modify
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
7 // it under the terms of the GNU General Public License as published by
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
8 // the Free Software Foundation; either version 3 of the License, or
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
9 // (at your option) any later version.
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
10 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
11 // This program is distributed in the hope that it will be useful,
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
14 // GNU General Public License for more details.
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
15 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
16 // You should have received a copy of the GNU General Public License
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
17 // along with this program; if not, see <http://www.gnu.org/licenses/>.
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
18 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
19 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
20 // sumskipnan: sums all non-NaN values
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
21 // usage:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
22 // [o,count,SSQ] = sumskipnan_mex(x,DIM,flag,W);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
23 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
24 // SUMSKIPNAN uses two techniques to reduce errors:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
25 // 1) long double (80bit) instead of 64-bit double is used internally
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
26 // 2) The Kahan Summation formula is used to reduce the error margin from N*eps to 2*eps
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
27 // The latter is only implemented in case of stride=1 (column vectors only, summation along 1st dimension).
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
28 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
29 // Input:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
30 // - x data array
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
31 // - DIM (optional) dimension to sum
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
32 // - flag (optional) is actually an output argument telling whether some NaN was observed
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
33 // - W (optional) weight vector to compute weighted sum (default 1)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
34 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
35 // Output:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
36 // - o (weighted) sum along dimension DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
37 // - count of valid elements
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
38 // - sums of squares
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
39 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
40 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
41 // $Id$
8037
6a419bec96bb update contact e-mail and www address
schloegl
parents: 8035
diff changeset
42 // Copyright (C) 2009,2010,2011 Alois Schloegl <alois.schloegl@gmail.com>
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
43 // This function is part of the NaN-toolbox
7889
c101c486d80a fix web address
schloegl
parents: 7888
diff changeset
44 // http://pub.ist.ac.at/~schloegl/matlab/NaN/
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
45 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
46 //-------------------------------------------------------------------
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
47
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
48
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
49
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
50
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
51 #include <math.h>
7888
b9f35668b55e replace <inttypes.h> with <stdint.h>
schloegl
parents: 7301
diff changeset
52 #include <stdint.h>
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
53 #include "mex.h"
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
54
12691
6d6285a2a633 use macro ISNAN() instead of C++'s isnan() - because it supports all floating point formats not just double
schloegl
parents: 12685
diff changeset
55 /*
6d6285a2a633 use macro ISNAN() instead of C++'s isnan() - because it supports all floating point formats not just double
schloegl
parents: 12685
diff changeset
56 math.h has isnan() defined for all sizes of floating point numbers,
6d6285a2a633 use macro ISNAN() instead of C++'s isnan() - because it supports all floating point formats not just double
schloegl
parents: 12685
diff changeset
57 but c++ assumes isnan(double), causing possible conversions for float and long double
6d6285a2a633 use macro ISNAN() instead of C++'s isnan() - because it supports all floating point formats not just double
schloegl
parents: 12685
diff changeset
58 */
6d6285a2a633 use macro ISNAN() instead of C++'s isnan() - because it supports all floating point formats not just double
schloegl
parents: 12685
diff changeset
59 #define ISNAN(a) (a!=a)
6d6285a2a633 use macro ISNAN() instead of C++'s isnan() - because it supports all floating point formats not just double
schloegl
parents: 12685
diff changeset
60
6d6285a2a633 use macro ISNAN() instead of C++'s isnan() - because it supports all floating point formats not just double
schloegl
parents: 12685
diff changeset
61
12555
6dea7ba11ac7 address compiler warnings on functions with no return argument
schloegl
parents: 12306
diff changeset
62 inline void __sumskipnan2w__(double *data, size_t Ni, double *s, double *No, char *flag_anyISNAN, double *W);
6dea7ba11ac7 address compiler warnings on functions with no return argument
schloegl
parents: 12306
diff changeset
63 inline void __sumskipnan3w__(double *data, size_t Ni, double *s, double *s2, double *No, char *flag_anyISNAN, double *W);
6dea7ba11ac7 address compiler warnings on functions with no return argument
schloegl
parents: 12306
diff changeset
64 inline void __sumskipnan2wr__(double *data, size_t Ni, double *s, double *No, char *flag_anyISNAN, double *W);
6dea7ba11ac7 address compiler warnings on functions with no return argument
schloegl
parents: 12306
diff changeset
65 inline void __sumskipnan3wr__(double *data, size_t Ni, double *s, double *s2, double *No, char *flag_anyISNAN, double *W);
6dea7ba11ac7 address compiler warnings on functions with no return argument
schloegl
parents: 12306
diff changeset
66 inline void __sumskipnan2we__(double *data, size_t Ni, double *s, double *No, char *flag_anyISNAN, double *W);
6dea7ba11ac7 address compiler warnings on functions with no return argument
schloegl
parents: 12306
diff changeset
67 inline void __sumskipnan3we__(double *data, size_t Ni, double *s, double *s2, double *No, char *flag_anyISNAN, double *W);
6dea7ba11ac7 address compiler warnings on functions with no return argument
schloegl
parents: 12306
diff changeset
68 inline void __sumskipnan2wer__(double *data, size_t Ni, double *s, double *No, char *flag_anyISNAN, double *W);
6dea7ba11ac7 address compiler warnings on functions with no return argument
schloegl
parents: 12306
diff changeset
69 inline void __sumskipnan3wer__(double *data, size_t Ni, double *s, double *s2, double *No, char *flag_anyISNAN, double *W);
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
70
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
71 //#define NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
72
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
73 #ifdef tmwtypes_h
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
74 #if (MX_API_VER<=0x07020000)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
75 typedef int mwSize;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
76 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
77 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
78
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
79
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
80 void mexFunction(int POutputCount, mxArray* POutput[], int PInputCount, const mxArray *PInputs[])
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
81 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
82 const mwSize *SZ;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
83 double* LInput;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
84 double* LOutputSum;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
85 double* LOutputCount;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
86 double* LOutputSum2;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
87 long double* LongOutputSum = NULL;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
88 long double* LongOutputCount = NULL;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
89 long double* LongOutputSum2 = NULL;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
90 double x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
91 double* W = NULL; // weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
92
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
93 mwSize DIM = 0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
94 mwSize D1, D2, D3; // NN; //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
95 mwSize ND, ND2; // number of dimensions: input, output
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
96 mwSize ix0, ix1, ix2; // index to input and output
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
97 mwSize j, l; // running indices
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
98 mwSize *SZ2; // size of output
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
99 char flag_isNaN = 0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
100
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
101 // check for proper number of input and output arguments
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
102 if ((PInputCount <= 0) || (PInputCount > 4))
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
103 mexErrMsgTxt("SUMSKIPNAN.MEX requires between 1 and 4 arguments.");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
104 if (POutputCount > 4)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
105 mexErrMsgTxt("SUMSKIPNAN.MEX has 1 to 3 output arguments.");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
106
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
107 // get 1st argument
12640
de98e4cb9248 check for sparse matrices and and convert to full if needed
schloegl
parents: 12555
diff changeset
108 if(mxIsDouble(PInputs[0]) && !mxIsComplex(PInputs[0]) && !mxIsSparse(PInputs[0]) )
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
109 LInput = mxGetPr(PInputs[0]);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
110 else
12640
de98e4cb9248 check for sparse matrices and and convert to full if needed
schloegl
parents: 12555
diff changeset
111 mexErrMsgTxt("First argument must be and not sparse REAL/DOUBLE.");
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
112
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
113 // get 2nd argument
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
114 if (PInputCount > 1) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
115 switch (mxGetNumberOfElements(PInputs[1])) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
116 case 0: x = 0.0; // accept empty element
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
117 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
118 case 1: x = (mxIsNumeric(PInputs[1]) ? mxGetScalar(PInputs[1]) : -1.0);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
119 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
120 default:x = -1.0; // invalid
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
121 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
122 if ((x < 0) || (x > 65535) || (x != floor(x)))
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
123 mexErrMsgTxt("Error SUMSKIPNAN.MEX: DIM-argument must be a positive integer scalar");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
124
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
125 DIM = (unsigned)floor(x);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
126 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
127
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
128 // get size
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
129 ND = mxGetNumberOfDimensions(PInputs[0]);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
130 // NN = mxGetNumberOfElements(PInputs[0]);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
131 SZ = mxGetDimensions(PInputs[0]);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
132
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
133 // if DIM==0 (undefined), look for first dimension with more than 1 element.
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
134 for (j = 0; (DIM < 1) && (j < ND); j++)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
135 if (SZ[j]>1) DIM = j+1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
136
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
137 if (DIM < 1) DIM=1; // in case DIM is still undefined
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
138
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
139 ND2 = (ND>DIM ? ND : DIM); // number of dimensions of output
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
140
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
141 SZ2 = (mwSize*)mxCalloc(ND2, sizeof(mwSize)); // allocate memory for output size
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
142
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
143 for (j=0; j<ND; j++) // copy size of input;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
144 SZ2[j] = SZ[j];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
145 for (j=ND; j<ND2; j++) // in case DIM > ND, add extra elements 1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
146 SZ2[j] = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
147
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
148 for (j=0, D1=1; j<DIM-1; D1=D1*SZ2[j++]); // D1 is the number of elements between two elements along dimension DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
149 D2 = SZ2[DIM-1]; // D2 contains the size along dimension DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
150 for (j=DIM, D3=1; j<ND; D3=D3*SZ2[j++]); // D3 is the number of blocks containing D1*D2 elements
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
151
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
152 SZ2[DIM-1] = 1; // size of output is same as size of input but SZ(DIM)=1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
153
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
154 // get weight vector for weighted sumskipnan
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
155 if (PInputCount > 3) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
156 if (!mxGetNumberOfElements(PInputs[3]))
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
157 ; // empty weight vector - no weighting
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
158 else if (mxGetNumberOfElements(PInputs[3])==D2)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
159 W = mxGetPr(PInputs[3]);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
160 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
161 mexErrMsgTxt("Error SUMSKIPNAN.MEX: length of weight vector does not match size of dimension");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
162 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
163
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
164 int ACC_LEVEL = 0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
165 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
166 mxArray *LEVEL = NULL;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
167 int s = mexCallMATLAB(1, &LEVEL, 0, NULL, "flag_accuracy_level");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
168 if (!s) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
169 ACC_LEVEL = (int) mxGetScalar(LEVEL);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
170 if ((D1>1) && (ACC_LEVEL>2))
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
171 mexWarnMsgTxt("Warning: Kahan summation not supported with stride > 1 !");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
172 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
173 mxDestroyArray(LEVEL);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
174 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
175 // mexPrintf("Accuracy Level=%i\n",ACC_LEVEL);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
176
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
177 // create outputs
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
178 #define TYP mxDOUBLE_CLASS
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
179
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
180 POutput[0] = mxCreateNumericArray(ND2, SZ2, TYP, mxREAL);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
181 LOutputSum = mxGetPr(POutput[0]);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
182 if (D1!=1 && D2>0) LongOutputSum = (long double*) mxCalloc(D1*D3,sizeof(long double));
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
183 if (POutputCount >= 2) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
184 POutput[1] = mxCreateNumericArray(ND2, SZ2, TYP, mxREAL);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
185 LOutputCount = mxGetPr(POutput[1]);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
186 if (D1!=1 && D2>0) LongOutputCount = (long double*) mxCalloc(D1*D3,sizeof(long double));
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
187 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
188 if (POutputCount >= 3) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
189 POutput[2] = mxCreateNumericArray(ND2, SZ2, TYP, mxREAL);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
190 LOutputSum2 = mxGetPr(POutput[2]);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
191 if (D1!=1 && D2>0) LongOutputSum2 = (long double*) mxCalloc(D1*D3,sizeof(long double));
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
192 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
193 mxFree(SZ2);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
194
7992
db5092052107 OpenMP support included
schloegl
parents: 7889
diff changeset
195
8232
49c1c23128ea use only single test for (D1==1)
schloegl
parents: 8037
diff changeset
196 if (!D1 || !D2 || !D3) // zero size array
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
197 ; // do nothing
8232
49c1c23128ea use only single test for (D1==1)
schloegl
parents: 8037
diff changeset
198 else if (D1==1) {
49c1c23128ea use only single test for (D1==1)
schloegl
parents: 8037
diff changeset
199 if (ACC_LEVEL<1) {
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
200 // double accuray, naive summation, error = N*2^-52
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
201 switch (POutputCount) {
8875
903fbabbd5fe return output for case of no output variable
schloegl
parents: 8232
diff changeset
202 case 0:
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
203 case 1:
8035
f5040c012714 move #PRAGMA OMP PARALLEL into #PRAGMA OMP PARALLEL FOR - this fixes SegFaulting on ML7.11
schloegl
parents: 7992
diff changeset
204 #pragma omp parallel for schedule(dynamic)
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
205 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
206 double count;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
207 __sumskipnan2wr__(LInput+l*D2, D2, LOutputSum+l, &count, &flag_isNaN, W);
7992
db5092052107 OpenMP support included
schloegl
parents: 7889
diff changeset
208 }
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
209 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
210 case 2:
8035
f5040c012714 move #PRAGMA OMP PARALLEL into #PRAGMA OMP PARALLEL FOR - this fixes SegFaulting on ML7.11
schloegl
parents: 7992
diff changeset
211 #pragma omp parallel for schedule(dynamic)
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
212 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
213 __sumskipnan2wr__(LInput+l*D2, D2, LOutputSum+l, LOutputCount+l, &flag_isNaN, W);
7992
db5092052107 OpenMP support included
schloegl
parents: 7889
diff changeset
214 }
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
215 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
216 case 3:
8035
f5040c012714 move #PRAGMA OMP PARALLEL into #PRAGMA OMP PARALLEL FOR - this fixes SegFaulting on ML7.11
schloegl
parents: 7992
diff changeset
217 #pragma omp parallel for schedule(dynamic)
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
218 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
219 __sumskipnan3wr__(LInput+l*D2, D2, LOutputSum+l, LOutputSum2+l, LOutputCount+l, &flag_isNaN, W);
7992
db5092052107 OpenMP support included
schloegl
parents: 7889
diff changeset
220 }
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
221 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
222 }
8232
49c1c23128ea use only single test for (D1==1)
schloegl
parents: 8037
diff changeset
223 }
49c1c23128ea use only single test for (D1==1)
schloegl
parents: 8037
diff changeset
224 else if (ACC_LEVEL==1) {
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
225 // extended accuray, naive summation, error = N*2^-64
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
226 switch (POutputCount) {
8875
903fbabbd5fe return output for case of no output variable
schloegl
parents: 8232
diff changeset
227 case 0:
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
228 case 1:
8035
f5040c012714 move #PRAGMA OMP PARALLEL into #PRAGMA OMP PARALLEL FOR - this fixes SegFaulting on ML7.11
schloegl
parents: 7992
diff changeset
229 #pragma omp parallel for schedule(dynamic)
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
230 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
231 double count;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
232 __sumskipnan2w__(LInput+l*D2, D2, LOutputSum+l, &count, &flag_isNaN, W);
7992
db5092052107 OpenMP support included
schloegl
parents: 7889
diff changeset
233 }
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
234 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
235 case 2:
8035
f5040c012714 move #PRAGMA OMP PARALLEL into #PRAGMA OMP PARALLEL FOR - this fixes SegFaulting on ML7.11
schloegl
parents: 7992
diff changeset
236 #pragma omp parallel for schedule(dynamic)
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
237 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
238 __sumskipnan2w__(LInput+l*D2, D2, LOutputSum+l, LOutputCount+l, &flag_isNaN, W);
7992
db5092052107 OpenMP support included
schloegl
parents: 7889
diff changeset
239 }
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
240 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
241 case 3:
8035
f5040c012714 move #PRAGMA OMP PARALLEL into #PRAGMA OMP PARALLEL FOR - this fixes SegFaulting on ML7.11
schloegl
parents: 7992
diff changeset
242 #pragma omp parallel for schedule(dynamic)
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
243 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
244 __sumskipnan3w__(LInput+l*D2, D2, LOutputSum+l, LOutputSum2+l, LOutputCount+l, &flag_isNaN, W);
7992
db5092052107 OpenMP support included
schloegl
parents: 7889
diff changeset
245 }
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
246 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
247 }
8232
49c1c23128ea use only single test for (D1==1)
schloegl
parents: 8037
diff changeset
248 }
49c1c23128ea use only single test for (D1==1)
schloegl
parents: 8037
diff changeset
249 else if (ACC_LEVEL==3) {
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
250 // ACC_LEVEL==3: extended accuracy and Kahan Summation, error = 2^-64
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
251 switch (POutputCount) {
8875
903fbabbd5fe return output for case of no output variable
schloegl
parents: 8232
diff changeset
252 case 0:
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
253 case 1:
8035
f5040c012714 move #PRAGMA OMP PARALLEL into #PRAGMA OMP PARALLEL FOR - this fixes SegFaulting on ML7.11
schloegl
parents: 7992
diff changeset
254 #pragma omp parallel for schedule(dynamic)
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
255 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
256 double count;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
257 __sumskipnan2we__(LInput+l*D2, D2, LOutputSum+l, &count, &flag_isNaN, W);
7992
db5092052107 OpenMP support included
schloegl
parents: 7889
diff changeset
258 }
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
259 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
260 case 2:
8035
f5040c012714 move #PRAGMA OMP PARALLEL into #PRAGMA OMP PARALLEL FOR - this fixes SegFaulting on ML7.11
schloegl
parents: 7992
diff changeset
261 #pragma omp parallel for schedule(dynamic)
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
262 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
263 __sumskipnan2we__(LInput+l*D2, D2, LOutputSum+l, LOutputCount+l, &flag_isNaN, W);
7992
db5092052107 OpenMP support included
schloegl
parents: 7889
diff changeset
264 }
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
265 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
266 case 3:
8035
f5040c012714 move #PRAGMA OMP PARALLEL into #PRAGMA OMP PARALLEL FOR - this fixes SegFaulting on ML7.11
schloegl
parents: 7992
diff changeset
267 #pragma omp parallel for schedule(dynamic)
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
268 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
269 __sumskipnan3we__(LInput+l*D2, D2, LOutputSum+l, LOutputSum2+l, LOutputCount+l, &flag_isNaN, W);
7992
db5092052107 OpenMP support included
schloegl
parents: 7889
diff changeset
270 }
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
271 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
272 }
8232
49c1c23128ea use only single test for (D1==1)
schloegl
parents: 8037
diff changeset
273 }
49c1c23128ea use only single test for (D1==1)
schloegl
parents: 8037
diff changeset
274 else if (ACC_LEVEL==2) {
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
275 // ACC_LEVEL==2: double accuracy and Kahan Summation, error = 2^-52
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
276 switch (POutputCount) {
8875
903fbabbd5fe return output for case of no output variable
schloegl
parents: 8232
diff changeset
277 case 0:
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
278 case 1:
8035
f5040c012714 move #PRAGMA OMP PARALLEL into #PRAGMA OMP PARALLEL FOR - this fixes SegFaulting on ML7.11
schloegl
parents: 7992
diff changeset
279 #pragma omp parallel for schedule(dynamic)
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
280 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
281 double count;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
282 __sumskipnan2wer__(LInput+l*D2, D2, LOutputSum+l, &count, &flag_isNaN, W);
7992
db5092052107 OpenMP support included
schloegl
parents: 7889
diff changeset
283 }
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
284 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
285 case 2:
8035
f5040c012714 move #PRAGMA OMP PARALLEL into #PRAGMA OMP PARALLEL FOR - this fixes SegFaulting on ML7.11
schloegl
parents: 7992
diff changeset
286 #pragma omp parallel for schedule(dynamic)
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
287 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
288 __sumskipnan2wer__(LInput+l*D2, D2, LOutputSum+l, LOutputCount+l, &flag_isNaN, W);
7992
db5092052107 OpenMP support included
schloegl
parents: 7889
diff changeset
289 }
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
290 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
291 case 3:
8035
f5040c012714 move #PRAGMA OMP PARALLEL into #PRAGMA OMP PARALLEL FOR - this fixes SegFaulting on ML7.11
schloegl
parents: 7992
diff changeset
292 #pragma omp parallel for schedule(dynamic)
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
293 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
294 __sumskipnan3wer__(LInput+l*D2, D2, LOutputSum+l, LOutputSum2+l, LOutputCount+l, &flag_isNaN, W);
7992
db5092052107 OpenMP support included
schloegl
parents: 7889
diff changeset
295 }
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
296 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
297 }
8232
49c1c23128ea use only single test for (D1==1)
schloegl
parents: 8037
diff changeset
298 }
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
299 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
300 else if (POutputCount <= 1) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
301 // OUTER LOOP: along dimensions > DIM
7992
db5092052107 OpenMP support included
schloegl
parents: 7889
diff changeset
302 for (l = 0; l<D3; l++) {
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
303 ix0 = l*D1; // index for output
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
304 ix1 = ix0*D2; // index for input
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
305 for (j=0; j<D2; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
306 // minimize cache misses
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
307 ix2 = ix0; // index for output
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
308 // Inner LOOP: along dimensions < DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
309 if (W) do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
310 long double x = *LInput;
12691
6d6285a2a633 use macro ISNAN() instead of C++'s isnan() - because it supports all floating point formats not just double
schloegl
parents: 12685
diff changeset
311 if (!ISNAN(x)) {
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
312 LongOutputSum[ix2] += W[j]*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
313 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
314 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
315 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
316 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
317 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
318 LInput++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
319 ix2++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
320 } while (ix2 != (l+1)*D1);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
321 else do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
322 long double x = *LInput;
12691
6d6285a2a633 use macro ISNAN() instead of C++'s isnan() - because it supports all floating point formats not just double
schloegl
parents: 12685
diff changeset
323 if (!ISNAN(x)) {
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
324 LongOutputSum[ix2] += x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
325 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
326 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
327 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
328 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
329 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
330 LInput++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
331 ix2++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
332 } while (ix2 != (l+1)*D1);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
333 } // end for (j=
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
334
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
335 /* copy to output */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
336 for (j=0; j<D1; j++) {
12685
f26b1170ea90 resulting values should be really converted to output data type
schloegl
parents: 12640
diff changeset
337 LOutputSum[ix0+j] = (typeof(*LOutputSum))LongOutputSum[ix0+j];
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
338 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
339 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
340 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
341
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
342 else if (POutputCount == 2) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
343 // OUTER LOOP: along dimensions > DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
344 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
345 ix0 = l*D1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
346 ix1 = ix0*D2; // index for input
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
347 for (j=0; j<D2; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
348 // minimize cache misses
12306
6eb642de94a3 update docu of TSA and NaN package: push version to 4.2.7 and 2.5.8 for TSA and NaN package, resp.
schloegl
parents: 8875
diff changeset
349 ix2 = ix0; // index for output
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
350 // Inner LOOP: along dimensions < DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
351 if (W) do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
352 long double x = *LInput;
12691
6d6285a2a633 use macro ISNAN() instead of C++'s isnan() - because it supports all floating point formats not just double
schloegl
parents: 12685
diff changeset
353 if (!ISNAN(x)) {
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
354 LongOutputCount[ix2] += W[j];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
355 LongOutputSum[ix2] += W[j]*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
356 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
357 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
358 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
359 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
360 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
361 LInput++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
362 ix2++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
363 } while (ix2 != (l+1)*D1);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
364 else do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
365 long double x = *LInput;
12691
6d6285a2a633 use macro ISNAN() instead of C++'s isnan() - because it supports all floating point formats not just double
schloegl
parents: 12685
diff changeset
366 if (!ISNAN(x)) {
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
367 LongOutputCount[ix2] += 1.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
368 LongOutputSum[ix2] += x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
369 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
370 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
371 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
372 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
373 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
374 LInput++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
375 ix2++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
376 } while (ix2 != (l+1)*D1);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
377 } // end for (j=
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
378
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
379 /* copy to output */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
380 for (j=0; j<D1; j++) {
12685
f26b1170ea90 resulting values should be really converted to output data type
schloegl
parents: 12640
diff changeset
381 LOutputSum[ix0+j] = (typeof(*LOutputSum))LongOutputSum[ix0+j];
f26b1170ea90 resulting values should be really converted to output data type
schloegl
parents: 12640
diff changeset
382 LOutputCount[ix0+j] = (typeof(*LOutputCount))LongOutputCount[ix0+j];
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
383 } // end else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
384 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
385 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
386
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
387 else if (POutputCount == 3) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
388 // OUTER LOOP: along dimensions > DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
389 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
390 ix0 = l*D1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
391 ix1 = ix0*D2; // index for input
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
392 for (j=0; j<D2; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
393 // minimize cache misses
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
394 ix2 = ix0; // index for output
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
395 // Inner LOOP: along dimensions < DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
396 if (W) do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
397 long double x = *LInput;
12691
6d6285a2a633 use macro ISNAN() instead of C++'s isnan() - because it supports all floating point formats not just double
schloegl
parents: 12685
diff changeset
398 if (!ISNAN(x)) {
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
399 LongOutputCount[ix2] += W[j];
12685
f26b1170ea90 resulting values should be really converted to output data type
schloegl
parents: 12640
diff changeset
400 long double t = W[j]*x;
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
401 LongOutputSum[ix2] += t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
402 LongOutputSum2[ix2] += x*t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
403 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
404 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
405 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
406 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
407 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
408 LInput++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
409 ix2++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
410 } while (ix2 != (l+1)*D1);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
411 else do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
412 long double x = *LInput;
12691
6d6285a2a633 use macro ISNAN() instead of C++'s isnan() - because it supports all floating point formats not just double
schloegl
parents: 12685
diff changeset
413 if (!ISNAN(x)) {
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
414 LongOutputCount[ix2] += 1.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
415 LongOutputSum[ix2] += x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
416 LongOutputSum2[ix2] += x*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
417 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
418 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
419 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
420 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
421 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
422 LInput++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
423 ix2++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
424 } while (ix2 != (l+1)*D1);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
425 } // end for (j=
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
426
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
427 /* copy to output */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
428 for (j=0; j<D1; j++) {
12685
f26b1170ea90 resulting values should be really converted to output data type
schloegl
parents: 12640
diff changeset
429 LOutputSum[ix0+j] = (typeof(*LOutputSum))LongOutputSum[ix0+j];
f26b1170ea90 resulting values should be really converted to output data type
schloegl
parents: 12640
diff changeset
430 LOutputCount[ix0+j] = (typeof(*LOutputCount))LongOutputCount[ix0+j];
f26b1170ea90 resulting values should be really converted to output data type
schloegl
parents: 12640
diff changeset
431 LOutputSum2[ix0+j] = (typeof(*LOutputSum2))LongOutputSum2[ix0+j];
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
432 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
433 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
434 }
7992
db5092052107 OpenMP support included
schloegl
parents: 7889
diff changeset
435
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
436 if (LongOutputSum) mxFree(LongOutputSum);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
437 if (LongOutputCount) mxFree(LongOutputCount);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
438 if (LongOutputSum2) mxFree(LongOutputSum2);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
439
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
440 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
441 //mexPrintf("Third argument must be not empty - otherwise status whether a NaN occured or not cannot be returned.");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
442 /* this is a hack, the third input argument is used to return whether a NaN occured or not.
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
443 this requires that the input argument is a non-empty variable
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
444 */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
445 if (flag_isNaN && (PInputCount > 2) && mxGetNumberOfElements(PInputs[2])) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
446 // set FLAG_NANS_OCCURED
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
447 switch (mxGetClassID(PInputs[2])) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
448 case mxLOGICAL_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
449 case mxCHAR_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
450 case mxINT8_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
451 case mxUINT8_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
452 *(uint8_t*)mxGetData(PInputs[2]) = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
453 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
454 case mxDOUBLE_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
455 *(double*)mxGetData(PInputs[2]) = 1.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
456 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
457 case mxSINGLE_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
458 *(float*)mxGetData(PInputs[2]) = 1.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
459 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
460 case mxINT16_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
461 case mxUINT16_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
462 *(uint16_t*)mxGetData(PInputs[2]) = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
463 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
464 case mxINT32_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
465 case mxUINT32_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
466 *(uint32_t*)mxGetData(PInputs[2])= 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
467 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
468 case mxINT64_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
469 case mxUINT64_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
470 *(uint64_t*)mxGetData(PInputs[2]) = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
471 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
472 case mxFUNCTION_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
473 case mxUNKNOWN_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
474 case mxCELL_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
475 case mxSTRUCT_CLASS:
7301
485d1594d155 addresses undesired side-effect off in-place sorting of data
schloegl
parents: 6549
diff changeset
476 default:
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
477 mexPrintf("Type of 3rd input argument not supported.");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
478 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
479 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
480 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
481 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
482
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
483 #define stride 1
12555
6dea7ba11ac7 address compiler warnings on functions with no return argument
schloegl
parents: 12306
diff changeset
484 inline void __sumskipnan2w__(double *data, size_t Ni, double *s, double *No, char *flag_anyISNAN, double *W)
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
485 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
486 long double sum=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
487 char flag=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
488 // LOOP along dimension DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
489
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
490 double *end = data + stride*Ni;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
491 if (W) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
492 // with weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
493 long double count = 0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
494 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
495 long double x = *data;
12691
6d6285a2a633 use macro ISNAN() instead of C++'s isnan() - because it supports all floating point formats not just double
schloegl
parents: 12685
diff changeset
496 if (!ISNAN(x))
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
497 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
498 count += *W;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
499 sum += *W*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
500 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
501 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
502 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
503 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
504 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
505
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
506 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
507 W++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
508 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
509 while (data < end);
12685
f26b1170ea90 resulting values should be really converted to output data type
schloegl
parents: 12640
diff changeset
510 *No = (typeof(*No))count;
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
511 } else {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
512 // w/o weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
513 size_t countI = 0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
514 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
515 long double x = *data;
12691
6d6285a2a633 use macro ISNAN() instead of C++'s isnan() - because it supports all floating point formats not just double
schloegl
parents: 12685
diff changeset
516 if (!ISNAN(x))
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
517 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
518 countI++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
519 sum += x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
520 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
521 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
522 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
523 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
524 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
525 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
526 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
527 while (data < end);
12685
f26b1170ea90 resulting values should be really converted to output data type
schloegl
parents: 12640
diff changeset
528 *No = (typeof(*No))countI;
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
529 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
530
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
531 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
532 if (flag && (flag_anyISNAN != NULL)) *flag_anyISNAN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
533 #endif
12685
f26b1170ea90 resulting values should be really converted to output data type
schloegl
parents: 12640
diff changeset
534 *s = (typeof(*s))sum;
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
535
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
536 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
537
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
538
12555
6dea7ba11ac7 address compiler warnings on functions with no return argument
schloegl
parents: 12306
diff changeset
539 inline void __sumskipnan3w__(double *data, size_t Ni, double *s, double *s2, double *No, char *flag_anyISNAN, double *W)
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
540 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
541 long double sum=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
542 long double msq=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
543 char flag=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
544 // LOOP along dimension DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
545
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
546 double *end = data + stride*Ni;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
547 if (W) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
548 // with weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
549 long double count = 0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
550 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
551 long double x = *data;
12691
6d6285a2a633 use macro ISNAN() instead of C++'s isnan() - because it supports all floating point formats not just double
schloegl
parents: 12685
diff changeset
552 if (!ISNAN(x)) {
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
553 count += *W;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
554 long double t = *W*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
555 sum += t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
556 msq += x*t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
557 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
558 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
559 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
560 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
561 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
562 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
563 W++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
564 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
565 while (data < end);
12685
f26b1170ea90 resulting values should be really converted to output data type
schloegl
parents: 12640
diff changeset
566 *No = (typeof(*No))count;
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
567 } else {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
568 // w/o weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
569 size_t countI = 0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
570 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
571 long double x = *data;
12691
6d6285a2a633 use macro ISNAN() instead of C++'s isnan() - because it supports all floating point formats not just double
schloegl
parents: 12685
diff changeset
572 if (!ISNAN(x)) {
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
573 countI++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
574 sum += x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
575 msq += x*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
576 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
577 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
578 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
579 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
580 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
581 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
582 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
583 while (data < end);
12685
f26b1170ea90 resulting values should be really converted to output data type
schloegl
parents: 12640
diff changeset
584 *No = (typeof(*No))countI;
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
585 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
586
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
587 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
588 if (flag && (flag_anyISNAN != NULL)) *flag_anyISNAN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
589 #endif
12685
f26b1170ea90 resulting values should be really converted to output data type
schloegl
parents: 12640
diff changeset
590 *s = (typeof(*s))sum;
f26b1170ea90 resulting values should be really converted to output data type
schloegl
parents: 12640
diff changeset
591 *s2 = (typeof(*s2))msq;
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
592 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
593
12555
6dea7ba11ac7 address compiler warnings on functions with no return argument
schloegl
parents: 12306
diff changeset
594 inline void __sumskipnan2wr__(double *data, size_t Ni, double *s, double *No, char *flag_anyISNAN, double *W)
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
595 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
596 double sum=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
597 char flag=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
598 // LOOP along dimension DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
599
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
600 double *end = data + stride*Ni;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
601 if (W) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
602 // with weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
603 double count = 0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
604 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
605 double x = *data;
12691
6d6285a2a633 use macro ISNAN() instead of C++'s isnan() - because it supports all floating point formats not just double
schloegl
parents: 12685
diff changeset
606 if (!ISNAN(x))
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
607 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
608 count += *W;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
609 sum += *W*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
610 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
611 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
612 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
613 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
614 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
615
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
616 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
617 W++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
618 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
619 while (data < end);
12685
f26b1170ea90 resulting values should be really converted to output data type
schloegl
parents: 12640
diff changeset
620 *No = (typeof(*No))count;
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
621 } else {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
622 // w/o weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
623 size_t countI = 0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
624 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
625 double x = *data;
12691
6d6285a2a633 use macro ISNAN() instead of C++'s isnan() - because it supports all floating point formats not just double
schloegl
parents: 12685
diff changeset
626 if (!ISNAN(x))
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
627 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
628 countI++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
629 sum += x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
630 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
631 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
632 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
633 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
634 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
635 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
636 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
637 while (data < end);
12685
f26b1170ea90 resulting values should be really converted to output data type
schloegl
parents: 12640
diff changeset
638 *No = (typeof(*No))countI;
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
639 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
640
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
641 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
642 if (flag && (flag_anyISNAN != NULL)) *flag_anyISNAN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
643 #endif
12685
f26b1170ea90 resulting values should be really converted to output data type
schloegl
parents: 12640
diff changeset
644 *s = (typeof(*s))sum;
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
645
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
646 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
647
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
648
12555
6dea7ba11ac7 address compiler warnings on functions with no return argument
schloegl
parents: 12306
diff changeset
649 inline void __sumskipnan3wr__(double *data, size_t Ni, double *s, double *s2, double *No, char *flag_anyISNAN, double *W)
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
650 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
651 double sum=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
652 double msq=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
653 char flag=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
654 // LOOP along dimension DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
655
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
656 double *end = data + stride*Ni;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
657 if (W) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
658 // with weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
659 double count = 0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
660 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
661 double x = *data;
12691
6d6285a2a633 use macro ISNAN() instead of C++'s isnan() - because it supports all floating point formats not just double
schloegl
parents: 12685
diff changeset
662 if (!ISNAN(x)) {
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
663 count += *W;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
664 double t = *W*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
665 sum += t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
666 msq += x*t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
667 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
668 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
669 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
670 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
671 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
672 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
673 W++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
674 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
675 while (data < end);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
676 *No = count;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
677 } else {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
678 // w/o weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
679 size_t countI = 0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
680 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
681 double x = *data;
12691
6d6285a2a633 use macro ISNAN() instead of C++'s isnan() - because it supports all floating point formats not just double
schloegl
parents: 12685
diff changeset
682 if (!ISNAN(x)) {
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
683 countI++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
684 sum += x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
685 msq += x*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
686 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
687 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
688 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
689 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
690 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
691 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
692 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
693 while (data < end);
12685
f26b1170ea90 resulting values should be really converted to output data type
schloegl
parents: 12640
diff changeset
694 *No = (typeof(*No))countI;
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
695 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
696
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
697 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
698 if (flag && (flag_anyISNAN != NULL)) *flag_anyISNAN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
699 #endif
12685
f26b1170ea90 resulting values should be really converted to output data type
schloegl
parents: 12640
diff changeset
700 *s = (typeof(*s))sum;
f26b1170ea90 resulting values should be really converted to output data type
schloegl
parents: 12640
diff changeset
701 *s2 = (typeof(*s2))msq;
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
702 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
703
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
704
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
705
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
706 /***************************************
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
707 using Kahan's summation formula [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
708 this gives more accurate results while the computational effort within the loop is about 4x as high
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
709 First tests show a penalty of about 40% in terms of computational time.
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
710
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
711 [1] David Goldberg,
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
712 What Every Computer Scientist Should Know About Floating-Point Arithmetic
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
713 ACM Computing Surveys, Vol 23, No 1, March 1991.
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
714 ****************************************/
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
715
12555
6dea7ba11ac7 address compiler warnings on functions with no return argument
schloegl
parents: 12306
diff changeset
716 inline void __sumskipnan2we__(double *data, size_t Ni, double *s, double *No, char *flag_anyISNAN, double *W)
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
717 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
718 long double sum=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
719 char flag=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
720 // LOOP along dimension DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
721
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
722 double *end = data + stride*Ni;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
723 if (W) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
724 // with weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
725 long double count = 0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
726 long double rc=0.0, rn=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
727 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
728 long double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
729 long double t,y;
12691
6d6285a2a633 use macro ISNAN() instead of C++'s isnan() - because it supports all floating point formats not just double
schloegl
parents: 12685
diff changeset
730 if (!ISNAN(x))
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
731 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
732 //count += *W; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
733 y = *W-rn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
734 t = count+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
735 rn= (t-count)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
736 count= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
737
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
738 //sum += *W*x; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
739 y = *W*x-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
740 t = sum+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
741 rc= (t-sum)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
742 sum= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
743 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
744 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
745 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
746 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
747 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
748
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
749 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
750 W++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
751 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
752 while (data < end);
12685
f26b1170ea90 resulting values should be really converted to output data type
schloegl
parents: 12640
diff changeset
753 *No = (typeof(*No))count;
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
754 } else {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
755 // w/o weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
756 size_t countI = 0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
757 long double rc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
758 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
759 long double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
760 long double t,y;
12691
6d6285a2a633 use macro ISNAN() instead of C++'s isnan() - because it supports all floating point formats not just double
schloegl
parents: 12685
diff changeset
761 if (!ISNAN(x))
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
762 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
763 countI++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
764 // sum += x; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
765 y = x-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
766 t = sum+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
767 rc= (t-sum)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
768 sum= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
769 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
770 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
771 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
772 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
773 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
774 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
775 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
776 while (data < end);
12685
f26b1170ea90 resulting values should be really converted to output data type
schloegl
parents: 12640
diff changeset
777 *No = (typeof(*No))countI;
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
778 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
779
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
780 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
781 if (flag && (flag_anyISNAN != NULL)) *flag_anyISNAN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
782 #endif
12685
f26b1170ea90 resulting values should be really converted to output data type
schloegl
parents: 12640
diff changeset
783 *s = (typeof(*s))sum;
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
784
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
785 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
786
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
787
12555
6dea7ba11ac7 address compiler warnings on functions with no return argument
schloegl
parents: 12306
diff changeset
788 inline void __sumskipnan3we__(double *data, size_t Ni, double *s, double *s2, double *No, char *flag_anyISNAN, double *W)
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
789 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
790 long double sum=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
791 long double msq=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
792 char flag=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
793 // LOOP along dimension DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
794
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
795 double *end = data + stride*Ni;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
796 if (W) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
797 // with weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
798 long double count = 0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
799 long double rc=0.0, rn=0.0, rq=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
800 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
801 long double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
802 long double t,y;
12691
6d6285a2a633 use macro ISNAN() instead of C++'s isnan() - because it supports all floating point formats not just double
schloegl
parents: 12685
diff changeset
803 if (!ISNAN(x)) {
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
804 //count += *W; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
805 y = *W-rn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
806 t = count+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
807 rn= (t-count)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
808 count= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
809
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
810 long double w = *W*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
811 //sum += *W*x; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
812 y = *W*x-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
813 t = sum+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
814 rc= (t-sum)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
815 sum= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
816
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
817 // msq += x*w;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
818 y = w*x-rq;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
819 t = msq+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
820 rq= (t-msq)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
821 msq= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
822 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
823 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
824 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
825 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
826 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
827 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
828 W++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
829 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
830 while (data < end);
12685
f26b1170ea90 resulting values should be really converted to output data type
schloegl
parents: 12640
diff changeset
831 *No = (typeof(*No))count;
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
832 } else {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
833 // w/o weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
834 size_t countI = 0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
835 long double rc=0.0, rq=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
836 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
837 long double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
838 long double t,y;
12691
6d6285a2a633 use macro ISNAN() instead of C++'s isnan() - because it supports all floating point formats not just double
schloegl
parents: 12685
diff changeset
839 if (!ISNAN(x)) {
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
840 countI++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
841 //sum += x; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
842 y = x-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
843 t = sum+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
844 rc= (t-sum)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
845 sum= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
846
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
847 // msq += x*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
848 y = x*x-rq;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
849 t = msq+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
850 rq= (t-msq)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
851 msq= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
852 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
853 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
854 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
855 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
856 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
857 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
858 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
859 while (data < end);
12685
f26b1170ea90 resulting values should be really converted to output data type
schloegl
parents: 12640
diff changeset
860 *No = (typeof(*No))countI;
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
861 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
862
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
863 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
864 if (flag && (flag_anyISNAN != NULL)) *flag_anyISNAN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
865 #endif
12685
f26b1170ea90 resulting values should be really converted to output data type
schloegl
parents: 12640
diff changeset
866 *s = (typeof(*s))sum;
f26b1170ea90 resulting values should be really converted to output data type
schloegl
parents: 12640
diff changeset
867 *s2 = (typeof(*s))msq;
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
868 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
869
12555
6dea7ba11ac7 address compiler warnings on functions with no return argument
schloegl
parents: 12306
diff changeset
870 inline void __sumskipnan2wer__(double *data, size_t Ni, double *s, double *No, char *flag_anyISNAN, double *W)
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
871 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
872 double sum=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
873 char flag=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
874 // LOOP along dimension DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
875
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
876 double *end = data + stride*Ni;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
877 if (W) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
878 // with weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
879 double count = 0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
880 double rc=0.0, rn=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
881 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
882 double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
883 double t,y;
12691
6d6285a2a633 use macro ISNAN() instead of C++'s isnan() - because it supports all floating point formats not just double
schloegl
parents: 12685
diff changeset
884 if (!ISNAN(x))
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
885 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
886 //count += *W; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
887 y = *W-rn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
888 t = count+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
889 rn= (t-count)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
890 count= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
891
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
892 //sum += *W*x; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
893 y = *W*x-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
894 t = sum+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
895 rc= (t-sum)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
896 sum= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
897 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
898 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
899 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
900 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
901 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
902
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
903 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
904 W++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
905 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
906 while (data < end);
12685
f26b1170ea90 resulting values should be really converted to output data type
schloegl
parents: 12640
diff changeset
907 *No = (typeof(*No))count;
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
908 } else {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
909 // w/o weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
910 size_t countI = 0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
911 double rc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
912 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
913 double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
914 double t,y;
12691
6d6285a2a633 use macro ISNAN() instead of C++'s isnan() - because it supports all floating point formats not just double
schloegl
parents: 12685
diff changeset
915 if (!ISNAN(x))
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
916 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
917 countI++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
918 // sum += x; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
919 y = x-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
920 t = sum+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
921 rc= (t-sum)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
922 sum= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
923 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
924 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
925 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
926 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
927 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
928 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
929 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
930 while (data < end);
12685
f26b1170ea90 resulting values should be really converted to output data type
schloegl
parents: 12640
diff changeset
931 *No = (typeof(*No))countI;
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
932 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
933
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
934 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
935 if (flag && (flag_anyISNAN != NULL)) *flag_anyISNAN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
936 #endif
12685
f26b1170ea90 resulting values should be really converted to output data type
schloegl
parents: 12640
diff changeset
937 *s = (typeof(*s))sum;
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
938
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
939 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
940
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
941
12555
6dea7ba11ac7 address compiler warnings on functions with no return argument
schloegl
parents: 12306
diff changeset
942 inline void __sumskipnan3wer__(double *data, size_t Ni, double *s, double *s2, double *No, char *flag_anyISNAN, double *W)
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
943 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
944 double sum=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
945 double msq=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
946 char flag=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
947 // LOOP along dimension DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
948
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
949 double *end = data + stride*Ni;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
950 if (W) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
951 // with weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
952 double count = 0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
953 double rc=0.0, rn=0.0, rq=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
954 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
955 double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
956 double t,y;
12691
6d6285a2a633 use macro ISNAN() instead of C++'s isnan() - because it supports all floating point formats not just double
schloegl
parents: 12685
diff changeset
957 if (!ISNAN(x)) {
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
958 //count += *W; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
959 y = *W-rn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
960 t = count+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
961 rn= (t-count)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
962 count= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
963
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
964 double w = *W*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
965 //sum += *W*x; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
966 y = *W*x-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
967 t = sum+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
968 rc= (t-sum)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
969 sum= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
970
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
971 // msq += x*w;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
972 y = w*x-rq;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
973 t = msq+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
974 rq= (t-msq)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
975 msq= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
976 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
977 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
978 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
979 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
980 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
981 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
982 W++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
983 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
984 while (data < end);
12685
f26b1170ea90 resulting values should be really converted to output data type
schloegl
parents: 12640
diff changeset
985 *No = (typeof(*No))count;
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
986 } else {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
987 // w/o weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
988 size_t countI = 0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
989 double rc=0.0, rq=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
990 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
991 double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
992 double t,y;
12691
6d6285a2a633 use macro ISNAN() instead of C++'s isnan() - because it supports all floating point formats not just double
schloegl
parents: 12685
diff changeset
993 if (!ISNAN(x)) {
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
994 countI++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
995 //sum += x; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
996 y = x-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
997 t = sum+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
998 rc= (t-sum)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
999 sum= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1000
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1001 // msq += x*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1002 y = x*x-rq;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1003 t = msq+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1004 rq= (t-msq)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1005 msq= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1006 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1007 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1008 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1009 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1010 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1011 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1012 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1013 while (data < end);
12685
f26b1170ea90 resulting values should be really converted to output data type
schloegl
parents: 12640
diff changeset
1014 *No = (typeof(*No))countI;
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1015 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1016
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1017 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1018 if (flag && (flag_anyISNAN != NULL)) *flag_anyISNAN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1019 #endif
12685
f26b1170ea90 resulting values should be really converted to output data type
schloegl
parents: 12640
diff changeset
1020 *s = (typeof(*s))sum;
f26b1170ea90 resulting values should be really converted to output data type
schloegl
parents: 12640
diff changeset
1021 *s2 = (typeof(*s))msq;
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1022 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1023