diff options
Diffstat (limited to 'thirdparty/bullet/LinearMath/btMatrixX.h')
-rw-r--r-- | thirdparty/bullet/LinearMath/btMatrixX.h | 554 |
1 files changed, 554 insertions, 0 deletions
diff --git a/thirdparty/bullet/LinearMath/btMatrixX.h b/thirdparty/bullet/LinearMath/btMatrixX.h new file mode 100644 index 0000000000..42caed42ef --- /dev/null +++ b/thirdparty/bullet/LinearMath/btMatrixX.h @@ -0,0 +1,554 @@ +/* +Bullet Continuous Collision Detection and Physics Library +Copyright (c) 2003-2013 Erwin Coumans http://bulletphysics.org + +This software is provided 'as-is', without any express or implied warranty. +In no event will the authors be held liable for any damages arising from the use of this software. +Permission is granted to anyone to use this software for any purpose, +including commercial applications, and to alter it and redistribute it freely, +subject to the following restrictions: + +1. The origin of this software must not be misrepresented; you must not claim that you wrote the original software. If you use this software in a product, an acknowledgment in the product documentation would be appreciated but is not required. +2. Altered source versions must be plainly marked as such, and must not be misrepresented as being the original software. +3. This notice may not be removed or altered from any source distribution. +*/ +///original version written by Erwin Coumans, October 2013 + +#ifndef BT_MATRIX_X_H +#define BT_MATRIX_X_H + +#include "LinearMath/btQuickprof.h" +#include "LinearMath/btAlignedObjectArray.h" +#include <stdio.h> + +//#define BT_DEBUG_OSTREAM +#ifdef BT_DEBUG_OSTREAM +#include <iostream> +#include <iomanip> // std::setw +#endif //BT_DEBUG_OSTREAM + +class btIntSortPredicate +{ + public: + bool operator() ( const int& a, const int& b ) const + { + return a < b; + } +}; + + +template <typename T> +struct btVectorX +{ + btAlignedObjectArray<T> m_storage; + + btVectorX() + { + } + btVectorX(int numRows) + { + m_storage.resize(numRows); + } + + void resize(int rows) + { + m_storage.resize(rows); + } + int cols() const + { + return 1; + } + int rows() const + { + return m_storage.size(); + } + int size() const + { + return rows(); + } + + T nrm2() const + { + T norm = T(0); + + int nn = rows(); + + { + if (nn == 1) + { + norm = btFabs((*this)[0]); + } + else + { + T scale = 0.0; + T ssq = 1.0; + + /* The following loop is equivalent to this call to the LAPACK + auxiliary routine: CALL SLASSQ( N, X, INCX, SCALE, SSQ ) */ + + for (int ix=0;ix<nn;ix++) + { + if ((*this)[ix] != 0.0) + { + T absxi = btFabs((*this)[ix]); + if (scale < absxi) + { + T temp; + temp = scale / absxi; + ssq = ssq * (temp * temp) + BT_ONE; + scale = absxi; + } + else + { + T temp; + temp = absxi / scale; + ssq += temp * temp; + } + } + } + norm = scale * sqrt(ssq); + } + } + return norm; + + } + void setZero() + { + if (m_storage.size()) + { + // for (int i=0;i<m_storage.size();i++) + // m_storage[i]=0; + //memset(&m_storage[0],0,sizeof(T)*m_storage.size()); + btSetZero(&m_storage[0],m_storage.size()); + } + } + const T& operator[] (int index) const + { + return m_storage[index]; + } + + T& operator[] (int index) + { + return m_storage[index]; + } + + T* getBufferPointerWritable() + { + return m_storage.size() ? &m_storage[0] : 0; + } + + const T* getBufferPointer() const + { + return m_storage.size() ? &m_storage[0] : 0; + } + +}; +/* + template <typename T> + void setElem(btMatrixX<T>& mat, int row, int col, T val) + { + mat.setElem(row,col,val); + } + */ + + +template <typename T> +struct btMatrixX +{ + int m_rows; + int m_cols; + int m_operations; + int m_resizeOperations; + int m_setElemOperations; + + btAlignedObjectArray<T> m_storage; + mutable btAlignedObjectArray< btAlignedObjectArray<int> > m_rowNonZeroElements1; + + T* getBufferPointerWritable() + { + return m_storage.size() ? &m_storage[0] : 0; + } + + const T* getBufferPointer() const + { + return m_storage.size() ? &m_storage[0] : 0; + } + btMatrixX() + :m_rows(0), + m_cols(0), + m_operations(0), + m_resizeOperations(0), + m_setElemOperations(0) + { + } + btMatrixX(int rows,int cols) + :m_rows(rows), + m_cols(cols), + m_operations(0), + m_resizeOperations(0), + m_setElemOperations(0) + { + resize(rows,cols); + } + void resize(int rows, int cols) + { + m_resizeOperations++; + m_rows = rows; + m_cols = cols; + { + BT_PROFILE("m_storage.resize"); + m_storage.resize(rows*cols); + } + } + int cols() const + { + return m_cols; + } + int rows() const + { + return m_rows; + } + ///we don't want this read/write operator(), because we cannot keep track of non-zero elements, use setElem instead + /*T& operator() (int row,int col) + { + return m_storage[col*m_rows+row]; + } + */ + + void addElem(int row,int col, T val) + { + if (val) + { + if (m_storage[col+row*m_cols]==0.f) + { + setElem(row,col,val); + } else + { + m_storage[row*m_cols+col] += val; + } + } + } + + + void setElem(int row,int col, T val) + { + m_setElemOperations++; + m_storage[row*m_cols+col] = val; + } + + void mulElem(int row,int col, T val) + { + m_setElemOperations++; + //mul doesn't change sparsity info + + m_storage[row*m_cols+col] *= val; + } + + + + + void copyLowerToUpperTriangle() + { + int count=0; + for (int row=0;row<rows();row++) + { + for (int col=0;col<row;col++) + { + setElem(col,row, (*this)(row,col)); + count++; + + } + } + //printf("copyLowerToUpperTriangle copied %d elements out of %dx%d=%d\n", count,rows(),cols(),cols()*rows()); + } + + const T& operator() (int row,int col) const + { + return m_storage[col+row*m_cols]; + } + + + void setZero() + { + { + BT_PROFILE("storage=0"); + btSetZero(&m_storage[0],m_storage.size()); + //memset(&m_storage[0],0,sizeof(T)*m_storage.size()); + //for (int i=0;i<m_storage.size();i++) + // m_storage[i]=0; + } + } + + void setIdentity() + { + btAssert(rows() == cols()); + + setZero(); + for (int row=0;row<rows();row++) + { + setElem(row,row,1); + } + } + + + + void printMatrix(const char* msg) + { + printf("%s ---------------------\n",msg); + for (int i=0;i<rows();i++) + { + printf("\n"); + for (int j=0;j<cols();j++) + { + printf("%2.1f\t",(*this)(i,j)); + } + } + printf("\n---------------------\n"); + + } + + + void rowComputeNonZeroElements() const + { + m_rowNonZeroElements1.resize(rows()); + for (int i=0;i<rows();i++) + { + m_rowNonZeroElements1[i].resize(0); + for (int j=0;j<cols();j++) + { + if ((*this)(i,j)!=0.f) + { + m_rowNonZeroElements1[i].push_back(j); + } + } + } + } + btMatrixX transpose() const + { + //transpose is optimized for sparse matrices + btMatrixX tr(m_cols,m_rows); + tr.setZero(); + for (int i=0;i<m_cols;i++) + for (int j=0;j<m_rows;j++) + { + T v = (*this)(j,i); + if (v) + { + tr.setElem(i,j,v); + } + } + return tr; + } + + + btMatrixX operator*(const btMatrixX& other) + { + //btMatrixX*btMatrixX implementation, brute force + btAssert(cols() == other.rows()); + + btMatrixX res(rows(),other.cols()); + res.setZero(); +// BT_PROFILE("btMatrixX mul"); + for (int j=0; j < res.cols(); ++j) + { + { + for (int i=0; i < res.rows(); ++i) + { + T dotProd=0; +// T dotProd2=0; + //int waste=0,waste2=0; + + { +// bool useOtherCol = true; + { + for (int v=0;v<rows();v++) + { + T w = (*this)(i,v); + if (other(v,j)!=0.f) + { + dotProd+=w*other(v,j); + } + + } + } + } + if (dotProd) + res.setElem(i,j,dotProd); + } + } + } + return res; + } + + // this assumes the 4th and 8th rows of B and C are zero. + void multiplyAdd2_p8r (const btScalar *B, const btScalar *C, int numRows, int numRowsOther ,int row, int col) + { + const btScalar *bb = B; + for ( int i = 0;i<numRows;i++) + { + const btScalar *cc = C; + for ( int j = 0;j<numRowsOther;j++) + { + btScalar sum; + sum = bb[0]*cc[0]; + sum += bb[1]*cc[1]; + sum += bb[2]*cc[2]; + sum += bb[4]*cc[4]; + sum += bb[5]*cc[5]; + sum += bb[6]*cc[6]; + addElem(row+i,col+j,sum); + cc += 8; + } + bb += 8; + } + } + + void multiply2_p8r (const btScalar *B, const btScalar *C, int numRows, int numRowsOther, int row, int col) + { + btAssert (numRows>0 && numRowsOther>0 && B && C); + const btScalar *bb = B; + for ( int i = 0;i<numRows;i++) + { + const btScalar *cc = C; + for ( int j = 0;j<numRowsOther;j++) + { + btScalar sum; + sum = bb[0]*cc[0]; + sum += bb[1]*cc[1]; + sum += bb[2]*cc[2]; + sum += bb[4]*cc[4]; + sum += bb[5]*cc[5]; + sum += bb[6]*cc[6]; + setElem(row+i,col+j,sum); + cc += 8; + } + bb += 8; + } + } + + void setSubMatrix(int rowstart,int colstart,int rowend,int colend,const T value) + { + int numRows = rowend+1-rowstart; + int numCols = colend+1-colstart; + + for (int row=0;row<numRows;row++) + { + for (int col=0;col<numCols;col++) + { + setElem(rowstart+row,colstart+col,value); + } + } + } + + void setSubMatrix(int rowstart,int colstart,int rowend,int colend,const btMatrixX& block) + { + btAssert(rowend+1-rowstart == block.rows()); + btAssert(colend+1-colstart == block.cols()); + for (int row=0;row<block.rows();row++) + { + for (int col=0;col<block.cols();col++) + { + setElem(rowstart+row,colstart+col,block(row,col)); + } + } + } + void setSubMatrix(int rowstart,int colstart,int rowend,int colend,const btVectorX<T>& block) + { + btAssert(rowend+1-rowstart == block.rows()); + btAssert(colend+1-colstart == block.cols()); + for (int row=0;row<block.rows();row++) + { + for (int col=0;col<block.cols();col++) + { + setElem(rowstart+row,colstart+col,block[row]); + } + } + } + + + btMatrixX negative() + { + btMatrixX neg(rows(),cols()); + for (int i=0;i<rows();i++) + for (int j=0;j<cols();j++) + { + T v = (*this)(i,j); + neg.setElem(i,j,-v); + } + return neg; + } + +}; + + + +typedef btMatrixX<float> btMatrixXf; +typedef btVectorX<float> btVectorXf; + +typedef btMatrixX<double> btMatrixXd; +typedef btVectorX<double> btVectorXd; + + +#ifdef BT_DEBUG_OSTREAM +template <typename T> +std::ostream& operator<< (std::ostream& os, const btMatrixX<T>& mat) + { + + os << " ["; + //printf("%s ---------------------\n",msg); + for (int i=0;i<mat.rows();i++) + { + for (int j=0;j<mat.cols();j++) + { + os << std::setw(12) << mat(i,j); + } + if (i!=mat.rows()-1) + os << std::endl << " "; + } + os << " ]"; + //printf("\n---------------------\n"); + + return os; + } +template <typename T> +std::ostream& operator<< (std::ostream& os, const btVectorX<T>& mat) + { + + os << " ["; + //printf("%s ---------------------\n",msg); + for (int i=0;i<mat.rows();i++) + { + os << std::setw(12) << mat[i]; + if (i!=mat.rows()-1) + os << std::endl << " "; + } + os << " ]"; + //printf("\n---------------------\n"); + + return os; + } + +#endif //BT_DEBUG_OSTREAM + + +inline void setElem(btMatrixXd& mat, int row, int col, double val) +{ + mat.setElem(row,col,val); +} + +inline void setElem(btMatrixXf& mat, int row, int col, float val) +{ + mat.setElem(row,col,val); +} + +#ifdef BT_USE_DOUBLE_PRECISION + #define btVectorXu btVectorXd + #define btMatrixXu btMatrixXd +#else + #define btVectorXu btVectorXf + #define btMatrixXu btMatrixXf +#endif //BT_USE_DOUBLE_PRECISION + + + +#endif//BT_MATRIX_H_H |