00001 #ifndef Impala_Core_Matrix_MatNorm2Dist_h
00002 #define Impala_Core_Matrix_MatNorm2Dist_h
00003
00004 #include "Core/Matrix/MatFunc.h"
00005 #include "Core/Matrix/MatSum.h"
00006 #include "Core/Matrix/MatMul.h"
00007 #include "Core/Array/Abs.h"
00008 #include "Core/Array/Add.h"
00009 #include "Core/Array/Mul.h"
00010 #include "Core/Array/MulVal.h"
00011 #include "Core/Array/Sqrt.h"
00012 #include "Basis/Timer.h"
00013 #ifdef CUDA
00014 #include "Link/Cuda/Cuda.h"
00015 #include "Link/Cuda/MatNorm2Dist.h"
00016 #endif
00017 #ifdef SSE_USED
00018 #include "Core/Matrix/MatNorm2DistSSE.h"
00019 #endif
00020
00021
00022 namespace Impala
00023 {
00024 namespace Core
00025 {
00026 namespace Matrix
00027 {
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081 template<class ArrayT>
00082 inline ArrayT*
00083 MatNorm2DistInternal(ArrayT* aT, ArrayT* b)
00084 {
00085 #ifdef SSE_USED
00086 return MatNorm2DistSSE(aT, b);
00087 #endif
00088 ILOG_VAR(Core.Matrix.MatNorm2Dist);
00089 if (MatNrCol(aT) != MatNrRow(b)) {
00090 ILOG_ERROR("MatNorm2DistInternal operands: dimensionality problem");
00091 }
00092
00093 Timer timer;
00094
00095
00096 ArrayT* tempAT = 0;
00097 Mul(tempAT, aT, aT);
00098 ArrayT* aaT = MatSumAxis1(tempAT);
00099 delete tempAT;
00100
00101 ILOG_DEBUG("aaT: " << timer.SplitTimeStr());
00102
00103
00104 ArrayT* tempB = 0;
00105 Mul(tempB, b, b);
00106 ArrayT* bb = MatSumAxis0(tempB);
00107 delete tempB;
00108
00109 ILOG_DEBUG("bb: " << timer.SplitTimeStr());
00110
00111
00112
00113
00114
00115 ILOG_DEBUG("aT " << MatNrRow(aT) << " " << MatNrCol(aT));
00116 ILOG_DEBUG("b " << MatNrRow(b) << " " << MatNrCol(b));
00117
00118 ArrayT* ab = MatMul(aT, b);
00119
00120 ILOG_DEBUG("matmul: " << timer.SplitTimeStr());
00121
00122
00123
00124
00125
00126
00127
00128 #pragma omp parallel for
00129 for (int i=0 ; i<MatNrRow(aaT) ; i++)
00130 {
00131 Real64 tmp = *MatE(aaT, i, 0);
00132 for (int j=0 ; j<MatNrCol(bb) ; j++)
00133 {
00134 *MatE(ab, i, j) = sqrt(fabs(*MatE(bb, 0, j) + tmp - 2.0 * *MatE(ab, i, j)));
00135 }
00136 }
00137 ILOG_DEBUG("sqrt(abs(aa + bb - 2ab)): " << timer.SplitTimeStr());
00138 delete bb;
00139 delete aaT;
00140
00141 ILOG_DEBUG("delete: " << timer.SplitTimeStr());
00142
00143
00144
00145
00146 ILOG_DEBUG(timer.SplitTime() << " (cpu-matnorm2dist-total)");
00147 return ab;
00148 }
00149
00150
00151 template<class ArrayT>
00152 inline ArrayT*
00153 MatNorm2Dist(ArrayT* a, ArrayT* b)
00154 {
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166 ILOG_VAR(Core.Matrix.MatNorm2Dist);
00167 if (MatNrRow(a) != MatNrRow(b)) {
00168 ILOG_ERROR("MatNorm2Dist operands: A and B should be of same dimensionality");
00169 }
00170 Mat* aT = MatTranspose(a);
00171 ArrayT* result = MatNorm2DistInternal(aT, b);
00172 delete aT;
00173 return result;
00174 }
00175
00176
00177 template<class ArrayT>
00178 inline ArrayT*
00179 MatNorm2DistTransposed(ArrayT* aT, ArrayT* bT)
00180 {
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192 ILOG_VAR(Core.Matrix.MatNorm2DistTransposed);
00193 if (MatNrCol(aT) != MatNrCol(bT)) {
00194 ILOG_ERROR("MatNorm2DistTransposed operands: A and B should be of same dimensionality");
00195 }
00196 #ifdef CUDA
00197 if(Link::Cuda::CudaUsed())
00198 {
00199 return Link::Cuda::MatNorm2DistTransposed(aT, bT);
00200 }
00201 #endif
00202 Mat* b = MatTranspose(bT);
00203 Mat* result = MatNorm2DistInternal(aT, b);
00204 delete b;
00205 return result;
00206 }
00207
00208 }
00209 }
00210 }
00211
00212 #endif