latmpl.h

Go to the documentation of this file.

/*-*-c++-*-****************************************************************
 *                     latmpl.h C++ Templates for lapackpp classes
                       -------------------
 begin                : 2005-12-29
 copyright            : (C) 2005 by Christian Stimming
 email                : stimming@tuhh.de
***************************************************************************/

// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 2, or (at
// your option) any later version.

// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public License along
// with this library; see the file COPYING.  If not, write to the Free
// Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307,
// USA.

#ifndef _LATMPL_H
#define _LATMPL_H

#include <cstdlib>
#include <cmath>
#include <lafnames.h>
#include LA_EXCEPTION_H

// Watch out: Only some compilers correctly implement the C++
// standard specifications, where type names from inside template
// classes need to be preceded by the "typename" keyword. Other
// compilers do not accept the keyword "typename". For those we
// define it as empty here.
#ifndef __GNUC__
   // Non-gcc compiler
#  ifdef _MSC_VER
     // This is Microsoft Visual C++.
     //  Microsoft Visual C++ 7.1  _MSC_VER = 1310
     //  Microsoft Visual C++ 7.0  _MSC_VER = 1300
     //  Microsoft Visual C++ 6.0  _MSC_VER = 1200
     //  Microsoft Visual C++ 5.0  _MSC_VER = 1100
#    if _MSC_VER <= 1300
       // This is Microsoft Visual C++ 7.0 or older, so define the
       // keyword as empty
#      define typename
#    else
       // Microsoft Visual C++ 7.1 or newer. typename should be ok.
#    endif
#  else
     // Compiler unknown
#  endif
#endif

namespace la {

template <class matT>
bool is_zero(const matT& mat)
{
  int i, j,  M=mat.rows(), N=mat.cols();

  // If your compiler gives an error in this line, please see the
  // note above on "typename".
  typename matT::value_type zero(0);

  for (j=0;j<N;j++)
    for (i=0;i<M; i++)
      if (mat(i, j) != zero)
        return false;
  return true;
}

template <class matT>
bool equal(const matT& mat1, const matT& mat2)
{
  int i, j,  M=mat1.rows(), N=mat1.cols();
  if (mat1.rows() != mat2.rows() || mat1.cols() != mat2.cols())
    throw LaException("equal<matT>(const matT&, const matT)",
                      "Both matrices have unequal sizes");
  for (j=0;j<N;j++)
    for (i=0;i<M; i++)
      if (mat1(i, j) != mat2(i, j))
        return false;
  return true;
}


template <class matT> 
void zeros(matT& mat, int N, int M=0)
{
  mat.resize(N, M == 0 ? N : M);
  mat = typename matT::value_type(0);
}

template <class matT> 
void ones(matT& mat, int N, int M=0)
{
  mat.resize(N, M == 0 ? N : M);
  mat = typename matT::value_type(1);
}

template <class matT> 
void eye(matT& mat, int N, int M=0)
{
  int nmin = (M == 0 ? N : (M < N ? M : N));
  mat.resize(N, M == 0 ? N : M);
  mat = typename matT::value_type(0);
  typename matT::value_type one(1);
  for (int k = 0; k < nmin; ++k)
    mat(k, k) = one;
}

template <class matT>
void rand(matT &A, typename matT::value_type low = 0,
          typename matT::value_type high = 1)
{
   int i, j,  M = A.rows(), N = A.cols();
   typename matT::value_type scale = high - low;
   for (j=0; j<N; ++j)
      for (i=0; i<M; ++i)
         A(i,j) = low + 
            scale * double(std::rand()) / double(RAND_MAX);
}

template <class matT>
void int_rand(matT &A, typename matT::value_type low = 0,
              typename matT::value_type high = 1)
{
   int i, j,  M = A.rows(), N = A.cols();
   double bins = high - low + 1;
   for (j=0; j<N; ++j)
      for (i=0; i<M; ++i)
         A(i,j) = low + 
            typename matT::value_type(
               std::floor(bins * double(std::rand()) / 
                          double(RAND_MAX)));
}

template <class matT>
void from_diag(matT& mat, const matT& vect)
{
  int nmin = (mat.rows() < mat.cols() ? mat.rows() : mat.cols());
  mat = typename matT::value_type(0);
  if (vect.rows() != 1 && vect.cols() != 1)
    throw LaException("diag<matT>(const matT& vect, matT& mat)",
                      "The argument 'vect' is not a vector: "
                      "neither dimension is equal to one");
  if (vect.rows() * vect.cols() != nmin)
    throw LaException("diag<matT>(const matT& vect, matT& mat)",
                      "The size of the vector is unequal to the matrix diagonal");
  if (vect.rows() == 1)
    for (int k = 0; k < nmin; ++k)
      mat(k, k) = vect(0, k);
  else
    for (int k = 0; k < nmin; ++k)
      mat(k, k) = vect(k, 0);
}



template <class matT>
matT zeros(int N, int M=0)
{
  matT mat;
  zeros(mat, N, M);
  return mat.shallow_assign();
}

template <class matT>
matT ones(int N, int M=0)
{
  matT mat;
  ones(mat, N, M);
  return mat.shallow_assign();
}

template <class matT>
matT eye(int N, int M=0)
{
  matT mat;
  eye(mat, N, M);
  return mat.shallow_assign();
}

template <class matT>
matT rand(int N, int M,
          typename matT::value_type low = 0,
          typename matT::value_type high = 1)
{
   matT mat(N, M);
   rand(mat, low, high);
   return mat.shallow_assign();
}

template <class matT>
matT int_rand(int N, int M,
          typename matT::value_type low = 0,
          typename matT::value_type high = 1)
{
   matT mat(N, M);
   int_rand(mat, low, high);
   return mat.shallow_assign();
}

template <class matT>
matT from_diag(const matT& vect)
{
  if (vect.rows() != 1 && vect.cols() != 1)
    throw LaException("diag<matT>(const matT& vect, matT& mat)",
                      "The argument 'vect' is not a vector: "
                      "neither dimension is equal to one");
  int nmax(vect.rows() > vect.cols() ? vect.rows() : vect.cols());
  matT mat(nmax, nmax);
  from_diag(mat, vect);
  return mat.shallow_assign();
}

template<class destT, class srcT>
destT convert_mat(const srcT& src)
{
   destT res(src.rows(), src.cols());
   // optimize later; for now use the correct but slow implementation
   int i, j,  M=src.rows(), N=src.cols();
   for (j=0; j<N; ++j)
      for (i=0; i<M; ++i)
         res(i, j) = typename destT::value_type ( src(i, j) );
   return res.shallow_assign();
}

template <class matT>
matT linspace(typename matT::value_type start, 
              typename matT::value_type end,
              int nr_points)
{
   matT mat(nr_points, 1);
   typename matT::value_type stepsize = 
      (end - start) / typename matT::value_type (nr_points - 1);
   for (int k = 0; k < nr_points; ++k)
   {
      mat(k, 0) = start;
      start += stepsize;
   }
   return mat.shallow_assign();
}

template <class matT>
matT repmat(const matT& A, int M, int N)
{
   int i, j, origM = A.rows(), origN = A.cols();
   matT mat(origM * M, origN * N);
   for (j = 0; j < N; ++j)
      for (i = 0; i < M; ++i)
         mat(LaIndex(i * origM, (i+1) * origM - 1),
             LaIndex(j * origN, (j+1) * origN - 1))
            .inject(A);

   return mat.shallow_assign();
}

template <class matT>
typename matT::value_type trace(const matT& mat)
{
  int nmin = (mat.rows() < mat.cols() ? mat.rows() : mat.cols());
  typename matT::value_type result(0);
  for (int k = 0; k < nmin; ++k)
    result += mat(k, k);
  return result;
}

template <class matT>
matT diag(const matT& mat)
{
  int nmin = (mat.rows() < mat.cols() ? mat.rows() : mat.cols());
  matT vect(nmin, 1);
  for (int k = 0; k < nmin; ++k)
      vect(k, 0) = mat(k, k);
  return vect.shallow_assign();
}


} // namespace la

#endif // _LATMPL_H

---