/**
* @file Matrix.cc
*
* Implementation of specialized Matrix template methods
*/
#include "Matrix.h"
#if defined(HAVE_CLAPACK)
#include "CLAPACK-3.1.1.1/INCLUDE/f2c.h"
extern "C" {
#include "CLAPACK-3.1.1.1/INCLUDE/clapack.h"
}
// These are some stupid clapack things that we want to get rid of
#ifdef min
#undef min
#endif
#ifdef max
#undef max
#endif
#endif
namespace TNet
{
//***************************************************************************
//***************************************************************************
#ifdef HAVE_ATLAS
//***************************************************************************
//***************************************************************************
template<>
Matrix<float> &
Matrix<float>::
Invert(float *LogDet, float *DetSign, bool inverse_needed)
{
assert(Rows() == Cols());
#if defined(HAVE_CLAPACK)
integer* pivot = new integer[mMRows];
integer M = Rows();
integer N = Cols();
integer LDA = mStride;
integer result;
integer l_work = std::max<integer>(1, N);
float* p_work = new float[l_work];
sgetrf_(&M, &N, mpData, &LDA, pivot, &result);
const int pivot_offset=1;
#else
int* pivot = new int[mMRows];
int result = clapack_sgetrf(CblasColMajor, Rows(), Cols(), mpData, mStride, pivot);
const int pivot_offset=0;
#endif
assert(result >= 0 && "Call to CLAPACK sgetrf_ or ATLAS clapack_sgetrf called with wrong arguments");
if(result != 0) {
Error("Matrix is singular");
}
if(DetSign!=NULL){ *DetSign=1.0; for(size_t i=0;i<mMRows;i++) if(pivot[i]!=(int)i+pivot_offset) *DetSign *= -1.0; }
if(LogDet!=NULL||DetSign!=NULL){ // Compute log determinant...
assert(mMRows==mMCols); // Can't take determinant of non-square matrix.
*LogDet = 0.0; float prod = 1.0;
for(size_t i=0;i<mMRows;i++){
prod *= (*this)(i,i);
if(i==mMRows-1 || fabs(prod)<1.0e-10 || fabs(prod)>1.0e+10){
if(LogDet!=NULL) *LogDet += log(fabs(prod));
if(DetSign!=NULL) *DetSign *= (prod>0?1.0:-1.0);
prod=1.0;
}
}
}
#if defined(HAVE_CLAPACK)
if(inverse_needed) sgetri_(&M, mpData, &LDA, pivot, p_work, &l_work, &result);
delete [] pivot;
#else
if(inverse_needed) result = clapack_sgetri(CblasColMajor, Rows(), mpData, mStride, pivot);
delete [] pivot;
#endif
assert(result == 0 && "Call to CLAPACK sgetri_ or ATLAS clapack_sgetri called with wrong arguments");
return *this;
}
//***************************************************************************
//***************************************************************************
template<>
Matrix<double> &
Matrix<double>::
Invert(double *LogDet, double *DetSign, bool inverse_needed)
{
assert(Rows() == Cols());
#if defined(HAVE_CLAPACK)
integer* pivot = new integer[mMRows];
integer M = Rows();
integer N = Cols();
integer LDA = mStride;
integer result;
integer l_work = std::max<integer>(1, N);
double* p_work = new double[l_work];
dgetrf_(&M, &N, mpData, &LDA, pivot, &result);
const int pivot_offset=1;
#else
int* pivot = new int[mMRows];
int result = clapack_dgetrf(CblasColMajor, Rows(), Cols(), mpData, mStride, pivot);
const int pivot_offset=0;
#endif
assert(result >= 0 && "Call to CLAPACK dgetrf_ or ATLAS clapack_dgetrf called with wrong arguments");
if(result != 0) {
Error("Matrix is singular");
}
if(DetSign!=NULL){ *DetSign=1.0; for(size_t i=0;i<mMRows;i++) if(pivot[i]!=(int)i+pivot_offset) *DetSign *= -1.0; }
if(LogDet!=NULL||DetSign!=NULL){ // Compute log determinant...
assert(mMRows==mMCols); // Can't take determinant of non-square matrix.
*LogDet = 0.0; double prod = 1.0;
for(size_t i=0;i<mMRows;i++){
prod *= (*this)(i,i);
if(i==mMRows-1 || fabs(prod)<1.0e-10 || fabs(prod)>1.0e+10){
if(LogDet!=NULL) *LogDet += log(fabs(prod));
if(DetSign!=NULL) *DetSign *= (prod>0?1.0:-1.0);
prod=1.0;
}
}
}
#if defined(HAVE_CLAPACK)
if(inverse_needed) dgetri_(&M, mpData, &LDA, pivot, p_work, &l_work, &result);
delete [] pivot;
#else
if(inverse_needed) result = clapack_dgetri(CblasColMajor, Rows(), mpData, mStride, pivot);
delete [] pivot;
#endif
assert(result == 0 && "Call to CLAPACK dgetri_ or ATLAS clapack_dgetri called with wrong arguments");
return *this;
}
template<>
Matrix<float> &
Matrix<float>::
BlasGer(const float alpha, const Vector<float>& rA, const Vector<float>& rB)
{
assert(rA.Dim() == mMRows && rB.Dim() == mMCols);
cblas_sger(CblasRowMajor, rA.Dim(), rB.Dim(), alpha, rA.pData(), 1, rB.pData(), 1, mpData, mStride);
return *this;
}
template<>
Matrix<double> &
Matrix<double>::
BlasGer(const double alpha, const Vector<double>& rA, const Vector<double>& rB)
{
assert(rA.Dim() == mMRows && rB.Dim() == mMCols);
cblas_dger(CblasRowMajor, rA.Dim(), rB.Dim(), alpha, rA.pData(), 1, rB.pData(), 1, mpData, mStride);
return *this;
}
template<>
Matrix<float>&
Matrix<float>::
BlasGemm(const float alpha,
const Matrix<float>& rA, MatrixTrasposeType transA,
const Matrix<float>& rB, MatrixTrasposeType transB,
const float beta)
{
assert((transA == NO_TRANS && transB == NO_TRANS && rA.Cols() == rB.Rows() && rA.Rows() == Rows() && rB.Cols() == Cols())
|| (transA == TRANS && transB == NO_TRANS && rA.Rows() == rB.Rows() && rA.Cols() == Rows() && rB.Cols() == Cols())
|| (transA == NO_TRANS && transB == TRANS && rA.Cols() == rB.Cols() && rA.Rows() == Rows() && rB.Rows() == Cols())
|| (transA == TRANS && transB == TRANS && rA.Rows() == rB.Cols() && rA.Cols() == Rows() && rB.Rows() == Cols()));
cblas_sgemm(CblasRowMajor, static_cast<CBLAS_TRANSPOSE>(transA), static_cast<CBLAS_TRANSPOSE>(transB),
Rows(), Cols(), transA == NO_TRANS ? rA.Cols() : rA.Rows(),
alpha, rA.mpData, rA.mStride, rB.mpData, rB.mStride,
beta, mpData, mStride);
return *this;
}
template<>
Matrix<double>&
Matrix<double>::
BlasGemm(const double alpha,
const Matrix<double>& rA, MatrixTrasposeType transA,
const Matrix<double>& rB, MatrixTrasposeType transB,
const double beta)
{
assert((transA == NO_TRANS && transB == NO_TRANS && rA.Cols() == rB.Rows() && rA.Rows() == Rows() && rB.Cols() == Cols())
|| (transA == TRANS && transB == NO_TRANS && rA.Rows() == rB.Rows() && rA.Cols() == Rows() && rB.Cols() == Cols())
|| (transA == NO_TRANS && transB == TRANS && rA.Cols() == rB.Cols() && rA.Rows() == Rows() && rB.Rows() == Cols())
|| (transA == TRANS && transB == TRANS && rA.Rows() == rB.Cols() && rA.Cols() == Rows() && rB.Rows() == Cols()));
cblas_dgemm(CblasRowMajor, static_cast<CBLAS_TRANSPOSE>(transA), static_cast<CBLAS_TRANSPOSE>(transB),
Rows(), Cols(), transA == NO_TRANS ? rA.Cols() : rA.Rows(),
alpha, rA.mpData, rA.mStride, rB.mpData, rB.mStride,
beta, mpData, mStride);
return *this;
}
template<>
Matrix<float>&
Matrix<float>::
Axpy(const float alpha,
const Matrix<float>& rA, MatrixTrasposeType transA){
int aStride = (int)rA.mStride, stride = mStride;
float *adata=rA.mpData, *data=mpData;
if(transA == NO_TRANS){
assert(rA.Rows()==Rows() && rA.Cols()==Cols());
for(size_t row=0;row<mMRows;row++,adata+=aStride,data+=stride)
cblas_saxpy(mMCols, alpha, adata, 1, data, 1);
} else {
assert(rA.Cols()==Rows() && rA.Rows()==Cols());
for(size_t row=0;row<mMRows;row++,adata++,data+=stride)
cblas_saxpy(mMCols, alpha, adata, aStride, data, 1);
}
return *this;
}
template<>
Matrix<double>&
Matrix<double>::
Axpy(const double alpha,
const Matrix<double>& rA, MatrixTrasposeType transA){
int aStride = (int)rA.mStride, stride = mStride;
double *adata=rA.mpData, *data=mpData;
if(transA == NO_TRANS){
assert(rA.Rows()==Rows() && rA.Cols()==Cols());
for(size_t row=0;row<mMRows;row++,adata+=aStride,data+=stride)
cblas_daxpy(mMCols, alpha, adata, 1, data, 1);
} else {
assert(rA.Cols()==Rows() && rA.Rows()==Cols());
for(size_t row=0;row<mMRows;row++,adata++,data+=stride)
cblas_daxpy(mMCols, alpha, adata, aStride, data, 1);
}
return *this;
}
template <> //non-member but friend!
double TraceOfProduct(const Matrix<double> &A, const Matrix<double> &B){ // tr(A B), equivalent to sum of each element of A times same element in B'
size_t aStride = A.mStride, bStride = B.mStride;
assert(A.Rows()==B.Cols() && A.Cols()==B.Rows());
double ans = 0.0;
double *adata=A.mpData, *bdata=B.mpData;
size_t arows=A.Rows(), acols=A.Cols();
for(size_t row=0;row<arows;row++,adata+=aStride,bdata++)
ans += cblas_ddot(acols, adata, 1, bdata, bStride);
return ans;
}
template <> //non-member but friend!
double TraceOfProductT(const Matrix<double> &A, const Matrix<double> &B){ // tr(A B), equivalent to sum of each element of A times same element in B'
size_t aStride = A.mStride, bStride = B.mStride;
assert(A.Rows()==B.Rows() && A.Cols()==B.Cols());
double ans = 0.0;
double *adata=A.mpData, *bdata=B.mpData;
size_t arows=A.Rows(), acols=A.Cols();
for(size_t row=0;row<arows;row++,adata+=aStride,bdata+=bStride)
ans += cblas_ddot(acols, adata, 1, bdata, 1);
return ans;
}
template <> //non-member but friend!
float TraceOfProduct(const Matrix<float> &A, const Matrix<float> &B){ // tr(A B), equivalent to sum of each element of A times same element in B'
size_t aStride = A.mStride, bStride = B.mStride;
assert(A.Rows()==B.Cols() && A.Cols()==B.Rows());
float ans = 0.0;
float *adata=A.mpData, *bdata=B.mpData;
size_t arows=A.Rows(), acols=A.Cols();
for(size_t row=0;row<arows;row++,adata+=aStride,bdata++)
ans += cblas_sdot(acols, adata, 1, bdata, bStride);
return ans;
}
template <> //non-member but friend!
float TraceOfProductT(const Matrix<float> &A, const Matrix<float> &B){ // tr(A B), equivalent to sum of each element of A times same element in B'
size_t aStride = A.mStride, bStride = B.mStride;
assert(A.Rows()==B.Rows() && A.Cols()==B.Cols());
float ans = 0.0;
float *adata=A.mpData, *bdata=B.mpData;
size_t arows=A.Rows(), acols=A.Cols();
for(size_t row=0;row<arows;row++,adata+=aStride,bdata+=bStride)
ans += cblas_sdot(acols, adata, 1, bdata, 1);
return ans;
}
#endif //HAVE_ATLAS
} //namespace STK