| /* |
| * Copyright (C) 2011 The Android Open Source Project |
| * |
| * Licensed under the Apache License, Version 2.0 (the "License"); |
| * you may not use this file except in compliance with the License. |
| * You may obtain a copy of the License at |
| * |
| * http://www.apache.org/licenses/LICENSE-2.0 |
| * |
| * Unless required by applicable law or agreed to in writing, software |
| * distributed under the License is distributed on an "AS IS" BASIS, |
| * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| * See the License for the specific language governing permissions and |
| * limitations under the License. |
| */ |
| |
| /* $Id: db_utilities_linalg.h,v 1.5 2011/06/17 14:03:31 mbansal Exp $ */ |
| |
| #ifndef DB_UTILITIES_LINALG |
| #define DB_UTILITIES_LINALG |
| |
| #include "db_utilities.h" |
| |
| |
| |
| /***************************************************************** |
| * Lean and mean begins here * |
| *****************************************************************/ |
| /*! |
| * \defgroup LMLinAlg (LM) Linear Algebra Utilities (QR factorization, orthogonal basis, etc.) |
| */ |
| |
| /*! |
| \ingroup LMBasicUtilities |
| */ |
| inline void db_MultiplyScalar6(double A[6],double mult) |
| { |
| (*A++) *= mult; (*A++) *= mult; (*A++) *= mult; (*A++) *= mult; (*A++) *= mult; |
| (*A++) *= mult; |
| } |
| /*! |
| \ingroup LMBasicUtilities |
| */ |
| inline void db_MultiplyScalar7(double A[7],double mult) |
| { |
| (*A++) *= mult; (*A++) *= mult; (*A++) *= mult; (*A++) *= mult; (*A++) *= mult; |
| (*A++) *= mult; (*A++) *= mult; |
| } |
| /*! |
| \ingroup LMBasicUtilities |
| */ |
| inline void db_MultiplyScalar9(double A[9],double mult) |
| { |
| (*A++) *= mult; (*A++) *= mult; (*A++) *= mult; (*A++) *= mult; (*A++) *= mult; |
| (*A++) *= mult; (*A++) *= mult; (*A++) *= mult; (*A++) *= mult; |
| } |
| |
| /*! |
| \ingroup LMBasicUtilities |
| */ |
| inline double db_SquareSum6Stride7(const double *x) |
| { |
| return(db_sqr(x[0])+db_sqr(x[7])+db_sqr(x[14])+ |
| db_sqr(x[21])+db_sqr(x[28])+db_sqr(x[35])); |
| } |
| |
| /*! |
| \ingroup LMBasicUtilities |
| */ |
| inline double db_SquareSum8Stride9(const double *x) |
| { |
| return(db_sqr(x[0])+db_sqr(x[9])+db_sqr(x[18])+ |
| db_sqr(x[27])+db_sqr(x[36])+db_sqr(x[45])+ |
| db_sqr(x[54])+db_sqr(x[63])); |
| } |
| |
| /*! |
| \ingroup LMLinAlg |
| Cholesky-factorize symmetric positive definite 6 x 6 matrix A. Upper |
| part of A is used from the input. The Cholesky factor is output as |
| subdiagonal part in A and diagonal in d, which is 6-dimensional |
| 1.9 microseconds on 450MHz*/ |
| DB_API void db_CholeskyDecomp6x6(double A[36],double d[6]); |
| |
| /*! |
| \ingroup LMLinAlg |
| Backsubstitute L%transpose(L)*x=b for x given the Cholesky decomposition |
| of a 6 x 6 matrix and the right hand side b. The vector b is unchanged |
| 1.3 microseconds on 450MHz*/ |
| DB_API void db_CholeskyBacksub6x6(double x[6],const double A[36],const double d[6],const double b[6]); |
| |
| /*! |
| \ingroup LMLinAlg |
| Cholesky-factorize symmetric positive definite n x n matrix A.Part |
| above diagonal of A is used from the input, diagonal of A is assumed to |
| be stored in d. The Cholesky factor is output as |
| subdiagonal part in A and diagonal in d, which is n-dimensional*/ |
| DB_API void db_CholeskyDecompSeparateDiagonal(double **A,double *d,int n); |
| |
| /*! |
| \ingroup LMLinAlg |
| Backsubstitute L%transpose(L)*x=b for x given the Cholesky decomposition |
| of an n x n matrix and the right hand side b. The vector b is unchanged*/ |
| DB_API void db_CholeskyBacksub(double *x,const double * const *A,const double *d,int n,const double *b); |
| |
| /*! |
| \ingroup LMLinAlg |
| Cholesky-factorize symmetric positive definite 3 x 3 matrix A. Part |
| above diagonal of A is used from the input, diagonal of A is assumed to |
| be stored in d. The Cholesky factor is output as subdiagonal part in A |
| and diagonal in d, which is 3-dimensional*/ |
| DB_API void db_CholeskyDecomp3x3SeparateDiagonal(double A[9],double d[3]); |
| |
| /*! |
| \ingroup LMLinAlg |
| Backsubstitute L%transpose(L)*x=b for x given the Cholesky decomposition |
| of a 3 x 3 matrix and the right hand side b. The vector b is unchanged*/ |
| DB_API void db_CholeskyBacksub3x3(double x[3],const double A[9],const double d[3],const double b[3]); |
| |
| /*! |
| \ingroup LMLinAlg |
| perform A-=B*mult*/ |
| inline void db_RowOperation3(double A[3],const double B[3],double mult) |
| { |
| *A++ -= mult*(*B++); *A++ -= mult*(*B++); *A++ -= mult*(*B++); |
| } |
| |
| /*! |
| \ingroup LMLinAlg |
| */ |
| inline void db_RowOperation7(double A[7],const double B[7],double mult) |
| { |
| *A++ -= mult*(*B++); *A++ -= mult*(*B++); *A++ -= mult*(*B++); *A++ -= mult*(*B++); *A++ -= mult*(*B++); |
| *A++ -= mult*(*B++); *A++ -= mult*(*B++); |
| } |
| |
| /*! |
| \ingroup LMLinAlg |
| */ |
| inline void db_RowOperation9(double A[9],const double B[9],double mult) |
| { |
| *A++ -= mult*(*B++); *A++ -= mult*(*B++); *A++ -= mult*(*B++); *A++ -= mult*(*B++); *A++ -= mult*(*B++); |
| *A++ -= mult*(*B++); *A++ -= mult*(*B++); *A++ -= mult*(*B++); *A++ -= mult*(*B++); |
| } |
| |
| /*! |
| \ingroup LMBasicUtilities |
| Swap values of A[7] and B[7] |
| */ |
| inline void db_Swap7(double A[7],double B[7]) |
| { |
| double temp; |
| temp= *A; *A++ = *B; *B++ =temp; temp= *A; *A++ = *B; *B++ =temp; temp= *A; *A++ = *B; *B++ =temp; |
| temp= *A; *A++ = *B; *B++ =temp; temp= *A; *A++ = *B; *B++ =temp; temp= *A; *A++ = *B; *B++ =temp; |
| temp= *A; *A++ = *B; *B++ =temp; |
| } |
| |
| /*! |
| \ingroup LMBasicUtilities |
| Swap values of A[9] and B[9] |
| */ |
| inline void db_Swap9(double A[9],double B[9]) |
| { |
| double temp; |
| temp= *A; *A++ = *B; *B++ =temp; temp= *A; *A++ = *B; *B++ =temp; temp= *A; *A++ = *B; *B++ =temp; |
| temp= *A; *A++ = *B; *B++ =temp; temp= *A; *A++ = *B; *B++ =temp; temp= *A; *A++ = *B; *B++ =temp; |
| temp= *A; *A++ = *B; *B++ =temp; temp= *A; *A++ = *B; *B++ =temp; temp= *A; *A++ = *B; *B++ =temp; |
| } |
| |
| |
| /*! |
| \ingroup LMLinAlg |
| */ |
| DB_API void db_Orthogonalize6x7(double A[42],int orthonormalize=0); |
| |
| /*! |
| \ingroup LMLinAlg |
| */ |
| DB_API void db_Orthogonalize8x9(double A[72],int orthonormalize=0); |
| |
| /*! |
| \ingroup LMLinAlg |
| */ |
| inline double db_OrthogonalizePair7(double *x,const double *v,double ssv) |
| { |
| double m,sp,sp_m; |
| |
| m=db_SafeReciprocal(ssv); |
| sp=db_ScalarProduct7(x,v); |
| sp_m=sp*m; |
| db_RowOperation7(x,v,sp_m); |
| return(sp*sp_m); |
| } |
| |
| /*! |
| \ingroup LMLinAlg |
| */ |
| inline double db_OrthogonalizePair9(double *x,const double *v,double ssv) |
| { |
| double m,sp,sp_m; |
| |
| m=db_SafeReciprocal(ssv); |
| sp=db_ScalarProduct9(x,v); |
| sp_m=sp*m; |
| db_RowOperation9(x,v,sp_m); |
| return(sp*sp_m); |
| } |
| |
| /*! |
| \ingroup LMLinAlg |
| */ |
| inline void db_OrthogonalizationSwap7(double *A,int i,double *ss) |
| { |
| double temp; |
| |
| db_Swap7(A,A+7*i); |
| temp=ss[0]; ss[0]=ss[i]; ss[i]=temp; |
| } |
| /*! |
| \ingroup LMLinAlg |
| */ |
| inline void db_OrthogonalizationSwap9(double *A,int i,double *ss) |
| { |
| double temp; |
| |
| db_Swap9(A,A+9*i); |
| temp=ss[0]; ss[0]=ss[i]; ss[i]=temp; |
| } |
| |
| /*! |
| \ingroup LMLinAlg |
| */ |
| DB_API void db_NullVectorOrthonormal6x7(double x[7],const double A[42]); |
| /*! |
| \ingroup LMLinAlg |
| */ |
| DB_API void db_NullVectorOrthonormal8x9(double x[9],const double A[72]); |
| |
| /*! |
| \ingroup LMLinAlg |
| */ |
| inline void db_NullVector6x7Destructive(double x[7],double A[42]) |
| { |
| db_Orthogonalize6x7(A,1); |
| db_NullVectorOrthonormal6x7(x,A); |
| } |
| |
| /*! |
| \ingroup LMLinAlg |
| */ |
| inline void db_NullVector8x9Destructive(double x[9],double A[72]) |
| { |
| db_Orthogonalize8x9(A,1); |
| db_NullVectorOrthonormal8x9(x,A); |
| } |
| |
| inline int db_ScalarProduct512_s(const short *f,const short *g) |
| { |
| #ifndef DB_USE_MMX |
| int back; |
| back=0; |
| for(int i=1; i<=512; i++) |
| back+=(*f++)*(*g++); |
| |
| return(back); |
| #endif |
| } |
| |
| |
| inline int db_ScalarProduct32_s(const short *f,const short *g) |
| { |
| #ifndef DB_USE_MMX |
| int back; |
| back=0; |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| |
| return(back); |
| #endif |
| } |
| |
| /*! |
| \ingroup LMLinAlg |
| Scalar product of 128-vectors (short) |
| Compile-time control: MMX, SSE2 or standard C |
| */ |
| inline int db_ScalarProduct128_s(const short *f,const short *g) |
| { |
| #ifndef DB_USE_MMX |
| int back; |
| back=0; |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| |
| return(back); |
| #else |
| #ifdef DB_USE_SSE2 |
| int back; |
| |
| _asm |
| { |
| mov eax,f |
| mov ecx,g |
| /*First iteration************************************/ |
| movdqa xmm0,[eax] |
| pxor xmm7,xmm7 /*Set xmm7 to zero*/ |
| pmaddwd xmm0,[ecx] |
| /*Stall*/ |
| /*Standard iteration************************************/ |
| movdqa xmm2,[eax+16] |
| paddd xmm7,xmm0 |
| pmaddwd xmm2,[ecx+16] |
| /*Stall*/ |
| movdqa xmm1,[eax+32] |
| paddd xmm7,xmm2 |
| pmaddwd xmm1,[ecx+32] |
| /*Stall*/ |
| /*Standard iteration************************************/ |
| movdqa xmm0,[eax+48] |
| paddd xmm7,xmm1 |
| pmaddwd xmm0,[ecx+48] |
| /*Stall*/ |
| /*Standard iteration************************************/ |
| movdqa xmm2,[eax+64] |
| paddd xmm7,xmm0 |
| pmaddwd xmm2,[ecx+64] |
| /*Stall*/ |
| movdqa xmm1,[eax+80] |
| paddd xmm7,xmm2 |
| pmaddwd xmm1,[ecx+80] |
| /*Stall*/ |
| /*Standard iteration************************************/ |
| movdqa xmm0,[eax+96] |
| paddd xmm7,xmm1 |
| pmaddwd xmm0,[ecx+96] |
| /*Stall*/ |
| /*Standard iteration************************************/ |
| movdqa xmm2,[eax+112] |
| paddd xmm7,xmm0 |
| pmaddwd xmm2,[ecx+112] |
| /*Stall*/ |
| movdqa xmm1,[eax+128] |
| paddd xmm7,xmm2 |
| pmaddwd xmm1,[ecx+128] |
| /*Stall*/ |
| /*Standard iteration************************************/ |
| movdqa xmm0,[eax+144] |
| paddd xmm7,xmm1 |
| pmaddwd xmm0,[ecx+144] |
| /*Stall*/ |
| /*Standard iteration************************************/ |
| movdqa xmm2,[eax+160] |
| paddd xmm7,xmm0 |
| pmaddwd xmm2,[ecx+160] |
| /*Stall*/ |
| movdqa xmm1,[eax+176] |
| paddd xmm7,xmm2 |
| pmaddwd xmm1,[ecx+176] |
| /*Stall*/ |
| /*Standard iteration************************************/ |
| movdqa xmm0,[eax+192] |
| paddd xmm7,xmm1 |
| pmaddwd xmm0,[ecx+192] |
| /*Stall*/ |
| /*Standard iteration************************************/ |
| movdqa xmm2,[eax+208] |
| paddd xmm7,xmm0 |
| pmaddwd xmm2,[ecx+208] |
| /*Stall*/ |
| movdqa xmm1,[eax+224] |
| paddd xmm7,xmm2 |
| pmaddwd xmm1,[ecx+224] |
| /*Stall*/ |
| /*Standard iteration************************************/ |
| movdqa xmm0,[eax+240] |
| paddd xmm7,xmm1 |
| pmaddwd xmm0,[ecx+240] |
| /*Stall*/ |
| /*Rest iteration************************************/ |
| paddd xmm7,xmm0 |
| |
| /* add up the sum squares */ |
| movhlps xmm0,xmm7 /* high half to low half */ |
| paddd xmm7,xmm0 /* add high to low */ |
| pshuflw xmm0,xmm7, 0xE /* reshuffle */ |
| paddd xmm7,xmm0 /* add remaining */ |
| movd back,xmm7 |
| |
| emms |
| } |
| |
| return(back); |
| #else |
| int back; |
| |
| _asm |
| { |
| mov eax,f |
| mov ecx,g |
| /*First iteration************************************/ |
| movq mm0,[eax] |
| pxor mm7,mm7 /*Set mm7 to zero*/ |
| pmaddwd mm0,[ecx] |
| /*Stall*/ |
| movq mm1,[eax+8] |
| /*Stall*/ |
| pmaddwd mm1,[ecx+8] |
| /*Stall*/ |
| /*Standard iteration************************************/ |
| movq mm2,[eax+16] |
| paddd mm7,mm0 |
| pmaddwd mm2,[ecx+16] |
| /*Stall*/ |
| movq mm0,[eax+24] |
| paddd mm7,mm1 |
| pmaddwd mm0,[ecx+24] |
| /*Stall*/ |
| movq mm1,[eax+32] |
| paddd mm7,mm2 |
| pmaddwd mm1,[ecx+32] |
| /*Stall*/ |
| /*Standard iteration************************************/ |
| movq mm2,[eax+40] |
| paddd mm7,mm0 |
| pmaddwd mm2,[ecx+40] |
| /*Stall*/ |
| movq mm0,[eax+48] |
| paddd mm7,mm1 |
| pmaddwd mm0,[ecx+48] |
| /*Stall*/ |
| movq mm1,[eax+56] |
| paddd mm7,mm2 |
| pmaddwd mm1,[ecx+56] |
| /*Stall*/ |
| /*Standard iteration************************************/ |
| movq mm2,[eax+64] |
| paddd mm7,mm0 |
| pmaddwd mm2,[ecx+64] |
| /*Stall*/ |
| movq mm0,[eax+72] |
| paddd mm7,mm1 |
| pmaddwd mm0,[ecx+72] |
| /*Stall*/ |
| movq mm1,[eax+80] |
| paddd mm7,mm2 |
| pmaddwd mm1,[ecx+80] |
| /*Stall*/ |
| /*Standard iteration************************************/ |
| movq mm2,[eax+88] |
| paddd mm7,mm0 |
| pmaddwd mm2,[ecx+88] |
| /*Stall*/ |
| movq mm0,[eax+96] |
| paddd mm7,mm1 |
| pmaddwd mm0,[ecx+96] |
| /*Stall*/ |
| movq mm1,[eax+104] |
| paddd mm7,mm2 |
| pmaddwd mm1,[ecx+104] |
| /*Stall*/ |
| /*Standard iteration************************************/ |
| movq mm2,[eax+112] |
| paddd mm7,mm0 |
| pmaddwd mm2,[ecx+112] |
| /*Stall*/ |
| movq mm0,[eax+120] |
| paddd mm7,mm1 |
| pmaddwd mm0,[ecx+120] |
| /*Stall*/ |
| movq mm1,[eax+128] |
| paddd mm7,mm2 |
| pmaddwd mm1,[ecx+128] |
| /*Stall*/ |
| /*Standard iteration************************************/ |
| movq mm2,[eax+136] |
| paddd mm7,mm0 |
| pmaddwd mm2,[ecx+136] |
| /*Stall*/ |
| movq mm0,[eax+144] |
| paddd mm7,mm1 |
| pmaddwd mm0,[ecx+144] |
| /*Stall*/ |
| movq mm1,[eax+152] |
| paddd mm7,mm2 |
| pmaddwd mm1,[ecx+152] |
| /*Stall*/ |
| /*Standard iteration************************************/ |
| movq mm2,[eax+160] |
| paddd mm7,mm0 |
| pmaddwd mm2,[ecx+160] |
| /*Stall*/ |
| movq mm0,[eax+168] |
| paddd mm7,mm1 |
| pmaddwd mm0,[ecx+168] |
| /*Stall*/ |
| movq mm1,[eax+176] |
| paddd mm7,mm2 |
| pmaddwd mm1,[ecx+176] |
| /*Stall*/ |
| /*Standard iteration************************************/ |
| movq mm2,[eax+184] |
| paddd mm7,mm0 |
| pmaddwd mm2,[ecx+184] |
| /*Stall*/ |
| movq mm0,[eax+192] |
| paddd mm7,mm1 |
| pmaddwd mm0,[ecx+192] |
| /*Stall*/ |
| movq mm1,[eax+200] |
| paddd mm7,mm2 |
| pmaddwd mm1,[ecx+200] |
| /*Stall*/ |
| /*Standard iteration************************************/ |
| movq mm2,[eax+208] |
| paddd mm7,mm0 |
| pmaddwd mm2,[ecx+208] |
| /*Stall*/ |
| movq mm0,[eax+216] |
| paddd mm7,mm1 |
| pmaddwd mm0,[ecx+216] |
| /*Stall*/ |
| movq mm1,[eax+224] |
| paddd mm7,mm2 |
| pmaddwd mm1,[ecx+224] |
| /*Stall*/ |
| /*Standard iteration************************************/ |
| movq mm2,[eax+232] |
| paddd mm7,mm0 |
| pmaddwd mm2,[ecx+232] |
| /*Stall*/ |
| movq mm0,[eax+240] |
| paddd mm7,mm1 |
| pmaddwd mm0,[ecx+240] |
| /*Stall*/ |
| movq mm1,[eax+248] |
| paddd mm7,mm2 |
| pmaddwd mm1,[ecx+248] |
| /*Stall*/ |
| /*Rest iteration************************************/ |
| paddd mm7,mm0 |
| /*Stall*/ |
| /*Stall*/ |
| /*Stall*/ |
| paddd mm7,mm1 |
| /*Stall*/ |
| movq mm0,mm7 |
| psrlq mm7,32 |
| paddd mm0,mm7 |
| /*Stall*/ |
| /*Stall*/ |
| /*Stall*/ |
| movd back,mm0 |
| emms |
| } |
| |
| return(back); |
| #endif |
| #endif /*DB_USE_MMX*/ |
| } |
| |
| /*! |
| \ingroup LMLinAlg |
| Scalar product of 16 byte aligned 128-vectors (float). |
| Compile-time control: SIMD (SSE) or standard C. |
| */ |
| inline float db_ScalarProduct128Aligned16_f(const float *f,const float *g) |
| { |
| #ifndef DB_USE_SIMD |
| float back; |
| back=0.0; |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| back+=(*f++)*(*g++); back+=(*f++)*(*g++); back+=(*f++)*(*g++); |
| |
| return(back); |
| #else |
| float back; |
| |
| _asm |
| { |
| mov eax,f |
| mov ecx,g |
| /*First iteration************************************/ |
| movaps xmm0,[eax] |
| xorps xmm7,xmm7 /*Set mm7 to zero*/ |
| mulps xmm0,[ecx] |
| /*Stall*/ |
| movaps xmm1,[eax+16] |
| /*Stall*/ |
| mulps xmm1,[ecx+16] |
| /*Stall*/ |
| /*Standard iteration************************************/ |
| movaps xmm2,[eax+32] |
| addps xmm7,xmm0 |
| mulps xmm2,[ecx+32] |
| /*Stall*/ |
| movaps xmm0,[eax+48] |
| addps xmm7,xmm1 |
| mulps xmm0,[ecx+48] |
| /*Stall*/ |
| movaps xmm1,[eax+64] |
| addps xmm7,xmm2 |
| mulps xmm1,[ecx+64] |
| /*Stall*/ |
| /*Standard iteration************************************/ |
| movaps xmm2,[eax+80] |
| addps xmm7,xmm0 |
| mulps xmm2,[ecx+80] |
| /*Stall*/ |
| movaps xmm0,[eax+96] |
| addps xmm7,xmm1 |
| mulps xmm0,[ecx+96] |
| /*Stall*/ |
| movaps xmm1,[eax+112] |
| addps xmm7,xmm2 |
| mulps xmm1,[ecx+112] |
| /*Stall*/ |
| /*Standard iteration************************************/ |
| movaps xmm2,[eax+128] |
| addps xmm7,xmm0 |
| mulps xmm2,[ecx+128] |
| /*Stall*/ |
| movaps xmm0,[eax+144] |
| addps xmm7,xmm1 |
| mulps xmm0,[ecx+144] |
| /*Stall*/ |
| movaps xmm1,[eax+160] |
| addps xmm7,xmm2 |
| mulps xmm1,[ecx+160] |
| /*Stall*/ |
| /*Standard iteration************************************/ |
| movaps xmm2,[eax+176] |
| addps xmm7,xmm0 |
| mulps xmm2,[ecx+176] |
| /*Stall*/ |
| movaps xmm0,[eax+192] |
| addps xmm7,xmm1 |
| mulps xmm0,[ecx+192] |
| /*Stall*/ |
| movaps xmm1,[eax+208] |
| addps xmm7,xmm2 |
| mulps xmm1,[ecx+208] |
| /*Stall*/ |
| /*Standard iteration************************************/ |
| movaps xmm2,[eax+224] |
| addps xmm7,xmm0 |
| mulps xmm2,[ecx+224] |
| /*Stall*/ |
| movaps xmm0,[eax+240] |
| addps xmm7,xmm1 |
| mulps xmm0,[ecx+240] |
| /*Stall*/ |
| movaps xmm1,[eax+256] |
| addps xmm7,xmm2 |
| mulps xmm1,[ecx+256] |
| /*Stall*/ |
| /*Standard iteration************************************/ |
| movaps xmm2,[eax+272] |
| addps xmm7,xmm0 |
| mulps xmm2,[ecx+272] |
| /*Stall*/ |
| movaps xmm0,[eax+288] |
| addps xmm7,xmm1 |
| mulps xmm0,[ecx+288] |
| /*Stall*/ |
| movaps xmm1,[eax+304] |
| addps xmm7,xmm2 |
| mulps xmm1,[ecx+304] |
| /*Stall*/ |
| /*Standard iteration************************************/ |
| movaps xmm2,[eax+320] |
| addps xmm7,xmm0 |
| mulps xmm2,[ecx+320] |
| /*Stall*/ |
| movaps xmm0,[eax+336] |
| addps xmm7,xmm1 |
| mulps xmm0,[ecx+336] |
| /*Stall*/ |
| movaps xmm1,[eax+352] |
| addps xmm7,xmm2 |
| mulps xmm1,[ecx+352] |
| /*Stall*/ |
| /*Standard iteration************************************/ |
| movaps xmm2,[eax+368] |
| addps xmm7,xmm0 |
| mulps xmm2,[ecx+368] |
| /*Stall*/ |
| movaps xmm0,[eax+384] |
| addps xmm7,xmm1 |
| mulps xmm0,[ecx+384] |
| /*Stall*/ |
| movaps xmm1,[eax+400] |
| addps xmm7,xmm2 |
| mulps xmm1,[ecx+400] |
| /*Stall*/ |
| /*Standard iteration************************************/ |
| movaps xmm2,[eax+416] |
| addps xmm7,xmm0 |
| mulps xmm2,[ecx+416] |
| /*Stall*/ |
| movaps xmm0,[eax+432] |
| addps xmm7,xmm1 |
| mulps xmm0,[ecx+432] |
| /*Stall*/ |
| movaps xmm1,[eax+448] |
| addps xmm7,xmm2 |
| mulps xmm1,[ecx+448] |
| /*Stall*/ |
| /*Standard iteration************************************/ |
| movaps xmm2,[eax+464] |
| addps xmm7,xmm0 |
| mulps xmm2,[ecx+464] |
| /*Stall*/ |
| movaps xmm0,[eax+480] |
| addps xmm7,xmm1 |
| mulps xmm0,[ecx+480] |
| /*Stall*/ |
| movaps xmm1,[eax+496] |
| addps xmm7,xmm2 |
| mulps xmm1,[ecx+496] |
| /*Stall*/ |
| /*Rest iteration************************************/ |
| addps xmm7,xmm0 |
| /*Stall*/ |
| addps xmm7,xmm1 |
| /*Stall*/ |
| movaps xmm6,xmm7 |
| /*Stall*/ |
| shufps xmm6,xmm6,4Eh |
| /*Stall*/ |
| addps xmm7,xmm6 |
| /*Stall*/ |
| movaps xmm6,xmm7 |
| /*Stall*/ |
| shufps xmm6,xmm6,11h |
| /*Stall*/ |
| addps xmm7,xmm6 |
| /*Stall*/ |
| movss back,xmm7 |
| } |
| |
| return(back); |
| #endif /*DB_USE_SIMD*/ |
| } |
| |
| #endif /* DB_UTILITIES_LINALG */ |
