Mercurial > forge

annotate extra/NaN/src/covm_mex.cpp @ 12691:6d6285a2a633 octave-forge

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