Mercurial > forge

annotate extra/NaN/src/covm_mex.cpp @ 7889:c101c486d80a octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
fix web address
author schloegl
date Thu, 27 Jan 2011 17:10:36 +0000
parents b9f35668b55e
children db5092052107
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 // C-MEX implementation of COVM - 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 // covm: in-product of matrices, NaN are skipped.
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
21 // usage:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
22 // [cc,nn] = covm_mex(X,Y,flag,W);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
23 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
24 // Input:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
25 // - X:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
26 // - Y: [optional], if empty, Y=X;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
27 // - flag: if not empty, it is set to 1 if some NaN was observed
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
28 // - W: weight vector to compute weighted correlation
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
29 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
30 // Output:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
31 // - CC = X' * sparse(diag(W)) * Y while NaN's are skipped
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
32 // - 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
33 // computed more efficiently
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
34 //
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
35 // $Id$
7889
c101c486d80a fix web address
schloegl
parents: 7888
diff changeset
36 // Copyright (C) 2009,2010,2011 Alois Schloegl <a.schloegl@ieee.org>
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
37 // This function is part of the NaN-toolbox
7889
c101c486d80a fix web address
schloegl
parents: 7888
diff changeset
38 // http://pub.ist.ac.at/~schloegl/matlab/NaN/
6549
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 */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
42
6585
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
43 #ifdef __GNUC__
7888
b9f35668b55e replace <inttypes.h> with <stdint.h>
schloegl
parents: 6585
diff changeset
44 #include <stdint.h>
6585
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
45 #endif
6549
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;
6585
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
55 double *CC,cc,nw;
6549
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;
6585
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
59 size_t i,j,k,nn;
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
60 char flag_isNaN = 0, flag_speed=0;
6585
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
61 int ACC_LEVEL;
6549
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
6585
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
118 #ifdef __GNUC__
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
119 ACC_LEVEL = 0;
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
120 {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
121 mxArray *LEVEL = NULL;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
122 int s = mexCallMATLAB(1, &LEVEL, 0, NULL, "flag_accuracy_level");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
123 if (!s) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
124 ACC_LEVEL = (int) mxGetScalar(LEVEL);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
125 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
126 mxDestroyArray(LEVEL);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
127 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
128 // mexPrintf("Accuracy Level=%i\n",ACC_LEVEL);
6585
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
129 #endif
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
130 if (Y0==NULL) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
131 Y0 = X0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
132 rY = rX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
133 cY = cX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
134 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
135 else {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
136 rY = mxGetM(PInputs[1]);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
137 cY = mxGetN(PInputs[1]);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
138 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
139 if (rX != rY)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
140 mexErrMsgTxt("number of rows in X and Y do not match");
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
141
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
142 /*********** create output arguments *****************/
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
143
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
144 POutput[0] = mxCreateDoubleMatrix(cX, cY, mxREAL);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
145 CC = mxGetPr(POutput[0]);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
146
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
147 if (POutputCount > 1) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
148 POutput[1] = mxCreateDoubleMatrix(cX, cY, mxREAL);
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
149 NN = mxGetPr(POutput[1]);
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
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
153 /*********** compute covariance *****************/
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
154
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
155 #if 0
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
156 /*------ version 1 ---------------------
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
157 this solution is slower than the alternative solution below
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
158 for transposed matrices, this might be faster.
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
159 */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
160 for (k=0; k<rX; k++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
161 double w;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
162 if (W) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
163 w = W[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
164 for (i=0; i<cX; i++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
165 double x = X0[k+i*rX];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
166 if (isnan(x)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
167 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
168 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
169 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
170 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
171 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
172 for (j=0; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
173 double y = Y0[k+j*rY];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
174 if (isnan(y)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
175 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
176 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
177 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
178 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
179 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
180 CC[i+j*cX] += x*y*w;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
181 if (NN != NULL)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
182 NN[i+j*cX] += w;
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 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
186 else for (i=0; i<cX; i++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
187 double x = X0[k+i*rX];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
188 if (isnan(x)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
189 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
190 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
191 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
192 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
193 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
194 for (j=0; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
195 double y = Y0[k+j*rY];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
196 if (isnan(y)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
197 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
198 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
199 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
200 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
201 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
202 CC[i+j*cX] += x*y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
203 if (NN != NULL)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
204 NN[i+j*cX] += 1.0;
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
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
209 #else
6585
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
210 #ifdef __GNUC__
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
211 if (ACC_LEVEL == 0)
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
212 #endif
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
213 {
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
214 /*------ version 2 ---------------------
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
215 using naive summation with double accuracy [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
216 */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
217 if ( (X0 != Y0) || (cX != cY) )
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
218 /******** X!=Y, output is not symetric *******/
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
219 if (W) /* weighted version */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
220 for (i=0; i<cX; i++)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
221 for (j=0; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
222 X = X0+i*rX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
223 Y = Y0+j*rY;
6585
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
224 cc=0.0;
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
225 nw=0.0;
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
226 for (k=0; k<rX; k++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
227 double z = X[k]*Y[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
228 if (isnan(z)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
229 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
230 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
231 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
232 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
233 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
234 cc += z*W[k];
6585
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
235 nw += W[k];
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
236 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
237 CC[i+j*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
238 if (NN != NULL)
6585
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
239 NN[i+j*cX] = nw;
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
240 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
241 else /* no weights, all weights are 1 */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
242 for (i=0; i<cX; i++)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
243 for (j=0; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
244 X = X0+i*rX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
245 Y = Y0+j*rY;
6585
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
246 cc=0.0;
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
247 nn=0;
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
248 for (k=0; k<rX; k++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
249 double z = X[k]*Y[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
250 if (isnan(z)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
251 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
252 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
253 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
254 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
255 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
256 cc += z;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
257 nn++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
258 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
259 CC[i+j*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
260 if (NN != NULL)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
261 NN[i+j*cX] = (double)nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
262 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
263 else // if (X0==Y0) && (cX==cY)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
264 /******** X==Y, output is symetric *******/
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
265 if (W) /* weighted version */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
266 for (i=0; i<cX; i++)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
267 for (j=i; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
268 X = X0+i*rX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
269 Y = Y0+j*rY;
6585
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
270 cc=0.0;
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
271 nw=0.0;
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
272 for (k=0; k<rX; k++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
273 double z = X[k]*Y[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
274 if (isnan(z)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
275 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
276 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
277 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
278 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
279 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
280 cc += z*W[k];
6585
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
281 nw += W[k];
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
282 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
283 CC[i+j*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
284 CC[j+i*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
285 if (NN != NULL) {
6585
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
286 NN[i+j*cX] = nw;
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
287 NN[j+i*cX] = nw;
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
288 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
289 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
290 else /* no weights, all weights are 1 */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
291 for (i=0; i<cX; i++)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
292 for (j=i; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
293 X = X0+i*rX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
294 Y = Y0+j*rY;
6585
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
295 cc=0.0;
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
296 nn=0;
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
297 for (k=0; k<rX; k++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
298 double z = X[k]*Y[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
299 if (isnan(z)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
300 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
301 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
302 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
303 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
304 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
305 cc += z;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
306 nn++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
307 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
308 CC[i+j*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
309 CC[j+i*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
310 if (NN != NULL) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
311 NN[i+j*cX] = (double)nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
312 NN[j+i*cX] = (double)nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
313 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
314 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
315
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
316 }
6585
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
317
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
318 #ifdef __GNUC__
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
319
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
320 else if (ACC_LEVEL == 1) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
321 /*------ version 2 ---------------------
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
322 using naive summation with extended accuracy [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
323 */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
324 if ( (X0 != Y0) || (cX != cY) )
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
325 /******** X!=Y, output is not symetric *******/
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
326 if (W) /* weighted version */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
327 for (i=0; i<cX; i++)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
328 for (j=0; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
329 X = X0+i*rX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
330 Y = Y0+j*rY;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
331 long double cc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
332 long double nn=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
333 for (k=0; k<rX; k++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
334 long double z = ((long double)X[k])*Y[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
335 if (isnan(z)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
336 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
337 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
338 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
339 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
340 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
341 cc += z*W[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
342 nn += W[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
343 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
344 CC[i+j*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
345 if (NN != NULL)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
346 NN[i+j*cX] = nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
347 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
348 else /* no weights, all weights are 1 */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
349 for (i=0; i<cX; i++)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
350 for (j=0; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
351 X = X0+i*rX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
352 Y = Y0+j*rY;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
353 long double cc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
354 size_t nn=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
355 for (k=0; k<rX; k++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
356 long double z = ((long double)X[k])*Y[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
357 if (isnan(z)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
358 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
359 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
360 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
361 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
362 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
363 cc += z;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
364 nn++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
365 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
366 CC[i+j*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
367 if (NN != NULL)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
368 NN[i+j*cX] = (double)nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
369 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
370 else // if (X0==Y0) && (cX==cY)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
371 /******** X==Y, output is symetric *******/
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
372 if (W) /* weighted version */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
373 for (i=0; i<cX; i++)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
374 for (j=i; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
375 X = X0+i*rX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
376 Y = Y0+j*rY;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
377 long double cc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
378 long double nn=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
379 for (k=0; k<rX; k++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
380 long double z = ((long double)X[k])*Y[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
381 if (isnan(z)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
382 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
383 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
384 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
385 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
386 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
387 cc += z*W[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
388 nn += W[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
389 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
390 CC[i+j*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
391 CC[j+i*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
392 if (NN != NULL) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
393 NN[i+j*cX] = nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
394 NN[j+i*cX] = nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
395 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
396 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
397 else /* no weights, all weights are 1 */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
398 for (i=0; i<cX; i++)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
399 for (j=i; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
400 X = X0+i*rX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
401 Y = Y0+j*rY;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
402 long double cc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
403 size_t nn=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
404 for (k=0; k<rX; k++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
405 long double z = ((long double)X[k])*Y[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
406 if (isnan(z)) {
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 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
409 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
410 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
411 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
412 cc += z;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
413 nn++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
414 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
415 CC[i+j*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
416 CC[j+i*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
417 if (NN != NULL) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
418 NN[i+j*cX] = (double)nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
419 NN[j+i*cX] = (double)nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
420 }
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 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
424 else if (ACC_LEVEL == 3) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
425 /*------ version 3 ---------------------
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
426 using Kahan's summation with extended (long double) accuracy [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
427 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
428 However, first test show an increase in computational time of only about 25 %.
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
429
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
430 [1] David Goldberg,
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
431 What Every Computer Scientist Should Know About Floating-Point Arithmetic
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
432 ACM Computing Surveys, Vol 23, No 1, March 1991
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
433 */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
434 if ( (X0 != Y0) || (cX != cY) )
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
435 /******** X!=Y, output is not symetric *******/
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
436 if (W) /* weighted version */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
437 for (i=0; i<cX; i++)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
438 for (j=0; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
439 X = X0+i*rX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
440 Y = Y0+j*rY;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
441 long double cc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
442 long double nn=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
443 long double rc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
444 long double rn=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
445 for (k=0; k<rX; k++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
446 long double t,y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
447 long double z = ((long double)X[k])*Y[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
448 if (isnan(z)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
449 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
450 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
451 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
452 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
453 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
454 // cc += z*W[k]; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
455 y = z*W[k]-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
456 t = cc+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
457 rc= (t-cc)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
458 cc= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
459
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
460 // nn += W[k]; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
461 y = z*W[k]-rn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
462 t = nn+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
463 rn= (t-nn)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
464 nn= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
465 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
466 CC[i+j*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
467 if (NN != NULL)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
468 NN[i+j*cX] = nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
469 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
470 else /* no weights, all weights are 1 */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
471 for (i=0; i<cX; i++)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
472 for (j=0; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
473 X = X0+i*rX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
474 Y = Y0+j*rY;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
475 long double cc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
476 long double rc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
477 size_t nn=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
478 for (k=0; k<rX; k++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
479 long double t,y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
480 long double z = ((long double)X[k])*Y[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
481 if (isnan(z)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
482 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
483 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
484 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
485 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
486 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
487 // cc += z; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
488 y = z-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
489 t = cc+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
490 rc= (t-cc)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
491 cc= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
492
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
493 nn++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
494 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
495 CC[i+j*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
496 if (NN != NULL)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
497 NN[i+j*cX] = (double)nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
498 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
499 else // if (X0==Y0) && (cX==cY)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
500 /******** X==Y, output is symetric *******/
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
501 if (W) /* weighted version */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
502 for (i=0; i<cX; i++)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
503 for (j=i; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
504 X = X0+i*rX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
505 Y = Y0+j*rY;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
506 long double cc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
507 long double nn=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
508 long double rc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
509 long double rn=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
510 for (k=0; k<rX; k++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
511 long double t,y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
512 long double z = ((long double)X[k])*Y[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
513 if (isnan(z)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
514 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
515 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
516 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
517 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
518 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
519 // cc += z*W[k]; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
520 y = z*W[k]-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
521 t = cc+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
522 rc= (t-cc)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
523 cc= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
524
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
525 // nn += W[k]; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
526 y = z*W[k]-rn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
527 t = nn+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
528 rn= (t-nn)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
529 nn= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
530 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
531 CC[i+j*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
532 CC[j+i*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
533 if (NN != NULL) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
534 NN[i+j*cX] = nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
535 NN[j+i*cX] = nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
536 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
537 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
538 else /* no weights, all weights are 1 */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
539 for (i=0; i<cX; i++)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
540 for (j=i; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
541 X = X0+i*rX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
542 Y = Y0+j*rY;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
543 long double cc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
544 long double rc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
545 size_t nn=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
546 for (k=0; k<rX; k++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
547 long double t,y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
548 long double z = ((long double)X[k])*Y[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
549 if (isnan(z)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
550 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
551 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
552 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
553 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
554 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
555 // cc += z; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
556 y = z-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
557 t = cc+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
558 rc= (t-cc)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
559 cc= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
560
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
561 nn++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
562 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
563 CC[i+j*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
564 CC[j+i*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
565 if (NN != NULL) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
566 NN[i+j*cX] = (double)nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
567 NN[j+i*cX] = (double)nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
568 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
569 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
570 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
571 else if (ACC_LEVEL == 2) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
572 /*------ version 3 ---------------------
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
573 using Kahan's summation with double accuracy [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
574 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
575 However, first test show an increase in computational time of only about 25 %.
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
576
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
577 [1] David Goldberg,
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
578 What Every Computer Scientist Should Know About Floating-Point Arithmetic
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
579 ACM Computing Surveys, Vol 23, No 1, March 1991
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
580 */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
581 if ( (X0 != Y0) || (cX != cY) )
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
582 /******** X!=Y, output is not symetric *******/
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
583 if (W) /* weighted version */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
584 for (i=0; i<cX; i++)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
585 for (j=0; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
586 X = X0+i*rX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
587 Y = Y0+j*rY;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
588 double cc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
589 double nn=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
590 double rc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
591 double rn=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
592 for (k=0; k<rX; k++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
593 double t,y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
594 double z = X[k]*Y[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
595 if (isnan(z)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
596 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
597 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
598 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
599 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
600 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
601 // cc += z*W[k]; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
602 y = z*W[k]-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
603 t = cc+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
604 rc= (t-cc)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
605 cc= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
606
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
607 // nn += W[k]; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
608 y = z*W[k]-rn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
609 t = nn+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
610 rn= (t-nn)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
611 nn= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
612 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
613 CC[i+j*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
614 if (NN != NULL)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
615 NN[i+j*cX] = nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
616 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
617 else /* no weights, all weights are 1 */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
618 for (i=0; i<cX; i++)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
619 for (j=0; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
620 X = X0+i*rX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
621 Y = Y0+j*rY;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
622 double cc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
623 double rc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
624 size_t nn=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
625 for (k=0; k<rX; k++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
626 double t,y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
627 double z = X[k]*Y[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
628 if (isnan(z)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
629 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
630 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
631 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
632 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
633 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
634 // cc += z; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
635 y = z-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
636 t = cc+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
637 rc= (t-cc)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
638 cc= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
639
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
640 nn++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
641 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
642 CC[i+j*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
643 if (NN != NULL)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
644 NN[i+j*cX] = (double)nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
645 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
646 else // if (X0==Y0) && (cX==cY)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
647 /******** X==Y, output is symetric *******/
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
648 if (W) /* weighted version */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
649 for (i=0; i<cX; i++)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
650 for (j=i; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
651 X = X0+i*rX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
652 Y = Y0+j*rY;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
653 double cc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
654 double nn=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
655 double rc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
656 double rn=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
657 for (k=0; k<rX; k++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
658 double t,y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
659 double z = X[k]*Y[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
660 if (isnan(z)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
661 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
662 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
663 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
664 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
665 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
666 // cc += z*W[k]; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
667 y = z*W[k]-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
668 t = cc+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
669 rc= (t-cc)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
670 cc= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
671
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
672 // nn += W[k]; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
673 y = z*W[k]-rn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
674 t = nn+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
675 rn= (t-nn)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
676 nn= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
677 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
678 CC[i+j*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
679 CC[j+i*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
680 if (NN != NULL) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
681 NN[i+j*cX] = nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
682 NN[j+i*cX] = nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
683 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
684 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
685 else /* no weights, all weights are 1 */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
686 for (i=0; i<cX; i++)
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
687 for (j=i; j<cY; j++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
688 X = X0+i*rX;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
689 Y = Y0+j*rY;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
690 double cc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
691 double rc=0.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
692 size_t nn=0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
693 for (k=0; k<rX; k++) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
694 double t,y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
695 double z = X[k]*Y[k];
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
696 if (isnan(z)) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
697 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
698 flag_isNaN = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
699 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
700 continue;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
701 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
702 // cc += z; [1]
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
703 y = z-rc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
704 t = cc+y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
705 rc= (t-cc)-y;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
706 cc= t;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
707
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
708 nn++;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
709 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
710 CC[i+j*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
711 CC[j+i*cX] = cc;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
712 if (NN != NULL) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
713 NN[i+j*cX] = (double)nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
714 NN[j+i*cX] = (double)nn;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
715 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
716 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
717 }
6585
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
718 #endif
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
719
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
720 #ifndef NO_FLAG
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
721 //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
722 /* 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
723 this requires that the input argument is a non-empty variable
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
724 */
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
725 if (flag_isNaN && (PInputCount > 2) && mxGetNumberOfElements(PInputs[2])) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
726 // set FLAG_NANS_OCCURED
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
727 switch (mxGetClassID(PInputs[2])) {
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
728 case mxDOUBLE_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
729 *(double*)mxGetData(PInputs[2]) = 1.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
730 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
731 case mxSINGLE_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
732 *(float*)mxGetData(PInputs[2]) = 1.0;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
733 break;
6585
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
734 case mxLOGICAL_CLASS:
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
735 case mxCHAR_CLASS:
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
736 case mxINT8_CLASS:
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
737 case mxUINT8_CLASS:
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
738 *(char*)mxGetData(PInputs[2]) = 1;
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
739 break;
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
740 #ifdef __GNUC__
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
741 case mxINT16_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
742 case mxUINT16_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
743 *(uint16_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 mxINT32_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
746 case mxUINT32_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
747 *(uint32_t*)mxGetData(PInputs[2])= 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
748 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
749 case mxINT64_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
750 case mxUINT64_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
751 *(uint64_t*)mxGetData(PInputs[2]) = 1;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
752 break;
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
753 case mxFUNCTION_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
754 case mxUNKNOWN_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
755 case mxCELL_CLASS:
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
756 case mxSTRUCT_CLASS:
6585
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
757 #endif
ae521dec5b54 partial support of lcc-win compiler
schloegl
parents: 6549
diff changeset
758 default:
6549
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
759 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
760 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
761 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
762 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
763 #endif
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
764 }
41e9854fe26d use *.cpp instead of *.c
schloegl
parents:
diff changeset
765