00001 #include "bnengine.h"
00002 #include "strctbse.h"
00003 #include "datahndl.h"
00004 #include "math.h"
00005 #include "blas.h"
00006
00007 void CAMbinaryEngine::doubleAequalToAlpha(CAMstructureBase &Astructure, CAMdataHandler &Adata,
00008 double alpha)
00009 {
00010 double* AdataP = (double*)Adata.getDataPointer();
00011
00012 #ifndef __NO_BLAS__
00013 long strideA;
00014 long Size;
00015 long longZero = 0;
00016 #endif
00017
00018 if(Astructure.isSubset() == 0)
00019 {
00020 #ifdef __NO_BLAS__
00021 double* aptr;
00022 for(aptr = AdataP; aptr < AdataP + Adata.getDataSize(); aptr++)
00023 { *aptr = alpha; };
00024 #else
00025 strideA = 1;
00026 Size = Adata.getDataSize();
00027 dcopy_(&Size,&alpha,&longZero, AdataP, &strideA);
00028 #endif
00029 }
00030 else
00031 {
00032
00033
00034
00035 long AloopDimension;
00036 long AOffset;
00037 MVAlongBase ACount;
00038 MVAlongBase AStride;
00039
00040 Astructure.getNormalizedLoops(AloopDimension,AOffset, ACount, AStride);
00041
00042 register long i1, i2, i3, i4, i5, i6, i7;
00043
00044 double* AOffptr = AdataP + AOffset;
00045
00046 register double* Ap1;
00047 register double* Ap2;
00048 register double* Ap3;
00049 register double* Ap4;
00050 register double* Ap5;
00051 register double* Ap6;
00052
00053 #ifndef __NO_BLAS__
00054 strideA = AStride[0];
00055 Size = ACount[0];
00056 #endif
00057
00058 switch(AloopDimension)
00059 {
00060 case 1 :
00061 #ifdef __NO_BLAS__
00062 for(i1=0; i1 < ACount[0]*AStride[0]; i1+=AStride[0])
00063 {*(AOffptr + i1) = alpha;};
00064 #else
00065 dcopy_(&Size,&alpha,&longZero, AOffptr, &strideA);
00066 #endif
00067
00068 break;
00069
00070 case 2 :
00071 for(i2=0; i2 < ACount[1]*AStride[1]; i2+=AStride[1])
00072 { Ap1 = AOffptr + i2;
00073 #ifdef __NO_BLAS__
00074 for(i1=0; i1 < ACount[0]*AStride[0]; i1+=AStride[0])
00075 {*(Ap1 + i1) = alpha;};
00076 #else
00077 dcopy_(&Size,&alpha,&longZero, Ap1, &strideA);
00078 #endif
00079 };
00080 break;
00081
00082 case 3 :
00083 for(i3=0; i3 < ACount[2]*AStride[2]; i3+=AStride[2])
00084 { Ap2 = AOffptr + i3;
00085 for(i2=0; i2 < ACount[1]*AStride[1]; i2+=AStride[1])
00086 { Ap1 = Ap2 + i2;
00087 #ifdef __NO_BLAS__
00088 for(i1=0; i1 < ACount[0]*AStride[0]; i1+=AStride[0])
00089 {*(Ap1 + i1) = alpha;};
00090 #else
00091 dcopy_(&Size,&alpha,&longZero, Ap1, &strideA);
00092 #endif
00093 }};
00094 break;
00095
00096 case 4 :
00097 for(i4=0; i4 < ACount[3]*AStride[3]; i4+=AStride[3])
00098 { Ap3 = AOffptr + i4;
00099 for(i3=0; i3 < ACount[2]*AStride[2]; i3+=AStride[2])
00100 { Ap2 = Ap3 + i3;
00101 for(i2=0; i2 < ACount[1]*AStride[1]; i2+=AStride[1])
00102 { Ap1 = Ap2 + i2;
00103 #ifdef __NO_BLAS__
00104 for(i1=0; i1 < ACount[0]*AStride[0]; i1+=AStride[0])
00105 {*(Ap1 + i1) = alpha;};
00106 #else
00107 dcopy_(&Size,&alpha,&longZero, Ap1, &strideA);
00108 #endif
00109 }}};
00110 break;
00111
00112 case 5 :
00113 for(i5=0; i5 < ACount[4]*AStride[4]; i5+=AStride[4])
00114 { Ap4 = AOffptr + i5;
00115 for(i4=0; i4 < ACount[3]*AStride[3]; i4+=AStride[3])
00116 { Ap3 = Ap4 + i4;
00117 for(i3=0; i3 < ACount[2]*AStride[2]; i3+=AStride[2])
00118 { Ap2 = Ap3 + i3;
00119 for(i2=0; i2 < ACount[1]*AStride[1]; i2+=AStride[1])
00120 { Ap1 = Ap2 + i2;
00121 #ifdef __NO_BLAS__
00122 for(i1=0; i1 < ACount[0]*AStride[0]; i1+=AStride[0])
00123 {*(Ap1 + i1) = alpha;};
00124 #else
00125 dcopy_(&Size,&alpha,&longZero, Ap1, &strideA);
00126 #endif
00127 }}}};
00128 break;
00129
00130 case 6 :
00131 for(i6=0; i6 < ACount[5]*AStride[5]; i6+=AStride[5])
00132 { Ap5 = AOffptr+ i6;
00133 for(i5=0; i5 < ACount[4]*AStride[4]; i5+=AStride[4])
00134 { Ap4 = Ap5 + i5;
00135 for(i4=0; i4 < ACount[3]*AStride[3]; i4+=AStride[3])
00136 { Ap3 = Ap4 + i4;
00137 for(i3=0; i3 < ACount[2]*AStride[2]; i3+=AStride[2])
00138 { Ap2 = Ap3 + i3;
00139 for(i2=0; i2 < ACount[1]*AStride[1]; i2+=AStride[1])
00140 { Ap1 = Ap2 + i2;
00141 #ifdef __NO_BLAS__
00142 for(i1=0; i1 < ACount[0]*AStride[0]; i1+=AStride[0])
00143 {*(Ap1 + i1) = alpha;};
00144 #else
00145 dcopy_(&Size,&alpha,&longZero, Ap1, &strideA);
00146 #endif
00147 }}}}};
00148 break;
00149
00150 case 7 :
00151 for(i7=0; i7 < ACount[6]*AStride[6]; i7+=AStride[6])
00152 { Ap6 = AOffptr + i7;
00153 for(i6=0; i6 < ACount[5]*AStride[5]; i6+=AStride[5])
00154 { Ap5 = Ap6 + i6;
00155 for(i5=0; i5 < ACount[4]*AStride[4]; i5+=AStride[4])
00156 { Ap4 = Ap5 + i5;
00157 for(i4=0; i4 < ACount[3]*AStride[3]; i4+=AStride[3])
00158 { Ap3 = Ap4 + i4;
00159 for(i3=0; i3 < ACount[2]*AStride[2]; i3+=AStride[2])
00160 { Ap2 = Ap3 + i3;
00161 for(i2=0; i2 < ACount[1]*AStride[1]; i2+=AStride[1])
00162 { Ap1 = Ap2 + i2;
00163 #ifdef __NO_BLAS__
00164 for(i1=0; i1 < ACount[0]*AStride[0]; i1+=AStride[0])
00165 {*(Ap1 + i1) = alpha;};
00166 #else
00167 dcopy_(&Size,&alpha,&longZero, Ap1, &strideA);
00168 #endif
00169 }}}}}};
00170 break;
00171
00172 }
00173 }
00174
00175 }
00176
00177
00178 void CAMbinaryEngine::doubleAplusEqualAlpha(CAMstructureBase &Astructure, CAMdataHandler &Adata,
00179 double alpha)
00180 {
00181 double* AdataP = (double*)Adata.getDataPointer();
00182
00183 #ifndef __NO_BLAS__
00184 long strideA;
00185 long Size;
00186 long longZero = 0;
00187 double dOne = 1.0;
00188 #endif
00189
00190 if(Astructure.isSubset() == 0)
00191 {
00192 #ifdef __NO_BLAS__
00193 double* aptr;
00194 for(aptr = AdataP; aptr < AdataP + Adata.getDataSize(); aptr++)
00195 { *aptr += alpha; };
00196 #else
00197 strideA = 1;
00198 Size = Adata.getDataSize();
00199 daxpy_(&Size,&dOne, &alpha,&longZero, AdataP, &strideA);
00200 #endif
00201 }
00202 else
00203 {
00204
00205
00206
00207 long AloopDimension;
00208 long AOffset;
00209 MVAlongBase ACount;
00210 MVAlongBase AStride;
00211
00212 Astructure.getNormalizedLoops(AloopDimension,AOffset, ACount, AStride);
00213
00214 register long i1, i2, i3, i4, i5, i6, i7;
00215
00216 double* AOffptr = AdataP + AOffset;
00217
00218 register double* Ap1;
00219 register double* Ap2;
00220 register double* Ap3;
00221 register double* Ap4;
00222 register double* Ap5;
00223 register double* Ap6;
00224
00225 #ifndef __NO_BLAS__
00226 strideA = AStride[0];
00227 Size = ACount[0];
00228 #endif
00229
00230 switch(AloopDimension)
00231 {
00232 case 1 :
00233 #ifdef __NO_BLAS__
00234 for(i1=0; i1 < ACount[0]*AStride[0]; i1+=AStride[0])
00235 {*(AOffptr + i1) += alpha;};
00236 #else
00237 daxpy_(&Size,&dOne,&alpha,&longZero, AOffptr, &strideA);
00238 #endif
00239
00240 break;
00241
00242 case 2 :
00243 for(i2=0; i2 < ACount[1]*AStride[1]; i2+=AStride[1])
00244 { Ap1 = AOffptr + i2;
00245 #ifdef __NO_BLAS__
00246 for(i1=0; i1 < ACount[0]*AStride[0]; i1+=AStride[0])
00247 {*(Ap1 + i1) += alpha;};
00248 #else
00249 daxpy_(&Size, &dOne, &alpha,&longZero, Ap1, &strideA);
00250 #endif
00251 };
00252 break;
00253
00254 case 3 :
00255 for(i3=0; i3 < ACount[2]*AStride[2]; i3+=AStride[2])
00256 { Ap2 = AOffptr + i3;
00257 for(i2=0; i2 < ACount[1]*AStride[1]; i2+=AStride[1])
00258 { Ap1 = Ap2 + i2;
00259 #ifdef __NO_BLAS__
00260 for(i1=0; i1 < ACount[0]*AStride[0]; i1+=AStride[0])
00261 {*(Ap1 + i1) += alpha;};
00262 #else
00263 daxpy_(&Size, &dOne, &alpha,&longZero, Ap1, &strideA);
00264 #endif
00265 }};
00266 break;
00267
00268 case 4 :
00269 for(i4=0; i4 < ACount[3]*AStride[3]; i4+=AStride[3])
00270 { Ap3 = AOffptr + i4;
00271 for(i3=0; i3 < ACount[2]*AStride[2]; i3+=AStride[2])
00272 { Ap2 = Ap3 + i3;
00273 for(i2=0; i2 < ACount[1]*AStride[1]; i2+=AStride[1])
00274 { Ap1 = Ap2 + i2;
00275 #ifdef __NO_BLAS__
00276 for(i1=0; i1 < ACount[0]*AStride[0]; i1+=AStride[0])
00277 {*(Ap1 + i1) += alpha;};
00278 #else
00279 daxpy_(&Size, &dOne, &alpha,&longZero, Ap1, &strideA);
00280 #endif
00281 }}};
00282 break;
00283
00284 case 5 :
00285 for(i5=0; i5 < ACount[4]*AStride[4]; i5+=AStride[4])
00286 { Ap4 = AOffptr + i5;
00287 for(i4=0; i4 < ACount[3]*AStride[3]; i4+=AStride[3])
00288 { Ap3 = Ap4 + i4;
00289 for(i3=0; i3 < ACount[2]*AStride[2]; i3+=AStride[2])
00290 { Ap2 = Ap3 + i3;
00291 for(i2=0; i2 < ACount[1]*AStride[1]; i2+=AStride[1])
00292 { Ap1 = Ap2 + i2;
00293 #ifdef __NO_BLAS__
00294 for(i1=0; i1 < ACount[0]*AStride[0]; i1+=AStride[0])
00295 {*(Ap1 + i1) += alpha;};
00296 #else
00297 daxpy_(&Size, &dOne, &alpha,&longZero, Ap1, &strideA);
00298 #endif
00299 }}}};
00300 break;
00301
00302 case 6 :
00303 for(i6=0; i6 < ACount[5]*AStride[5]; i6+=AStride[5])
00304 { Ap5 = AOffptr+ i6;
00305 for(i5=0; i5 < ACount[4]*AStride[4]; i5+=AStride[4])
00306 { Ap4 = Ap5 + i5;
00307 for(i4=0; i4 < ACount[3]*AStride[3]; i4+=AStride[3])
00308 { Ap3 = Ap4 + i4;
00309 for(i3=0; i3 < ACount[2]*AStride[2]; i3+=AStride[2])
00310 { Ap2 = Ap3 + i3;
00311 for(i2=0; i2 < ACount[1]*AStride[1]; i2+=AStride[1])
00312 { Ap1 = Ap2 + i2;
00313 #ifdef __NO_BLAS__
00314 for(i1=0; i1 < ACount[0]*AStride[0]; i1+=AStride[0])
00315 {*(Ap1 + i1) += alpha;};
00316 #else
00317 daxpy_(&Size, &dOne, &alpha,&longZero, Ap1, &strideA);
00318 #endif
00319 }}}}};
00320 break;
00321
00322 case 7 :
00323 for(i7=0; i7 < ACount[6]*AStride[6]; i7+=AStride[6])
00324 { Ap6 = AOffptr + i7;
00325 for(i6=0; i6 < ACount[5]*AStride[5]; i6+=AStride[5])
00326 { Ap5 = Ap6 + i6;
00327 for(i5=0; i5 < ACount[4]*AStride[4]; i5+=AStride[4])
00328 { Ap4 = Ap5 + i5;
00329 for(i4=0; i4 < ACount[3]*AStride[3]; i4+=AStride[3])
00330 { Ap3 = Ap4 + i4;
00331 for(i3=0; i3 < ACount[2]*AStride[2]; i3+=AStride[2])
00332 { Ap2 = Ap3 + i3;
00333 for(i2=0; i2 < ACount[1]*AStride[1]; i2+=AStride[1])
00334 { Ap1 = Ap2 + i2;
00335 #ifdef __NO_BLAS__
00336 for(i1=0; i1 < ACount[0]*AStride[0]; i1+=AStride[0])
00337 {*(Ap1 + i1 ) += alpha;};
00338 #else
00339 daxpy_(&Size, &dOne, &alpha,&longZero, Ap1, &strideA);
00340 #endif
00341 }}}}}};
00342 break;
00343
00344 }
00345 }
00346
00347 }
00348
00349 void CAMbinaryEngine::doubleAtimesEqualAlpha(CAMstructureBase &Astructure, CAMdataHandler &Adata,
00350 double alpha)
00351 {
00352 double* AdataP = (double*)Adata.getDataPointer();
00353
00354 if(Astructure.isSubset() == 0)
00355 {
00356 double* aptr;
00357 for(aptr = AdataP; aptr < AdataP + Adata.getDataSize(); aptr++)
00358 { *aptr *= alpha; };
00359 }
00360 else
00361 {
00362
00363
00364
00365 long AloopDimension;
00366 long AOffset;
00367 MVAlongBase ACount;
00368 MVAlongBase AStride;
00369
00370 Astructure.getNormalizedLoops(AloopDimension,AOffset, ACount, AStride);
00371
00372 register long i1, i2, i3, i4, i5, i6, i7;
00373
00374 double* AOffptr = AdataP + AOffset;
00375
00376 register double* Ap1;
00377 register double* Ap2;
00378 register double* Ap3;
00379 register double* Ap4;
00380 register double* Ap5;
00381 register double* Ap6;
00382
00383 switch(AloopDimension)
00384 {
00385 case 1 :
00386 for(i1=0; i1 < ACount[0]*AStride[0]; i1+=AStride[0])
00387 {*(AOffptr + i1) *= alpha;};
00388 break;
00389
00390 case 2 :
00391 for(i2=0; i2 < ACount[1]*AStride[1]; i2+=AStride[1])
00392 { Ap1 = AOffptr + i2;
00393 for(i1=0; i1 < ACount[0]*AStride[0]; i1+=AStride[0])
00394 {*(Ap1 + i1 ) *= alpha;};
00395 };
00396 break;
00397
00398 case 3 :
00399 for(i3=0; i3 < ACount[2]*AStride[2]; i3+=AStride[2])
00400 { Ap2 = AOffptr + i3;
00401 for(i2=0; i2 < ACount[1]*AStride[1]; i2+=AStride[1])
00402 { Ap1 = Ap2 + i2;
00403 for(i1=0; i1 < ACount[0]*AStride[0]; i1+=AStride[0])
00404 {*(Ap1 + i1 ) *= alpha;};
00405 }};
00406 break;
00407
00408 case 4 :
00409 for(i4=0; i4 < ACount[3]*AStride[3]; i4+=AStride[3])
00410 { Ap3 = AOffptr + i4;
00411 for(i3=0; i3 < ACount[2]*AStride[2]; i3+=AStride[2])
00412 { Ap2 = Ap3 + i3;
00413 for(i2=0; i2 < ACount[1]*AStride[1]; i2+=AStride[1])
00414 { Ap1 = Ap2 + i2;
00415 for(i1=0; i1 < ACount[0]*AStride[0]; i1+=AStride[0])
00416 {*(Ap1 + i1 ) *= alpha;};
00417 }}};
00418 break;
00419
00420 case 5 :
00421 for(i5=0; i5 < ACount[4]*AStride[4]; i5+=AStride[4])
00422 { Ap4 = AOffptr + i5;
00423 for(i4=0; i4 < ACount[3]*AStride[3]; i4+=AStride[3])
00424 { Ap3 = Ap4 + i4;
00425 for(i3=0; i3 < ACount[2]*AStride[2]; i3+=AStride[2])
00426 { Ap2 = Ap3 + i3;
00427 for(i2=0; i2 < ACount[1]*AStride[1]; i2+=AStride[1])
00428 { Ap1 = Ap2 + i2;
00429 for(i1=0; i1 < ACount[0]*AStride[0]; i1+=AStride[0])
00430 {*(Ap1 + i1 ) *= alpha;};
00431 }}}};
00432 break;
00433
00434 case 6 :
00435 for(i6=0; i6 < ACount[5]*AStride[5]; i6+=AStride[5])
00436 { Ap5 = AOffptr+ i6;
00437 for(i5=0; i5 < ACount[4]*AStride[4]; i5+=AStride[4])
00438 { Ap4 = Ap5 + i5;
00439 for(i4=0; i4 < ACount[3]*AStride[3]; i4+=AStride[3])
00440 { Ap3 = Ap4 + i4;
00441 for(i3=0; i3 < ACount[2]*AStride[2]; i3+=AStride[2])
00442 { Ap2 = Ap3 + i3;
00443 for(i2=0; i2 < ACount[1]*AStride[1]; i2+=AStride[1])
00444 { Ap1 = Ap2 + i2;
00445 for(i1=0; i1 < ACount[0]*AStride[0]; i1+=AStride[0])
00446 {*(Ap1 + i1 ) *= alpha;};
00447 }}}}};
00448 break;
00449
00450 case 7 :
00451 for(i7=0; i7 < ACount[6]*AStride[6]; i7+=AStride[6])
00452 { Ap6 = AOffptr + i7;
00453 for(i6=0; i6 < ACount[5]*AStride[5]; i6+=AStride[5])
00454 { Ap5 = Ap6 + i6;
00455 for(i5=0; i5 < ACount[4]*AStride[4]; i5+=AStride[4])
00456 { Ap4 = Ap5 + i5;
00457 for(i4=0; i4 < ACount[3]*AStride[3]; i4+=AStride[3])
00458 { Ap3 = Ap4 + i4;
00459 for(i3=0; i3 < ACount[2]*AStride[2]; i3+=AStride[2])
00460 { Ap2 = Ap3 + i3;
00461 for(i2=0; i2 < ACount[1]*AStride[1]; i2+=AStride[1])
00462 { Ap1 = Ap2 + i2;
00463 for(i1=0; i1 < ACount[0]*AStride[0]; i1+=AStride[0])
00464 {*(Ap1 + i1 ) *= alpha;};
00465 }}}}}};
00466 break;
00467
00468 }
00469 }
00470
00471 }
00472
00473 void CAMbinaryEngine::doubleAminusEqualAlpha(CAMstructureBase &Astructure, CAMdataHandler &Adata,
00474 double alpha)
00475 {CAMbinaryEngine::doubleAplusEqualAlpha(Astructure, Adata, -alpha);}
00476
00477 void CAMbinaryEngine::doubleAdivideEqualAlpha(CAMstructureBase &Astructure, CAMdataHandler &Adata,
00478 double alpha)
00479 {CAMbinaryEngine::doubleAtimesEqualAlpha(Astructure, Adata, 1.0/alpha);}
00480
00481
00482
00483
00484 void CAMbinaryEngine::doubleMaxValue(double* data, long n, double& maxValue)
00485 {
00486 long i;
00487 maxValue = *data;
00488 for(i = 1; i < n; i++)
00489 if(maxValue < *(data + i)) maxValue = *(data + i);
00490 }
00491 void CAMbinaryEngine::doubleMinValue(double* data, long n, double& minValue)
00492 {
00493 long i;
00494 minValue = *data;
00495 for(i = 1; i < n; i++)
00496 if(minValue > *(data + i)) minValue = *(data + i);
00497 }
00498 void CAMbinaryEngine::doubleMaxAbsValue(double* data, long n, double& maxAbsValue)
00499 {
00500 long i;
00501 maxAbsValue = fabs(*data);
00502 for(i = 1; i < n; i++)
00503 if(maxAbsValue < fabs(*(data + i))) maxAbsValue = fabs(*(data + i));
00504 }
00505 void CAMbinaryEngine::doubleMinAbsValue(double* data, long n, double& minAbsValue)
00506 {
00507 long i;
00508 minAbsValue = fabs(*data);
00509 for(i = 1; i < n; i++)
00510 if(minAbsValue > fabs(*(data + i))) minAbsValue = fabs(*(data + i));
00511 }
00512
00513 void CAMbinaryEngine::doublepNorm(double* data, long n, double p, double& pNormValue)
00514 {
00515 pNormValue = 0.0;
00516 double* dptr;
00517 for(dptr = (double*)data; dptr < (double*)data + n; dptr++)
00518 pNormValue += ::pow(fabs(*(dptr)),p);
00519
00520 pNormValue = ::pow(pNormValue,1.0/p);
00521 }
00522
