/* * 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 */