Mercurial > forge

annotate extra/NaN/src/sumskipnan_mex.cpp @ 8037:6a419bec96bb octave-forge

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