utlVNLBlas.h

Go to the documentation of this file.

#ifndef __utlVNLBlas_h
#define __utlVNLBlas_h

#include "utlCore.h"
#include "utlVNL.h"
#include "utlBlas.h"

#ifdef UTL_USE_MKL
#include "utlMKL.h"
#endif

namespace utl 
{

template <class T> inline void 
MatrixCopy(const vnl_matrix<T>& mat, vnl_matrix<T>& matOut, const T alpha, const char trans='N')
{
#ifdef UTL_USE_MKL
  if (trans=='N' || trans=='n' || trans=='R' || trans=='r')
    matOut.set_size(mat.rows(), mat.cols());
  else if (trans=='T' || trans=='t' || trans=='C' || trans=='c')
    matOut.set_size(mat.cols(), mat.rows());
  else
    utlException(true, "wrong trans");
  if ((trans=='N' || trans=='n') && std::fabs(alpha-1.0)<1e-10 )
    utl::cblas_copy(mat.rows()*mat.cols(), mat.data_block(), 1, matOut.data_block(), 1);
  else
    utl::mkl_omatcopy<T>('R', trans, mat.rows(), mat.cols(), alpha, mat.data_block(), mat.cols(), matOut.data_block(), matOut.cols()); 
#else
  if (trans=='N' || trans=='n')
    matOut = mat;
  else if (trans=='T' || trans=='t')
    matOut = mat.transpose();
  else if (trans=='C' || trans=='c')
    matOut = mat.conjugate_transpose();
  else if (trans=='R' || trans=='r')
    {
    matOut = mat;
    vnl_c_vector<T>::conjugate(matOut.begin(),matOut.begin(),matOut.size());  
    }

  if (std::fabs(alpha-1)>1e-10)
    matOut *= alpha;
#endif
}


__utl_gemm_MatrixTimesMatrix(T, gemm_VnlMatrixTimesMatrix, Vnl, vnl_matrix<T>, rows, cols, data_block, set_size);


__utl_gemv_MatrixTimesVector(T, gemv_VnlMatrixTimesVector, Vnl, vnl_matrix<T>, rows, cols, data_block, vnl_vector<T>, size, data_block, set_size);


__utl_gevm_MatrixTimesVector(T, gemm_VnlVectorTimesMatrix, Vnl, vnl_matrix<T>, rows, cols, data_block, vnl_vector<T>, size, data_block, set_size);

__utl_syrk_Matrix(T, syrk_VnlMatrix, Vnl, vnl_matrix<T>, rows, cols, data_block, set_size); 

template <class T> 
inline T 
InnerProduct(const vnl_vector<T>& v1, const vnl_vector<T>& v2)
{
  utlSAException(v1.size() != v2.size())(v1.size())(v2.size()).msg("vector sizes mismatch");
  return utl::cblas_dot<T>(v1.size(), v1.data_block(), 1, v2.data_block(), 1);
}

template <class T> 
inline void 
OuterProduct(const vnl_vector<T>& v1, const vnl_vector<T>& v2, vnl_matrix<T>& mat, const double alpha=1.0)
{
  int M = v1.size(), N = v2.size();
  if (M!=mat.rows() || N!=mat.cols())
    {
    mat.set_size(M, N);
    mat.fill(0.0);
    }
  utl::cblas_ger<T>(CblasRowMajor, M, N, alpha, v1.data_block(), 1, v2.data_block(), 1, mat.data_block(), mat.cols());
}

template <class T> 
inline void 
OuterProduct(const vnl_vector<T>& v1, vnl_matrix<T>& mat, const double alpha=1.0)
{
  int M = v1.size();
  if (M!=mat.rows() || M!=mat.cols())
    {
    mat.set_size(M, M);
    mat.fill(0.0);
    }
  utl::cblas_syr<T>(CblasRowMajor, CblasUpper, M, alpha, v1.data_block(), 1, mat.data_block(), mat.cols());
  T* data = mat.data_block();
  for ( int i = 0; i < M; ++i ) 
    for ( int j = 0; j < i; ++j ) 
      data[i*M+j] = data[j*M+i];
}

template <class T> 
inline void 
GetRow(const vnl_matrix<T>& mat, const int index, vnl_vector<T>& v1)
{
  if (v1.size()!=mat.cols())
    v1.set_size(mat.cols());
  utl::cblas_copy<T>(v1.size(), mat.data_block()+mat.cols()*index, 1, v1.data_block(), 1);
}

template <class T> 
inline void 
GetColumn(const vnl_matrix<T>& mat, const int index, vnl_vector<T>& v1)
{
  if (v1.size()!=mat.rows())
    v1.set_size(mat.rows());
  utl::cblas_copy<T>(v1.size(), mat.data_block()+index, mat.cols(), v1.data_block(), 1);
}


}


#endif