Mercurial > forge

annotate extra/NaN/src/sumskipnan_mex.cpp @ 8232:49c1c23128ea octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
use only single test for (D1==1)
author schloegl
date Wed, 27 Jul 2011 22:20:09 +0000
parents 6a419bec96bb
children 903fbabbd5fe
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
8232
49c1c23128ea use only single test for (D1==1)
schloegl
parents: 8037
diff changeset
189 if (!D1 || !D2 || !D3) // zero size array
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
190 ; // do nothing
8232
49c1c23128ea use only single test for (D1==1)
schloegl
parents: 8037
diff changeset
191 else if (D1==1) {
49c1c23128ea use only single test for (D1==1)
schloegl
parents: 8037
diff changeset
192 if (ACC_LEVEL<1) {
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
193 // double accuray, naive summation, error = N*2^-52
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
194 switch (POutputCount) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
195 case 1:
8035
f5040c012714 move #PRAGMA OMP PARALLEL into #PRAGMA OMP PARALLEL FOR - this fixes SegFaulting on ML7.11
schloegl
parents: 7992
diff changeset
196 #pragma omp parallel for schedule(dynamic)
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
197 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
198 double count;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
199 __sumskipnan2wr__(LInput+l*D2, D2, LOutputSum+l, &count, &flag_isNaN, W);
7992
db5092052107 OpenMP support included
schloegl
parents: 7889
diff changeset
200 }
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
201 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
202 case 2:
8035
f5040c012714 move #PRAGMA OMP PARALLEL into #PRAGMA OMP PARALLEL FOR - this fixes SegFaulting on ML7.11
schloegl
parents: 7992
diff changeset
203 #pragma omp parallel for schedule(dynamic)
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
204 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
205 __sumskipnan2wr__(LInput+l*D2, D2, LOutputSum+l, LOutputCount+l, &flag_isNaN, W);
7992
db5092052107 OpenMP support included
schloegl
parents: 7889
diff changeset
206 }
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
207 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
208 case 3:
8035
f5040c012714 move #PRAGMA OMP PARALLEL into #PRAGMA OMP PARALLEL FOR - this fixes SegFaulting on ML7.11
schloegl
parents: 7992
diff changeset
209 #pragma omp parallel for schedule(dynamic)
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
210 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
211 __sumskipnan3wr__(LInput+l*D2, D2, LOutputSum+l, LOutputSum2+l, LOutputCount+l, &flag_isNaN, W);
7992
db5092052107 OpenMP support included
schloegl
parents: 7889
diff changeset
212 }
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
213 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
214 }
8232
49c1c23128ea use only single test for (D1==1)
schloegl
parents: 8037
diff changeset
215 }
49c1c23128ea use only single test for (D1==1)
schloegl
parents: 8037
diff changeset
216 else if (ACC_LEVEL==1) {
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
217 // extended accuray, naive summation, error = N*2^-64
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
218 switch (POutputCount) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
219 case 1:
8035
f5040c012714 move #PRAGMA OMP PARALLEL into #PRAGMA OMP PARALLEL FOR - this fixes SegFaulting on ML7.11
schloegl
parents: 7992
diff changeset
220 #pragma omp parallel for schedule(dynamic)
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
221 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
222 double count;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
223 __sumskipnan2w__(LInput+l*D2, D2, LOutputSum+l, &count, &flag_isNaN, W);
7992
db5092052107 OpenMP support included
schloegl
parents: 7889
diff changeset
224 }
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
225 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
226 case 2:
8035
f5040c012714 move #PRAGMA OMP PARALLEL into #PRAGMA OMP PARALLEL FOR - this fixes SegFaulting on ML7.11
schloegl
parents: 7992
diff changeset
227 #pragma omp parallel for schedule(dynamic)
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
228 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
229 __sumskipnan2w__(LInput+l*D2, D2, LOutputSum+l, LOutputCount+l, &flag_isNaN, W);
7992
db5092052107 OpenMP support included
schloegl
parents: 7889
diff changeset
230 }
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
231 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
232 case 3:
8035
f5040c012714 move #PRAGMA OMP PARALLEL into #PRAGMA OMP PARALLEL FOR - this fixes SegFaulting on ML7.11
schloegl
parents: 7992
diff changeset
233 #pragma omp parallel for schedule(dynamic)
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
234 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
235 __sumskipnan3w__(LInput+l*D2, D2, LOutputSum+l, LOutputSum2+l, LOutputCount+l, &flag_isNaN, W);
7992
db5092052107 OpenMP support included
schloegl
parents: 7889
diff changeset
236 }
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
237 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
238 }
8232
49c1c23128ea use only single test for (D1==1)
schloegl
parents: 8037
diff changeset
239 }
49c1c23128ea use only single test for (D1==1)
schloegl
parents: 8037
diff changeset
240 else if (ACC_LEVEL==3) {
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
241 // ACC_LEVEL==3: extended accuracy and Kahan Summation, error = 2^-64
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
242 switch (POutputCount) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
243 case 1:
8035
f5040c012714 move #PRAGMA OMP PARALLEL into #PRAGMA OMP PARALLEL FOR - this fixes SegFaulting on ML7.11
schloegl
parents: 7992
diff changeset
244 #pragma omp parallel for schedule(dynamic)
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
245 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
246 double count;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
247 __sumskipnan2we__(LInput+l*D2, D2, LOutputSum+l, &count, &flag_isNaN, W);
7992
db5092052107 OpenMP support included
schloegl
parents: 7889
diff changeset
248 }
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
249 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
250 case 2:
8035
f5040c012714 move #PRAGMA OMP PARALLEL into #PRAGMA OMP PARALLEL FOR - this fixes SegFaulting on ML7.11
schloegl
parents: 7992
diff changeset
251 #pragma omp parallel for schedule(dynamic)
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
252 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
253 __sumskipnan2we__(LInput+l*D2, D2, LOutputSum+l, LOutputCount+l, &flag_isNaN, W);
7992
db5092052107 OpenMP support included
schloegl
parents: 7889
diff changeset
254 }
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
255 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
256 case 3:
8035
f5040c012714 move #PRAGMA OMP PARALLEL into #PRAGMA OMP PARALLEL FOR - this fixes SegFaulting on ML7.11
schloegl
parents: 7992
diff changeset
257 #pragma omp parallel for schedule(dynamic)
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
258 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
259 __sumskipnan3we__(LInput+l*D2, D2, LOutputSum+l, LOutputSum2+l, LOutputCount+l, &flag_isNaN, W);
7992
db5092052107 OpenMP support included
schloegl
parents: 7889
diff changeset
260 }
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
261 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
262 }
8232
49c1c23128ea use only single test for (D1==1)
schloegl
parents: 8037
diff changeset
263 }
49c1c23128ea use only single test for (D1==1)
schloegl
parents: 8037
diff changeset
264 else if (ACC_LEVEL==2) {
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
265 // ACC_LEVEL==2: double accuracy and Kahan Summation, error = 2^-52
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
266 switch (POutputCount) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
267 case 1:
8035
f5040c012714 move #PRAGMA OMP PARALLEL into #PRAGMA OMP PARALLEL FOR - this fixes SegFaulting on ML7.11
schloegl
parents: 7992
diff changeset
268 #pragma omp parallel for schedule(dynamic)
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
269 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
270 double count;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
271 __sumskipnan2wer__(LInput+l*D2, D2, LOutputSum+l, &count, &flag_isNaN, W);
7992
db5092052107 OpenMP support included
schloegl
parents: 7889
diff changeset
272 }
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
273 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
274 case 2:
8035
f5040c012714 move #PRAGMA OMP PARALLEL into #PRAGMA OMP PARALLEL FOR - this fixes SegFaulting on ML7.11
schloegl
parents: 7992
diff changeset
275 #pragma omp parallel for schedule(dynamic)
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
276 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
277 __sumskipnan2wer__(LInput+l*D2, D2, LOutputSum+l, LOutputCount+l, &flag_isNaN, W);
7992
db5092052107 OpenMP support included
schloegl
parents: 7889
diff changeset
278 }
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
279 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
280 case 3:
8035
f5040c012714 move #PRAGMA OMP PARALLEL into #PRAGMA OMP PARALLEL FOR - this fixes SegFaulting on ML7.11
schloegl
parents: 7992
diff changeset
281 #pragma omp parallel for schedule(dynamic)
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
282 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
283 __sumskipnan3wer__(LInput+l*D2, D2, LOutputSum+l, LOutputSum2+l, LOutputCount+l, &flag_isNaN, W);
7992
db5092052107 OpenMP support included
schloegl
parents: 7889
diff changeset
284 }
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
285 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
286 }
8232
49c1c23128ea use only single test for (D1==1)
schloegl
parents: 8037
diff changeset
287 }
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
288 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
289 else if (POutputCount <= 1) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
290 // OUTER LOOP: along dimensions > DIM
7992
db5092052107 OpenMP support included
schloegl
parents: 7889
diff changeset
291 for (l = 0; l<D3; l++) {
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
292 ix0 = l*D1; // index for output
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
293 ix1 = ix0*D2; // index for input
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
294 for (j=0; j<D2; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
295 // minimize cache misses
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
296 ix2 = ix0; // index for output
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
297 // Inner LOOP: along dimensions < DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
298 if (W) do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
299 long double x = *LInput;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
300 if (!isnan(x)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
301 LongOutputSum[ix2] += W[j]*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
302 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
303 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
304 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
305 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
306 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
307 LInput++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
308 ix2++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
309 } while (ix2 != (l+1)*D1);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
310 else do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
311 long double x = *LInput;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
312 if (!isnan(x)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
313 LongOutputSum[ix2] += x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
314 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
315 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
316 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
317 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
318 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
319 LInput++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
320 ix2++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
321 } while (ix2 != (l+1)*D1);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
322 } // end for (j=
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
323
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
324 /* copy to output */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
325 for (j=0; j<D1; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
326 LOutputSum[ix0+j] = LongOutputSum[ix0+j];
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 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
330
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
331 else if (POutputCount == 2) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
332 // OUTER LOOP: along dimensions > DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
333 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
334 ix0 = l*D1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
335 ix1 = ix0*D2; // index for input
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
336 for (j=0; j<D2; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
337 // minimize cache misses
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
338 ix2 = ix0; // index for output
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
339 // Inner LOOP: along dimensions < DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
340 if (W) do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
341 long double x = *LInput;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
342 if (!isnan(x)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
343 LongOutputCount[ix2] += W[j];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
344 LongOutputSum[ix2] += W[j]*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
345 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
346 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
347 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
348 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
349 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
350 LInput++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
351 ix2++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
352 } while (ix2 != (l+1)*D1);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
353 else do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
354 long double x = *LInput;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
355 if (!isnan(x)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
356 LongOutputCount[ix2] += 1.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
357 LongOutputSum[ix2] += x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
358 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
359 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
360 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
361 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
362 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
363 LInput++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
364 ix2++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
365 } while (ix2 != (l+1)*D1);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
366 } // end for (j=
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
367
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
368 /* copy to output */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
369 for (j=0; j<D1; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
370 LOutputSum[ix0+j] = LongOutputSum[ix0+j];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
371 LOutputCount[ix0+j] = LongOutputCount[ix0+j];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
372 } // end else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
373 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
374 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
375
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
376 else if (POutputCount == 3) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
377 // OUTER LOOP: along dimensions > DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
378 for (l = 0; l<D3; l++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
379 ix0 = l*D1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
380 ix1 = ix0*D2; // index for input
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
381 for (j=0; j<D2; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
382 // minimize cache misses
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
383 ix2 = ix0; // index for output
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
384 // Inner LOOP: along dimensions < DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
385 if (W) do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
386 long double x = *LInput;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
387 if (!isnan(x)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
388 LongOutputCount[ix2] += W[j];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
389 double t = W[j]*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
390 LongOutputSum[ix2] += t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
391 LongOutputSum2[ix2] += x*t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
392 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
393 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
394 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
395 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
396 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
397 LInput++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
398 ix2++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
399 } while (ix2 != (l+1)*D1);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
400 else do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
401 long double x = *LInput;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
402 if (!isnan(x)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
403 LongOutputCount[ix2] += 1.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
404 LongOutputSum[ix2] += x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
405 LongOutputSum2[ix2] += x*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
406 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
407 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
408 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
409 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
410 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
411 LInput++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
412 ix2++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
413 } while (ix2 != (l+1)*D1);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
414 } // end for (j=
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
415
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
416 /* copy to output */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
417 for (j=0; j<D1; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
418 LOutputSum[ix0+j] = LongOutputSum[ix0+j];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
419 LOutputCount[ix0+j] = LongOutputCount[ix0+j];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
420 LOutputSum2[ix0+j] = LongOutputSum2[ix0+j];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
421 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
422 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
423 }
7992
db5092052107 OpenMP support included
schloegl
parents: 7889
diff changeset
424
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
425 if (LongOutputSum) mxFree(LongOutputSum);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
426 if (LongOutputCount) mxFree(LongOutputCount);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
427 if (LongOutputSum2) mxFree(LongOutputSum2);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
428
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
429 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
430 //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
431 /* 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
432 this requires that the input argument is a non-empty variable
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
433 */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
434 if (flag_isNaN && (PInputCount > 2) && mxGetNumberOfElements(PInputs[2])) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
435 // set FLAG_NANS_OCCURED
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
436 switch (mxGetClassID(PInputs[2])) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
437 case mxLOGICAL_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
438 case mxCHAR_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
439 case mxINT8_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
440 case mxUINT8_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
441 *(uint8_t*)mxGetData(PInputs[2]) = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
442 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
443 case mxDOUBLE_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
444 *(double*)mxGetData(PInputs[2]) = 1.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
445 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
446 case mxSINGLE_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
447 *(float*)mxGetData(PInputs[2]) = 1.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
448 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
449 case mxINT16_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
450 case mxUINT16_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
451 *(uint16_t*)mxGetData(PInputs[2]) = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
452 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
453 case mxINT32_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
454 case mxUINT32_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
455 *(uint32_t*)mxGetData(PInputs[2])= 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
456 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
457 case mxINT64_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
458 case mxUINT64_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
459 *(uint64_t*)mxGetData(PInputs[2]) = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
460 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
461 case mxFUNCTION_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
462 case mxUNKNOWN_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
463 case mxCELL_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
464 case mxSTRUCT_CLASS:
7301
485d1594d155 addresses undesired side-effect off in-place sorting of data
schloegl
parents: 6549
diff changeset
465 default:
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
466 mexPrintf("Type of 3rd input argument not supported.");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
467 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
468 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
469 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
470 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
471
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
472 #define stride 1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
473 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
474 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
475 long double sum=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
476 char flag=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
477 // LOOP along dimension DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
478
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
479 double *end = data + stride*Ni;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
480 if (W) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
481 // with weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
482 long double count = 0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
483 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
484 long double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
485 if (!isnan(x))
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
486 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
487 count += *W;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
488 sum += *W*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
489 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
490 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
491 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
492 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
493 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
494
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
495 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
496 W++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
497 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
498 while (data < end);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
499 *No = count;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
500 } else {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
501 // w/o weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
502 size_t countI = 0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
503 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
504 long double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
505 if (!isnan(x))
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
506 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
507 countI++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
508 sum += x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
509 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
510 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
511 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
512 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
513 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
514 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
515 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
516 while (data < end);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
517 *No = (double)countI;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
518 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
519
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
520 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
521 if (flag && (flag_anyISNAN != NULL)) *flag_anyISNAN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
522 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
523 *s = sum;
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
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
527
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
528 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
529 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
530 long double sum=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
531 long double msq=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
532 char flag=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
533 // LOOP along dimension DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
534
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
535 double *end = data + stride*Ni;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
536 if (W) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
537 // with weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
538 long double count = 0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
539 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
540 long double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
541 if (!isnan(x)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
542 count += *W;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
543 long double t = *W*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
544 sum += t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
545 msq += x*t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
546 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
547 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
548 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
549 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
550 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
551 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
552 W++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
553 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
554 while (data < end);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
555 *No = count;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
556 } else {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
557 // w/o weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
558 size_t countI = 0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
559 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
560 long double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
561 if (!isnan(x)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
562 countI++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
563 sum += x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
564 msq += x*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
565 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
566 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
567 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
568 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
569 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
570 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
571 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
572 while (data < end);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
573 *No = (double)countI;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
574 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
575
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
576 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
577 if (flag && (flag_anyISNAN != NULL)) *flag_anyISNAN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
578 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
579 *s = sum;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
580 *s2 = msq;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
581 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
582
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
583 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
584 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
585 double sum=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
586 char flag=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
587 // LOOP along dimension DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
588
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
589 double *end = data + stride*Ni;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
590 if (W) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
591 // with weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
592 double count = 0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
593 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
594 double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
595 if (!isnan(x))
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
596 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
597 count += *W;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
598 sum += *W*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
599 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
600 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
601 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
602 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
603 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
604
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
605 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
606 W++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
607 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
608 while (data < end);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
609 *No = count;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
610 } else {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
611 // w/o weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
612 size_t countI = 0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
613 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
614 double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
615 if (!isnan(x))
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
616 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
617 countI++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
618 sum += x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
619 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
620 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
621 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
622 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
623 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
624 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
625 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
626 while (data < end);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
627 *No = (double)countI;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
628 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
629
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
630 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
631 if (flag && (flag_anyISNAN != NULL)) *flag_anyISNAN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
632 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
633 *s = sum;
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
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
637
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
638 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
639 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
640 double sum=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
641 double msq=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
642 char flag=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
643 // LOOP along dimension DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
644
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
645 double *end = data + stride*Ni;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
646 if (W) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
647 // with weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
648 double count = 0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
649 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
650 double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
651 if (!isnan(x)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
652 count += *W;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
653 double t = *W*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
654 sum += t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
655 msq += x*t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
656 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
657 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
658 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
659 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
660 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
661 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
662 W++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
663 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
664 while (data < end);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
665 *No = count;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
666 } else {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
667 // w/o weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
668 size_t countI = 0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
669 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
670 double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
671 if (!isnan(x)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
672 countI++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
673 sum += x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
674 msq += x*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
675 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
676 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
677 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
678 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
679 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
680 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
681 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
682 while (data < end);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
683 *No = (double)countI;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
684 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
685
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
686 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
687 if (flag && (flag_anyISNAN != NULL)) *flag_anyISNAN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
688 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
689 *s = sum;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
690 *s2 = msq;
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
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
695 /***************************************
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
696 using Kahan's summation formula [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
697 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
698 First tests show a penalty of about 40% in terms of computational time.
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
699
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
700 [1] David Goldberg,
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
701 What Every Computer Scientist Should Know About Floating-Point Arithmetic
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
702 ACM Computing Surveys, Vol 23, No 1, March 1991.
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
703 ****************************************/
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
704
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
705 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
706 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
707 long double sum=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
708 char flag=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
709 // LOOP along dimension DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
710
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
711 double *end = data + stride*Ni;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
712 if (W) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
713 // with weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
714 long double count = 0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
715 long double rc=0.0, rn=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
716 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
717 long double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
718 long double t,y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
719 if (!isnan(x))
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
720 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
721 //count += *W; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
722 y = *W-rn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
723 t = count+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
724 rn= (t-count)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
725 count= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
726
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
727 //sum += *W*x; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
728 y = *W*x-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
729 t = sum+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
730 rc= (t-sum)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
731 sum= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
732 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
733 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
734 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
735 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
736 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
737
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
738 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
739 W++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
740 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
741 while (data < end);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
742 *No = count;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
743 } else {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
744 // w/o weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
745 size_t countI = 0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
746 long double rc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
747 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
748 long double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
749 long double t,y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
750 if (!isnan(x))
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
751 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
752 countI++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
753 // sum += x; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
754 y = x-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
755 t = sum+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
756 rc= (t-sum)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
757 sum= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
758 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
759 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
760 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
761 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
762 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
763 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
764 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
765 while (data < end);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
766 *No = (double)countI;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
767 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
768
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
769 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
770 if (flag && (flag_anyISNAN != NULL)) *flag_anyISNAN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
771 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
772 *s = sum;
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
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
776
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
777 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
778 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
779 long double sum=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
780 long double msq=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
781 char flag=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
782 // LOOP along dimension DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
783
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
784 double *end = data + stride*Ni;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
785 if (W) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
786 // with weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
787 long double count = 0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
788 long double rc=0.0, rn=0.0, rq=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
789 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
790 long double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
791 long double t,y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
792 if (!isnan(x)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
793 //count += *W; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
794 y = *W-rn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
795 t = count+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
796 rn= (t-count)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
797 count= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
798
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
799 long double w = *W*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
800 //sum += *W*x; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
801 y = *W*x-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
802 t = sum+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
803 rc= (t-sum)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
804 sum= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
805
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
806 // msq += x*w;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
807 y = w*x-rq;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
808 t = msq+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
809 rq= (t-msq)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
810 msq= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
811 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
812 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
813 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
814 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
815 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
816 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
817 W++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
818 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
819 while (data < end);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
820 *No = count;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
821 } else {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
822 // w/o weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
823 size_t countI = 0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
824 long double rc=0.0, rq=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
825 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
826 long double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
827 long double t,y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
828 if (!isnan(x)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
829 countI++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
830 //sum += x; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
831 y = x-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
832 t = sum+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
833 rc= (t-sum)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
834 sum= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
835
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
836 // msq += x*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
837 y = x*x-rq;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
838 t = msq+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
839 rq= (t-msq)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
840 msq= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
841 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
842 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
843 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
844 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
845 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
846 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
847 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
848 while (data < end);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
849 *No = (double)countI;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
850 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
851
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
852 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
853 if (flag && (flag_anyISNAN != NULL)) *flag_anyISNAN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
854 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
855 *s = sum;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
856 *s2 = msq;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
857 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
858
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
859 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
860 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
861 double sum=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
862 char flag=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
863 // LOOP along dimension DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
864
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
865 double *end = data + stride*Ni;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
866 if (W) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
867 // with weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
868 double count = 0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
869 double rc=0.0, rn=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
870 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
871 double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
872 double t,y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
873 if (!isnan(x))
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
874 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
875 //count += *W; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
876 y = *W-rn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
877 t = count+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
878 rn= (t-count)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
879 count= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
880
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
881 //sum += *W*x; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
882 y = *W*x-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
883 t = sum+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
884 rc= (t-sum)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
885 sum= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
886 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
887 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
888 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
889 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
890 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
891
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
892 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
893 W++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
894 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
895 while (data < end);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
896 *No = count;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
897 } else {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
898 // w/o weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
899 size_t countI = 0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
900 double rc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
901 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
902 double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
903 double t,y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
904 if (!isnan(x))
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
905 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
906 countI++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
907 // sum += x; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
908 y = x-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
909 t = sum+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
910 rc= (t-sum)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
911 sum= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
912 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
913 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
914 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
915 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
916 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
917 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
918 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
919 while (data < end);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
920 *No = (double)countI;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
921 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
922
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
923 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
924 if (flag && (flag_anyISNAN != NULL)) *flag_anyISNAN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
925 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
926 *s = sum;
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
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
930
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
931 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
932 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
933 double sum=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
934 double msq=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
935 char flag=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
936 // LOOP along dimension DIM
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
937
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
938 double *end = data + stride*Ni;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
939 if (W) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
940 // with weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
941 double count = 0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
942 double rc=0.0, rn=0.0, rq=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
943 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
944 double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
945 double t,y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
946 if (!isnan(x)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
947 //count += *W; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
948 y = *W-rn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
949 t = count+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
950 rn= (t-count)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
951 count= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
952
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
953 double w = *W*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
954 //sum += *W*x; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
955 y = *W*x-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
956 t = sum+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
957 rc= (t-sum)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
958 sum= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
959
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
960 // msq += x*w;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
961 y = w*x-rq;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
962 t = msq+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
963 rq= (t-msq)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
964 msq= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
965 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
966 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
967 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
968 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
969 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
970 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
971 W++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
972 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
973 while (data < end);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
974 *No = count;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
975 } else {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
976 // w/o weight vector
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
977 size_t countI = 0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
978 double rc=0.0, rq=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
979 do {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
980 double x = *data;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
981 double t,y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
982 if (!isnan(x)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
983 countI++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
984 //sum += x; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
985 y = x-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
986 t = sum+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
987 rc= (t-sum)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
988 sum= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
989
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
990 // msq += x*x;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
991 y = x*x-rq;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
992 t = msq+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
993 rq= (t-msq)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
994 msq= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
995 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
996 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
997 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
998 flag = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
999 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1000 data++; // stride=1
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1001 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1002 while (data < end);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1003 *No = (double)countI;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1004 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1005
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1006 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1007 if (flag && (flag_anyISNAN != NULL)) *flag_anyISNAN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1008 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1009 *s = sum;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1010 *s2 = msq;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1011 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1012