Mercurial > forge

annotate extra/NaN/src/covm_mex.cpp @ 6549:41e9854fe26d octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
use *.cpp instead of *.c
author schloegl
date Sun, 10 Jan 2010 22:05:59 +0000
parents
children ae521dec5b54
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
1 /*
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
2 //-------------------------------------------------------------------
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
3 #pragma hdrstop
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
4 //-------------------------------------------------------------------
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
5 // C-MEX implementation of COVM - this function is part of the NaN-toolbox.
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
6 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
7 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
8 // This program is free software; you can redistribute it and/or modify
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
9 // it under the terms of the GNU General Public License as published by
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
10 // the Free Software Foundation; either version 3 of the License, or
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
11 // (at your option) any later version.
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
12 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
13 // This program is distributed in the hope that it will be useful,
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
16 // GNU General Public License for more details.
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
17 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
18 // You should have received a copy of the GNU General Public License
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
19 // along with this program; if not, see <http://www.gnu.org/licenses/>.
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
20 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
21 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
22 // covm: in-product of matrices, NaN are skipped.
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
23 // usage:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
24 // [cc,nn] = covm_mex(X,Y,flag,W);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
25 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
26 // Input:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
27 // - X:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
28 // - Y: [optional], if empty, Y=X;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
29 // - flag: if not empty, it is set to 1 if some NaN was observed
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
30 // - W: weight vector to compute weighted correlation
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
31 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
32 // Output:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
33 // - CC = X' * sparse(diag(W)) * Y while NaN's are skipped
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
34 // - NN = real(~isnan(X)')*sparse(diag(W))*real(~isnan(Y)) count of valid (non-NaN) elements
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
35 // computed more efficiently
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
36 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
37 // $Id$
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
38 // Copyright (C) 2009 Alois Schloegl <a.schloegl@ieee.org>
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
39 // This function is part of the NaN-toolbox
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
40 // http://hci.tugraz.at/~schloegl/matlab/NaN/
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
41 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
42 //-------------------------------------------------------------------
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
43 */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
44
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
45 #include <inttypes.h>
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
46 #include <math.h>
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
47 #include "mex.h"
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
48
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
49 /*#define NO_FLAG*/
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
50
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
51
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
52 void mexFunction(int POutputCount, mxArray* POutput[], int PInputCount, const mxArray *PInputs[])
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
53 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
54 double *X0,*Y0=NULL,*X,*Y,*W=NULL;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
55 double *CC;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
56 double *NN=NULL;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
57
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
58 size_t rX,cX,rY,cY,nW = 0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
59 size_t i,j,k; // running indices
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
60 char flag_isNaN = 0, flag_speed=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
61
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
62
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
63 /*********** check input arguments *****************/
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
64
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
65 // check for proper number of input and output arguments
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
66 if ((PInputCount <= 0) || (PInputCount > 5)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
67 mexPrintf("usage: [CC,NN] = covm_mex(X [,Y [,flag [,W [,'E']]]])\n\n");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
68 mexPrintf("Do not use COVM_MEX directly, use COVM instead. \n");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
69 /*
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
70 mexPrintf("\nCOVM_MEX computes the covariance matrix of real matrices and skips NaN's\n");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
71 mexPrintf("\t[CC,NN] = covm_mex(...)\n\t\t computes CC=X'*Y, NN contains the number of not-NaN elements\n");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
72 mexPrintf("\t\t CC./NN is the unbiased covariance matrix\n");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
73 mexPrintf("\t... = covm_mex(X,Y,...)\n\t\t computes CC=X'*sparse(diag(W))*Y, number of rows of X and Y must match\n");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
74 mexPrintf("\t... = covm_mex(X,[], ...)\n\t\t computes CC=X'*sparse(diag(W))*X\n");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
75 mexPrintf("\t... = covm_mex(...,flag,...)\n\t\t if flag is not empty, it is set to 1 if some NaN occured in X or Y\n");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
76 mexPrintf("\t... = covm_mex(...,W)\n\t\t W to compute weighted covariance, number of elements must match the number of rows of X\n");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
77 mexPrintf("\t\t if isempty(W), all weights are 1\n");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
78 mexPrintf("\t[CC,NN]=covm_mex(X,Y,flag,W)\n");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
79 */ return;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
80 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
81 if (POutputCount > 2)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
82 mexErrMsgTxt("covm.MEX has 1 to 2 output arguments.");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
83
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
84
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
85 // get 1st argument
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
86 if(mxIsDouble(PInputs[0]) && !mxIsComplex(PInputs[0]))
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
87 X0 = mxGetPr(PInputs[0]);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
88 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
89 mexErrMsgTxt("First argument must be REAL/DOUBLE.");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
90 rX = mxGetM(PInputs[0]);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
91 cX = mxGetN(PInputs[0]);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
92
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
93 // get 2nd argument
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
94 if (PInputCount > 1) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
95 if (!mxGetNumberOfElements(PInputs[1]))
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
96 ; // Y0 = NULL;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
97
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
98 else if (mxIsDouble(PInputs[1]) && !mxIsComplex(PInputs[1]))
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
99 Y0 = mxGetPr(PInputs[1]);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
100
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
101 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
102 mexErrMsgTxt("Second argument must be REAL/DOUBLE.");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
103 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
104
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
105
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
106 // get weight vector for weighted sumskipnan
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
107 if (PInputCount > 3) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
108 // get 4th argument
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
109 nW = mxGetNumberOfElements(PInputs[3]);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
110 if (!nW)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
111 ;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
112 else if (nW == rX)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
113 W = mxGetPr(PInputs[3]);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
114 else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
115 mexErrMsgTxt("number of elements in W must match numbers of rows in X");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
116 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
117
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
118 int ACC_LEVEL = 0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
119 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
120 mxArray *LEVEL = NULL;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
121 int s = mexCallMATLAB(1, &LEVEL, 0, NULL, "flag_accuracy_level");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
122 if (!s) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
123 ACC_LEVEL = (int) mxGetScalar(LEVEL);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
124 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
125 mxDestroyArray(LEVEL);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
126 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
127 // mexPrintf("Accuracy Level=%i\n",ACC_LEVEL);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
128
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
129 if (Y0==NULL) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
130 Y0 = X0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
131 rY = rX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
132 cY = cX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
133 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
134 else {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
135 rY = mxGetM(PInputs[1]);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
136 cY = mxGetN(PInputs[1]);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
137 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
138 if (rX != rY)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
139 mexErrMsgTxt("number of rows in X and Y do not match");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
140
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
141 /*********** create output arguments *****************/
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
142
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
143 POutput[0] = mxCreateDoubleMatrix(cX, cY, mxREAL);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
144 CC = mxGetPr(POutput[0]);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
145
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
146 if (POutputCount > 1) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
147 POutput[1] = mxCreateDoubleMatrix(cX, cY, mxREAL);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
148 NN = mxGetPr(POutput[1]);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
149 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
150
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
151
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
152 /*********** compute covariance *****************/
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
153
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
154 #if 0
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
155 /*------ version 1 ---------------------
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
156 this solution is slower than the alternative solution below
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
157 for transposed matrices, this might be faster.
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
158 */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
159 for (k=0; k<rX; k++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
160 double w;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
161 if (W) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
162 w = W[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
163 for (i=0; i<cX; i++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
164 double x = X0[k+i*rX];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
165 if (isnan(x)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
166 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
167 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
168 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
169 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
170 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
171 for (j=0; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
172 double y = Y0[k+j*rY];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
173 if (isnan(y)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
174 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
175 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
176 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
177 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
178 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
179 CC[i+j*cX] += x*y*w;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
180 if (NN != NULL)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
181 NN[i+j*cX] += w;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
182 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
183 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
184 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
185 else for (i=0; i<cX; i++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
186 double x = X0[k+i*rX];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
187 if (isnan(x)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
188 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
189 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
190 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
191 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
192 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
193 for (j=0; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
194 double y = Y0[k+j*rY];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
195 if (isnan(y)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
196 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
197 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
198 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
199 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
200 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
201 CC[i+j*cX] += x*y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
202 if (NN != NULL)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
203 NN[i+j*cX] += 1.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
204 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
205 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
206 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
207
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
208 #else
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
209 if (ACC_LEVEL == 0) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
210 /*------ version 2 ---------------------
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
211 using naive summation with double accuracy [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
212 */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
213 if ( (X0 != Y0) || (cX != cY) )
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
214 /******** X!=Y, output is not symetric *******/
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
215 if (W) /* weighted version */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
216 for (i=0; i<cX; i++)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
217 for (j=0; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
218 X = X0+i*rX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
219 Y = Y0+j*rY;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
220 double cc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
221 double nn=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
222 for (k=0; k<rX; k++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
223 double z = X[k]*Y[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
224 if (isnan(z)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
225 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
226 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
227 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
228 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
229 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
230 cc += z*W[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
231 nn += W[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
232 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
233 CC[i+j*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
234 if (NN != NULL)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
235 NN[i+j*cX] = nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
236 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
237 else /* no weights, all weights are 1 */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
238 for (i=0; i<cX; i++)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
239 for (j=0; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
240 X = X0+i*rX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
241 Y = Y0+j*rY;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
242 double cc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
243 size_t nn=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
244 for (k=0; k<rX; k++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
245 double z = X[k]*Y[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
246 if (isnan(z)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
247 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
248 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
249 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
250 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
251 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
252 cc += z;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
253 nn++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
254 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
255 CC[i+j*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
256 if (NN != NULL)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
257 NN[i+j*cX] = (double)nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
258 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
259 else // if (X0==Y0) && (cX==cY)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
260 /******** X==Y, output is symetric *******/
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
261 if (W) /* weighted version */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
262 for (i=0; i<cX; i++)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
263 for (j=i; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
264 X = X0+i*rX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
265 Y = Y0+j*rY;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
266 double cc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
267 double nn=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
268 for (k=0; k<rX; k++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
269 double z = X[k]*Y[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
270 if (isnan(z)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
271 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
272 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
273 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
274 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
275 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
276 cc += z*W[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
277 nn += W[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
278 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
279 CC[i+j*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
280 CC[j+i*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
281 if (NN != NULL) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
282 NN[i+j*cX] = nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
283 NN[j+i*cX] = nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
284 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
285 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
286 else /* no weights, all weights are 1 */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
287 for (i=0; i<cX; i++)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
288 for (j=i; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
289 X = X0+i*rX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
290 Y = Y0+j*rY;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
291 double cc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
292 size_t nn=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
293 for (k=0; k<rX; k++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
294 double z = X[k]*Y[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
295 if (isnan(z)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
296 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
297 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
298 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
299 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
300 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
301 cc += z;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
302 nn++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
303 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
304 CC[i+j*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
305 CC[j+i*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
306 if (NN != NULL) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
307 NN[i+j*cX] = (double)nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
308 NN[j+i*cX] = (double)nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
309 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
310 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
311
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
312 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
313 else if (ACC_LEVEL == 1) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
314 /*------ version 2 ---------------------
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
315 using naive summation with extended accuracy [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
316 */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
317 if ( (X0 != Y0) || (cX != cY) )
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
318 /******** X!=Y, output is not symetric *******/
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
319 if (W) /* weighted version */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
320 for (i=0; i<cX; i++)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
321 for (j=0; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
322 X = X0+i*rX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
323 Y = Y0+j*rY;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
324 long double cc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
325 long double nn=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
326 for (k=0; k<rX; k++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
327 long double z = ((long double)X[k])*Y[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
328 if (isnan(z)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
329 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
330 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
331 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
332 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
333 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
334 cc += z*W[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
335 nn += W[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
336 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
337 CC[i+j*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
338 if (NN != NULL)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
339 NN[i+j*cX] = nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
340 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
341 else /* no weights, all weights are 1 */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
342 for (i=0; i<cX; i++)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
343 for (j=0; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
344 X = X0+i*rX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
345 Y = Y0+j*rY;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
346 long double cc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
347 size_t nn=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
348 for (k=0; k<rX; k++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
349 long double z = ((long double)X[k])*Y[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
350 if (isnan(z)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
351 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
352 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
353 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
354 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
355 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
356 cc += z;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
357 nn++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
358 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
359 CC[i+j*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
360 if (NN != NULL)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
361 NN[i+j*cX] = (double)nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
362 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
363 else // if (X0==Y0) && (cX==cY)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
364 /******** X==Y, output is symetric *******/
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
365 if (W) /* weighted version */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
366 for (i=0; i<cX; i++)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
367 for (j=i; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
368 X = X0+i*rX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
369 Y = Y0+j*rY;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
370 long double cc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
371 long double nn=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
372 for (k=0; k<rX; k++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
373 long double z = ((long double)X[k])*Y[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
374 if (isnan(z)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
375 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
376 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
377 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
378 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
379 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
380 cc += z*W[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
381 nn += W[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
382 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
383 CC[i+j*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
384 CC[j+i*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
385 if (NN != NULL) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
386 NN[i+j*cX] = nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
387 NN[j+i*cX] = nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
388 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
389 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
390 else /* no weights, all weights are 1 */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
391 for (i=0; i<cX; i++)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
392 for (j=i; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
393 X = X0+i*rX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
394 Y = Y0+j*rY;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
395 long double cc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
396 size_t nn=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
397 for (k=0; k<rX; k++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
398 long double z = ((long double)X[k])*Y[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
399 if (isnan(z)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
400 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
401 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
402 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
403 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
404 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
405 cc += z;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
406 nn++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
407 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
408 CC[i+j*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
409 CC[j+i*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
410 if (NN != NULL) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
411 NN[i+j*cX] = (double)nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
412 NN[j+i*cX] = (double)nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
413 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
414 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
415
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
416 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
417 else if (ACC_LEVEL == 3) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
418 /*------ version 3 ---------------------
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
419 using Kahan's summation with extended (long double) accuracy [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
420 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
421 However, first test show an increase in computational time of only about 25 %.
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
422
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
423 [1] David Goldberg,
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
424 What Every Computer Scientist Should Know About Floating-Point Arithmetic
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
425 ACM Computing Surveys, Vol 23, No 1, March 1991
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
426 */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
427 if ( (X0 != Y0) || (cX != cY) )
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
428 /******** X!=Y, output is not symetric *******/
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
429 if (W) /* weighted version */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
430 for (i=0; i<cX; i++)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
431 for (j=0; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
432 X = X0+i*rX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
433 Y = Y0+j*rY;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
434 long double cc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
435 long double nn=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
436 long double rc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
437 long double rn=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
438 for (k=0; k<rX; k++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
439 long double t,y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
440 long double z = ((long double)X[k])*Y[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
441 if (isnan(z)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
442 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
443 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
444 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
445 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
446 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
447 // cc += z*W[k]; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
448 y = z*W[k]-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
449 t = cc+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
450 rc= (t-cc)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
451 cc= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
452
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
453 // nn += W[k]; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
454 y = z*W[k]-rn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
455 t = nn+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
456 rn= (t-nn)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
457 nn= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
458 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
459 CC[i+j*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
460 if (NN != NULL)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
461 NN[i+j*cX] = nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
462 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
463 else /* no weights, all weights are 1 */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
464 for (i=0; i<cX; i++)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
465 for (j=0; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
466 X = X0+i*rX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
467 Y = Y0+j*rY;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
468 long double cc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
469 long double rc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
470 size_t nn=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
471 for (k=0; k<rX; k++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
472 long double t,y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
473 long double z = ((long double)X[k])*Y[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
474 if (isnan(z)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
475 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
476 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
477 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
478 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
479 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
480 // cc += z; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
481 y = z-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
482 t = cc+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
483 rc= (t-cc)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
484 cc= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
485
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
486 nn++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
487 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
488 CC[i+j*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
489 if (NN != NULL)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
490 NN[i+j*cX] = (double)nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
491 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
492 else // if (X0==Y0) && (cX==cY)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
493 /******** X==Y, output is symetric *******/
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
494 if (W) /* weighted version */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
495 for (i=0; i<cX; i++)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
496 for (j=i; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
497 X = X0+i*rX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
498 Y = Y0+j*rY;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
499 long double cc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
500 long double nn=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
501 long double rc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
502 long double rn=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
503 for (k=0; k<rX; k++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
504 long double t,y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
505 long double z = ((long double)X[k])*Y[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
506 if (isnan(z)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
507 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
508 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
509 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
510 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
511 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
512 // cc += z*W[k]; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
513 y = z*W[k]-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
514 t = cc+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
515 rc= (t-cc)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
516 cc= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
517
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
518 // nn += W[k]; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
519 y = z*W[k]-rn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
520 t = nn+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
521 rn= (t-nn)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
522 nn= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
523 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
524 CC[i+j*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
525 CC[j+i*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
526 if (NN != NULL) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
527 NN[i+j*cX] = nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
528 NN[j+i*cX] = nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
529 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
530 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
531 else /* no weights, all weights are 1 */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
532 for (i=0; i<cX; i++)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
533 for (j=i; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
534 X = X0+i*rX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
535 Y = Y0+j*rY;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
536 long double cc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
537 long double rc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
538 size_t nn=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
539 for (k=0; k<rX; k++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
540 long double t,y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
541 long double z = ((long double)X[k])*Y[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
542 if (isnan(z)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
543 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
544 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
545 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
546 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
547 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
548 // cc += z; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
549 y = z-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
550 t = cc+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
551 rc= (t-cc)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
552 cc= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
553
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
554 nn++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
555 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
556 CC[i+j*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
557 CC[j+i*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
558 if (NN != NULL) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
559 NN[i+j*cX] = (double)nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
560 NN[j+i*cX] = (double)nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
561 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
562 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
563 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
564 else if (ACC_LEVEL == 2) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
565 /*------ version 3 ---------------------
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
566 using Kahan's summation with double accuracy [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
567 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
568 However, first test show an increase in computational time of only about 25 %.
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
569
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
570 [1] David Goldberg,
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
571 What Every Computer Scientist Should Know About Floating-Point Arithmetic
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
572 ACM Computing Surveys, Vol 23, No 1, March 1991
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
573 */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
574 if ( (X0 != Y0) || (cX != cY) )
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
575 /******** X!=Y, output is not symetric *******/
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
576 if (W) /* weighted version */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
577 for (i=0; i<cX; i++)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
578 for (j=0; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
579 X = X0+i*rX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
580 Y = Y0+j*rY;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
581 double cc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
582 double nn=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
583 double rc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
584 double rn=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
585 for (k=0; k<rX; k++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
586 double t,y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
587 double z = X[k]*Y[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
588 if (isnan(z)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
589 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
590 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
591 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
592 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
593 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
594 // cc += z*W[k]; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
595 y = z*W[k]-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
596 t = cc+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
597 rc= (t-cc)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
598 cc= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
599
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
600 // nn += W[k]; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
601 y = z*W[k]-rn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
602 t = nn+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
603 rn= (t-nn)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
604 nn= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
605 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
606 CC[i+j*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
607 if (NN != NULL)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
608 NN[i+j*cX] = nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
609 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
610 else /* no weights, all weights are 1 */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
611 for (i=0; i<cX; i++)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
612 for (j=0; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
613 X = X0+i*rX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
614 Y = Y0+j*rY;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
615 double cc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
616 double rc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
617 size_t nn=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
618 for (k=0; k<rX; k++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
619 double t,y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
620 double z = X[k]*Y[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
621 if (isnan(z)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
622 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
623 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
624 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
625 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
626 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
627 // cc += z; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
628 y = z-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
629 t = cc+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
630 rc= (t-cc)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
631 cc= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
632
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
633 nn++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
634 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
635 CC[i+j*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
636 if (NN != NULL)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
637 NN[i+j*cX] = (double)nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
638 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
639 else // if (X0==Y0) && (cX==cY)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
640 /******** X==Y, output is symetric *******/
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
641 if (W) /* weighted version */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
642 for (i=0; i<cX; i++)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
643 for (j=i; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
644 X = X0+i*rX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
645 Y = Y0+j*rY;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
646 double cc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
647 double nn=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
648 double rc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
649 double rn=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
650 for (k=0; k<rX; k++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
651 double t,y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
652 double z = X[k]*Y[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
653 if (isnan(z)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
654 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
655 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
656 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
657 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
658 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
659 // cc += z*W[k]; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
660 y = z*W[k]-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
661 t = cc+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
662 rc= (t-cc)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
663 cc= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
664
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
665 // nn += W[k]; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
666 y = z*W[k]-rn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
667 t = nn+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
668 rn= (t-nn)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
669 nn= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
670 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
671 CC[i+j*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
672 CC[j+i*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
673 if (NN != NULL) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
674 NN[i+j*cX] = nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
675 NN[j+i*cX] = nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
676 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
677 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
678 else /* no weights, all weights are 1 */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
679 for (i=0; i<cX; i++)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
680 for (j=i; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
681 X = X0+i*rX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
682 Y = Y0+j*rY;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
683 double cc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
684 double rc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
685 size_t nn=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
686 for (k=0; k<rX; k++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
687 double t,y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
688 double z = X[k]*Y[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
689 if (isnan(z)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
690 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
691 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
692 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
693 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
694 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
695 // cc += z; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
696 y = z-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
697 t = cc+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
698 rc= (t-cc)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
699 cc= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
700
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
701 nn++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
702 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
703 CC[i+j*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
704 CC[j+i*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
705 if (NN != NULL) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
706 NN[i+j*cX] = (double)nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
707 NN[j+i*cX] = (double)nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
708 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
709 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
710 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
711
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
712
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
713 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
714 //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
715 /* 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
716 this requires that the input argument is a non-empty variable
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
717 */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
718 if (flag_isNaN && (PInputCount > 2) && mxGetNumberOfElements(PInputs[2])) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
719 // set FLAG_NANS_OCCURED
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
720 switch (mxGetClassID(PInputs[2])) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
721 case mxLOGICAL_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
722 case mxCHAR_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
723 case mxINT8_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
724 case mxUINT8_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
725 *(uint8_t*)mxGetData(PInputs[2]) = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
726 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
727 case mxDOUBLE_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
728 *(double*)mxGetData(PInputs[2]) = 1.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
729 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
730 case mxSINGLE_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
731 *(float*)mxGetData(PInputs[2]) = 1.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
732 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
733 case mxINT16_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
734 case mxUINT16_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
735 *(uint16_t*)mxGetData(PInputs[2]) = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
736 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
737 case mxINT32_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
738 case mxUINT32_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
739 *(uint32_t*)mxGetData(PInputs[2])= 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
740 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
741 case mxINT64_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
742 case mxUINT64_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
743 *(uint64_t*)mxGetData(PInputs[2]) = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
744 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
745 case mxFUNCTION_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
746 case mxUNKNOWN_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
747 case mxCELL_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
748 case mxSTRUCT_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
749 mexPrintf("Type of 3rd input argument cannot be used to return status of NaN occurence.");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
750 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
751 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
752 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
753 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
754 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
755