Mercurial > forge

annotate extra/NaN/src/covm_mex.cpp @ 6585:ae521dec5b54 octave-forge

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