/*-*-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
/** @file
* @brief Template functions for matrices
*
* The definition of template functions for easier usage of lapackpp's
* matrix classes.
*/
#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 {
/** @name Matrix Predicates */
//@{
/** Returns true if this is an all-zero matrix. (New in lapackpp-2.4.4) */
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;
}
/** Returns true if both matrices are exactly equal. (New in
* lapackpp-2.4.4) */
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;
}
//@}
/** @name Create elementary matrices */
//@{
/** Sets the given matrix \c M to an all-zero matrix of dimension \c
* NxN, if \c M is not given, or \c NxM if \c M is given.
* (New in lapackpp-2.4.4) */
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);
}
/** Sets the given matrix \c M to an all-one matrix of dimension \c
* NxN, if \c M is not given, or \c NxM if \c M is given. (New in
* lapackpp-2.4.4) */
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);
}
/** Sets the given matrix \c M to an identity matrix (a diagonal of
* ones and all other elements zeros) of square dimension \c NxN, if
* \c M is not given, or of rectangular dimension \c NxM if \c M is
* given. (New in lapackpp-2.4.4) */
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;
}
/** Fills the given matrix \c A with pseudo-random values. The
* values are uniformly distributed in the interval \c (0,1) or,
* if specified, \c (low,high).
*
* Note: Since this uses the system's \c rand() call, the
* randomness of the values might be questionable -- don't use
* this if you need really strong random numbers. */
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);
}
/** Fills the given matrix \c A with pseudo-random values, where
* the value type of the matrix is an integer type. The values are
* uniformly distributed in the interval \c [0,1] or, if
* specified, \c (low,high), in both cases including the interval
* edges.
*
* Note: Since this uses the system's \c rand() call, the
* randomness of the values might be questionable -- don't use
* this if you need really strong random numbers. */
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)));
}
/** Sets the given matrix \c mat to a diagonal matrix with the
* given vector \c vect on the diagonal and zeros elsewhere. The
* matrix \c mat is allowed to be non-square, only the length of
* the diagonal has to fit to the vector's length. The vector \c
* vect is allowed to be either a row vector (dimension 1xN) or a
* column vector (dimension Nx1). (New in lapackpp-2.4.4) */
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);
}
//@}
/** @name Constructors for elementary and special matrices */
//@{
/** Returns a newly allocated all-zero matrix of dimension \c NxN, if
* \c M is not given, or \c NxM if \c M is given.
* (New in lapackpp-2.4.4) */
template <class matT>
matT zeros(int N, int M=0)
{
matT mat;
zeros(mat, N, M);
return mat.shallow_assign();
}
/** Returns a newly allocated all-one matrix of dimension \c NxN, if
* \c M is not given, or \c NxM if \c M is given.
* (New in lapackpp-2.4.4) */
template <class matT>
matT ones(int N, int M=0)
{
matT mat;
ones(mat, N, M);
return mat.shallow_assign();
}
/** Returns a newly allocated identity matrix of dimension \c NxN, if
* \c M is not given, or a rectangular matrix \c NxM if \c M is given.
* (New in lapackpp-2.4.4) */
template <class matT>
matT eye(int N, int M=0)
{
matT mat;
eye(mat, N, M);
return mat.shallow_assign();
}
/** Returns a newly allocated matrix of dimension \c NxM with
* pseudo-random values. The values are uniformly distributed in
* the interval \c (0,1) or, if specified, \c (low,high). (New in
* lapackpp-2.4.5)
*
* Note: Since this uses the system's \c rand() call, the
* randomness of the values might be questionable -- don't use
* this if you need really strong random numbers. */
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();
}
/** Returns a newly allocated matrix of dimension \c NxM with
* pseudo-random values, where the matrix element type is an
* integer type. The values are uniformly distributed in the
* interval \c [0,1] or, if specified, \c (low,high), in both
* cases including the interval edges. (New in lapackpp-2.4.5)
*
* Note: Since this uses the system's \c rand() call, the
* randomness of the values might be questionable -- don't use
* this if you need really strong random numbers. */
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();
}
/** Returns a newly allocated diagonal matrix of dimension \c NxN
* that has the vector \c vect of length \c N on the diagonal.
* The vector \c vect is allowed to be either a row vector
* (dimension 1xN) or a column vector (dimension Nx1). (New in
* lapackpp-2.4.5) */
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();
}
/** Returns a newly allocated matrix of type \c destT, containing
* the element-by-element converted values of the matrix \c src
* which was of type \c srcT.
*
* The template argument \c destT must be specified; the template
* argument \c srcT is deduced automatically. To convert a
* LaGenMatDouble into a LaGenMatInt:
\verbatim
LaGenMatDouble foo(5, 6);
LaGenMatInt bar(la::convert_to<LaGenMatInt>(foo));
\endverbatim
*
* Note: This conversion should even work from and to various
* matrices of the IT++ library, http://itpp.sourceforge.net (New
* in lapackpp-2.4.5.)
*/
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();
}
/** Returns a newly allocated linarly spaced column vector with \c
* nr_points elements, between and including \c start and \c
* end. (New in lapackpp-2.4.5.) */
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();
}
/** Returns a newly allocated large matrix that consists of \c
* M-by-N copies of the given matrix \c A. (New in
* lapackpp-2.4.5.) */
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();
}
//@}
/** @name Calculate some matrix measures */
//@{
/** Returns the trace, i.e. the sum of all diagonal elements of
* the matrix. The matrix \c mat does not have to be square. (New
* in lapackpp-2.4.4) */
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;
}
/** Returns a newly allocated column vector of dimension \c Nx1 that
* contains the diagonal of the given matrix \c mat. (New in
* lapackpp-2.4.5) */
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
syntax highlighted by Code2HTML, v. 0.9.1