RE: MTL on VS2003 (_Ptrit Problems)
"Laws, Aaron D" <[email protected]> Fri, 11 Jun 2004 09:24:22 -0500
| Newsgroups | gmane.comp.lib.mtl.devel |
|---|---|
| Message-ID | <A9D1D3D3ED1BF24D9BF22ADDF0EFAE48053CBD96@xch-mw-11.mw.nos.boeing.com> |
I made the attached hacks to get it to compile under .NET2003, in lieu of a real fix, mostly by following the compile errors and commenting a lot of the workarounds put in for .NET2002. I haven't done much testing to verify results, so use at your own risk. Hope this helps, Aaron -----Original Message----- From: David Horner [mailto:[email protected]] Sent: Friday, June 11, 2004 3:21 AM To: [email protected] Subject: MTL: MTL on VS2003 (_Ptrit Problems) Hey everyone, I've downloaded the mtl-2.1.2-21 version of MTL. However, I seem to be having the same problem as many others with the _Ptrit light1D.h(87): error C2039: '_Ptrit' : is not a member of 'std' I even went as far as to apply the contributed zip mtl_fixes_for_msvc7, but that didn't work either. Is anyone out there using VS2003 (7.1.3088) with success? I read in some of the other messages of a CVS server. Is this public and does it contain the fix for this problem? I didn't see a link to the CVS server on the website. Any ideas? Thanks, Dave _______________________________________________ This list is archived at http://www.osl.iu.edu/MailArchives/mtl-devel/
mtl.h
(application/octet-stream, 91.2 KB)
// -*- c++ -*-
//
// Copyright 1997, 1998, 1999 University of Notre Dame.
// Authors: Andrew Lumsdaine, Jeremy G. Siek, Lie-Quan Lee
//
// This file is part of the Matrix Template Library
//
// You should have received a copy of the License Agreement for the
// Matrix Template Library along with the software; see the
// file LICENSE. If not, contact Office of Research, University of Notre
// Dame, Notre Dame, IN 46556.
//
// Permission to modify the code and to distribute modified code is
// granted, provided the text of this NOTICE is retained, a notice that
// the code was modified is included with the above COPYRIGHT NOTICE and
// with the COPYRIGHT NOTICE in the LICENSE file, and that the LICENSE
// file is distributed with the modified code.
//
// LICENSOR MAKES NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED.
// By way of example, but not limitation, Licensor MAKES NO
// REPRESENTATIONS OR WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY
// PARTICULAR PURPOSE OR THAT THE USE OF THE LICENSED SOFTWARE COMPONENTS
// OR DOCUMENTATION WILL NOT INFRINGE ANY PATENTS, COPYRIGHTS, TRADEMARKS
// OR OTHER RIGHTS.
//
//===========================================================================
#ifndef _MTL_MTL_H_
#define _MTL_MTL_H_
#include <functional>
#include <iostream>
#include "mtl/mtl_limits.h"
#include "mtl/mtl_complex.h"
#include "mtl/fast.h"
#include "mtl/dense1D.h"
#include "mtl/mtl_exception.h"
#include "mtl/matrix_traits.h"
#include "mtl/transform_iterator.h"
#include "mtl/scaled1D.h"
#include "mtl/abs.h"
#ifdef USE_DOUBLE_DOUBLE
#include "contrib/double_double/double_double.h"
#endif
#if USE_BLAIS
#include "mtl/blais.h"
#endif
#include "mtl/matrix.h"
/*
This is a nasty hack necessitated by several things:
1. C++ does not allow temporaries to be passed
into a reference argument.
2. Many MTL expressions result in temporaries
3. Some MTL matrix classes (static matrix) can
not be passed by value for the output argument
since they are not handles.
*/
#define MTL_OUT(X) const X&
namespace mtl {
template <class T>
inline T sign(const T& x) { return (x < 0) ? T(-1) : T(1); }
template <class T>
inline T xfer_sign(const T& x, const T& y)
{
return (y < 0) ? -MTL_ABS(x) : MTL_ABS(x);
}
//: for tri_solve and others
//!noindex:
class right_side { };
//: for tri_solve and others
//!noindex:
class left_side { };
#include "mtl/dim_calc.h"
template <class Vector> inline
typename linalg_traits<Vector>::value_type
sum__(const Vector& x, fast::count<0>)
{
typedef typename linalg_traits<Vector>::value_type vt;
return mtl_algo::accumulate(x.begin(), x.end(), vt());
}
#if USE_BLAIS
template <class Vector, int N> inline
typename linalg_traits<Vector>::value_type
sum__(const Vector& x, fast::count<N>)
{
typedef typename linalg_traits<Vector>::value_type vt;
return fast::accumulate(x.begin(), fast::count<N>(), vt());
}
#endif
//: Sum: <tt>s <- sum_i(x(i))</tt>
//!category: algorithms
//!component: function
//!definition: mtl.h
//!example: vec_sum.cc
//!complexity: linear
//!typereqs: The addition operator must be defined for <TT>Vector::value_type</TT>.
// The sum of all of the elements in the container.
template <class Vector> inline
typename linalg_traits<Vector>::value_type
sum(const Vector& x)
{
return sum__(x, dim_n<Vector>::RET());
}
#include "mtl/mtl_set.h"
template <class S, class T, class R>
struct mtl_multiplies : std::binary_function<S, T, R> {
typedef S first_argument_type;
typedef T second_argument_type;
typedef R result_type;
R operator () (const S& x, const T& y) const { return x * y; }
};
template <class Vector, class T> inline
void
oned_scale(Vector& x, const T& alpha, fast::count<0>)
{
typedef typename Vector::value_type VT;
mtl_algo::transform(x.begin(), x.end(), x.begin(),
std::bind1st(mtl_multiplies<T,VT,VT>(), alpha));
}
#if USE_BLAIS
template <class Vector, class T, int N> inline
void
oned_scale(Vector& x, const T& alpha, fast::count<N>)
{
typedef typename Vector::value_type VT;
fast::transform(x.begin(), fast::count<N>(), x.begin(),
std::bind1st(mtl_multiplies<T,VT,VT>(), alpha));
}
#endif
template <class Vector, class T> inline
void
scale_dim(Vector& x, const T& alpha, oned_tag)
{
oned_scale(x, alpha, dim_n<Vector>::RET());
}
template <class Matrix, class T>
inline void
scale_dim(Matrix& A, const T& alpha, twod_tag)
{
typename Matrix::iterator i;
typename Matrix::OneD::iterator j, jend;
for (i = A.begin(); i != A.end(); ++i) {
j = (*i).begin(); jend = (*i).end();
for (; j != jend; ++j)
*j *= alpha;
}
}
//: Scale: <tt>A <- alpha*A or x <- alpha x</tt>
//
// Multiply all the elements in <tt>A</tt> (or <tt>x</tt>) by
// <tt>alpha</tt>.
//
//!category: algorithms
//!component: function
//!example: vec_scale_algo.cc
//!complexity: O(n)
//!definition: mtl.h
//!typereqs: <TT>Vector</TT> must be mutable
//!typereqs: <TT>T</TT> is convertible to <TT>Vector</TT>'s <TT>value_type</TT>
//!typereqs: The multiplication operator must be defined for <TT>Vector::value_type</TT> and <tt>T</tt>
template <class LinalgObj, class T>
inline void
scale(MTL_OUT(LinalgObj) A, const T& alpha)
{
typedef typename linalg_traits<LinalgObj>::dimension Dim;
scale_dim(const_cast<LinalgObj&>(A), alpha, Dim());
}
//: Set Diagonal: <tt>A(i,i) <- alpha</tt>
//
// Set the value of the elements on the main diagonal of A to alpha.
//
//!category: algorithms
//!component: function
//!definition: mtl.h
//!example: tri_pack_sol.cc
//!typereqs: <tt>T</tt> must be convertible to <tt>Matrix::value_type</tt>.
//!complexity: O(min(m,n)) for dense matrices, O(nnz) for sparse matrices (except envelope, which is O(m))
template <class Matrix, class T>
inline void
set_diagonal(MTL_OUT(Matrix) A_, const T& alpha)
{
Matrix& A = const_cast<Matrix&>(A_);
typedef typename mtl::matrix_traits<Matrix>::size_type Int;
if (! A.is_unit())
for (Int i = 0; i < A.nrows() && i < A.ncols(); ++i)
A(i,i) = alpha;
}
//: add absolute value
//!noindex:
struct abs_add {
template <class T, class U>
T operator()(const T& a, const U& b) {
return a + MTL_ABS(b);
}
};
template <class Vector>
inline typename linalg_traits<Vector>::magnitude_type
oned_one_norm(const Vector& x, fast::count<0>)
{
typedef typename linalg_traits<Vector>::magnitude_type T;
return mtl_algo::accumulate(x.begin(), x.end(), T(), abs_add());
}
#if USE_BLAIS
template <class Vector, int N>
inline typename linalg_traits<Vector>::magnitude_type
oned_one_norm(const Vector& x, fast::count<N>)
{
typedef typename
number_traits<typename Vector::value_type>::magnitude_type T;
return fast::accumulate(x.begin(), fast::count<N>(), T(), abs_add());
}
#endif
template <class Vector>
inline typename linalg_traits<Vector>::magnitude_type
one_norm(const Vector& x, oned_tag)
{
return oned_one_norm(x, dim_n<Vector>::RET());
}
//: add square
//!noindex:
struct sqr_add {
template <class T, class U>
T operator()(const T& a, const U& b) {
return a + MTL_ABS(b * b);
}
};
template <class Vector>
inline typename linalg_traits<Vector>::magnitude_type
oned_two_norm(const Vector& x, fast::count<0>)
{
typedef typename Vector::value_type T;
typedef typename number_traits<T>::magnitude_type M;
using std::sqrt;
return ::sqrt(mtl_algo::accumulate(x.begin(), x.end(), M(), sqr_add()));
}
#if USE_BLAIS
template <class Vector, int N>
inline typename linalg_traits<Vector>::magnitude_type
oned_two_norm(const Vector& x, fast::count<N>)
{
typedef typename Vector::value_type T;
typedef typename number_traits<T>::magnitude_type M;
using std::sqrt;
return ::sqrt(fast::accumulate(x.begin(), fast::count<N>(), M(), sqr_add()));
}
#endif
//: Two Norm: <tt>s <- sqrt(sum_i(|x(i)^2|))</tt>
//
// The square root of the sum of the squares of the elements of the container.
//
//!category: algorithms
//!component: function
//!definition: mtl.h
//!example: vec_two_norm.cc
//!complexity: O(n)
//!typereqs: <tt>Vector</tt> must have an associated magnitude_type that is the type of the absolute value of <tt>Vector::value_type</tt>.
//!typereqs: There must be <tt>abs()</tt> defined for <tt>Vector::value_type</tt>.
//!typereqs: The addition must be defined for magnitude_type.
//!typereqs: <tt>sqrt()</tt> must be defined for magnitude_type.
template <class Vector>
inline typename linalg_traits<Vector>::magnitude_type
two_norm(const Vector& x)
{
return oned_two_norm(x, dim_n<Vector>::RET());
}
//: add square
//!noindex:
struct sqr_ {
template <class T, class U>
T operator()(const T& a, const U& b) {
return a + MTL_ABS(b * b);
}
};
//: Sum of the Squares
//
//!category: algorithms
//!component: function
//!definition: mtl.h
//!complexity: O(n)
template <class Vector>
inline typename linalg_traits<Vector>::value_type
sum_squares(const Vector& x)
{
typedef typename linalg_traits<Vector>::value_type T;
return mtl_algo::accumulate(x.begin(), x.end(), T(), sqr_add());
}
//: compare absolute values
//!noindex:
struct abs_cmp { template <class T>
bool operator()(const T& a, const T& b) {
return MTL_ABS(a) < MTL_ABS(b);
}};
template <class Vec>
inline typename linalg_traits<Vec>::magnitude_type
infinity_norm(const Vec& x, oned_tag)
{
return MTL_ABS(*mtl_algo::max_element(x.begin(), x.end(), abs_cmp()));
}
//: use by one and inf norm
//!noindex:
template <class Matrix>
inline typename linalg_traits<Matrix>::magnitude_type
major_norm__(const Matrix& A)
{
typedef typename linalg_traits<Matrix>::magnitude_type T;
typedef typename matrix_traits<Matrix>::size_type Int;
T norm = 0;
T sum = 0;
typename Matrix::const_iterator i;
typename Matrix::OneD::const_iterator j;
i = A.begin();
/* get the first sum */
if (i != A.end()) {
j = (*i).begin();
sum = T(0);
for (; j != (*i).end(); ++j)
sum = sum + MTL_ABS(*j);
norm = sum;
++i;
}
for (; i != A.end(); ++i) {
j = (*i).begin();
if (A.is_unit() && Int(i.index()) < MTL_MIN(A.nrows(), A.ncols()))
sum = T(1);
else sum = T(0);
for (; j != (*i).end(); ++j)
sum = sum + MTL_ABS(*j);
norm = MTL_MAX(MTL_ABS(norm), MTL_ABS(sum));
}
return norm;
}
//: used by one and inf norm
//!noindex:
template <class Matrix>
inline typename linalg_traits<Matrix>::magnitude_type
minor_norm__(const Matrix& A)
{
typedef typename linalg_traits<Matrix>::magnitude_type T;
typedef typename matrix_traits<Matrix>::size_type Int;
typename Matrix::const_iterator i;
typename Matrix::OneD::const_iterator j, jend;
dense1D<T> sums(A.minor(), T());
if (A.is_unit()) {
for (Int x = 0; x < MTL_MIN(A.nrows(), A.ncols()); ++x)
sums[x] = T(1);
}
for (i = A.begin(); i != A.end(); ++i) {
j = (*i).begin(); jend = (*i).end();
for (; j != jend; ++j)
sums[j.index()] += MTL_ABS(*j);
}
return infinity_norm(sums, oned_tag());
}
/* this handles both the major and minor norm
for symmetric matrices */
template <class Matrix>
inline typename linalg_traits<Matrix>::magnitude_type
symmetric_norm(const Matrix& A, row_tag)
{
typedef typename linalg_traits<Matrix>::magnitude_type T;
typename Matrix::const_iterator i;
typename Matrix::OneD::const_iterator j, jend;
dense1D<T> sums(A.minor(), T(0));
for (i = A.begin(); i != A.end(); ++i) {
j = (*i).begin();
jend = (*i).end();
if (A.is_upper()) { /* handle the diagonal elements */
sums[j.row()] += MTL_ABS(*j);
++j;
} else
--jend;
for (; j != jend; ++j) {
sums[j.row()] += MTL_ABS(*j);
sums[j.column()] += MTL_ABS(*j);
}
if (A.is_lower())
sums[j.row()] += MTL_ABS(*j);
}
return infinity_norm(sums, oned_tag());
}
template <class Matrix>
inline typename linalg_traits<Matrix>::magnitude_type
symmetric_norm(const Matrix& A, column_tag)
{
typedef typename linalg_traits<Matrix>::magnitude_type T;
typename Matrix::const_iterator i;
typename Matrix::OneD::const_iterator j, jend;
dense1D<T> sums(A.minor(), T(0));
for (i = A.begin(); i != A.end(); ++i) {
j = (*i).begin();
jend = (*i).end();
if (A.is_lower()) { /* handle the diagonal elements */
sums[j.row()] += MTL_ABS(*j);
++j;
} else
--jend;
for (; j != jend; ++j) {
sums[j.row()] += MTL_ABS(*j);
sums[j.column()] += MTL_ABS(*j);
}
if (A.is_upper())
sums[j.row()] += MTL_ABS(*j);
}
return infinity_norm(sums, oned_tag());
}
template <class Matrix>
inline typename linalg_traits<Matrix>::magnitude_type
symmetric_norm(const Matrix& A)
{
typedef typename matrix_traits<Matrix>::orientation Orien;
return symmetric_norm(A, Orien());
}
template <class Matrix>
inline typename linalg_traits<Matrix>::magnitude_type
diagonal_one_norm(const Matrix& A)
{
typedef typename linalg_traits<Matrix>::magnitude_type T;
typename Matrix::const_iterator i;
typename Matrix::OneD::const_iterator j, jend;
dense1D<T> sums(A.ncols(), T(0));
for (i = A.begin(); i != A.end(); ++i) {
j = (*i).begin(); jend = (*i).end();
for (; j != jend; ++j)
sums[j.column()] += MTL_ABS(*j);
}
return infinity_norm(sums);
}
template <class Matrix>
inline typename linalg_traits<Matrix>::magnitude_type
diagonal_infinity_norm(const Matrix& A)
{
typedef typename linalg_traits<Matrix>::magnitude_type T;
typename Matrix::const_iterator i;
typename Matrix::OneD::const_iterator j, jend;
dense1D<T> sums(A.nrows(), T(0));
for (i = A.begin(); i != A.end(); ++i) {
j = (*i).begin(); jend = (*i).end();
for (; j != jend; ++j)
sums[j.row()] += MTL_ABS(*j);
}
return infinity_norm(sums);
}
//: dispatch function
//!noindex:
template <class Matrix>
inline typename linalg_traits<Matrix>::magnitude_type
one_norm__(const Matrix& A, column_tag)
{
return major_norm__(A);
}
//: dispatch function
//!noindex:
template <class Matrix>
inline typename linalg_traits<Matrix>::magnitude_type
one_norm__(const Matrix& A, row_tag)
{
return minor_norm__(A);
}
template <class Matrix, class Shape>
inline typename linalg_traits<Matrix>::magnitude_type
twod_one_norm(const Matrix& A, Shape)
{
typedef typename Matrix::orientation Orien;
return one_norm__(A, Orien());
}
template <class Matrix>
inline typename linalg_traits<Matrix>::magnitude_type
twod_one_norm(const Matrix& A, symmetric_tag)
{
return symmetric_norm(A);
}
template <class Matrix>
inline typename linalg_traits<Matrix>::magnitude_type
twod_one_norm(const Matrix& A, diagonal_tag)
{
return diagonal_one_norm(A);
}
template <class Linalg>
inline typename linalg_traits<Linalg>::magnitude_type
one_norm(const Linalg& A, twod_tag)
{
typedef typename matrix_traits<Linalg>::shape Shape;
return twod_one_norm(A, Shape());
}
//: One Norm: <tt>s <- sum(|x_i|) or s <- max_i(sum_j(|A(i,j)|))</tt>
//
// For vectors, the sum of the absolute values of the elements.
// For matrices, the maximum of the column sums.
// Note: not implemented yet for unit triangle matrices.
//
//!category: algorithms
//!component: function
//!definition: mtl.h
//!example: vec_one_norm.cc
//!complexity: O(n)
//!typereqs: The vector or matrix must have an associated magnitude_type that
// is the type of the absolute value of its <tt>value_type</tt>.
//!typereqs: There must be <tt>abs()</tt> defined for <tt>Vector::value_type</tt>.
//!typereqs: The addition must be defined for magnitude_type.
template <class LinalgObj>
inline typename linalg_traits<LinalgObj>::magnitude_type
one_norm(const LinalgObj& A)
{
typedef typename linalg_traits<LinalgObj>::dimension Dim;
return one_norm(A, Dim());
}
//: dispatch function
//!noindex:
template <class Matrix>
inline typename linalg_traits<Matrix>::magnitude_type
infinity_norm__(const Matrix& A, row_tag)
{
return major_norm__(A);
}
//: dispatch function
//!noindex:
template <class Matrix>
inline typename linalg_traits<Matrix>::magnitude_type
infinity_norm__(const Matrix& A, column_tag)
{
return minor_norm__(A);
}
template <class Matrix, class Shape>
inline typename linalg_traits<Matrix>::magnitude_type
twod_infinity_norm(const Matrix& A, Shape)
{
typedef typename Matrix::orientation Orien;
return infinity_norm__(A, Orien());
}
template <class Matrix>
inline typename linalg_traits<Matrix>::magnitude_type
twod_infinity_norm(const Matrix& A, symmetric_tag)
{
return symmetric_norm(A);
}
template <class Matrix>
inline typename linalg_traits<Matrix>::magnitude_type
twod_infinity_norm(const Matrix& A, diagonal_tag)
{
return diagonal_infinity_norm(A);
}
template <class Matrix>
inline typename linalg_traits<Matrix>::magnitude_type
infinity_norm(const Matrix& A, twod_tag)
{
typedef typename matrix_traits<Matrix>::shape Shape;
return twod_infinity_norm(A, Shape());
}
//: Infinity Norm: <tt>s <- max_j(sum_i(|A(i,j)|)) or s <- max_i(|x(i)|)</tt>
//
// For matrices, the maximum of the row sums.
// For vectors, the maximum absolute value of any of its element.
//
//!category: algorithms
//!component: function
//!definition: mtl.h
//!complexity: O(n) for vectors, O(m*n) for dense matrices, O(nnz) for sparse
//!example: vec_inf_norm.cc
//!typereqs: The vector or matrix must have an associated magnitude_type that is the type of the absolute value of its <tt>value_type</tt>.
//!typereqs: There must be <tt>abs()</tt> defined for <tt>Vector::value_type</tt>.
//!typereqs: The addition must be defined for magnitude_type.
template <class LinalgObj>
inline typename linalg_traits<LinalgObj>::magnitude_type
infinity_norm(const LinalgObj& A)
{
typedef typename linalg_traits<LinalgObj>::dimension Dim;
return infinity_norm(A, Dim());
}
//: Max Index: <tt>i <- index of max(|x(i)|)</tt>
//!category: algorithms
//!component: function
//!definition: mtl.h
//!complexity: O(n)
// The location (index) of the element with the maximum absolute value.
//!example: max_index.cc
//!typereqs: <tt>Vec::value_type</tt> must be LessThanComparible.
template <class Vec>
inline typename Vec::size_type
max_index(const Vec& x)
{
typename Vec::const_iterator maxi =
mtl_algo::max_element(x.begin(), x.end(), abs_cmp());
return maxi.index();
}
//: Maximum Absolute Index: <tt>i <- index of max(|x(i)|)</tt>
//!category: algorithms
//!component: function
//!definition: mtl.h
//!complexity: O(n)
// The location (index) of the element with the maximum absolute value.
//!example: max_abs_index.cc
//!typereqs: The vector or matrix must have an associated magnitude_type that
// is the type of the absolute value of its <tt>value_type</tt>.
//!typereqs: There must be <tt>abs()</tt> defined for <tt>Vector::value_type</tt>.
//!typereqs: The magnitude type must be LessThanComparible.
template <class Vec>
inline typename Vec::size_type
max_abs_index(const Vec& x)
{
typename Vec::const_iterator maxi =
mtl_algo::max_element(x.begin(), x.end(), abs_cmp());
return maxi.index();
}
//: Minimum Index: <tt>i <- index of min(x(i))</tt>
//!category: algorithms
//!component: function
//!definition: mtl.h
//!complexity: O(n)
// The location (index) of the element with the minimum value.
//!example: min_abs_index.cc
//!typereqs: <tt>Vec::value_type</tt> must be LessThanComparible.
template<class Vec>
inline typename Vec::size_type
min_index(const Vec& x)
{
typename Vec::const_iterator mini =
mtl_algo::min_element(x.begin(), x.end());
return mini.index();
}
//: Minimum Absolute Index: <tt>i <- index of min(|x(i)|)</tt>
//!category: algorithms
//!component: function
//!definition: mtl.h
//!complexity: O(n)
// The location (index) of the element with the minimum absolute value.
//!example: max_index.cc
//!typereqs: The vector or matrix must have an associated magnitude_type that
// is the type of the absolute value of its <tt>value_type</tt>.
//!typereqs: There must be <tt>abs()</tt> defined for <tt>Vector::value_type</tt>.
//!typereqs: The magnitude type must be LessThanComparible.
template<class Vec>
inline typename Vec::size_type
min_abs_index(const Vec& x)
{
typename Vec::const_iterator mini =
mtl_algo::min_element(x.begin(), x.end(), abs_cmp());
return mini.index();
}
//: Max Value: <tt>s <- max(x(i))</tt>
//!category: algorithms
//!component: function
//!definition: mtl.h
//!example: vec_max.cc
//!complexity: O(n)
//!typereqs: <tt>Vec::value_type</tt> must be LessThanComparible.
// Returns the value of the element with the maximum value
template <class VectorT>
inline typename VectorT::value_type
max(const VectorT& x)
{
return *mtl_algo::max_element(x.begin(), x.end());
}
//: Min Value: <tt>s <- min(x_i)</tt>
//!category: algorithms
//!component: function
//!complexity: O(n)
//!definition: mtl.h
//!typereqs: <tt>Vec::value_type</tt> must be LessThanComparible.
template <class VectorT>
inline typename VectorT::value_type
min(const VectorT& x)
{
return *mtl_algo::min_element(x.begin(), x.end());
}
#define MTL_BLAS_GROT
//use blas version always, since there is a bug in the lapack verions
// of givens_rotation according to Andy's email
//: Givens Plane Rotation
//!category: functors
//!component: type
//!definition: mtl.h
//!example: apply_givens.cc
//
// Input a and b to the constructor to create a givens plane rotation
// object. Then apply the rotation to two vectors. There is a
// specialization of the givens rotation for complex numbers.
//
// <codeblock>
// [ c s ] [ a ] = [ r ]
// [ -s c ] [ b ] [ 0 ]
// </codeblock>
//
//!typereqs: the addition operator must be defined for <tt>T</tt>
//!typereqs: the multiplication operator must be defined for <tt>T</tt>
//!typereqs: the division operator must be defined for <tt>T</tt>
//!typereqs: the abs() function must be defined for <tt>T</tt>
template <class T>
class givens_rotation {
public:
//: Default constructor
inline givens_rotation()
:
#ifdef MTL_BLAS_GROT
a_(0), b_(0),
#endif
c_(0), s_(0)
#ifndef MTL_BLAS_GROT
, r_(0)
#endif
{ }
//: Givens Plane Rotation Constructor
inline givens_rotation(T a_in, T b_in) {
#ifdef MTL_BLAS_GROT // old BLAS version
T roe;
if (MTL_ABS(a_in) > MTL_ABS(b_in))
roe = a_in;
else
roe = b_in;
T scal = MTL_ABS(a_in) + MTL_ABS(b_in);
T r, z;
if (scal != T(0)) {
T a_scl = a_in / scal;
T b_scl = b_in / scal;
r = scal * sqrt(a_scl * a_scl + b_scl * b_scl);
if (roe < T(0)) r *= -1;
c_ = a_in / r;
s_ = b_in / r;
z = 1;
if (MTL_ABS(a_in) > MTL_ABS(b_in))
z = s_;
else if (MTL_ABS(b_in) >= MTL_ABS(a_in) && c_ != T(0))
z = T(1) / c_;
} else {
c_ = 1; s_ = 0; r = 0; z = 0;
}
a_ = r;
b_ = z;
#else // similar LAPACK slartg version, modified to the NEW BLAS proposal
T a = a_in, b = b_in;
if (b == T(0)) {
c_ = T(1);
s_ = T(0);
r_ = a;
} else if (a == T(0)) {
c_ = T(0);
s_ = sign(b);
r_ = b;
} else {
// cs = |a| / sqrt(|a|^2 + |b|^2)
// sn = sign(a) * b / sqrt(|a|^2 + |b|^2)
T abs_a = MTL_ABS(a);
T abs_b = MTL_ABS(b);
if (abs_a > abs_b) {
// 1/cs = sqrt( 1 + |b|^2 / |a|^2 )
T t = abs_b / abs_a;
T tt = sqrt(T(1) + t * t);
c_ = T(1) / tt;
s_ = t * c_;
r_ = a * tt;
} else {
// 1/sn = sign(a) * sqrt( 1 + |a|^2/|b|^2 )
T t = abs_a / abs_b;
T tt = sqrt(T(1) + t * t);
s_ = sign(a) / tt;
c_ = t * s_;
r_ = b * tt;
}
}
#endif
}
inline void set_cs(T cin, T sin) { c_ = cin; s_ = sin; }
//: Apply plane rotation to two real scalars. (name change a VC++ workaround)
inline void scalar_apply(T& x, T& y) {
T tmp = c_ * x + s_ * y;
y = c_ * y - s_ * x;
x = tmp;
}
//: Apply plane rotation to two vectors.
template <class VecX, class VecY>
inline void apply(MTL_OUT(VecX) x_, MTL_OUT(VecY) y_) MTL_THROW_ASSERTION {
VecX& x = const_cast<VecX&>(x_);
VecY& y = const_cast<VecY&>(y_);
MTL_ASSERT(x.size() <= y.size(), "mtl::givens_rotation::apply()");
typename VecX::iterator xi = x.begin();
typename VecX::iterator xend = x.end();
typename VecY::iterator yi = y.begin();
while (mtl::not_at(xi, xend)) {
scalar_apply(*xi, *yi);
++xi; ++yi;
}
}
#ifdef MTL_BLAS_GROT
inline T a() { return a_; }
inline T b() { return b_; }
#endif
inline T c() { return c_; }
inline T s() { return s_; }
#ifndef MTL_BLAS_GROT
inline T r() { return r_; }
#endif
protected:
#ifdef MTL_BLAS_GROT
T a_, b_;
#endif
T c_, s_;
#ifndef MTL_BLAS_GROT
T r_;
#endif
};
using std::real;
using std::imag;
#if MTL_PARTIAL_SPEC
//: The specialization for complex numbers.
//!category: functors
//!component: type
template <class T>
class givens_rotation < std::complex<T> > {
typedef std::complex<T> C;
public:
//:
inline givens_rotation() : cs(0), sn(0)
#ifndef MTL_BLAS_GROT
, r_(0)
#endif
{ }
inline T abs_sq(C t) { return real(t) * real(t) + imag(t) * imag(t); }
inline T abs1(C t) { return MTL_ABS(real(t)) + MTL_ABS(imag(t)); }
//:
inline givens_rotation(C a_in, C b_in) {
#ifdef MTL_BLAS_GROT
T a = std::abs(a_in), b = std::abs(b_in);
if ( a == T(0) ) {
cs = T(0);
sn = C(1.);
//in zrotg there is an assignment for ca, what is that for?
} else {
T scale = a + b;
T norm = std::sqrt(abs_sq(a_in/scale)+abs_sq(b_in/scale)) * scale;
cs = a / norm;
sn = a_in/a * std::conj(b_in)/norm;
//in zrotg there is an assignment for ca, what is that for?
}
#else // LAPACK version, clartg
C f(a_in), g(b_in);
if (g == C(0)) {
cs = T(1);
sn = C(0);
r_ = f;
} else if (f == C(0)) {
cs = T(0);
sn = MTL_CONJ(g) / MTL_ABS(g);
r_ = MTL_ABS(g);
} else {
C fs, gs, ss, t;
T d, di, f1, f2, fa, g1, g2, ga;
f1 = abs1(f);
g1 = abs1(g);
if (f1 >= g1) {
gs = g / f1;
g2 = abs_sq(gs);
fs = f / f1;
f2 = abs_sq(fs);
d = sqrt(T(1) + g2 / f2);
cs = T(1) / d;
sn = MTL_CONJ(gs) * fs * (cs / f2);
r_ = f * d;
} else {
fs = f / g1;
f2 = abs_sq(fs);
fa = sqrt(f2);
gs = g / g1;
g2 = abs_sq(gs);
ga = sqrt(g2);
d = sqrt(T(1) + f2 / g2);
di = T(1) / d;
cs = (fa / ga ) * di;
ss = (MTL_CONJ(gs) * fs) / (fa * ga);
sn = ss * di;
r_ = g * ss * d;
}
}
#endif
}
//: Apply plane rotation to two vectors.
template <class VecX, class VecY>
inline void apply(MTL_OUT(VecX) x_, MTL_OUT(VecY) y_) MTL_THROW_ASSERTION {
VecX& x = const_cast<VecX&>(x_);
VecY& y = const_cast<VecY&>(y_);
MTL_ASSERT(x.size() <= y.size(), "mtl::givens_rotation::apply()");
typename VecX::iterator xi = x.begin();
typename VecX::iterator xend = x.end();
typename VecY::iterator yi = y.begin();
while (mtl::not_at(xi, xend)) {
scalar_apply(*xi, *yi);
++xi; ++yi;
}
}
//: Apply plane rotation to two complex scalars.
inline void scalar_apply(C& x, C& y) {
complex<T> temp = MTL_CONJ(cs) * x + MTL_CONJ(sn) * y;
y = cs * y - sn * x;
x = temp;
}
inline void set_cs(const T& cs_, const C& sn_) {
cs = cs_; sn = sn_;
}
inline T c() { return cs; }
inline C s() { return sn; }
#ifndef MTL_BLAS_GROT
inline C r() { return r_; }
#endif
protected:
T cs;
C sn;
#ifndef MTL_BLAS_GROT
C r_;
#endif
};
#else
//: The specialization for complex numbers.
//!category: functors
//!component: type
template<>
class givens_rotation < std::complex<double> > {
typedef double T;
typedef std::complex<T> C;
public:
//:
inline givens_rotation() : cs(0), sn(0)
#ifndef MTL_BLAS_GROT
, r_(0)
#endif
{ }
inline T abs_sq(C t) { return real(t) * real(t) + imag(t) * imag(t); }
inline T abs1(C t) { return MTL_ABS(real(t)) + MTL_ABS(imag(t)); }
//:
inline givens_rotation(C a_in, C b_in) {
#ifdef MTL_BLAS_GROT
T a = std::abs(a_in), b = std::abs(b_in);
if ( a == T(0) ) {
cs = T(0);
sn = C(1.);
//in zrotg there is an assignment for ca, what is that for?
} else {
T scale = a + b;
T norm = std::sqrt(abs_sq(a_in/scale)+abs_sq(b_in/scale)) * scale;
cs = a / norm;
sn = a_in/a * std::conj(b_in)/norm;
//in zrotg there is an assignment for ca, what is that for?
}
#else // LAPACK version, clartg
C f(a_in), g(b_in);
if (g == C(0)) {
cs = T(1);
sn = C(0);
r_ = f;
} else if (f == C(0)) {
cs = T(0);
sn = MTL_CONJ(g) / MTL_ABS(g);
r_ = MTL_ABS(g);
} else {
C fs, gs, ss, t;
T d, di, f1, f2, fa, g1, g2, ga;
f1 = abs1(f);
g1 = abs1(g);
if (f1 >= g1) {
gs = g / f1;
g2 = abs_sq(gs);
fs = f / f1;
f2 = abs_sq(fs);
d = sqrt(T(1) + g2 / f2);
cs = T(1) / d;
sn = MTL_CONJ(gs) * fs * (cs / f2);
r_ = f * d;
} else {
fs = f / g1;
f2 = abs_sq(fs);
fa = sqrt(f2);
gs = g / g1;
g2 = abs_sq(gs);
ga = sqrt(g2);
d = sqrt(T(1) + f2 / g2);
di = T(1) / d;
cs = (fa / ga ) * di;
ss = (MTL_CONJ(gs) * fs) / (fa * ga);
sn = ss * di;
r_ = g * ss * d;
}
}
#endif
}
//: Apply plane rotation to two vectors.
template <class VecX, class VecY>
inline void apply(MTL_OUT(VecX) x_, MTL_OUT(VecY) y_) MTL_THROW_ASSERTION {
VecX& x = const_cast<VecX&>(x_);
VecY& y = const_cast<VecY&>(y_);
MTL_ASSERT(x.size() <= y.size(), "mtl::givens_rotation::apply()");
typename VecX::iterator xi = x.begin();
typename VecX::iterator xend = x.end();
typename VecY::iterator yi = y.begin();
while (mtl::not_at(xi, xend)) {
scalar_apply(*xi, *yi);
++xi; ++yi;
}
}
//: Apply plane rotation to two complex scalars.
inline void scalar_apply(C& x, C& y) {
complex<T> temp = MTL_CONJ(cs) * x + MTL_CONJ(sn) * y;
y = cs * y - sn * x;
x = temp;
}
T c() { return cs; }
C s() { return sn; }
#ifndef MTL_BLAS_GROT
inline C r() { return r_; }
#endif
protected:
T cs;
C sn;
#ifndef MTL_BLAS_GROT
C r_;
#endif
};
//: The specialization for complex numbers.
//!category: functors
//!component: type
template<>
class givens_rotation < std::complex<float> > {
typedef float T;
typedef std::complex<T> C;
public:
//:
inline givens_rotation() : cs(0), sn(0)
#ifndef MTL_BLAS_GROT
, r_(0)
#endif
{ }
inline T abs_sq(C t) { return real(t) * real(t) + imag(t) * imag(t); }
inline T abs1(C t) { return MTL_ABS(real(t)) + MTL_ABS(imag(t)); }
//:
inline givens_rotation(C a_in, C b_in) {
#ifdef MTL_BLAS_GROT
T a = std::abs(a_in), b = std::abs(b_in);
if ( a == T(0) ) {
cs = T(0);
sn = C(1.);
//in zrotg there is an assignment for ca, what is that for?
} else {
T scale = a + b;
T norm = std::sqrt(abs_sq(a_in/scale)+abs_sq(b_in/scale)) * scale;
cs = a / norm;
sn = a_in/a * std::conj(b_in)/norm;
//in zrotg there is an assignment for ca, what is that for?
}
#else // LAPACK version, clartg
C f(a_in), g(b_in);
if (g == C(0)) {
cs = T(1);
sn = C(0);
r_ = f;
} else if (f == C(0)) {
cs = T(0);
sn = MTL_CONJ(g) / MTL_ABS(g);
r_ = MTL_ABS(g);
} else {
C fs, gs, ss, t;
T d, di, f1, f2, fa, g1, g2, ga;
f1 = abs1(f);
g1 = abs1(g);
if (f1 >= g1) {
gs = g / f1;
g2 = abs_sq(gs);
fs = f / f1;
f2 = abs_sq(fs);
d = sqrt(T(1) + g2 / f2);
cs = T(1) / d;
sn = MTL_CONJ(gs) * fs * (cs / f2);
r_ = f * d;
} else {
fs = f / g1;
f2 = abs_sq(fs);
fa = sqrt(f2);
gs = g / g1;
g2 = abs_sq(gs);
ga = sqrt(g2);
d = sqrt(T(1) + f2 / g2);
di = T(1) / d;
cs = (fa / ga ) * di;
ss = (MTL_CONJ(gs) * fs) / (fa * ga);
sn = ss * di;
r_ = g * ss * d;
}
}
#endif
}
//: Apply plane rotation to two vectors.
template <class VecX, class VecY>
inline void apply(MTL_OUT(VecX) x_, MTL_OUT(VecY) y_) MTL_THROW_ASSERTION {
VecX& x = const_cast<VecX&>(x_);
VecY& y = const_cast<VecY&>(y_);
MTL_ASSERT(x.size() <= y.size(), "mtl::givens_rotation::apply()");
typename VecX::iterator xi = x.begin();
typename VecX::iterator xend = x.end();
typename VecY::iterator yi = y.begin();
while (mtl::not_at(xi, xend)) {
scalar_apply(*xi, *yi);
++xi; ++yi;
}
}
//: Apply plane rotation to two complex scalars.
inline void scalar_apply(C& x, C& y) {
complex<T> temp = MTL_CONJ(cs) * x + MTL_CONJ(sn) * y;
y = cs * y - sn * x;
x = temp;
}
T c() { return cs; }
C s() { return sn; }
#ifndef MTL_BLAS_GROT
inline C r() { return r_; }
#endif
protected:
T cs;
C sn;
#ifndef MTL_BLAS_GROT
C r_;
#endif
};
#endif
#undef MTL_BLAS_GROT
//do allow internal macro escape out of the scope
//: Modified Givens Transformation
//!category: functors
//!component: type
//
// This class is under construction. Like the givens rotation class,
// there will be a real and complex class.
template <class T>
class modified_givens {
};
template <class T>
inline T two_norm3(const T& x, const T& y, const T& z) {
return sqrt(x*x + y*y + z*z);
}
//: Generate Householder Transform
//
// Ok to alias x and v to the same vector.
// T can be real or complex.
// Equivalent to LAPACK's xLARFG
template <class T, class Vec>
inline void generate_householder(T& alpha, const Vec& x,
Vec& v, T& tau) MTL_THROW_ASSERTION
{
MTL_ASSERT(x.size() == v.size(), "mtl::generate_householder");
typedef typename number_traits<T>::magnitude_type Real;
typename Vec::subrange_type subx = x(0, x.size() - 1);
typename Vec::subrange_type subv = v(0, v.size() - 1);
v[v.size() - 1] = x[x.size() - 1];
Real xnorm = two_norm(x);
Real alpha_r = real(alpha);
Real alpha_i = imag(alpha);
if (xnorm == Real(0) && alpha_i == Real(0))
tau = T(0); // H = I
else {
Real beta = -xfer_sign(two_norm3(alpha_r, alpha_i, xnorm), alpha_r);
Real safe_min = std::numeric_limits<Real>::min();
Real r_safe_min = Real(1) / safe_min;
int count = 0;
while (MTL_ABS(beta) < safe_min) { // xnorm and beta may be inaccurate
if (count == 0) // so scale x and recompute them
copy(mtl::scaled(subx, r_safe_min), subv);
else
scale(subv, r_safe_min);
beta *= r_safe_min;
alpha *= r_safe_min;
++count;
}
if (count != 0) {
alpha_r = real(alpha);
alpha_i = imag(alpha);
xnorm = two_norm(x);
beta = -xfer_sign(two_norm3(alpha_r, alpha_i, xnorm), alpha_r);
}
tau = beta - (alpha / beta);
alpha = T(1) / (alpha - beta);
scale(subv, alpha);
alpha = beta;
for (int j = 0; j < count; ++j)
alpha *= safe_min;
}
}
//: Householder Transform
//
// Constructor does the generation, then call apply
//
template <class T>
class householder_transform {
typedef typename number_traits<T>::magnitude_type Real;
typedef dense1D<T> Vec;
public:
template <class Vec>
inline householder_transform(const T& alpha, Vec& x, const T& tau)
: _alpha(alpha), _v(x.size()), _tau(tau) {
generate_householder(_alpha, x, _v, _tau);
}
#if 0
// JGS conj ???
// Equivalent to LAPACK xLARF
template <class MatrixC>
inline void apply(MatrixC& C, right_side) {
if (_tau != T(0)) {
Vec w(C.nrows());
mult(C, _v, w); // w <- C * v
rank_one_update(mtl::caled(w, -_tau),// C <- C - w * v'
conj(v), C);
}
}
template <class MatrixC>
inline void apply(MatrixC& C, left_side) {
if (_tau != T(0)) {
Vec w(C.nrows());
mult(conj(C), _v, w); // w <- C' * v
rank_one_update(mtl::scaled(w,-_tau), // C <- C - w * v'
conj(v), C);
}
}
#endif
protected:
T _alpha;
Vec _v;
T _tau;
};
//: Transpose in Place: <tt>A <- A^T</tt>
// Currently this algorithm only applies to square dense matrices
// Plan to include all rectangular dense matrices..
//!category: algorithms
//!component: function
//!definition: mtl.h
template <class Matrix>
inline void
transpose(MTL_OUT(Matrix) A_) MTL_THROW_ASSERTION
{
Matrix& A = const_cast<Matrix&>(A_);
MTL_ASSERT(A.nrows() == A.ncols(), "mat::transpose()");
typedef typename matrix_traits<Matrix>::value_type T;
typedef typename mtl::matrix_traits<Matrix>::size_type Int;
for (Int i = 0; i < A.nrows(); ++i)
for (Int j = i; j < A.ncols(); ++j) {
T tmp = A(i, j);
A(i, j) = A(j, i);
A(j, i) = tmp;
}
}
//: Transpose: <tt>B <- A^T</tt>
//!precond: <tt> B(i,j) = 0 & B = A^T </tt>
//
// When matrix B is banded, it is up to the user to ensure
// that the bandwidth is sufficient to contain the elements
// from A^T. If there are elements of A^T that do not
// fall within the bandwidth, an exception will be thrown.
// (exception not implemented yet).
//
//!category: algorithms
//!component: function
//!definition: mtl.h
//!complexity: O(n^2)
template <class MatA, class MatB>
inline void
transpose(const MatA& A, MTL_OUT(MatB) B_) MTL_THROW_ASSERTION
{
MatB& B = const_cast<MatB&>(B_);
MTL_ASSERT(A.nrows() <= B.ncols(), "matmat::transpose()");
MTL_ASSERT(A.ncols() <= B.nrows(), "matmat::transpose()");
typename MatA::const_iterator i;
typename MatA::OneD::const_iterator j, jend;
for (i = A.begin(); i != A.end(); ++i) {
j = (*i).begin(); jend = (*i).end();
for (; j != jend; ++j)
B(j.column(), j.row()) = *j;
}
}
/*
This version of the algorithm depends on the compiler
hoisting the reference of z[j.row()] out of the inner loop
(for the row major case)
KCC doesn't do this, and niether does the underlying Sun C
compiler.
In order to hoist the reference by hand, I'll have to write
specializations for column major matrix and for row major matrices.
While I'm at it I'll the the unrolling stuff too.
*/
/* this is generic
*/
template <class Matrix, class VecX, class VecZ>
inline void
mult_generic__(const Matrix& A, const VecX& xx, VecZ& zz) MTL_THROW_ASSERTION
{
MTL_ASSERT(A.nrows() <= zz.size(), "mtl::mult()");
MTL_ASSERT(A.ncols() <= xx.size(), "mtl::mult()");
typedef typename matrix_traits<Matrix>::value_type T;
typename Matrix::const_iterator i;
typename Matrix::OneD::const_iterator j, jend;
typename VecX::const_iterator x = xx.begin();
typename VecZ::iterator z = zz.begin();
for (i = A.begin(); i != A.end(); ++i) {
j = (*i).begin(); jend = (*i).end();
for (; j != jend; ++j)
z[j.row()] += *j * x[j.column()];
}
}
template <class Matrix, class VecX, class VecZ>
inline void
mult_shape__(const Matrix& A, const VecX& x, VecZ& z,
banded_tag) MTL_THROW_ASSERTION
{
mult_generic__(A, x, z);
}
/* this is fast
*/
template <class Matrix, class VecX, class VecZ>
inline void
rect_mult(const Matrix& A, const VecX& xx, VecZ& zz,
row_tag, dense_tag) MTL_THROW_ASSERTION
{
MTL_ASSERT(A.nrows() <= zz.size(), "mtl::mult()");
MTL_ASSERT(A.ncols() <= xx.size(), "mtl::mult()");
typedef typename matrix_traits<Matrix>::value_type T;
typename Matrix::const_iterator i, iend;
typename Matrix::OneD::const_iterator j, jend;
typename VecX::const_iterator x = xx.begin();
typename VecZ::iterator z = zz.begin();
i = A.begin();
iend = A.end();
for (; i != iend; ++i) {
j = (*i).begin(); jend = (*i).end();
T tmp = z[j.row()];
for (; j != jend; ++j)
tmp += *j * x[j.column()];
z[j.row()] = tmp;
}
}
/*
This is slow
*/
template <class Matrix, class VecX, class VecZ>
inline void
rect_mult(const Matrix& A, const VecX& xx, VecZ& zz,
column_tag, dense_tag) MTL_THROW_ASSERTION
{
MTL_ASSERT(A.nrows() <= zz.size(), "mtl::mult()");
MTL_ASSERT(A.ncols() <= xx.size(), "mtl::mult()");
typedef typename matrix_traits<Matrix>::value_type T;
typedef typename matrix_traits<Matrix>::size_type Int;
typename VecX::const_iterator x = xx.begin();
typename VecZ::iterator z = zz.begin();
typename Matrix::const_iterator i, iend;
typename Matrix::OneD::const_iterator j, jend;
i = A.begin();
iend = A.end();
for (; i != iend; ++i) {
j = (*i).begin(); jend = (*i).end();
for (; j != jend; ++j)
z[j.row()] += *j * x[j.column()];
}
}
// x is sparse
// A is column oriented
template <class Matrix, class VecX, class VecY>
void
rect_mult(const Matrix& A, const VecX& x, VecY& y,
column_tag, sparse_tag)
{
typename VecX::const_iterator xi = x.begin();
for (; xi != x.end(); ++xi) {
mtl::add(mtl::scaled(A[xi.index()], *xi), y);
}
}
// x is sparse
// A is row oriented
template <class Matrix, class VecX, class VecY>
void
rect_mult(const Matrix& A, const VecX& x, VecY& y,
row_tag, sparse_tag)
{
typename Matrix::const_iterator Ai;
for (Ai = A.begin(); Ai != A.end(); ++Ai) {
y[Ai.index()] = mtl::dot(*Ai, x); // this is a sparse dot
}
}
template <class Matrix, class VecX, class VecZ>
inline void
mult_shape__(const Matrix& A, const VecX& x, VecZ z,
rectangle_tag) MTL_THROW_ASSERTION
{
typedef typename matrix_traits<Matrix>::orientation Orien;
typedef typename linalg_traits<VecX>::sparsity SparseX;
rect_mult(A, x, z, Orien(), SparseX());
}
template <class Matrix, class VecX, class VecZ>
inline void
mult_shape__(const Matrix& A, const VecX& x, VecZ& z,
triangle_tag)
{
mult_shape__(A, x, z, rectangle_tag());
if (A.is_unit()) {
/* actually, this still isn't quite right,
should do
add_n(x, z, z, MTL_MIN(A.nrows(), A.ncols()));
instead
*/
if (z.size() <= x.size())
mtl::add(z, x, z);
else
mtl::add(x, z, z);
}
}
template <class Matrix, class VecX, class VecZ>
inline void
mult_symm__(const Matrix& A, const VecX& x, VecZ& z, row_tag)
{
typedef typename matrix_traits<Matrix>::value_type T;
typename Matrix::const_iterator i;
typename Matrix::OneD::const_iterator j, jend;
for (i = A.begin(); i != A.end(); ++i) {
T tmp = z[i.index()];
j = (*i).begin();
jend = (*i).end();
if (A.is_upper()) {
tmp += *j * x[j.column()];
++j;
} else
--jend;
for (; j != jend; ++j) {
/* normal side */
tmp += *j * x[j.column()];
/* symmetric side */
z[j.column()] += *j * x[j.row()];
}
if (A.is_lower())
tmp += *j * x[j.column()];
z[i.index()] = tmp;
}
}
template <class Matrix, class VecX, class VecZ>
inline void
mult_symm__(const Matrix& A, const VecX& x, VecZ& z, column_tag)
{
typedef typename matrix_traits<Matrix>::value_type T;
typename Matrix::const_iterator i;
typename Matrix::OneD::const_iterator j, jend;
for (i = A.begin(); i != A.end(); ++i) {
T tmp = T(0);
j = (*i).begin();
jend = (*i).end();
if (A.is_lower()) {
z[j.column()] += *j * x[j.column()];
++j;
} else
--jend;
for (; j != jend; ++j) {
/* normal side */
z[j.row()] += *j * x[j.column()];
/* symmetric side */
tmp += *j * x[j.row()];
}
if (A.is_upper())
tmp += *j * x[j.row()];
z[i.index()] += tmp;
}
}
template <class Matrix, class VecX, class VecZ>
inline void
mult_shape__(const Matrix& A, const VecX& x, VecZ& z,
symmetric_tag)
{
typedef typename matrix_traits<Matrix>::orientation Orien;
mult_symm__(A, x, z, Orien());
}
//: Multiplication: <tt>z <- A x + y</tt>
//!category: algorithms
//!component: function
//!definition: mtl.h
//!precond: <TT>A.nrows() <= y.size()</TT>
//!precond: <TT>A.nrows() <= z.size()</TT>
//!precond: <TT>A.ncols() <= x.size()</TT>
//!precond: no aliasing in the arguments
//!example: symm_sparse_vec_prod.cc
//!typereqs: <tt>Matrix::value_type</tt>, <tt>VecX::value_type</tt>, <tt>VecY::value_type</tt>, and <tt>VecZ::value_type</tt> must be the same type
//!typereqs: the multiplication operator must be defined for <tt>Matrix::value_type</tt>
//!typereqs: the addition operator must be defined for <tt>Matrix::value_type</tt>
template <class Matrix, class VecX, class VecY, class VecZ>
inline void
mult(const Matrix& A, const VecX& x, const VecY& y, MTL_OUT(VecZ) z_)
MTL_THROW_ASSERTION
{
VecZ& z = const_cast<VecZ&>(z_);
mtl::copy(y, z);
typedef typename matrix_traits<Matrix>::shape Shape;
mult_shape__(A, x, z, Shape());
}
//: Matrix Vector Multiplication: <tt>y <- A x</tt>
//
// Multiplies matrix A times vector x and stores the result in vector y.
// <p>
// Note: ignore the <tt>oned_tag</tt> parameter and the underscores in
// the name of this function.
//
//!category: algorithms
//!component: function
//!definition: mtl.h
//!example: general_matvec_mult.cc, banded_matvec_mult.cc, symm_matvec_mult.cc
//!precond: <TT>A.nrows() <= y.size()</TT>
//!precond: <TT>A.ncols() <= x.size()</TT>
//!precond: x and y not same vector
//!example: symm_matvec_mult.cc
//!typereqs: <tt>Matrix::value_type</tt>, <tt>VecX::value_type</tt>, and <tt>VecY::value_type</tt> must be the same type
//!typereqs: the multiplication operator must be defined for <tt>Matrix::value_type</tt>
//!typereqs: the addition operator must be defined for <tt>Matrix::value_type</tt>
template <class Matrix, class VecX, class VecY>
inline void
mult_dim__(const Matrix& A, const VecX& x, VecY& y, oned_tag) MTL_THROW_ASSERTION
{
mtl::mult(A, x, mtl::scaled(y, 0), y);
#if 0
typedef typename matrix_traits<Matrix>::shape Shape;
mult_shape__(A, x, y, Shape());
#endif
}
template <class Matrix, class VecX, class VecY>
inline void
mult_add(const Matrix& A, const VecX& x, MTL_OUT(VecY) y_) MTL_THROW_ASSERTION
{
VecY& y = const_cast<VecY&>(y_);
typedef typename matrix_traits<Matrix>::shape Shape;
mult_shape__(A, x, y, Shape());
}
//: simple 3 loop version of matmat mult
//!noindex:
template <class MatA, class MatB, class MatC, class Orien>
inline void
simple_mult(const MatA& A, const MatB& B, MatC& C, dense_tag, Orien)
{
typedef typename matrix_traits<MatA>::size_type Int;
typename MatA::const_iterator A_k;
typename MatA::OneD::const_iterator A_ki;
A_k = A.begin();
while (not_at(A_k, A.end())) {
for (Int j = 0; j < B.ncols(); ++j) {
A_ki = (*A_k).begin();
while (not_at(A_ki, (*A_k).end())) {
Int k = A_ki.column();
Int i = A_ki.row();
C(i,j) += *A_ki * B(k,j);
++A_ki;
}
}
++A_k;
}
}
/* Assumes A and B are also row oriented */
template <class MatrixA, class MatrixB, class MatrixC>
inline void
simple_mult(const MatrixA& A, const MatrixB& B, MatrixC& C,
sparse_tag, row_tag)
{
typedef typename matrix_traits<MatrixA>::value_type T;
typedef typename matrix_traits<MatrixA>::size_type Int;
T scal;
Int len = 0;
Int jj, k;
Int nzmax = C.capacity();
Int M = A.nrows();
Int N = B.ncols();
dense1D<Int> ic(M + 1, 0);
dense1D<Int> jc(nzmax);
dense1D<T> c(nzmax);
typedef typename dense1D<Int>::iterator di_iter;
typedef typename dense1D<T>::iterator dt_iter;
compressed1D<T> tmp1(N), tmp2(N), tmp3(N);
tmp1.reserve(N);
tmp2.reserve(N);
typedef typename compressed1D<T>::iterator tmpiter;
typename MatrixA::const_iterator Ai;
typename MatrixA::Row::const_iterator Aij;
for (Ai = A.begin(); Ai != A.end(); ++Ai) {
copy(C[Ai.index()], tmp1);
for (Aij = (*Ai).begin(); Aij != (*Ai).end(); ++Aij) {
scal = *Aij;
jj = Aij.column();
// add B[jj] and tmp1 into tmp2
add(mtl::scaled(B[jj], scal), tmp1, tmp2);
tmp1.clear();
// swap tmp1 and tmp2
tmp3 = tmp1; tmp1 = tmp2; tmp2 = tmp3;
}
// copy tmp1 into C[ii]
k = len;
if (k + tmp1.nnz() > nzmax) {
std::cerr << "Not enough work space, increase capacity of C" << std::endl;
return;
}
for (tmpiter t = tmp1.begin(); t != tmp1.end(); ++t, ++k) {
c[k] = *t;
jc[k] = t.index();
}
len += tmp1.nnz();
ic[Ai.index() + 1] = len;
}
typedef typename matrix<T, rectangle<>,
compressed<Int, external>,
row_major>::type SpMat;
SpMat CC(M, N, len, c.data(), ic.data(), jc.data());
copy(CC, C);
}
/* Assumes A and B are also column oriented */
template <class MatrixA, class MatrixB, class MatrixC>
inline void
simple_mult(const MatrixA& A, const MatrixB& B, MatrixC& C,
sparse_tag, column_tag)
{
typedef typename matrix_traits<MatrixA>::value_type T;
typedef typename matrix_traits<MatrixA>::size_type Int;
T scal;
Int len = 0;
Int kk, k;
Int nzmax = C.capacity();
Int M = A.nrows();
Int N = B.ncols();
dense1D<Int> ic(N + 1, 0);
dense1D<Int> jc(nzmax);
dense1D<T> c(nzmax);
typedef typename dense1D<Int>::iterator di_iter;
typedef typename dense1D<T>::iterator dt_iter;
compressed1D<T> tmp1(M), tmp2(M), tmp3(M);
tmp1.reserve(M);
tmp2.reserve(M);
typedef typename compressed1D<T>::iterator tmpiter;
typename MatrixB::const_iterator Bj;
typename MatrixB::Column::const_iterator Bjk;
for (Bj = B.begin(); Bj != B.end(); ++Bj) {
copy(C[Bj.index()], tmp1);
for (Bjk = (*Bj).begin(); Bjk != (*Bj).end(); ++Bjk) {
scal = *Bjk;
kk = Bjk.row();
// add A[kk] and tmp1 into tmp2
add(mtl::scaled(A[kk], scal), tmp1, tmp2);
tmp1.clear();
// swap tmp1 and tmp2
tmp3 = tmp1; tmp1 = tmp2; tmp2 = tmp3;
}
// copy tmp1 into C[ii]
k = len;
if (k + tmp1.nnz() > nzmax) {
std::cerr << "Not enough work space, increase capacity of C" << std::endl;
return;
}
for (tmpiter t = tmp1.begin(); t != tmp1.end(); ++t, ++k) {
c[k] = *t;
jc[k] = t.index();
}
len += tmp1.nnz();
ic[Bj.index() + 1] = len;
}
typedef typename matrix<T, rectangle<>,
compressed<Int, external>,
column_major>::type SpMat;
SpMat CC(M, N, len, c.data(), ic.data(), jc.data());
copy(CC, C);
}
//: Symmetric version, row-major
//!noindex:
template <class MatA, class MatB, class MatC>
inline void
symm_simple_mult(const MatA& A, const MatB& B, MatC& C, row_tag)
{
typedef typename matrix_traits<MatA>::size_type Int;
typename MatA::const_iterator A_k;
typename MatA::OneD::const_iterator A_ki, A_kiend;
A_k = A.begin();
while (not_at(A_k, A.end())) {
for (Int j = 0; j < B.ncols(); ++j) {
A_ki = (*A_k).begin();
A_kiend = (*A_k).end();
Int k = A_ki.column();
Int i = A_ki.row();
if (A.is_upper()) { /* handle the diagonal elements */
C(i,j) += *A_ki * B(k,j);
++A_ki;
} else
--A_kiend;
while (not_at(A_ki, A_kiend)) {
k = A_ki.column();
i = A_ki.row();
C(i,j) += *A_ki * B(k,j);
C(k,j) += *A_ki * B(i,j);
++A_ki;
}
k = A_ki.column();
i = A_ki.row();
if (A.is_lower())
C(i,j) += *A_ki * B(k,j);
}
++A_k;
}
}
//: Symmetric version, column-major
//!noindex:
template <class MatA, class MatB, class MatC>
inline void
symm_simple_mult(const MatA& A, const MatB& B, MatC& C, column_tag)
{
typedef typename matrix_traits<MatA>::size_type Int;
typename MatA::const_iterator A_k;
typename MatA::OneD::const_iterator A_ki, A_kiend;
A_k = A.begin();
while (not_at(A_k, A.end())) {
for (Int j = 0; j < B.ncols(); ++j) {
A_ki = (*A_k).begin();
A_kiend = (*A_k).end();
Int k = A_ki.column();
Int i = A_ki.row();
if (A.is_lower()) { /* handle the diagonal elements */
C(i,j) += *A_ki * B(k,j);
++A_ki;
} else
--A_kiend;
while (not_at(A_ki, A_kiend)) {
k = A_ki.column();
i = A_ki.row();
C(i,j) += *A_ki * B(k,j);
C(k,j) += *A_ki * B(i,j);
++A_ki;
}
k = A_ki.column();
i = A_ki.row();
if (A.is_upper())
C(i,j) += *A_ki * B(k,j);
}
++A_k;
}
}
//: Specialization for triangular matrices
//!noindex:
template <class MatA, class MatB, class MatC>
inline void
matmat_mult(const MatA& A, const MatB& B, MatC& C, symmetric_tag)
{
typedef typename matrix_traits<MatA>::orientation Orien;
symm_simple_mult(A, B, C, Orien());
}
//: Specialization for triangular matrices
//!noindex
template <class MatA, class MatB, class MatC>
inline void
matmat_mult(const MatA& A, const MatB& B, MatC& C, triangle_tag)
{
typedef typename matrix_traits<MatA>::size_type Int;
typedef typename matrix_traits<MatA>::orientation Orien;
if (A.is_unit()) {
Int M = MTL_MIN(A.nrows(), A.ncols());
Int N = B.ncols();
for (Int i = 0; i < M; ++i)
for (Int j = 0; j < N; ++j)
C(i,j) += B(i,j);
}
simple_mult(A, B, C, mtl::dense_tag(), Orien());
}
//: Dispatch to row/column general and banded matrices
//!noindex:
template <class MatA, class MatB, class MatC>
inline void
matmat_mult(const MatA& A, const MatB& B, MatC& C, rectangle_tag)
{
typedef typename matrix_traits<MatA>::sparsity Sparsity;
typedef typename matrix_traits<MatA>::orientation Orien;
simple_mult(A, B, C, Sparsity(), Orien());
}
template <class MatA, class MatB, class MatC>
inline void
matmat_mult(const MatA& A, const MatB& B, MatC& C, banded_tag)
{
typedef typename matrix_traits<MatC>::sparsity Sparsity;
typedef typename matrix_traits<MatA>::orientation Orien;
simple_mult(A, B, C, Sparsity(), Orien());
}
//: Matrix multiplication C <- C + A * B
//
// The actual specialization of the algorithm used depends of the
// types of matrices used. If all the matrices are dense and
// rectangular the blocked algorithm is used (when --with-blais is
// specified in the configure). Otherwise the traversal depends on
// matrix A. Therefore if one is multiplying a sparse matrix by a
// dense, one would want the sparse matrix as the A
// argument. Typically, for performance reasons, one would not want
// to use a sparse matrix for C.
// <p>
// Note: ignore the <tt>twod_tag</tt> argument and the underscores in
// the name of this function.
//
//!precond: <tt>A.nrows() == C.nrows()</tt>
//!precond: <tt>A.ncols() == B.nrows()</tt>
//!precond: <tt>B.ncols() == C.ncols()</tt>
//!category: algorithms
//!component: function
//!definition: mtl.h
//!typereqs: the value types for each of the matrices must be compatible
//!typereqs: the multiplication operator must be defined for <tt>MatA::value_type</tt>
//!typereqs: the addition operator must be defined for <tt>MatA::value_type</tt>
template <class MatA, class MatB, class MatC>
inline void
mult_dim__(const MatA& A, const MatB& B, MatC& C, twod_tag)
{
typedef typename MatA::shape Shape;
matmat_mult(A, B, C, Shape());
}
//: Dispatch between matrix matrix and matrix vector mult.
//!noindex:
template <class LinalgA, class LinalgB, class LinalgC>
inline void
mult(const LinalgA& A, const LinalgB& B, MTL_OUT(LinalgC) C_)
{
LinalgC& C = const_cast<LinalgC&>(C_);
typedef typename linalg_traits<LinalgB>::dimension Dim;
mult_dim__(A, B, C, Dim());
}
//: for column oriented
//!noindex:
template <class TriMatrix, class VecX>
inline void
tri_solve__(const TriMatrix& T, VecX& x, column_tag)
{
typedef typename matrix_traits<TriMatrix>::size_type Int;
typedef typename matrix_traits<TriMatrix>::value_type VT;
typename VecX::value_type x_j;
if (T.is_upper()) {
typename TriMatrix::const_reverse_iterator T_j;
typename TriMatrix::Column::const_reverse_iterator T_ji, T_jrend;
for (T_j = T.rbegin(); T_j != T.rend(); ++T_j) {
T_ji = (*T_j).rbegin();
T_jrend = (*T_j).rend();
//Int j = T_ji.column();
Int j = T_j.index();
//Paul C. Leopardi <[email protected]> reported the fix
//for for a sparse matrix to have a completely empty row (or column)
if ( (T_ji != T_jrend) && ! T.is_unit()) {
x[j] /= *T_ji; /* the diagonal */
++T_ji;
}
x_j = x[j];
while (T_ji != T_jrend) {
Int i = T_ji.row();
x[i] -= x_j * *T_ji;
++T_ji;
}
}
} else { /* T is lower */
typename TriMatrix::const_iterator T_j;
typename TriMatrix::Column::const_iterator T_ji, T_jend;
for (T_j = T.begin(); T_j != T.end(); ++T_j) {
T_ji = (*T_j).begin();
T_jend = (*T_j).end();
//Int j = T_ji.column(); //T_ji could be T_jend
Int j = T_j.index();
if ( (T_ji != T_jend) && ! T.is_unit()) {
x[j] /= *T_ji; /* the diagonal */
++T_ji;
}
x_j = x[j];
while (T_ji != T_jend) {
Int i = T_ji.row();
x[i] -= x_j * *T_ji;
++T_ji;
}
}
}
}
//: for row major
//!noindex:
template <class TriMatrix, class VecX>
inline void
tri_solve__(const TriMatrix& T, VecX& x, row_tag)
{
typedef typename matrix_traits<TriMatrix>::value_type VT;
typedef typename matrix_traits<TriMatrix>::size_type Int;
if (T.is_upper()) {
typename TriMatrix::const_reverse_iterator T_i, T_iend;
typename TriMatrix::Row::const_reverse_iterator T_ij;
T_i = T.rbegin();
T_iend = T.rend();
if ( (T_i != T_iend) && ! T.is_unit()) {
T_ij = (*T_i).rbegin();
x[T_ij.row()] /= *T_ij;
++T_i;
}
while (T_i != T_iend) {
T_ij = (*T_i).rbegin();
//Int i = T_ij.row();
Int i = T_i.index();
VT t = x[i];
typename TriMatrix::Row::const_reverse_iterator T_iend;
T_iend = (*T_i).rend();
if ( (T_ij != T_iend) && ! T.is_unit())
--T_iend;
Int j;
while (T_ij != T_iend) {
j = T_ij.column();
t -= (*T_ij) * x[j];
++T_ij;
}
if ( (*T_i).rbegin() != (*T_i).rend() && !T.is_unit()) //T_i is not empty
t /= *T_ij;
x[i] = t;
++T_i;
}
} else { /* T is lower */
typename TriMatrix::const_iterator T_i;
typename TriMatrix::Row::const_iterator T_ij;
T_i = T.begin();
if (T_i != T.end() && ! T.is_unit()) {
T_ij = (*T_i).begin();
x[T_ij.row()] *= VT(1) / *T_ij;
++T_i;
}
while (T_i != T.end()) {
T_ij = (*T_i).begin();
//Int i = T_ij.row(); //T_ij could be bad
Int i = T_i.index();
VT t = x[i];
typename TriMatrix::Row::const_iterator T_iend;
T_iend = (*T_i).end();
if ( ( T_ij != T_iend ) && ! T.is_unit())
--T_iend;
Int j;
while (T_ij != T_iend) {
j = T_ij.column();
t -= (*T_ij) * x[j];
++T_ij;
}
if ( (*T_i).begin() !=(*T_i).end() && !T.is_unit())
t /= *T_ij;
x[i] = t;
++T_i;
}
}
}
//: Triangular Solve: <tt>x <- T^{-1} * x</tt>
// Use with trianguler matrixes only ie. use the <TT>triangle</TT>
// adaptor class.
//
// To use with a sparse matrix, the sparse matrix must be wrapped with
// a triangle adaptor. You must specify "packed" in the triangle
// adaptor. The sparse matrix must only have elements in the correct
// side.
//
//!category: algorithms
//!component: function
//!definition: mtl.h
//!example: tri_solve.cc
//!typereqs: <tt>Matrix::value_type</tt> and <tt>VecX::value_type</tt> must be the same type
//!typereqs: the multiplication operator must be defined for <tt>Matrix::value_type</tt>
//!typereqs: the division operator must be defined for <tt>Matrix::value_type</tt>
//!typereqs: the addition operator must be defined for <tt>Matrix::value_type</tt>
template <class TriMatrix, class VecX>
inline void
tri_solve(const TriMatrix& T, MTL_OUT(VecX) x_) MTL_THROW_ASSERTION
{
VecX& x = const_cast<VecX&>(x_);
MTL_ASSERT(T.nrows() <= x.size(), "mtl::tri_solve()");
MTL_ASSERT(T.ncols() <= x.size(), "mtl::tri_solve()");
MTL_ASSERT(T.ncols() == T.nrows(), "mtl::tri_solve()");
typedef typename TriMatrix::orientation orien;
tri_solve__(T, x, orien());
}
//: tri solve for left side
//!noindex:
template <class MatT, class MatB>
inline void
tri_solve__(const MatT& T, MatB& B, left_side)
{
/* const int M = B.nrows(); */
const int N = B.ncols();
/* unoptimized version */
for (int j = 0; j < B.ncols(); ++j)
mtl::tri_solve(T, columns(B)[j]);
/* JGS need to do an optimized version of this
if (T.is_upper()) {
for (int k = M-1; k > 0; --k) {
if (B(k,j) != 0) {
if (! T.is_unit())
B(k,j) /= T(k,k);
for (int i = 0; i < k; ++i)
B(i,j) -= B(k,j) * T(i,k);
}
}
} else {
for (int j = 0; j < N; ++j)
for (int k = 0; k < M; ++k) {
if (B(k,j) != 0) {
if (! T.is_unit())
B(k,j) /= T(k,k);
for (int i = k; i < M; ++i)
B(i,j) -= B(k,j) * T(i,k);
}
}
}
*/
}
/* JGS untested!!! */
//: tri solve for right side
//!noindex:
template <class MatT, class MatB>
inline void
tri_solve__(const MatT& T, MatB& B, right_side)
{
const int M = B.nrows();
const int N = B.ncols();
typedef typename MatT::PR PR;
if (T.is_upper()) {
for (int j = 0; j < N; ++j) {
for (int k = 0; k < j; ++k)
if (T(k,j) != PR(0))
for (int i = 0; i < M; ++i)
B(i,j) -= T(k,j) * B(i,k);
if (! T.is_unit()) {
PR tmp = PR(1) / T(j,j);
for (int i = 1; i < M; ++i)
B(i,j) = tmp * B(i,j);
}
}
} else { // T is lower
for (int j = N - 1; j > 0; --j) {
for (int k = j; k < N; ++k)
if (T(k,j) != PR(0))
for (int i = 0; i < M; ++i)
B(i,j) -= T(k,j) * B(i,k);
if (! T.is_unit()) {
PR tmp = PR(1) / T(j,j);
for (int i = 1; i < M; ++i)
B(i,j) = tmp * B(i,j);
}
}
}
}
//: Triangular Solve: <tt>B <- A^{-1} * B or B <- B * A^{-1}</tt>
//
// This solves the equation <tt>T*X = B</tt> or <tt>X*T = B</tt> where T
// is an upper or lower triangular matrix, and B is a general
// matrix. The resulting matrix X is written onto matrix B. The first
// equation is solved if <tt>left_side</tt> is specified. The second
// equation is solved if <tt>right_side</tt> is specified.
//
// Currently only works with dense storage format.
//
//!category: algorithms
//!component: function
//!definition: mtl.h
//!complexity: O(n^3)
//!example: matmat_trisolve.cc
//!typereqs: <tt>MatT::value_type</tt> and <tt>MatB::value_type</tt> must be the same type
//!typereqs: the multiplication operator must be defined for <tt>MatT::value_type</tt>
//!typereqs: the division operator must be defined for <tt>MatT::value_type</tt>
//!typereqs: the addition operator must be defined for <tt>MatT::value_type</tt>
template <class MatT, class MatB, class Side>
inline void
tri_solve(const MatT& T, MTL_OUT(MatB) B, Side s)
{
tri_solve__(T, const_cast<MatB&>(B), s);
}
//: Rank One Update: <tt>A <- A + x * y^T</tt>
//
// Also known as the outer product of two vectors.
// <codeblock>
// y = [ 1 2 3 ]
//
// [ 1 ] [ 1 2 3 ]
// x = [ 2 ] [ 2 4 6 ] => A
// [ 3 ] [ 3 6 9 ]
// [ 4 ] [ 4 8 12 ]
// </codeblock>
// <p>
// When using this algorithm with a symmetric matrix, x and y
// must be the same vector, or at least have the same values.
// Otherwise the resulting matrix is not symmetric.
//
//!precond: <TT>A.nrows() <= x.size()</TT>
//!precond: <TT>A.ncols() <= y.size()</TT>
//!precond: A has rectangle shape and is dense
//!category: algorithms
//!component: function
//!definition: mtl.h
//!example: rank_one.cc
//!typereqs: <tt>Matrix::value_type</tt>, <tt>VecX::value_type</tt>, and <tt>VecY::value_type</tt> must be the same type
//!typereqs: the multiplication operator must be defined for <tt>Matrix::value_type</tt>
//!typereqs: the addition operator must be defined for <tt>Matrix::value_type</tt>
template <class Matrix, class VecX, class VecY>
inline void
rank_one_update(MTL_OUT(Matrix) A_,
const VecX& x, const VecY& y) MTL_THROW_ASSERTION
{
Matrix& A = const_cast<Matrix&>(A_);
MTL_ASSERT(A.nrows() <= x.size(), "mtl::rank_one_update()");
MTL_ASSERT(A.ncols() <= y.size(), "mtl::rank_one_update()");
typename Matrix::iterator i;
typename Matrix::OneD::iterator j, jend;
for (i = A.begin(); i != A.end(); ++i) {
j = (*i).begin(); jend = (*i).end();
for (; j != jend; ++j)
*j += x[j.row()] * MTL_CONJ(y[j.column()]);
}
}
/* 1. how will the scaling by alpha work into this
* 2. is my placement of conj() ok with respect
* to both row and column oriented matrices
* 3. Perhaps split this in two, have diff version for complex
*/
//: Rank Two Update: <tt>A <- A + x * y^T + y * x^T</tt>
//
//
//!category: algorithms
//!component: function
//!precond: <TT>A.nrows() == A.ncols()</TT>
//!precond: <TT>A.nrows() == x.size()</TT>
//!precond: <TT>x.size() == y.size()</TT>
//!precond: A has rectangle shape and is dense.
//!definition: mtl.h
//!example: rank_2_symm_sparse.cc
//!typereqs: <tt>Matrix::value_type</tt>, <tt>VecX::value_type</tt>, and <tt>VecY::value_type</tt> must be the same type.
//!typereqs: The multiplication operator must be defined for <tt>Matrix::value_type</tt>.
//!typereqs: The addition operator must be defined for <tt>Matrix::value_type</tt>.
template <class Matrix, class VecX, class VecY>
inline void
rank_two_update(MTL_OUT(Matrix) A_,
const VecX& x, const VecY& y) MTL_THROW_ASSERTION
{
Matrix& A = const_cast<Matrix&>(A_);
MTL_ASSERT(A.nrows() == A.ncols(), "mtl::rank_two_update()");
MTL_ASSERT(A.nrows() <= x.size(), "mtl::rank_two_update()");
MTL_ASSERT(A.nrows() <= y.size(), "mtl::rank_two_update()");
typename Matrix::iterator i;
typename Matrix::OneD::iterator j, jend;
for (i = A.begin(); i != A.end(); ++i) {
j = (*i).begin(); jend = (*i).end();
for (; j != jend; ++j)
*j += x[j.row()] * MTL_CONJ(y[j.column()])
+ y[j.row()] * MTL_CONJ(x[j.column()]);
}
}
template <class VecX, class VecY>
inline void
copy__(const VecX& x, VecY& y, fast::count<0>)
{
mtl_algo::copy(x.begin(), x.end(), y.begin());
}
#if USE_BLAIS
template <class VecX, class VecY, int N>
inline void
copy__(const VecX& x, VecY& y, fast::count<N>)
{
fast::copy(x.begin(), fast::count<N>(), y.begin());
}
#endif
template <class VecX, class VecY>
inline void
oned_copy(const VecX& x, VecY& y, dense_tag, dense_tag) MTL_THROW_ASSERTION
{
MTL_ASSERT(x.size() <= y.size(), "mtl::copy()");
copy__(x, y, dim_n<VecX>::RET());
}
#if 0
/* perform a scatter */
template <class VecX, class VecY>
inline void
oned_copy(const VecX& x, VecY y, sparse_tag, dense_tag) MTL_THROW_ASSERTION
{
typename VecX::const_iterator xi;
for (xi = x.begin(); xi != x.end(); ++xi)
y[xi.index()] = *xi;
}
/* perform a gather JGS, does this really make sense? */
template <class VecX, class VecY>
inline void
oned_copy(const VecX& x, VecY y, dense_tag, sparse_tag) MTL_THROW_ASSERTION
{
typedef typename VecX::value_type T;
typename VecY::iterator yi;
for (yi = y.begin(); yi != y.end(); ++yi)
*yi = x[yi.index()];
}
#else
template <class VecX, class VecY>
inline void
oned_copy(const VecX& x, VecY& y, sparse_tag, dense_tag) MTL_THROW_ASSERTION
{
typedef typename linalg_traits<VecY>::value_type T;
mtl::set_value(y, T(0));
typename VecX::const_iterator xi;
for (xi = x.begin(); xi != x.end(); ++xi)
y[xi.index()] = *xi;
}
//: Scatter <tt>y <- x</tt>
//
// Scatters the elements of the sparse vector x into
// the dense vector y.
//
//!category: algorithms
//!component: function
//!definition: mtl.h
//!complexity: O(n) where n is the size of the sparse vector
template <class VecX, class VecY>
inline void
scatter(const VecX& x, MTL_OUT(VecY) y_) MTL_THROW_ASSERTION
{
VecY& y = const_cast<VecY&>(y_);
typename VecX::const_iterator xi;
for (xi = x.begin(); xi != x.end(); ++xi)
y[xi.index()] = *xi;
}
//: Gather <tt>y <- x</tt>
//
// Gathers the elements of the dense vector x into
// the sparse vector y, based on the non-zero structure of y.
//
//!category: algorithms
//!component: function
//!definition: mtl.h
//!complexity: O(n) where n is the size of the sparse vector
template <class VecX, class VecY>
inline void
gather(const VecX& x, MTL_OUT(VecY) y_) MTL_THROW_ASSERTION
{
VecY& y = const_cast<VecY&>(y_);
typedef typename VecX::value_type T;
typename VecY::iterator yi;
for (yi = y.begin(); yi != y.end(); ++yi)
*yi = x[yi.index()];
}
#endif
template <class VecX, class VecY, class Tag>
inline void
oned_copy(const VecX& x, VecY& y, Tag, sparse_tag) MTL_THROW_ASSERTION
{
MTL_ASSERT(x.size() <= y.size(), "mtl::copy()");
y.clear();
typename VecX::const_iterator i = x.begin(), iend = x.end();
for (; i != iend; ++i)
y.push_back(i.index(), *i);
}
template <class VecX, class VecY>
inline void
copy__(const VecX& x, VecY& y, oned_tag) MTL_THROW_ASSERTION
{
typedef typename linalg_traits<VecX>::sparsity SpX;
typedef typename linalg_traits<VecY>::sparsity SpY;
oned_copy(x, y, SpX(), SpY());
}
template <class MatA, class MatB>
inline void
twod_copy_default(const MatA& A, MatB& B) MTL_THROW_ASSERTION
{
typename MatA::const_iterator i;
typename MatA::OneD::const_iterator j, jend;
for (i = A.begin(); i != A.end(); ++i) {
j = (*i).begin(); jend = (*i).end();
for (; j != jend; ++j)
B(j.row(),j.column()) = *j;
}
}
template <class MatA, class MatB>
inline void
twod_copy(const MatA& A, MatB& B, rectangle_tag) MTL_THROW_ASSERTION
{
twod_copy_default(A, B);
}
template <class MatA, class MatB>
inline void
twod_copy(const MatA& A, MatB& B, banded_tag) MTL_THROW_ASSERTION
{
twod_copy_default(A, B);
}
template <class MatA, class MatB>
inline void
twod_copy(const MatA& A, MatB& B, symmetric_tag) MTL_THROW_ASSERTION
{
typename MatA::const_iterator i;
typename MatA::OneD::const_iterator j, jend;
for (i = A.begin(); i != A.end(); ++i) {
j = (*i).begin(); jend = (*i).end();
for (; j != jend; ++j) {
B(j.row(),j.column()) = *j;
B(j.column(),j.row()) = *j;
}
}
}
template <class MatA, class MatB>
inline void
twod_copy(const MatA& A, MatB& B, triangle_tag) MTL_THROW_ASSERTION
{
typedef typename matrix_traits<MatB>::value_type T;
if (A.is_unit())
set_diagonal(B, T(1));
twod_copy(A, B, rectangle_tag());
}
template <class MatA, class MatB>
inline void
twod_copy__(const MatA& A, MatB& B, dense_tag)
{
typedef typename matrix_traits<MatA>::shape Shape;
twod_copy(A, B, Shape());
}
/*
Sparse matrices have specialized copy functions since
they need to optimize the creation of the non-zero structure.
only good for same orientation!!!
*/
template <class MatA, class MatB>
inline void
twod_copy(const MatA& A, MatB& B, row_tag, row_tag)
{
B.fast_copy(A);
}
template <class MatA, class MatB>
inline void
twod_copy(const MatA& A, MatB& B, column_tag, column_tag)
{
B.fast_copy(A);
}
template <class MatA, class MatB>
inline void
twod_copy(const MatA& A, MatB& B, row_tag, column_tag)
{
twod_copy__(A, B, dense_tag());
}
template <class MatA, class MatB>
inline void
twod_copy(const MatA& A, MatB& B, column_tag, row_tag)
{
twod_copy__(A, B, dense_tag());
}
template <class MatA, class MatB>
inline void
twod_copy__(const MatA& A, MatB& B, sparse_tag)
{
typedef typename matrix_traits<MatA>::orientation OrienA;
typedef typename matrix_traits<MatB>::orientation OrienB;
twod_copy(A, B, OrienA(), OrienB());
}
template <class MatA, class MatB>
inline void
copy__(const MatA& A, MatB& B, twod_tag) MTL_THROW_ASSERTION
{
MTL_ASSERT(A.nrows() <= B.nrows(), "copy(A, B, twod_tag)");
MTL_ASSERT(A.ncols() <= B.ncols(), "copy(A, B, twod_tag)");
typedef typename matrix_traits<MatB>::sparsity Sparsity;
twod_copy__(A, B, Sparsity());
}
//: Copy: <tt>B <- A or y <- x</tt>
//
// Copy the elements of matrix A into matrix B, or copy the elements
// of vector x into vector y. For shaped and sparse matrices, this
// copies only the elements stored in A to B. If x is a sparse
// vector and y is dense, a "scatter" is performed. If y is sparse
// and x is dense, then a "gather" is performed. If both vectors
// are sparse, but of different structure the result is undefined.
//
//!category: algorithms
//!component: function
//!definition: mtl.h
//!complexity: O(m*n) for matrices. O(nnz) if either A or B are sparse and of the same orientation (otherwise it can be O(nnz^2). O(n) for vectors.
//!example: vecvec_copy.cc
template <class LinalgA, class LinalgB>
inline void
copy(const LinalgA& A, MTL_OUT(LinalgB) B_) MTL_THROW_ASSERTION
{
LinalgB& B = const_cast<LinalgB&>(B_);
typedef typename linalg_traits<LinalgA>::dimension Dim;
copy__(A, B, Dim());
}
template <class VecX, class VecY> inline
void
add__(const VecX& x, VecY& y, fast::count<0>)
{
typedef typename VecX::value_type T;
mtl_algo::transform_add(x.begin(), x.end(), y.begin());
}
#if USE_BLAIS
template <class VecX, class VecY, int N> inline
void
add__(const VecX& x, VecY& y, fast::count<N>)
{
typedef typename VecX::value_type T;
fast::transform(x.begin(), fast::count<N>(), y.begin(),
y.begin(), std::plus<T>());
}
#endif
template <class VecX, class VecY> inline
void
add__(const VecX& x, VecY& y, oned_tag) MTL_THROW_ASSERTION
{
MTL_ASSERT(x.size() <= y.size(), "mtl::add()");
add__(x, y, dim_n<VecX>::RET());
}
template <class VecX, class VecY, class VecZ> inline
void
oned_add(const VecX& x, const VecY& y, VecZ& z, fast::count<0>)
{
typedef typename VecX::value_type T;
mtl_algo::transform(x.begin(), x.end(), y.begin(), z.begin(), std::plus<T>());
}
#if USE_BLAIS
template <class VecX, class VecY, class VecZ, int N> inline
void
oned_add(const VecX& x, const VecY& y, VecZ& z, fast::count<N>)
{
typedef typename VecX::value_type T;
fast::transform(x.begin(), fast::count<N>(), y.begin(), z.begin(), std::plus<T>());
}
#endif
template <class VecX, class VecY, class VecZ>
inline void
oned_add(const VecX& x, const VecY& y, VecZ& z_, sparse_tag)
{
typedef typename VecZ::value_type T;
compressed1D<T> z;
typedef typename VecX::const_iterator xiter;
typedef typename VecY::const_iterator yiter;
xiter xi = x.begin();
xiter xiend = x.end();
yiter yi = y.begin();
yiter yiend = y.end();
while (xi != xiend && yi != yiend) {
if (yi.index() < xi.index()) {
z.push_back(yi.index(), *yi);
++yi;
} else if (xi.index() < yi.index()) {
z.push_back(xi.index(), *xi);
++xi;
} else {
z.push_back(xi.index(), *yi + *xi);
++xi; ++yi;
}
}
while (xi != xiend) {
z.push_back(xi.index(), *xi);
++xi;
}
while (yi != yiend) {
z.push_back(yi.index(), *yi);
++yi;
}
z_.clear();
mtl::copy(z, z_);
}
template <class VecX, class VecY, class VecZ>
inline void
oned_add(const VecX& x, const VecY& y, VecZ& z, dense_tag) MTL_THROW_ASSERTION
{
oned_add(x, y, z, dim_n<VecX>::RET());
}
//: Add: <tt>z <- x + y</tt>
//
// Add the elements of x and y and assign into z.
//
//!category: algorithms
//!component: function
//!definition: mtl.h
//!example: y_ax_y.cc, vecvec_add.cc
//!typereqs: <tt>VecX::value_type</tt>, <tt>VecY::value_type</tt>, and <tt>VecZ::value_type</tt> should be the same type
//!typereqs: The addition operator must be defined for the value_type.
//!complexity: linear time
template <class VecX, class VecY, class VecZ>
inline void
add(const VecX& x, const VecY& y, MTL_OUT(VecZ) z_) MTL_THROW_ASSERTION
{
VecZ& z = const_cast<VecZ&>(z_);
MTL_ASSERT(x.size() <= y.size(), "mtl::add()");
MTL_ASSERT(x.size() <= z.size(), "mtl::add()");
typedef typename linalg_traits<VecZ>::sparsity Sparsity;
oned_add(x, y, z, Sparsity());
}
//: Add: <tt>w <- x + y + z</tt>
//
// Add the elements of x, y, and z and assign into w.
// For now just dense vectors.
//
//!category: algorithms
//!component: function
//!definition: mtl.h
//!example: vecvec_add3.cc
//!typereqs: <tt>VecX::value_type</tt>, <tt>VecY::value_type</tt>, <tt>VecZ::value_type</tt>, and <tt>VecW::value_type</tt> should be the same type
//!typereqs: The addition operator must be defined for the value_type.
//!complexity: linear time
template <class VecW, class VecX, class VecY, class VecZ>
inline void
add(const VecX& x, const VecY& y, const VecZ& z, MTL_OUT(VecW) w_)
MTL_THROW_ASSERTION
{
VecW& w = const_cast<VecW&>(w_);
MTL_ASSERT(x.size() <= y.size(), "mtl::add()");
MTL_ASSERT(x.size() <= z.size(), "mtl::add()");
MTL_ASSERT(x.size() <= w.size(), "mtl::add()");
typename VecX::const_iterator x_i = x.begin();
typename VecY::const_iterator y_i = y.begin();
typename VecZ::const_iterator z_i = z.begin();
typename VecW::iterator w_i = w.begin();
while (not_at(x_i, x.end())) {
*w_i = *x_i + *y_i + *z_i;
++x_i; ++y_i; ++z_i; ++w_i;
}
}
template <class MatA, class MatB>
inline void
twod_add_default(const MatA& A, MatB& B)
{
typename MatA::const_iterator i;
typename MatA::OneD::const_iterator j, jend;
for (i = A.begin(); i != A.end(); ++i) {
j = (*i).begin(); jend = (*i).end();
for (; j != jend; ++j)
B(j.row(), j.column()) += *j;
}
}
template <class MatA, class MatB>
inline void
twod_add(const MatA& A, MatB& B, banded_tag)
{
twod_add_default(A, B);
}
template <class MatA, class MatB>
inline void
twod_add(const MatA& A, MatB& B, rectangle_tag)
{
twod_add_default(A, B);
}
template <class MatA, class MatB>
inline void
twod_add(const MatA& A, MatB& B, triangle_tag)
{
typedef typename matrix_traits<MatA>::size_type Int;
typedef typename matrix_traits<MatA>::value_type T;
if (A.is_unit())
for (Int i = 0; i < MTL_MIN(A.nrows(), A.ncols()); ++i)
B(i,i) += T(1);
twod_add(A, B, banded_tag());
}
/* perhaps I should add is_row() and is_column()
methods to the matrices
*/
template <class MatA, class MatB>
inline void
twod_symmetric_add(const MatA& A, MatB& B, row_tag)
{
typename MatA::const_iterator i;
typename MatA::Row::const_iterator j, jend;
for (i = A.begin(); i != A.end(); ++i) {
j = (*i).begin();
jend = (*i).end();
if (A.is_upper()) { /* handle the diagonal elements */
B(j.column(), j.row()) += *j;
++j;
} else
--jend;
for (; j != jend; ++j) {
B(j.row(), j.column()) += *j;
B(j.column(), j.row()) += *j;
}
if (A.is_lower())
B(j.column(), j.row()) += *j;
}
}
template <class MatA, class MatB>
inline void
twod_symmetric_add(const MatA& A, MatB& B, column_tag)
{
typename MatA::const_iterator i;
typename MatA::Column::const_iterator j, jend;
for (i = A.begin(); i != A.end(); ++i) {
j = (*i).begin();
jend = (*i).end();
if (A.is_lower()) { /* handle the diagonal elements */
B(j.column(), j.row()) += *j;
++j;
} else
--jend;
for (; j != jend; ++j) {
B(j.row(), j.column()) += *j;
B(j.column(), j.row()) += *j;
}
if (A.is_upper())
B(j.column(), j.row()) += *j;
}
}
template <class MatA, class MatB>
inline void
twod_add(const MatA& A, MatB& B, symmetric_tag)
{
typedef typename matrix_traits<MatA>::orientation Orien;
twod_symmetric_add(A, B, Orien());
}
template <class MatA, class MatB>
inline void
add__(const MatA& A, MatB& B, twod_tag) MTL_THROW_ASSERTION
{
MTL_ASSERT(A.nrows() <= B.nrows(), "matmat::add()");
MTL_ASSERT(A.ncols() <= B.ncols(), "matmat::add()");
typedef typename matrix_traits<MatA>::shape Shape;
twod_add(A, B, Shape());
}
//: Add: <tt>B <- A + B or y <- x + y</tt>
// The function adds the element of A to B, or the elements of x to y.
//
//!category: algorithms
//!component: function
//!definition: mtl.h
//!complexity: O(m*n) for a dense A, O(nnz) for a sparse A. O(n) for a vector.
template <class LinalgA, class LinalgB>
inline void
add(const LinalgA& A, MTL_OUT(LinalgB) B_) MTL_THROW_ASSERTION
{
LinalgB& B = const_cast<LinalgB&>(B_);
typedef typename linalg_traits<LinalgA>::dimension Dim;
add__(A, B, Dim());
}
template <class VecX, class VecY, class VecZ>
inline void
ele_mult(const VecX& x, const VecY& y, MTL_OUT(VecZ) z_, fast::count<0>)
{
VecZ& z = const_cast<VecZ&>(z_);
typedef typename VecX::value_type T;
mtl_algo::transform(x.begin(), x.end(), y.begin(), z.begin(),
std::multiplies<T>());
}
#if USE_BLAIS
template <class VecX, class VecY, class VecZ, int N>
inline void
ele_mult(const VecX& x, const VecY& y, MTL_OUT(VecZ) z_, fast::count<N>)
{
VecZ& z = const_cast<VecZ&>(z_);
typedef typename VecX::value_type T;
fast::transform(x.begin(), fast::count<N>(), y.begin(), z.begin(),
std::multiplies<T>());
}
#endif
//: Element-wise Multiplication: <tt>z <- x O* y</tt>
//!category: algorithms
//!component: function
//!definition: mtl.h
//!example: vecvec_ele_mult.cc
template <class VecX, class VecY, class VecZ>
inline void
ele_mult(const VecX& x, const VecY& y, MTL_OUT(VecZ) z_) MTL_THROW_ASSERTION
{
VecZ& z = const_cast<VecZ&>(z_);
MTL_ASSERT(x.size() <= y.size(), "mtl::ele_mult()");
MTL_ASSERT(x.size() <= z.size(), "mtl::ele_mult()");
ele_mult(x, y, z, dim_n<VecX>::RET());
}
//: Element-wise Multiply: <tt>B <- A O* B</tt>
//
// This function multiplies each of the elements
// of B by the corresponding element of A.
//
//!category: algorithms
//!component: function
//!definition: mtl.h
//!complexity: O(n^2)
template <class MatA, class MatB>
inline void
ele_mult(const MatA& A, MTL_OUT(MatB) B_) MTL_THROW_ASSERTION
{
MatB& B = const_cast<MatB&>(B_);
/* Note: have to iterator over B, since
* elements of B may get zeroed out,
* but zero elements of B stay zero
*/
//typename MatB::row_2Diterator B_i;
//typename MatB::RowVector::iterator j, jend;
typename MatB::iterator i;
typename MatB::OneD::iterator j, jend;
for (i = B.begin(); i != B.end(); ++i) {
j = (*i).begin(); jend = (*i).end();
for (; j != jend; ++j)
*j *= A(j.row(),j.column());
}
}
//: Element-wise Division: <tt>z <- x O/ y</tt>
//!category: algorithms
//!component: function
//!definition: mtl.h
//!example: vecvec_ele_div.cc
template <class VecX, class VecY, class VecZ>
inline void
ele_div(const VecX& x, const VecY& y, MTL_OUT(VecZ) z_) MTL_THROW_ASSERTION
{
VecZ& z = const_cast<VecZ&>(z_);
MTL_ASSERT(x.size() <= y.size(), "mtl::ele_div()");
MTL_ASSERT(x.size() <= z.size(), "mtl::ele_div()");
typedef typename VecX::value_type T;
mtl_algo::transform(x.begin(), x.end(), y.begin(), z.begin(),
std::divides<T>());
}
template <class VecX, class VecY>
inline void
swap(VecX& x, VecY& y, fast::count<0>)
{
mtl_algo::swap_ranges(x.begin(), x.end(), y.begin());
}
#if USE_BLAIS
template <class VecX, class VecY, int N>
inline void
swap(VecX& x, VecY& y, fast::count<N>)
{
fast::swap_ranges(x.begin(), fast::count<N>(), y.begin());
}
#endif
template <class VecX, class VecY>
inline void
swap(VecX& x, VecY& y, oned_tag) MTL_THROW_ASSERTION
{
MTL_ASSERT(x.size() <= y.size(), "mtl::swap()");
swap(x, y, dim_n<VecX>::RET());
}
template <class MatA, class MatB>
inline void
swap(MatA& A, MatB& B, twod_tag) MTL_THROW_ASSERTION
{
MTL_ASSERT(A.nrows() == B.nrows(), "matmat::swap()");
MTL_ASSERT(A.ncols() == B.ncols(), "matmat::swap()");
typename MatA::iterator A_i;
typename MatA::OneD::iterator A_ij, A_ijend;
typename MatB::iterator B_i;
typename MatB::Row::iterator B_ij;
A_i = A.begin(); B_i = B.begin();
while (A_i != A.end()) {
A_ij = (*A_i).begin(); B_ij = (*B_i).begin();
A_ijend = (*A_i).end();
while (A_ij != A_ijend) {
typename matrix_traits<MatA>::value_type tmp = *B_ij;
*B_ij = *A_ij;
*A_ij = tmp;
++A_ij; ++B_ij;
}
++A_i; ++B_i;
}
}
//: Swap: <tt>B <-> A or y <-> x</tt>
//
// Exchanges the elements of the containers.
// Not compatible with sparse matrices. For banded matrices
// and other shaped matrices, A and B must be the same shape.
// Also, the two matrices must be the same orientation.
//
//!category: algorithms
//!component: function
//!definition: mtl.h
//!complexity: O(n^2)
//!example: vecvec_swap.cc
template <class LinalgA, class LinalgB>
inline void
swap(MTL_OUT(LinalgA) A, MTL_OUT(LinalgB) B) MTL_THROW_ASSERTION
{
typedef typename linalg_traits<LinalgA>::dimension Dim;
swap(const_cast<LinalgA&>(A), const_cast<LinalgB&>(B), Dim());
}
template <class VecX, class VecY, class T>
inline T
dot(const VecX& x, const VecY& y, T s, fast::count<0>)
{
return mtl_algo::inner_product(x.begin(), x.end(), y.begin(), s);
}
#if USE_BLAIS
template <class VecX, class VecY, class T, int N>
inline T
dot(const VecX& x, const VecY& y, T s, fast::count<N>)
{
return fast::inner_product(x.begin(), fast::count<N>(), y.begin(), s);
}
#endif
template <class VecX, class VecY, class T>
inline T
dot(const VecX& x, const VecY& y, T s, dense_tag, dense_tag)
{
return dot(x, y, s, dim_n<VecX>::RET());
}
template <class InputIterator1, class InputIterator2, class T>
inline T
sparse_inner_product(InputIterator1 f1, InputIterator1 l1,
InputIterator2 f2, InputIterator2 l2, T init)
{
InputIterator1 first1 = f1;
InputIterator1 last1 = l1;
InputIterator2 first2 = f2;
InputIterator2 last2 = l2;
while (first1 != last1 && first2 != last2) {
if (first1.index() == first2.index())
init += (*first1++ * *first2++);
else if (first1.index() < first2.index())
++first1;
else
++first2;
}
return init;
}
template <class IndexedIterator, class RandomAccessIterator, class T>
inline T
sparse_dense_inner_product(IndexedIterator f1, IndexedIterator l1,
RandomAccessIterator f2, T init)
{
IndexedIterator first1 = f1, last1 = l1;
RandomAccessIterator first2 = f2;
while (first1 != last1) {
init += (*first1 * first2[first1.index()]);
++first1;
}
return init;
}
template <class VecX, class VecY, class T>
inline T
dot(const VecX& x, const VecY& y, T s, sparse_tag, sparse_tag)
{
if (x.nnz() < y.nnz())
return sparse_inner_product(x.begin(), x.end(), y.begin(), y.end(), s);
else
return sparse_inner_product(y.begin(), y.end(), x.begin(), x.end(), s);
}
template <class VecX, class VecY, class T>
inline T
dot(const VecX& x, const VecY& y, T s, dense_tag, sparse_tag)
{
return sparse_dense_inner_product(y.begin(), y.end(), x.begin(), s);
}
template <class VecX, class VecY, class T>
inline T
dot(const VecX& x, const VecY& y, T s, sparse_tag, dense_tag)
{
return sparse_dense_inner_product(x.begin(), x.end(), y.begin(), s);
}
//: Dot Product: <tt>s <- x . y + s</tt>
// The type used for argument s determines the
// type of the resulting product.
//!category: algorithms
//!component: function
//!definition: mtl.h
template <class VecX, class VecY, class T>
inline T
dot(const VecX& x, const VecY& y, T s) MTL_THROW_ASSERTION
{
MTL_ASSERT(x.size() == y.size(), "mtl::dot()");
typedef typename linalg_traits<VecX>::sparsity SparseX;
typedef typename linalg_traits<VecY>::sparsity SparseY;
return dot(x, y, s, SparseX(), SparseY());
}
//: Dot Product: <tt>s <- x . y</tt>
// The type of the resulting product is <TT>VecX::value_type</TT>.
//!category: algorithms
//!component: function
//!example: vecvec_dot.cc, dot_prod.cc
//!definition: mtl.h
template <class VecX, class VecY>
inline typename VecX::value_type
dot(const VecX& x, const VecY& y) MTL_THROW_ASSERTION
{
typedef typename VecX::value_type T;
return mtl::dot(x, y, T(0));
}
#ifdef USE_DOUBLE_DOUBLE
//: Dot Product (extended precision): <tt>s <- x . y + s</tt>
// The type of the resulting product is double_double
// Extended precision is used internally.
//!category: algorithms
//!component: function
//!definition: mtl.h
template <class VecX, class VecY>
inline double_double
dot(const VecX& x, const VecY& y, double_double s) MTL_THROW_ASSERTION
{
typedef typename VecX::value_type x_type;
typedef typename VecY::value_type y_type;
// x_type and y_type must be either float or double
multiply<double_double, x_type, y_type> m;
addition<double_double, double_double, double_double> a;
return mtl_algo::inner_product(x.begin(), x.end(), y.begin(), s, a, m);
}
#endif /* USE_DOUBLE_DOUBLE */
template <class T>
struct conj_func {
typedef T result_type;
inline T operator()(const T& x) const { return MTL_CONJ(x); }
};
template <class VecX, class VecY, class T>
inline T
dot_conj(const VecX& x, const VecY& y, T s, fast::count<0>)
{
return mtl_algo::inner_product(x.begin(), x.end(),
trans_iter(y.begin(), conj_func<T>()), s);
}
#if USE_BLAIS
template <class VecX, class VecY, class T, int N>
inline T
dot_conj(const VecX& x, const VecY& y, T s, fast::count<N>)
{
return fast::inner_product(x.begin(), x.end(),
trans_iter(y.begin(), conj_func<T>()), s);
}
#endif
//: Dot Conjugate: <tt>s <- x . conj(y) + s</tt>
// Similar to dot product. The complex conjugate of the elements of y
// is used. For real numbers, the conjugate is just that real number.
// Note that the type of parameter s is the return type of this
// function.
//!category: algorithms
//!component: function
//!definition: mtl.h
template <class VecX, class VecY, class T>
inline T
dot_conj(const VecX& x, const VecY& y, T s) MTL_THROW_ASSERTION
{
MTL_ASSERT(x.size() <= y.size(), "mtl::dot_conj()");
return dot_conj(x, y, s, dim_n<VecX>::RET());
}
//: Dot Conjugate: <tt>s <- x . conj(y)</tt>
// A slightly simpler version of the dot conjugate.
// The return type is the element type of vector x.
//!category: algorithms
//!component: function
//!definition: mtl.h
template <class VecX, class VecY>
inline typename VecX::value_type
dot_conj(const VecX& x, const VecY& y) MTL_THROW_ASSERTION
{
typedef typename VecX::value_type T;
return mtl::dot_conj(x, y, T(0));
}
} /* namespace mtl */
#endif /* _MTL_MTL_H_ */
dense2D.h
(application/octet-stream, 44.8 KB)
// -*- c++ -*-
//
// Copyright 1997, 1998, 1999 University of Notre Dame.
// Authors: Andrew Lumsdaine, Jeremy G. Siek, Lie-Quan Lee
//
// This file is part of the Matrix Template Library
//
// You should have received a copy of the License Agreement for the
// Matrix Template Library along with the software; see the
// file LICENSE. If not, contact Office of Research, University of Notre
// Dame, Notre Dame, IN 46556.
//
// Permission to modify the code and to distribute modified code is
// granted, provided the text of this NOTICE is retained, a notice that
// the code was modified is included with the above COPYRIGHT NOTICE and
// with the COPYRIGHT NOTICE in the LICENSE file, and that the LICENSE
// file is distributed with the modified code.
//
// LICENSOR MAKES NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED.
// By way of example, but not limitation, Licensor MAKES NO
// REPRESENTATIONS OR WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY
// PARTICULAR PURPOSE OR THAT THE USE OF THE LICENSED SOFTWARE COMPONENTS
// OR DOCUMENTATION WILL NOT INFRINGE ANY PATENTS, COPYRIGHTS, TRADEMARKS
// OR OTHER RIGHTS.
//
//===========================================================================
#ifndef MTL_DENSE2D_H
#define MTL_DENSE2D_H
#include "mtl/mtl_iterator.h"
#include <utility>
#include <assert.h>
#include <vector>
#include "mtl/mtl_config.h"
#include "mtl/linalg_vec.h"
#include "mtl/strided1D.h"
#include "mtl/initialize.h"
#include "mtl/reverse_iter.h"
#include "mtl/matrix_traits.h"
#include "mtl/dimension.h"
#ifndef MTL_DISABLE_BLOCKING
#include "mtl/block2D.h"
#endif
namespace mtl {
template <class size_t, int MM, int NN>
class strided_offset;
struct strided_tag { enum { id = 1 }; };
struct not_strided_tag { enum { id = 0 }; };
template <class size_t, int MM, int NN>
class band_view_offset;
template <class size_t, int MM, int NN>
class strided_band_view_offset;
template <int M, int N> struct gen_rect_offset;
template <int M, int N> struct gen_strided_offset;
template <int M, int N> struct gen_banded_offset;
template <int M, int N> struct gen_banded_view_offset;
template <int M, int N> struct gen_strided_band_view_offset;
template <int M, int N> struct gen_packed_offset;
//: Rectangular Offset Class
//!models: Offset
//!category: utilities
//!component: type
template <class size_t, int MM, int NN>
class rect_offset {
public:
#if !defined(_MSVCPP_)
template <class Vec>
struct bind_oned {
typedef Vec type;
};
#endif
typedef not_strided_tag is_strided;
typedef size_t size_type;
enum { M = MM, N = NN, IS_STRIDED = 0 };
typedef dimension<size_type, MM, NN> dim_type;
typedef dimension<int> band_type;
typedef strided_offset<size_type, MM, NN> transpose_type;
typedef strideable strideability;
// VC++ doesn't like this
//friend class transpose_type;
//what is that for? -- llee
//inline rect_offset() : dim(4444,4444), ld(4444) { }
inline rect_offset() : dim(0,0), ld(0) { }
inline rect_offset(const rect_offset& x) : dim(x.dim), ld(x.ld) { }
inline rect_offset(size_type m, size_type n, size_type ld_)
: dim(m, n), ld(ld_) { }
inline rect_offset(size_type m, size_type n, size_type ld_, band_type)
: dim(m, n), ld(ld_) { }
rect_offset(const transpose_type& x); /* see below strided_offset for def */
inline rect_offset& operator=(const rect_offset& x) {
dim = x.dim; ld = x.ld; return *this;
}
inline size_type elt(size_type i, size_type j) const { return i * ld + j; }
inline size_type oned_offset(size_type i) const { return i * ld; }
inline size_type oned_length(size_type) const { return dim.second(); }
inline size_type twod_length() const { return dim.first(); }
inline size_type stride() const { return 1; }
inline static size_type size(size_type m, size_type n,
size_type , size_type) { return m * n; }
inline size_type major() const { return dim.first(); }
inline size_type minor() const { return dim.second(); }
/* private: */
dim_type dim;
size_type ld;
};
//: blah
//!noindex:
template <int M, int N>
struct gen_rect_offset {
#if defined( _MSVCPP_ )
typedef rect_offset<unsigned int, M, N> type;
#else
template <class size_type>
struct bind {
typedef rect_offset<size_type, M, N> type;
};
#endif
typedef gen_strided_offset<M,N> transpose_type;
typedef gen_banded_view_offset<M,N> banded_view_type;
};
//: Strided Rectangular Offset Class
//!models: Offset
//!category: utilities
//!component: type
template <class size_t, int MM, int NN>
class strided_offset {
public:
#if !defined(_MSVCPP_)
template <class Vec>
struct bind_oned {
typedef strided1D<Vec> type;
};
#endif
/* typedef strided_band_view_offset<size_t,MM,NN> banded_view_type;
*/
typedef strided_tag is_strided;
enum { M = MM, N = NN, IS_STRIDED = 1 };
typedef size_t size_type;
typedef dimension<size_type, MM, NN> dim_type;
typedef dimension<int> band_type;
typedef rect_offset<size_type,MM,NN> transpose_type;
typedef strideable strideability;
// VC++ doesn't like this
//friend class transpose_type;
inline strided_offset() : dim(0,0), ld(0) { }
inline strided_offset(size_type m, size_type n, size_type ld_)
: dim(m, n), ld(ld_) { }
inline strided_offset(const transpose_type& x) : dim(x.dim), ld(x.ld) { }
inline strided_offset& operator=(const strided_offset& x) {
dim = x.dim; ld = x.ld; return *this;
}
inline size_type elt(size_type i, size_type j) const { return j * ld + i; }
inline size_type oned_offset(size_type i) const { return i; }
inline size_type oned_length(size_type) const { return dim.first() * ld; }
inline size_type twod_length() const { return dim.second(); }
inline size_type stride() const { return ld; }
inline static size_type size(size_type m, size_type n,
size_type , size_type) { return m * n; }
inline size_type major() const { return dim.first(); }
inline size_type minor() const { return dim.second(); }
/* private: */
dim_type dim;
size_type ld;
};
//: blah
//!noindex:
template <int M, int N>
struct gen_strided_offset {
#if defined( _MSVCPP_ )
typedef strided_offset<unsigned int, M, N> type;
#else
template <class size_type>
struct bind {
typedef strided_offset<size_type, M, N> type;
};
#endif
typedef gen_rect_offset<M,N> transpose_type;
typedef gen_strided_band_view_offset<M,N> banded_view_type;
};
/*
template <class size_t, int MM, int NN>
inline rect_offset<size_t,MM,NN>::rect_offset(const rect_offset<size_t,MM,NN>::transpose_type& x)
: dim(x.dim), ld(x.ld) { }
*/
//: Banded View Offset Class
// This creates a banded view into a full matrix.
//!models: Offset
//!category: utilities
//!component: type
template <class size_t, int MM, int NN>
class banded_view_offset {
public:
#if !defined(_MSVCPP_)
template <class Vec>
struct bind_oned {
typedef Vec type;
};
#endif
typedef not_strided_tag is_strided;
enum { M = MM, N = NN, IS_STRIDED = 0 };
typedef size_t size_type;
typedef dimension<size_type, MM, NN> dim_type;
typedef dimension<int, MM, NN> band_type;
typedef strided_band_view_offset<size_type, MM, NN> transpose_type;
typedef not_strideable strideability;
// VC++ doesn't like this
//friend class transpose_type;
inline banded_view_offset()
: dim(0,0), ld(0), bw(std::make_pair(0,0)) { }
inline banded_view_offset(size_type m, size_type n, size_type leading_dim,
band_type band)
: dim(m, n), ld(leading_dim), bw(band) { }
inline banded_view_offset(size_type m, size_type n, size_type leading_dim)
: dim(m, n), ld(leading_dim), bw(band_type(0,0)) { }
template <class Offset>
inline banded_view_offset(Offset os, band_type band)
: dim(os.dim), ld(os.ld), bw(band) { }
inline banded_view_offset& operator=(const banded_view_offset& x) {
dim = x.dim; ld = x.ld; bw = x.bw; return *this;
}
inline size_type elt(size_type i, size_type j) const {
size_type start = MTL_MAX(int(i) - bw.first(), 0);
return i * ld + j + start;
}
inline size_type oned_offset(size_type i) const {
size_type start = MTL_MAX(int(i) - bw.first(), 0);
return i * ld + start;
}
inline size_type oned_length(size_type i) const {
return MTL_MAX(0, MTL_MIN(int(dim.second()), int(i) + bw.second() + 1)
- MTL_MAX(0, int(i) - bw.first()));
}
inline size_type twod_length() const { return dim.first(); }
inline int stride() const { return 1; }
inline static size_type size(size_type m, size_type n,
size_type , size_type) {
return m * n;
}
inline size_type major() const { return dim.first(); }
inline size_type minor() const { return dim.second(); }
/* private: */
dim_type dim;
size_type ld;
band_type bw; /* bandwidth */
};
//: blah
//!noindex:
template <int M, int N>
struct gen_banded_view_offset {
#if defined( _MSVCPP_ )
typedef banded_view_offset<unsigned int, M, N> type;
#else
template <class size_type>
struct bind {
typedef banded_view_offset<size_type, M, N> type;
};
#endif
typedef gen_strided_band_view_offset<M,N> transpose_type;
typedef gen_banded_view_offset<M,N> banded_view_type; /* bogus */
};
//: Strided Band View Offset Class
//
// This creates a strided band view into a full matrix.
// This class is to banded_view as strided_offset is to rect_offset.
//
//!models: Offset
//!category: utilities
//!component: type
template <class size_t, int MM, int NN>
class strided_band_view_offset {
public:
#if !defined(_MSVCPP_)
template <class Vec>
struct bind_oned {
typedef strided1D<Vec> type;
};
#endif
typedef strided_tag is_strided;
enum { M = MM, N = NN, IS_STRIDED = 1 };
typedef size_t size_type;
typedef dimension<size_type, MM, NN> dim_type;
typedef dimension<int, MM, NN> band_type;
typedef banded_view_offset<size_type, MM, NN> transpose_type;
typedef not_strideable strideability;
// VC++ doesn't like this
//friend class transpose_type;
inline strided_band_view_offset()
: dim(0,0), ld(0), bw(std::make_pair(0,0)) { }
inline strided_band_view_offset(size_type m, size_type n,
size_type leading_dim,
band_type band)
: dim(m, n), ld(leading_dim), bw(band) { }
template <class Offset>
inline strided_band_view_offset(Offset os, band_type band)
: dim(os.dim), ld(os.ld), bw(band) { }
inline strided_band_view_offset&
operator=(const strided_band_view_offset& x) {
dim = x.dim; ld = x.ld; bw = x.bw; return *this;
}
inline size_type elt(size_type i, size_type j) const {
size_type start = MTL_MAX(int(i) - bw.first(), 0);
return (j + start) * ld + i;
}
inline size_type oned_offset(size_type i) const {
size_type start = MTL_MAX(int(i) - bw.first(), 0);
return start * ld + i;
}
inline size_type oned_length(size_type i) const {
/* use dim.first() here */
size_type len = MTL_MAX(0, MTL_MIN(int(dim.first()), int(i) + bw.second() + 1)
- MTL_MAX(0, int(i) - bw.first()));
return len * ld;
}
inline size_type twod_length() const { return dim.second(); }
inline int stride() const { return ld; }
inline static size_type size(size_type m, size_type n,
size_type , size_type) {
return m * n;
}
inline size_type major() const { return dim.first(); }
inline size_type minor() const { return dim.second(); }
/* private: */
dim_type dim;
size_type ld;
band_type bw; /* bandwidth */
};
//: blah
//!noindex:
template <int M, int N>
struct gen_strided_band_view_offset {
#if defined( _MSVCPP_ )
typedef strided_band_view_offset<unsigned int, M, N> type;
#else
template <class size_type>
struct bind {
typedef strided_band_view_offset<size_type, M, N> type;
};
#endif
typedef gen_banded_view_offset<M,N> transpose_type;
typedef gen_strided_band_view_offset<M,N> banded_view_type; /* bogus */
};
template <class size_t, int MM, int NN>
class packed_offset;
//: Banded Offset Class
// This cooresponds to lapack/blas banded storage format.
//!models: Offset
//!category: utilities
//!component: type
template <class size_t, int MM, int NN>
class banded_offset {
public:
#if !defined(_MSVCPP_)
template <class Vec>
struct bind_oned {
typedef Vec type;
};
#endif
typedef not_strided_tag is_strided;
enum { M = MM, N = NN, IS_STRIDED = 0 };
typedef size_t size_type;
typedef dimension<size_type, MM, NN> dim_type;
typedef dimension<int> band_type;
typedef packed_offset<size_type,MM,NN> transpose_type; /* bogus */
typedef not_strideable strideability;
inline banded_offset()
: dim(0,0), bw(band_type(0,0)), ndiag(0) { }
inline banded_offset(size_type m, size_type n, size_type /* lead */,
band_type band)
: dim(m,n), bw(band), ndiag(band.first() + band.second() + 1) { }
inline banded_offset(size_type m, size_type n, size_type /* lead */)
: dim(m,n), bw(band_type(0,0)), ndiag(0) { }
inline banded_offset& operator=(const banded_offset& x) {
dim = x.dim; ndiag = x.ndiag; bw = x.bw; return *this;
}
inline size_type elt(size_type i, size_type j) const {
return this->oned_offset(i) + j;
}
inline size_type oned_offset(size_type i) const {
return i * ndiag + MTL_MAX(0, bw.first() - int(i));
}
inline size_type oned_length(size_type i) const {
return MTL_MAX(0, MTL_MIN(int(dim.second()), int(i) + bw.second() + 1)
- MTL_MAX(0, int(i) - bw.first()));
}
inline size_type twod_length() const { return dim.first(); }
inline int stride() const { return 1; }
inline static size_type size(size_type m, size_type n,
size_type low, size_type up) {
/* M' = number of diagonals = low + up + 1
N' = min (m, n + low) */
return (low + up + 1) * MTL_MIN(m, n + low);
}
inline size_type major() const { return dim.first(); }
inline size_type minor() const { return dim.second(); }
private:
dim_type dim;
band_type bw; /* bandwidth */
size_type ndiag;
};
//: blah
//!noindex:
template <int M, int N>
struct gen_banded_offset {
#if defined( _MSVCPP_ )
typedef banded_offset<unsigned int, M, N> type;
#else
template <class size_type>
struct bind {
typedef banded_offset<size_type, M, N> type;
};
#endif
typedef gen_packed_offset<M,N> transpose_type; // bogus
typedef gen_banded_view_offset<M,N> banded_view_type; /* bogus */
};
//: Packed Offset Class
// This cooresponds to lapack/blas packed storage format
//!models: Offset
//!category: utilities
//!component: type
template <class size_t, int MM, int NN>
class packed_offset {
public:
#if !defined(_MSVCPP_)
template <class Vec>
struct bind_oned {
typedef Vec type;
};
#endif
typedef not_strided_tag is_strided;
enum { M = MM, N = NN, IS_STRIDED = 0 };
typedef size_t size_type;
typedef dimension<size_type, MM, NN> dim_type;
typedef dimension<int> band_type;
typedef banded_offset<size_type, MM, NN> transpose_type; /* bogus */
typedef not_strideable strideability;
inline packed_offset()
: dim(0,0), bw(band_type(0,0)) { }
inline packed_offset(size_type m, size_type n, size_type /* lead */,
band_type bandwidth)
: dim(m,n), bw(bandwidth) { }
inline packed_offset& operator=(const packed_offset& x) {
dim = x.dim; bw = x.bw; return *this;
}
inline int elt(size_type i, size_type j) const {
return this->oned_offset(i) + j;
}
inline int calc_low(int i, int low) const {
int l = MTL_MIN(low, int(i));
int lower_area = low * i;
lower_area -= ( - l*l + 2*low*l + l) / 2;
return lower_area;
}
inline int calc_up(int i, int up) const {
int upper_area = up * i;
int n = i + up - dim.second();
if (n > 0) {
int n1 = MTL_MAX(n - up, 0);
int n2 = n - n1;
upper_area -= n1 * up;
upper_area -= ((n2 + 1) * n2) / 2;
}
return upper_area;
}
inline int oned_offset(size_type i) const { /* the ith major container */
int low = bw.first();
int up = bw.second();
int upper_area, lower_area;
if (up < -1)
upper_area = - calc_low(i, - (up + 1));
else if (up > 0)
upper_area = calc_up(i, up);
else
upper_area = 0;
if (low < -1)
lower_area = - calc_up(i, - (low + 1));
else if (low > 0)
lower_area = calc_low(i, low);
else
lower_area = 0;
size_type diagonal_len;
if (up < 0 || low < 0)
diagonal_len = 0;
else
diagonal_len = MTL_MIN(MTL_MIN(i, dim.first()), dim.second());
size_type ret = upper_area + lower_area + diagonal_len;
return ret;
}
inline int stride() const { return 1; }
inline size_type oned_length(size_type i) const {
return MTL_MAX(0, MTL_MIN(int(dim.second()), int(i) + bw.second() + 1)
- MTL_MAX(0, int(i) - bw.first()));
}
inline size_type twod_length() const { return dim.first(); }
inline static size_type size(int m, int n, int low, int up) {
packed_offset offset(m, n, n, band_type(low, up));
return offset.oned_offset(m);
}
inline size_type major() const { return dim.first(); }
inline size_type minor() const { return dim.second(); }
private:
dim_type dim;
band_type bw; /* bandwidth */
};
//: blah
//!noindex:
template <int M, int N>
struct gen_packed_offset {
#if defined( _MSVCPP_ )
typedef packed_offset<unsigned int, M, N> type;
#else
template <class size_type>
struct bind {
typedef packed_offset<size_type, M, N> type;
};
#endif
typedef gen_banded_offset<M,N> transpose_type; /* bogus */
typedef gen_banded_view_offset<M,N> banded_view_type; /* bogus */
};
/* egcs doesn't "see" the friend functions
* when the dense2D_iterator class is defined inside of dense2D
*/
//: blah
//!noindex:
template <int isConst, class T, class Offset, class InnerOneD, class OneD>
class dense2D_iterator {
public:
typedef typename Offset::size_type size_type;
typedef std::pair<size_type,size_type> pair_type;
typedef typename IF<isConst, const T*,T*>::RET Iterator;
typedef dense2D_iterator self;
typedef int distance_type;
typedef int difference_type;
typedef std::random_access_iterator_tag iterator_category;
typedef OneD* pointer;
typedef OneD value_type;
typedef OneD reference;
typedef difference_type Distance;
typedef Iterator iterator_type;
protected:
Iterator start;
size_type pos;
size_type ld; /* leading dimension */
pair_type starts;
Offset offset;
public:
inline size_type index() const { return pos + starts.second; }
inline dense2D_iterator () {}
inline dense2D_iterator(const self& x)
: start(x.start), pos(x.pos),
ld(x.ld), starts(x.starts), offset(x.offset) { }
inline self& operator=(const self& x) {
start = x.start;
pos = x.pos;
ld = x.ld;
starts = x.starts;
offset = x.offset;
return *this;
}
inline explicit
dense2D_iterator(Iterator x, size_type ld_, size_type p, pair_type s,
Offset os)
: start(x), pos(p), ld(ld_), starts(s), offset(os) { }
inline Iterator base () const { return start + pos; }
inline reference deref(Distance pos, not_strided_tag) const {
return reference((T*)start + offset.oned_offset(pos),
offset.oned_length(pos),
starts.first);
}
inline reference deref(Distance pos, strided_tag) const {
InnerOneD vec((T*)start + offset.oned_offset(pos),
offset.oned_length(pos),
starts.first);
return strided(vec, offset.stride());
}
inline reference operator*() const {
typedef typename Offset::is_strided Strided;
return deref(pos, Strided());
}
inline reference operator[] (Distance n) const {
typedef typename Offset::is_strided Strided;
return deref(pos + n, Strided());
}
/* won't work, the OneD is temporary
pointer operator-> () const { return & (operator* ()); }
*/
inline self& operator++ () { ++pos; return *this; }
inline self operator++ (int) { self tmp = *this; ++pos; return tmp; }
inline self& operator-- () { --pos; return *this; }
inline self operator-- (int) { self tmp = *this; --pos; return tmp; }
inline self& operator+=(size_type n) { pos += n; return *this; }
inline self operator+(size_type n) const {
return self(start, ld, pos + n, starts);
}
inline self& operator-=(size_type n) { pos -= n; return *this; }
};
template <int isConst, class T, class Offset, class InnerOneD, class OneD>
inline typename dense2D_iterator<isConst,T,Offset,InnerOneD,OneD>::difference_type
operator-(const dense2D_iterator<isConst,T,Offset,InnerOneD,OneD>& x,
const dense2D_iterator<isConst,T,Offset,InnerOneD,OneD>& y)
{
return x.index() - y.index();
}
template <int isConst,class T, class Offset, class InnerOneD, class OneD>
inline bool
operator== (const dense2D_iterator<isConst,T,Offset,InnerOneD,OneD>& x,
const dense2D_iterator<isConst,T,Offset,InnerOneD,OneD>& y)
{
return x.index() == y.index();
}
template <int isConst, class T, class Offset, class InnerOneD, class OneD>
inline bool
operator!= (const dense2D_iterator<isConst, T,Offset,InnerOneD,OneD>& x,
const dense2D_iterator<isConst, T,Offset,InnerOneD,OneD>& y)
{
return x.index() != y.index();
}
template <int isConst,class T, class Offset, class InnerOneD, class OneD>
inline bool
operator< (const dense2D_iterator<isConst,T,Offset,InnerOneD,OneD>& x,
const dense2D_iterator<isConst,T,Offset,InnerOneD,OneD>& y)
{
return x.index() < y.index();
}
/*
Workaround (g++ 2.91) helper class
*/
template <class Strided>
struct __bracket { };
template <>
struct __bracket<strided_tag> {
template <class OneD, class InnerOneD, class elt_type, class size_type>
inline OneD
operator()(elt_type* d, size_type len, size_type f, size_type ld,
const OneD*, const InnerOneD*) {
InnerOneD vec(d , len, f);
return OneD(vec, ld);
}
};
template <>
struct __bracket<not_strided_tag> {
template <class OneD, class InnerOneD, class elt_type, class size_type>
inline OneD
operator()(elt_type* d, size_type len, size_type f, size_type,
const OneD*, const InnerOneD*) {
return OneD(d, len, f);
}
};
template<class T, class OffsetGen, int MM, int NN>
class dense2D;
template <class T, class OffsetGen, int MM, int NN>
class external2D;
//: Generic Dense 2-D Container
//!category: containers
//!component: type
//
// The generic_dense2D container implements sevaral of the MTL storage
// types. They include dense, packed, banded, and banded_view. The
// common theme here is that the matrix is stored in a contiguous
// piece of memory. The differences in these storage types has to do
// with where to find the OneD segements in the linear
// memory. Caclulating these offsets is the job of the Offset concept,
// which has a model to handle each of the different storage types:
// rect_offset, strided_offset, banded_offset, packed_offset, and
// banded_view_offset.
//
// There are two derived classes of generic_dense2D that specify the
// memory management, dense2D and external2D. The dense2D version owns
// its memory, while the external2D imports its memory from somewhere
// else through a pointer (which allows for interoperability with
// other codes -- even with Fortran!). <p>
//
//!definition: dense2D.h
//!tparam: RepType - The Container used to store the elements
//!tparam: RepPtr - The type used to reference to the container
//!tparam: OffsetGen - The generator that creates the Offset class
//!tparam: MM - For static sized matrix, the major dimension
//!tparam: NN - For static sized matrix, the minor dimension
//!models: TwoDStorage
template <class RepType, class RepPtr, class OffsetGen, int MM, int NN>
class generic_dense2D {
public:
//: Static sizes (0 if dynamic)
enum { M = MM, N = NN };
//: The type for dimensions and indices
typedef typename RepType::size_type size_type;
//: The type for differences between iterators
typedef typename RepType::difference_type difference_type;
protected:
typedef std::pair<size_type,size_type> pair_type;
typedef RepType reptype;
typedef RepPtr rep_ptr;
typedef typename RepType::value_type elt_type;
#if defined(_MSVCPP_)
//JGS Nasty VC++ workaround
typedef typename OffsetGen::type Offset;
#else
typedef typename OffsetGen:: MTL_TEMPLATE bind<size_type>::type Offset;
#endif
typedef dimension<elt_type> dyn_dim;
public:
//: A pair type for dimensions
typedef typename Offset::dim_type dim_type;
//: A pair type for bandwidth
typedef typename Offset::band_type band_type;
/* Type Definitions */
//: This is a dense matrix
typedef dense_tag sparsity;
typedef typename Offset::is_strided is_strided;
protected:
typedef external_vec<elt_type, N> InnerOneD;
#if defined(_MSVCPP_)
enum { offset_strided = Offset::IS_STRIDED };
typedef typename IF<offset_strided, strided1D<InnerOneD>, InnerOneD>::RET OneD;
#else
typedef typename Offset:: MTL_TEMPLATE bind_oned<InnerOneD>::type OneD;
#endif
typedef OneD OneDRef;
typedef OneD ConstOneDRef;
public:
//: The 1D container type
typedef OneD value_type;
//: The type for a reference to value_type
typedef value_type reference;
//: The type for a const reference to value_type
typedef value_type const_reference;
//: The iterator type
typedef dense2D_iterator<0,elt_type, Offset, InnerOneD, OneD> iterator;
//: The const iterator type
typedef dense2D_iterator<1,elt_type, Offset, InnerOneD, OneD> const_iterator;
//: The reverse iterator type
typedef reverse_iter<iterator> reverse_iterator;
//: The const reverse iterator type
typedef reverse_iter<const_iterator> const_reverse_iterator;
//: The type for the transpose of this container
typedef generic_dense2D<RepType, RepPtr,
typename OffsetGen::transpose_type, MM, NN> transpose_type;
//: The type for a banded view of this container
typedef generic_dense2D<RepType, RepPtr,
typename OffsetGen::banded_view_type, MM, NN> banded_view_type;
//: The type for a sub-section of this 2D container
typedef external2D<elt_type, OffsetGen, MM, NN> submatrix_type;
#ifndef MTL_DISABLE_BLOCKING
template <class Block>
struct blocked_view {
typedef block2D<Block, OffsetGen> type;
};
#endif
//: This is a stridable container, can use rows(A), columns(A)
typedef typename Offset::strideability strideability;
/* Constructors */
//: Default Constructor
inline generic_dense2D()
: ld_(0), data_(0), starts(std::make_pair(0,0)) { }
//: Normal Constructor
inline generic_dense2D(rep_ptr data, size_type m, size_type n, size_type ld)
: ld_(ld), data_(data),
starts(std::make_pair(0,0)), offset(m, n, ld) { }
//: Constructor with non-zero upper-left corner indices
inline generic_dense2D(rep_ptr data, size_type m, size_type n,
size_type ld, dyn_dim s, char)
: ld_(ld), data_(data),
starts(std::make_pair(s.first(),s.second())), offset(m, n, ld) { }
//: Static M, N constructor
inline generic_dense2D(rep_ptr data, size_type ld)
: ld_(ld), data_(data),
starts(std::make_pair(0,0)), offset(M, N, ld) { }
//: with bandwidth constructor
inline generic_dense2D(rep_ptr data, size_type m, size_type n, size_type ld,
band_type bw)
: ld_(ld), data_(data), starts(std::make_pair(0,0)),
offset(m, n, ld, bw) { }
//: Static M, N with bandwith?
//: Copy Constructor
inline generic_dense2D(const generic_dense2D& x)
: ld_(x.ld_), data_(x.data_), starts(x.starts),
offset(x.offset) { }
//: Assignment Operator
inline generic_dense2D& operator=(const generic_dense2D& x) {
ld_ = x.ld_; data_ = x.data_; starts = x.starts; offset = x.offset;
return *this;
}
//: Subclass Constructor
inline generic_dense2D(rep_ptr d, const generic_dense2D& x)
: ld_(x.ld_), data_(d), starts(x.starts), offset(x.offset) { }
//: Transpose Constructor
inline generic_dense2D(const transpose_type& x, do_transpose, do_transpose)
: ld_(x.ld_), data_(x.data_), starts(x.starts), offset(x.offset) { }
/* JGS, remove stream constructor, not very necessary
just have them call another constructor
*/
//: Matrix Stream Constructor
template <class MatrixStream, class Orien>
inline generic_dense2D(rep_ptr data, MatrixStream& s, Orien)
: ld_(Orien::map(dim_type(s.nrows(),s.ncols())).second()),
data_(data),
starts(std::make_pair(0,0)),
offset(Orien::map(dim_type(s.nrows(),s.ncols())).first(),
Orien::map(dim_type(s.nrows(),s.ncols())).second(),
Orien::map(dim_type(s.nrows(),s.ncols())).second()) { }
//: Banded Matrix Stream Constructor
template <class MatrixStream, class Orien>
inline generic_dense2D(rep_ptr data, MatrixStream& s,
Orien, band_type bw)
: ld_(Orien::map(dim_type(s.nrows(),s.ncols())).second()),
data_(data),
starts(std::make_pair(0,0)),
offset(Orien::map(dim_type(s.nrows(),s.ncols())).first(),
Orien::map(dim_type(s.nrows(),s.ncols())).second(),
Orien::map(dim_type(s.nrows(),s.ncols())).second(),
bw) { }
//: Banded View Constructor
template <class TwoD>
inline generic_dense2D(rep_ptr data, const TwoD& x, band_type bw, banded_tag)
: ld_(x.ld_),
data_(data),
starts(x.starts),
offset(x.offset, bw) { }
// VC++ doesn't like this
//friend class transpose_type;
//: The destructor.
inline ~generic_dense2D() { }
/* Access Methods */
/* Iterator Access Methods */
//: Return an iterator pointing to the first 1D container
inline iterator begin() {
return iterator(data(), ld_, 0, starts, offset);
}
//: Return an iterator pointing past the end of the 2D container
inline iterator end() {
return iterator(data(), ld_, offset.twod_length(), starts, offset);
}
//: Return a const iterator pointing to the first 1D container
inline const_iterator begin() const {
return const_iterator(data(), ld_, 0, starts, offset);
}
//: Return a const iterator pointing past the end of the 2D container
inline const_iterator end() const {
return const_iterator(data(), ld_, offset.twod_length(),
starts, offset);
}
/* reverse iterators */
//: Return a reverse iterator pointing to the last 1D container
inline reverse_iterator rbegin() {
return reverse_iterator(end());
}
//: Return a reverse iterator pointing past the start of the 2D container
inline reverse_iterator rend() {
return reverse_iterator(begin());
}
//: Return a const reverse iterator pointing to the last 1D container
inline const_reverse_iterator rbegin() const {
return const_reverse_iterator(end());
}
//: Return a const reverse iterator pointing past the start of the 2D container
inline const_reverse_iterator rend() const {
return const_reverse_iterator(begin());
}
/* Element Access Methods */
//: Return a reference to the (i,j) element, where (i,j) is in the 2D coordinate system
inline const elt_type& operator()(size_type i, size_type j) const {
return *(data() + offset.elt(i, j));
}
//: Return a const reference to the (i,j) element, where (i,j) is in the 2D coordinate system
inline elt_type& operator()(size_type i, size_type j) {
return *(data() + offset.elt(i, j));
}
/* Size Methods */
//: Number of non-zeroes
inline size_type nnz() const { return offset.major() * offset.minor(); }
//: Capacity
inline size_type capacity() const { return offset.major() * offset.minor(); }
//: Major axis size
inline size_type major() const { return offset.major(); }
//: Minor axis size
inline size_type minor() const { return offset.minor(); }
//: Leading Dimension
inline size_type ld() const { return ld_; }
//: Memory Access
inline const elt_type* data() const { return &(*data_)[0]; }
inline elt_type* data() { return &(*data_)[0]; }
/* obsolete
inline const elt_type* get_contiguous() const { return data_->data(); }
inline elt_type* get_contiguous() { return data_->data(); }
inline void set_contiguous(elt_type*) { }
*/
/* Vector Access Methods */
//: OneD Access
inline OneD operator[](size_type i) const {
typedef OneD* oned_ptr;
typedef InnerOneD* inner_oned_ptr;
return __bracket<is_strided>()((elt_type*)data() + offset.oned_offset(i),
offset.oned_length(i),
starts.first, ld_,
oned_ptr(),
inner_oned_ptr());
}
/* All the submatrix stuff is in matrix_implementation for now
inline submatrix_type sub_matrix(size_type m_start, size_type m_finish,
size_type n_start, size_type n_finish) {
return submatrix_type(data_->data() + m_start * ld_ + n_start,
dim_type(m_finish - m_start, n_finish - n_start), ld_);
}
inline submatrix_type sub_matrix(size_type m_start, size_type n_start,
size_type m, size_type n) {
return submatrix_type(data_->data() + m_start * ld_ + n_start,
dim_type(m, n), ld_);
}
inline submatrix_type section(size_type m_start, size_type n_start,
size_type m, size_type n) {
return submatrix_type(data_->data() + m_start * ld_ + n_start,
dim_type(m, n), ld_, dim_type(m_start, n_start));
}
typedef range<size_type> Range;
inline generic_dense2D operator()(Range m, Range n) {
return generic_dense2D(data_->data() + m.start * ld_ + n.start,
m.finish - m.start, n.finish - n.start, ld_);
}
inline OneD::subrange_type operator()(size_type i, Range n) {
return operator[i](n);
}
typedef strided1D< InnerOneD > MinorVector;
inline MinorVector minor_vector(size_type i) const {
InnerOneD vec((elt_type*)data_->data() + i,
offset.major() * ld_, starts.first);
return MinorVector(vec, ld_);
}
inline MinorVector::subrange_type operator()(Range m, size_type j) {
return minor_vector(j)(m);
}
*/
/*JGS friend not working for transpose constructor
protected:
*/
size_type ld_;/* JGS redundant */
rep_ptr data_;
pair_type starts;
Offset offset;
};
template <class T, class OffsetGen, int M, int N>
struct gen_dense2D;
/* why didn't I use std::vector here?
or perhaps I should just use plain old memory here?
will that work with the reference counting in terms
of deallocating?
*/
//: Dense2D Storage Type
//
// Inherits from generic_dense2D. The class "owns" its data.
//
//!category: containers
//!component: type
//!tparam: T - the element type
//!tparam: OffsetGen - the Offset class generator
//!tparam: MM - For static sized matrix, the major dimension
//!tparam: NN - For static sized matrix, the minor dimension
//!models: TwoDStorage
template<class T, class OffsetGen, int MM = 0, int NN = 0>
class dense2D
: public generic_dense2D< std::vector<T> ,
refcnt_ptr< std::vector<T> >, OffsetGen, MM, NN >
/* : public generic_dense2D< bare_bones_array<T> ,
refcnt_ptr< bare_bones_array<T> >, OffsetGen, MM, NN >
*/
{
public:
typedef generic_dense2D< std::vector<T> ,
refcnt_ptr< std::vector<T> >, OffsetGen, MM, NN> super;
/* typedef generic_dense2D< bare_bones_array<T> ,
refcnt_ptr< bare_bones_array<T> >, OffsetGen, MM, NN> super;
*/
typedef typename super::Offset Offset;
//: Pair type for dimension
typedef typename Offset::dim_type dim_type;
//: Pair type for bandwidth
typedef typename Offset::band_type band_type;
typedef typename super::reptype reptype;
typedef typename super::rep_ptr rep_ptr;
//: Unsigned integral type for dimensions and indices
typedef typename super::size_type size_type;
//: The transpose type
typedef dense2D<T, typename OffsetGen::transpose_type,
MM, NN> transpose_type;
// VC++ doesn't like this
//friend class transpose_type;
//: This has internal storage
typedef internal_tag storage_loc;
//: Default Constructor
inline dense2D() { }
//: Constructor from Dimension Pair
inline dense2D(dim_type dim)
: super(new reptype(Offset::size(dim.first(), dim.second(), 0, 0)),
dim.first(),
dim.second(),
dim.second()) { }
//: Constructor from Dimension and Bandwidth Pairs
inline dense2D(dim_type dim, band_type bw)
: super(new reptype(Offset::size(dim.first(),dim.second(),
bw.first(), bw.second())),
dim.first(),
dim.second(),
dim.second(),
bw) { }
//: Copy Constructor
inline dense2D(const dense2D& x)
: super(x) { }
//: Assignment Operator
inline dense2D& operator=(const dense2D& x) {
super::operator=(x);
return *this;
}
//: Transpose Constructor
inline dense2D(const transpose_type& x, do_transpose t, do_transpose)
: super(x, t, t) { }
#if !defined(_MSVCPP_)
// JGS, use actual stream types
//: Matrix Stream Constructor
template <class MatrixStream, class Orien>
inline dense2D(MatrixStream& s, Orien)
: super(new reptype(Offset::size(Orien::map(dim_type(s.nrows(),
s.ncols())).first(),
Orien::map(dim_type(s.nrows(),
s.ncols())).second(),
0, 0)),
s,
Orien()) { }
//: Matrix Stream Constructor with bandwidth
template <class MatrixStream, class Orien>
inline dense2D(MatrixStream& s, Orien, band_type bw)
: super(new reptype(Offset::size(Orien::map(dim_type(s.nrows(),
s.ncols())).first(),
Orien::map(dim_type(s.nrows(),
s.ncols())).second(),
bw.first(), bw.second())),
s,
Orien(),
bw) { }
#endif
#if 0
// deprecated
template <class SubMatrix>
struct partitioned {
typedef dense2D<SubMatrix, OffsetGen> type;
typedef gen_dense2D<SubMatrix, OffsetGen> generator;
};
#endif
#if 1 // This makes no sense. dense2D can not be a "view"
//: banded view constructor
template <class TwoD>
inline dense2D(const TwoD& x, band_type bw, banded_tag)
: super(x.data_, x, bw, banded_tag()) { }
#endif
//: Destructor
inline ~dense2D() { }
inline void resize(size_type m, size_type n) {
rep_ptr newdata = new reptype(Offset::size(m, n, 0, 0));
size_type i, j;
size_type M = MTL_MIN(m, offset.major());
size_type N = MTL_MIN(n, offset.minor());
for (i = 0; i < M; ++i)
for (j = 0; j < N; ++j)
(*newdata)[i * n + j] = (*this)(i,j);
for (; i < m; ++i)
for (; j < n; ++j)
(*newdata)[i * n + j] = T();
data_ = newdata;
ld_ = n;
offset.dim = dim_type(m, n);
offset.ld = n;
}
};
template <class T, class OffsetGen, int M, int N>
struct gen_external2D;
#ifndef MTL_DISABLE_BLOCKING
template <class Block, class OffsetGen, int M, int N>
struct gen_block2D;
#endif
//: blah
//!noindex:
template <class T, class OffsetGen, int M, int N>
struct gen_dense2D {
typedef gen_dense2D<T, typename OffsetGen::transpose_type,N,M> transpose_type;
typedef gen_external2D<T, OffsetGen,M,N> submatrix_type;
#ifndef MTL_DISABLE_BLOCKING
template <class Block>
struct blocked_view {
typedef gen_block2D<Block, OffsetGen, M, N> type;
};
#endif
typedef gen_dense2D<T, typename OffsetGen::banded_view_type,M,N>
banded_view_type;
typedef dense2D<T, OffsetGen, M, N> type;
};
//: External2D Storage Type
//
// Inherits from generic_dense2D. The class does not "own" its data.
//
//!category: containers
//!component: type
//!tparam: T - the element type
//!tparam: OffsetGen - the Offset class generator
//!tparam: MM - For static sized matrix, the major dimension
//!tparam: NN - For static sized matrix, the minor dimension
//!models: TwoDStorage
//
template <class T, class OffsetGen, int MM = 0, int NN = 0>
class external2D
: public generic_dense2D< external_vec<T,NN>,
external_vec<T,NN>*, OffsetGen, MM, NN >
{
typedef generic_dense2D< external_vec<T,NN>,
external_vec<T,NN>*, OffsetGen, MM, NN > super;
public:
external_vec<T,NN> rep;
typedef dimension<T> dyn_dim;
typedef typename super::Offset Offset;
//: Pair type for dimension
typedef typename Offset::dim_type dim_type;
//: Pair type for bandwidth
typedef typename Offset::band_type band_type;
typedef typename super::reptype reptype;
//: Unsigned integral type for dimensions and indices
typedef typename super::size_type size_type;
//: Type for the transpose
typedef external2D<T, typename OffsetGen::transpose_type,
MM, NN> transpose_type;
// VC++ doesn't like this
//friend class transpose_type;
//: This has external storage
typedef external_tag storage_loc;
//: Default Constructor
inline external2D() { }
//: Construct from pointer and dimensions
inline external2D(T* data, dim_type dim)
: super(&rep, dim.first(), dim.second(), dim.second()),
rep(data, dim.first() * dim.second())
{ }
//: Construct from pointer, dimensions, and leading dimension
inline external2D(T* data, dim_type dim, size_type ld)
: super(&rep, dim.first(), dim.second(), ld),
rep(data, dim.first() * ld)
{ }
//: non-zero indices in upper left corner
inline external2D(T* data, dim_type dim, size_type ld,
dyn_dim s, char)
: super(&rep, dim.first(), dim.second(), ld, s, char()),
rep(data, dim.first() * ld)
{ }
//: Constructor with bandwith
inline external2D(T* data, dim_type dim, band_type bw)
: super(&rep, dim.first(), dim.second(), dim.second(), bw),
rep(data, dim.first() * dim.second())
{ }
//: Constructor with leading dimension and bandwith
inline external2D(T* data, dim_type dim, size_type ld, band_type bw)
: super(&rep, dim.first(), dim.second(), ld, bw),
rep(data, dim.first() * ld)
{ }
//: Copy Constructor
inline external2D(const external2D& x)
: super(&rep, x), rep(x.rep)
{ }
//: Assignment Operator
inline external2D& operator=(const external2D& x) {
rep = x.rep;
super::operator=(x);
data_ = &rep;
return *this;
}
//: Transpose Constructor
inline external2D(const transpose_type& x, do_transpose t, do_transpose)
: super(x, t, t), rep(x.rep) { }
/* JGS This conflicts with the external2D(T* data, dim_type dim,
band_type bw) constructor, and I am not sure this is really
needed anyway.
//: Matrix Stream Constructor
template <class MatrixStream, class Orien>
inline external2D(T* data, MatrixStream& s, Orien o)
: rep(data, s.nrows() * s.ncols()), super(&rep, s, o) { }
template <class MatrixStream, class Orien>
inline external2D(T* data, MatrixStream& s, Orien o,
band_type bw)
: super(&rep, s, o, bw),
rep(data, s.nrows() * s.ncols())
{ }
*/
//: banded view constructor
template <class TwoD>
inline external2D(const TwoD& x, band_type bw, banded_tag)
: super(&rep, x, bw, banded_tag()), rep((T*)x.data(), x.major() * x.ld()) { }
inline ~external2D() { }
#if 0
// deprecated
template <class SubMatrix>
struct partitioned {
typedef dense2D<SubMatrix, OffsetGen> type;
typedef gen_dense2D<SubMatrix, OffsetGen> generator;
};
#endif
};
//: blah
//!noindex:
template <class T, class OffsetGen, int M, int N>
struct gen_external2D {
typedef gen_external2D<T, typename OffsetGen::transpose_type, N, M> transpose_type;
typedef gen_external2D<T, OffsetGen,M,N> submatrix_type;
typedef gen_external2D<T, typename OffsetGen::banded_view_type,M,N>
banded_view_type;
#ifndef MTL_DISABLE_BLOCKING
template <class Block>
struct blocked_view {
typedef gen_block2D<Block, OffsetGen, M, N> type;
};
#endif
typedef external2D<T, OffsetGen, M, N> type;
};
} /* namespace mtl */
#endif /* MTL_DENSE2D_H */
light_matrix.h
(application/octet-stream, 12.8 KB)
#ifndef MTL_LIGHT_MATRIX_H
#define MTL_LIGHT_MATRIX_H
#include "mtl/matrix_traits.h"
#include "mtl/dimension.h"
#include "mtl/meta_if.h"
#include "mtl/meta_equal.h"
namespace mtl {
template <int Orien>
struct TRANS {
enum { RET = 0 };
};
template<>
struct TRANS<ROW_MAJOR> {
enum { RET = COL_MAJOR };
};
template<>
struct TRANS<COL_MAJOR> {
enum { RET = ROW_MAJOR };
};
template <class T, class SizeType, int Orien, int Strided>
class light_matrix {
public:
typedef light_matrix self;
typedef light_matrix light_matrix_t; // VC++ workaround
typedef T* DataPtr;
typedef rectangle_tag shape;
typedef typename IF< EQUAL<Orien,ROW_MAJOR>::RET,
row_tag, column_tag>::RET orientation; // mostly wrong
typedef typename IF< EQUAL<Orien,ROW_MAJOR>::RET,
row_orien, column_orien>::RET orien;
typedef light_matrix<T, SizeType, TRANS<Orien>::RET, Strided> transpose_type;
typedef light_matrix<T, SizeType, Orien, !Strided> strided_type;
typedef light_matrix<T, SizeType, Orien, Strided> scaled_type;// wrong
typedef light_matrix<T, SizeType, Orien, Strided> submatrix_type;
typedef int DiffType;
//: The size type
typedef SizeType size_type;
//: The type for differences between iterators
typedef DiffType difference_type;
typedef T value_type;
typedef value_type& reference;
typedef const value_type& const_reference;
typedef value_type* pointer;
enum { M = 0, N = 0 };
protected:
static inline size_type& twod_pos(size_type& i, size_type& j) {
if (Orien == ROW_MAJOR)
return i;
else
return j;
}
static inline const size_type& twod_pos(const size_type& i,
const size_type& j) {
if (Orien == ROW_MAJOR)
return i;
else
return j;
}
static inline size_type& oned_pos(size_type& i, size_type& j) {
if (Orien == ROW_MAJOR)
return j;
else
return i;
}
static inline const size_type& oned_pos(const size_type& i,
const size_type& j) {
if (Orien == ROW_MAJOR)
return j;
else
return i;
}
// idea: completely separate stride/offset/positioning from indexing
// but encapsulate both somehow
public:
//: This is a dense 2D container
typedef dense_tag sparsity;
//: This has external storage
typedef external_tag storage_loc;
//: This is strideable
typedef strideable strideability;
class oned {
public:
typedef T& reference;
typedef const T& const_reference;
typedef T value_type;
typedef T* pointer;
typedef SizeType size_type;
typedef int difference_type;
enum { M = 0, N = 0 };
typedef oned subrange_type;
typedef dense_tag sparsity;
typedef oned IndexArray; /* bogus */
typedef oned IndexArrayRef; /* bogus */
typedef oned_tag dimension; /* bogus */
template <int isConst>
class __iterator {
typedef __iterator self;
public:
typedef typename oned::value_type value_type;
typedef oned::pointer pointer;
typedef typename oned::size_type size_type;
typedef oned::difference_type difference_type;
typedef typename IF<isConst, oned::const_reference, oned::reference>::RET reference;
typedef std::random_access_iterator_tag iterator_category;
inline __iterator(DataPtr d,
size_type ii, size_type jj,
size_type os, size_type s)
: data(d), i(ii), j(jj), offset(os), stride(s) { }
inline __iterator(const self& x)
: data(x.data), i(x.i), j(x.j), offset(x.offset), stride(x.stride) { }
inline self& operator=(const self& x) {
data = x.data; i = x.i; j = x.j; offset = x.offset; stride = x.stride;
return *this;
}
inline __iterator() : data(0), i(0), j(0), offset(0), stride(0) { }
inline reference operator*() const { return data[offset]; }
inline self& operator++() { ++pos(); offset += stride; return *this; }
inline self& operator+=(size_type n) {
pos() += n; offset += stride*n; return *this;
}
inline self operator++(int) { self t = *this; ++(*this); return t; }
inline self& operator--() { --pos(); offset -= stride; return *this; }
inline self& operator-=(size_type n) {
pos() -= n; offset -= stride*n; return *this; }
inline self operator--(int) { self t = *this; --(*this); return t; }
inline bool operator!=(const self& x) const { return pos() != x.pos(); }
inline bool operator==(const self& x) const { return pos() == x.pos(); }
inline bool operator<(const self& x) const { return pos() < x.pos(); }
inline size_type index() const { return pos(); }
inline size_type& pos() { return oned_pos(i,j); }
inline const size_type& pos() const { return oned_pos(i,j); }
inline size_type row() const { return i; }
inline size_type column() const { return j; }
protected:
DataPtr data;
size_type i, j;
size_type offset;
size_type stride;
};
typedef __iterator<0> iterator;
typedef __iterator<1> const_iterator;
inline oned(DataPtr d, size_type ii, size_type jj,
size_type ie, size_type je,
size_type os, size_type ld)
: data(d), i(ii), j(jj), iend(ie), jend(je),
offset(os), ldim(ld) { }
inline oned(const oned& x)
: data(x.data), i(x.i), j(x.j),
iend(x.iend), jend(x.jend),
offset(x.offset), ldim(x.ldim) { }
inline oned& operator=(const oned& x) {
data = x.data; i = x.i; j = x.j;
iend = x.iend; jend = x.jend;
offset = x.offset; ldim = x.ldim;
return *this;
}
inline oned()
: data(0), i(0), j(0), iend(0), jend(0), offset(0), ldim(0) { }
inline ~oned() { }
inline reference operator[](size_type n) {
return data[ Strided ? offset + n * ldim : offset + n];
}
inline const_reference operator[](size_type n) const {
return data[ Strided ? offset + n * ldim : offset + n];
}
inline iterator begin() {
return iterator(data, i, j, offset, Strided ? ldim : 1);
}
inline iterator end() {
size_type iiend, jjend;
if (Orien == ROW_MAJOR) { iiend = i; jjend = jend; }
else { iiend = iend; jjend = j; }
return iterator(data, iiend, jjend, offset, Strided ? ldim: 1);
}
inline const_iterator begin() const {
return const_iterator(data, i, j, offset, Strided ? ldim : 1);
}
inline const_iterator end() const {
size_type iiend, jjend;
if (Orien == ROW_MAJOR) { iiend = i; jjend = jend; }
else { iiend = iend; jjend = j; }
return const_iterator(data, iiend, jjend, offset, Strided ? ldim : 1);
}
protected:
DataPtr data;
size_type i, j;
size_type iend, jend;
size_type offset;
size_type ldim;
};
typedef oned OneD;
typedef OneD OneDRef;
typedef OneD Row;
typedef OneD RowRef;
typedef OneD Column;
typedef OneD ColumnRef;
//: The iterator type
template <int Const>
class __iterator {
typedef __iterator self;
public:
typedef std::random_access_iterator_tag iterator_category;
typedef oned value_type;
typedef value_type* pointer;
#if defined(_MSVCPP_)
typedef typename light_matrix_t::size_type size_type;
typedef typename light_matrix_t::difference_type difference_type;
#else
typedef SizeType size_type;
typedef DiffType difference_type;
#endif
typedef typename IF<Const, const oned, oned>::RET reference;
inline __iterator(DataPtr d, size_type ii, size_type jj,
size_type ie, size_type je, size_type ld)
: data(d), i(ii), j(jj), iend(ie), jend(je), offset(0), ldim(ld) {
if (Strided) stride = 1; else stride = ldim;
}
inline __iterator() : data(0), i(0), j(0), iend(0), jend(0),
offset(0), stride(0), ldim(0) { }
inline __iterator(const self& x)
: data(x.data), i(x.i), j(x.j),
iend(x.iend), jend(x.jend), offset(x.offset),
stride(x.stride), ldim(x.ldim) { }
inline self& operator=(const self& x) {
data = x.data; i = x.i; j = x.j;
iend = x.iend; jend = x.jend;
offset = x.offset; stride = x.stride; ldim = x.ldim;
return *this;
}
inline reference operator*() const {
return oned(data, i, j, iend, jend, offset, ldim);
}
inline self& operator++() { ++pos(); offset += stride; return *this; }
inline self& operator+=(size_type n) {
pos() += n; offset += stride*n; return *this; }
inline self operator++(int) { self t = *this; ++(*this); return t; }
inline self& operator--() { --pos(); offset -= stride; return *this; }
inline self& operator-=(size_type n) {
pos() -= n; offset -= stride*n; return *this; }
inline self operator--(int) { self t = *this; --(*this); return t; }
inline bool operator!=(const self& x) const { return pos() != x.pos(); }
inline bool operator==(const self& x) const { return pos() == x.pos(); }
inline bool operator<(const self& x) const { return pos() < x.pos(); }
inline size_type index() const { return pos(); }
inline size_type& pos() { return twod_pos(i,j); }
inline const size_type& pos() const { return twod_pos(i,j); }
inline size_type row() const { return i; }
inline size_type column() const { return j; }
protected:
DataPtr data;
size_type i, j;
size_type iend, jend;
size_type offset;
size_type stride;
size_type ldim;
};
typedef __iterator<0> iterator;
typedef __iterator<1> const_iterator;
//: Standard Constructor
inline light_matrix(DataPtr d, size_type m, size_type n, size_type ld)
: data_(d), nrows_(m), ncols_(n), ldim(ld) { }
inline light_matrix(DataPtr d, size_type m, size_type n)
: data_(d), nrows_(m), ncols_(n), ldim(Orien == ROW_MAJOR ? n : m) { }
//: Copy Constructor
inline light_matrix(const light_matrix& x)
: data_(x.data_), nrows_(x.nrows_), ncols_(x.ncols_), ldim(x.ldim) { }
//: Assignment Operator
inline const light_matrix& operator=(const light_matrix& x) {
data_ = x.data_; nrows_ = x.nrows_; ncols_ = x.ncols_; ldim = x.ldim;
return *this;
}
//: Default Constructor
inline light_matrix() : data_(0), nrows_(0), ncols_(0), ldim(0) { }
inline light_matrix(const strided_type& x, do_strided s)
: data_(x.data_), nrows_(x.nrows_), ncols_(x.ncols_), ldim(x.ldim) { }
template <class StridedType>
inline light_matrix(const StridedType& x, do_strided s)
: data_(x.data_), nrows_(x.nrows_), ncols_(x.ncols_), ldim(x.ldim) { }
//: Destructor
inline ~light_matrix() { }
//: Return an iterator pointing to the first 1D container
inline iterator begin() {
return iterator(data_, 0, 0, nrows_, ncols_, ldim);
}
//: Return an iterator pointing past the end of the 2D container
inline iterator end() {
return iterator(data_, nrows_, ncols_, nrows_, ncols_, ldim);
}
//: Return a const iterator pointing to the first 1D container
inline const_iterator begin() const {
return const_iterator(data_, 0, 0, nrows_, ncols_, ldim);
}
//: Return a const iterator pointing past the end of the 2D container
inline const_iterator end() const {
return const_iterator(data_, nrows_, ncols_, nrows_, ncols_, ldim);
}
//: Return a reference to the ith 1D container
inline oned operator[](size_type n) {
if (Orien == ROW_MAJOR)
return oned(data_, n, 0, nrows_, ncols_, Strided ? n : ldim * n, ldim);
else
return oned(data_, 0, n, nrows_, ncols_, Strided ? n : ldim * n, ldim);
}
inline const oned operator[](size_type n) const {
if (Orien == ROW_MAJOR)
return oned(data_, n, 0, nrows_, ncols_, Strided ? n : ldim * n, ldim);
else
return oned(data_, 0, n, nrows_, ncols_, Strided ? n : ldim * n, ldim);
}
//: Return a reference to the (i,j) element, where (i,j) is in the 2D coordinate system
inline reference operator()(size_type i, size_type j) {
return Orien == ROW_MAJOR ? operator[](i)[j] : operator[](j)[i];
}
//: Return a const reference to the (i,j) element, where (i,j) is in the 2D coordinate system
inline const_reference operator()(size_type i, size_type j) const {
return Orien == ROW_MAJOR ? operator[](i)[j] : operator[](j)[i];
}
inline size_type nrows() const { return nrows_; }
inline size_type ncols() const { return ncols_; }
inline submatrix_type sub_matrix(size_type i, size_type iend,
size_type j, size_type jend) const
{
if (Strided)
return submatrix_type(data_ + oned_pos(i,j) * ldim + twod_pos(iend,jend),
iend - i, jend - j, ldim);
else
return submatrix_type(data_ + twod_pos(i,j) * ldim + oned_pos(iend,jend),
iend - i, jend - j, ldim);
}
DataPtr data_;
size_type nrows_, ncols_;
size_type ldim;
};
} /* namespace mtl */
#endif /* MTL_LIGHT_MATRIX_H */
light1D.h
(application/octet-stream, 5.4 KB)
//
// Copyright 1997, 1998, 1999 University of Notre Dame.
// Authors: Andrew Lumsdaine, Jeremy G. Siek, Lie-Quan Lee
//
// This file is part of the Matrix Template Library
//
// You should have received a copy of the License Agreement for the
// Matrix Template Library along with the software; see the
// file LICENSE. If not, contact Office of Research, University of Notre
// Dame, Notre Dame, IN 46556.
//
// Permission to modify the code and to distribute modified code is
// granted, provided the text of this NOTICE is retained, a notice that
// the code was modified is included with the above COPYRIGHT NOTICE and
// with the COPYRIGHT NOTICE in the LICENSE file, and that the LICENSE
// file is distributed with the modified code.
//
// LICENSOR MAKES NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED.
// By way of example, but not limitation, Licensor MAKES NO
// REPRESENTATIONS OR WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY
// PARTICULAR PURPOSE OR THAT THE USE OF THE LICENSED SOFTWARE COMPONENTS
// OR DOCUMENTATION WILL NOT INFRINGE ANY PATENTS, COPYRIGHTS, TRADEMARKS
// OR OTHER RIGHTS.
//
#ifndef MTL_LIGHT1D_H
#define MTL_LIGHT1D_H
#include "mtl/mtl_iterator.h"
#include "mtl/mtl_config.h"
#include "mtl/dense_iterator.h"
#include "mtl/reverse_iter.h"
#include "mtl/matrix_traits.h"
#include "mtl/scaled1D.h"
#include <stdlib.h>
namespace mtl {
/**
This is a {\em light} version of {\tt dense1D}.
It does no memory management (or reference counting)
and can only be used with pre-existing memory.
The purpose of {\tt light1D} is to be used
in the high performance kernels.
@memo Light 1-D Container
*/
template <class T, int NN = 0, int IND_OFFSET = 0>
class light1D {
typedef light1D self;
public:
enum { N = NN };
typedef light1D<int> IndexArray; /* JGS */
/**@name Type Definitions */
//@{
///
typedef dense_tag sparsity;
///
typedef scaled1D< light1D<T> > scaled_type;
///
typedef T value_type;
///
typedef T& reference;
///
typedef T* pointer;
///
typedef const T& const_reference;
///
typedef const T* const_pointer;
typedef int size_type;
typedef ptrdiff_t difference_type;
#if defined(_MSVCPP_)
///
typedef dense_iterator<T, 0, IND_OFFSET> iterator;
///
typedef dense_iterator<T, 1, IND_OFFSET> const_iterator;
///
/*
#elif defined( _MSVCPP7_ )
/// used std::_Ptrit in order to support iterator_traits for
/// pointers masquerading as iterators as per std::vector and std::basic_string - BEL
//
typedef std::_Ptrit<value_type, difference_type, pointer, reference, pointer, reference> ptr_iterator;
typedef std::_Ptrit<value_type, difference_type, const_pointer, const_reference, pointer, reference> ptr_const_iterator;
///
typedef dense_iterator<ptr_iterator, IND_OFFSET> iterator;
///
typedef dense_iterator<ptr_const_iterator, IND_OFFSET> const_iterator;
///
*/
#else
///
typedef dense_iterator<T*, IND_OFFSET> iterator;
///
typedef dense_iterator<const T*, IND_OFFSET> const_iterator;
#endif
///
typedef reverse_iter<iterator> reverse_iterator;
///
typedef reverse_iter<const_iterator> const_reverse_iterator;
typedef self IndexArrayRef;
typedef self subrange_type;
typedef oned_tag dimension;
//@}
/**@name Constructors */
//@{
/// Default Constructor
inline light1D() : rep(0), size_(0), first(0) { }
/// Preallocated Memory Constructor with optional non-zero starting index
inline light1D(T* data, size_type n, size_type start = 0)
: rep(data), size_(n), first(start) { }
/// Copy Constructor
inline light1D(const self& x)
: rep(x.rep), size_(x.size_), first(x.first) { }
inline ~light1D() { }
//@}
/**@name Access Methods */
//@{
/**@name Iterator Access Methods */
//@{
///
inline iterator begin() { return iterator(rep, 0, first); }
///
inline iterator end() { return iterator(rep, size_, first); }
///
inline const_iterator begin() const {
return const_iterator(rep, 0, first);
}
///
inline const_iterator end() const{
return const_iterator(rep, size_, first);
}
///
inline reverse_iterator rbegin() {
return reverse_iterator(end());
}
///
inline reverse_iterator rend() { return reverse_iterator(begin()); }
///
inline const_reverse_iterator rbegin() const {
return const_reverse_iterator(end());
}
///
inline const_reverse_iterator rend() const{
return const_reverse_iterator(begin());
}
//@}
/**@name Element Access Methods */
//@{
///
inline reference operator[](size_type n) { return rep[n - first]; }
///
inline const_reference operator[](size_type n) const {
return rep[n - first];
}
inline subrange_type operator()(size_type s, size_type f) const {
return subrange_type(rep + s - first, f - s, 0);
}
//@}
/**@name Size Methods */
//@{
///
inline int size() const { return size_; }
///
inline int nnz() const { return size_; }
///
inline void resize(int n) {
if (rep) delete [] rep;
size_ = n;
rep = new T[size_];
}
inline self& adjust_index(size_type delta) {
first += delta;
return *this;
}
/// Memory Access
inline T* data() const { return rep; }
protected:
T* rep;
int size_;
int first;
};
} /* namespace mtl */
#endif // MTL_LIGHT1D_H
linalg_vec.h
(application/octet-stream, 19 KB)
//
// Copyright 1997, 1998, 1999 University of Notre Dame.
// Authors: Andrew Lumsdaine, Jeremy G. Siek, Lie-Quan Lee
//
// This file is part of the Matrix Template Library
//
// You should have received a copy of the License Agreement for the
// Matrix Template Library along with the software; see the
// file LICENSE. If not, contact Office of Research, University of Notre
// Dame, Notre Dame, IN 46556.
//
// Permission to modify the code and to distribute modified code is
// granted, provided the text of this NOTICE is retained, a notice that
// the code was modified is included with the above COPYRIGHT NOTICE and
// with the COPYRIGHT NOTICE in the LICENSE file, and that the LICENSE
// file is distributed with the modified code.
//
// LICENSOR MAKES NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED.
// By way of example, but not limitation, Licensor MAKES NO
// REPRESENTATIONS OR WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY
// PARTICULAR PURPOSE OR THAT THE USE OF THE LICENSED SOFTWARE COMPONENTS
// OR DOCUMENTATION WILL NOT INFRINGE ANY PATENTS, COPYRIGHTS, TRADEMARKS
// OR OTHER RIGHTS.
//
//===========================================================================
#ifndef MTL_LINALG_VECTOR_H
#define MTL_LINALG_VECTOR_H
#include <utility>
#include <vector>
#include "mtl/refcnt_ptr.h"
#include "mtl/dense_iterator.h"
#include "mtl/reverse_iter.h"
#include "mtl/light1D.h"
#include "mtl/mtl_config.h"
#include "mtl/matrix_traits.h"
#include "mtl/scaled1D.h"
#include "mtl/mtl_exception.h"
#include "mtl/external_vector.h"
namespace mtl {
//: Linalg Vector Adaptor
//!category: containers, adaptors
//!component: type
//
// This captures the main functionality of a dense MTL vector. The
// dense1D and external1D derive from this class, and specialize
// this class to use either internal or external storage.
//
//!definition: linalg_vector.h
//!tparam: RepType - the underlying representation
//!models: Linalg_Vector
template <class RepType, class RepPtr = RepType*, int NN = 0>
class linalg_vec {
public:
typedef linalg_vec self;
typedef RepType rep_type;
typedef RepPtr rep_ptr;
enum { N = NN };
/**@name Type Definitions */
/* enum { dimension = 1 }; */
typedef oned_tag dimension;
//: The sparsity tag
typedef dense_tag sparsity;
//: The scaled type of this container
//!wheredef: Scalable
typedef scaled1D< self > scaled_type;
//: The value type
//!wheredef: Container
typedef typename rep_type::value_type value_type;
//: The reference type
//!wheredef: Container
typedef typename rep_type::reference reference;
//: The const reference type
//!wheredef: Container
typedef typename rep_type::const_reference const_reference;
//: The pointer (to the value_type) type
//!wheredef: Container
typedef typename rep_type::pointer pointer;
//: The size type (non negative)
//!wheredef: Container
typedef typename rep_type::size_type size_type;
//: The difference type (an integral type)
//!wheredef: Container
typedef typename rep_type::difference_type difference_type;
#if !defined( _MSVCPP_ )
//: The iterator type
//!wheredef: Container
typedef dense_iterator<typename rep_type::iterator> iterator;
//: The const iterator type
//!wheredef: Container
typedef dense_iterator<typename rep_type::const_iterator> const_iterator;
#else
typedef dense_iterator<typename rep_type::value_type, 0, 0, size_type> iterator;
typedef dense_iterator<typename rep_type::value_type, 1, 0, size_type> const_iterator;
#endif
//: The reverse iterator type
//!wheredef: Reversible Container
typedef reverse_iter<iterator> reverse_iterator;
//: The const reverse iterator type
//!wheredef: Reversible Container
typedef reverse_iter<const_iterator> const_reverse_iterator;
/* skip over the zeros and report the indices
this implements the nonzero structure array
*/
typedef linalg_vec<RepType, RepPtr, NN> Vec;
typedef size_type Vec_size_type;
typedef difference_type Vec_difference_type;
typedef iterator Vec_iterator;
typedef const_iterator Vec_const_iterator;
class IndexArray {
public:
typedef Vec_size_type size_type;
typedef Vec_difference_type difference_type;
typedef Vec_size_type value_type;
class iterator {
public:
typedef size_type value_type;
typedef size_type reference;
typedef size_type* pointer;
typedef Vec_difference_type difference_type;
typedef typename std::iterator_traits<Vec_iterator>::iterator_category iterator_category;
iterator(Vec_iterator iter, Vec_iterator e) : i(iter), end(e) {
while (*i == Vec_value_type(0)) ++i;
}
reference operator*() const { return i.index(); }
iterator& operator++() {
++i; while (*i == Vec_value_type(0) && i != end) ++i;
return *this; }
iterator operator++(int) { iterator t = *this; ++(*this); return t; }
iterator& operator--() {
--i; while (*i == Vec_value_type(0) && i != end) --i;
return *this; }
iterator operator--(int) { iterator t = *this; --(*this); return t; }
difference_type operator-(const iterator& x) const { return i - x.i; }
bool operator==(const iterator& x) const { return i == x.i; }
bool operator!=(const iterator& x) const { return i != x.i; }
bool operator<(const iterator& x) const { return i < x.i; }
Vec_iterator i;
Vec_iterator end;
};
class const_iterator {
public:
typedef size_type value_type;
typedef size_type reference;
typedef size_type* pointer;
typedef Vec_difference_type difference_type;
typedef typename std::iterator_traits<Vec_iterator>::iterator_category iterator_category;
const_iterator(Vec_const_iterator iter, Vec_const_iterator e)
: i(iter), end(e) {
while (*i == Vec_value_type(0) && i != end) ++i;
}
reference operator*() const { return i.index(); }
const_iterator& operator++() {
++i; while (*i == Vec_value_type(0) && i != end) ++i;
return *this; }
const_iterator operator++(int) {
const_iterator t = *this; ++(*this); return t; }
const_iterator& operator--() {
--i; while (*i == Vec_value_type(0)) --i;
return *this; }
const_iterator operator--(int) {
const_iterator t = *this; --(*this); return t; }
difference_type operator-(const const_iterator& x) const {
return i - x.i; }
bool operator==(const const_iterator& x) const { return i == x.i; }
bool operator!=(const const_iterator& x) const { return i != x.i; }
bool operator<(const const_iterator& x) const { return i < x.i; }
Vec_const_iterator i;
Vec_const_iterator end;
};
inline IndexArray(const Vec& v) : vec((Vec*)&v) { }
inline iterator begin() { return iterator(vec->begin(), vec->end()); }
inline iterator end() { return iterator(vec->end(), vec->end()); }
inline const_iterator begin() const{
return const_iterator(((const Vec*)vec)->begin(),
((const Vec*)vec)->end());
}
inline const_iterator end() const {
return const_iterator(((const Vec*)vec)->end(),
((const Vec*)vec)->end()); }
size_type size() const {
size_type s = 0;
Vec_const_iterator i;
for (i = ((const Vec*)vec)->begin(); i != ((const Vec*)vec)->end(); ++i)
if (*i != Vec_value_type(0)) ++s;
return s;
}
Vec* vec;
};
//: The type for an array of the indices of the element in the vector
//!wheredef: Vector
typedef IndexArray IndexArrayRef;
//: The type for a subrange vector-view of the original vector
//!wheredef: Vector
typedef light1D<value_type> subrange_type;
typedef std::pair<size_type, size_type> range;
/**@name Constructors */
//: Default Constructor (allocates the container)
//!wheredef: Container
inline linalg_vec() : rep(0) { }
//: Normal Constructor
inline linalg_vec(rep_ptr x, size_type start_index)
: rep(x), first(start_index) { }
//: Copy Constructor (shallow copy)
//!wheredef: ContainerRef
inline linalg_vec(const self& x) : rep(x.rep), first(x.first) { }
//: The destructor.
//!wheredef: Container
inline ~linalg_vec() { }
//: Assignment Operator (shallow copy)
//!wheredef: AssignableRef
inline self& operator=(const self& x) {
rep = x.rep;
first = x.first;
return *this;
}
/**@name Access Methods */
/**@name Iterator Access Methods */
//: Return an iterator pointing to the beginning of the vector
//!wheredef: Container
inline iterator begin() { return iterator(rep->begin(), 0, first); }
//: Return an iterator pointing past the end of the vector
//!wheredef: Container
inline iterator end() { return iterator(rep->begin(), rep->size(), first); }
//: Return a const iterator pointing to the begining of the vector
//!wheredef: Container
inline const_iterator begin() const { return const_iterator(rep->begin(),
0, first); }
//: Return a const iterator pointing past the end of the vector
//!wheredef: Container
inline const_iterator end() const{ return const_iterator(rep->begin(),
rep->size(), first); }
//: Return a reverse iterator pointing to the last element of the vector
//!wheredef: Reversible Container
inline reverse_iterator rbegin() { return reverse_iterator(end()); }
//: Return a reverse iterator pointing past the end of the vector
//!wheredef: Reversible Container
inline reverse_iterator rend() { return reverse_iterator(begin()); }
//: Return a const reverse iterator pointing to the last element of the vector
//!wheredef: Reversible Container
inline const_reverse_iterator rbegin() const {
return reverse_iterator(end()); }
//: Return a const reverse iterator pointing past the end of the vector
//!wheredef: Reversible Container
inline const_reverse_iterator rend() const{
return reverse_iterator(begin()); }
/**@name Element Access Methods */
//: Return a reference to the element with the ith index
//!wheredef: Vector
inline reference operator[](size_type i) MTL_THROW_ASSERTION {
MTL_ASSERT(i < size(), "linalg_vec::operator[]");
return (*rep)[i - first];
}
inline subrange_type operator()(range r) MTL_THROW_ASSERTION {
return subrange_type(data() + r.first, r.second - r.first);
}
inline subrange_type operator()(size_type s, size_type f)
MTL_THROW_ASSERTION
{
return subrange_type(data() + s, f - s);
}
//: Return a const reference to the element with the ith index
//!wheredef: Vector
inline const_reference operator[](size_type i) const MTL_THROW_ASSERTION {
MTL_ASSERT(i < size(), "linalg_vec::operator[]");
return (*rep)[i - first];
}
/**@name Size Methods */
//: The size of the vector
//!wheredef: Container
inline size_type size() const { return rep->size(); }
//: The number of non-zeroes in the vector
inline size_type nnz() const { return rep->size(); }
//: Resize the vector to n
inline void resize(size_type n) { rep->resize(n); }
//: Resize the vector to n, and assign x to the new positions
inline void resize(size_type n, const value_type& x) { rep->resize(n, x); }
//: Return the total capacity of the vector
size_type capacity() const { return rep->capacity(); }
//: Reserve more space in the vector
void reserve(size_type n) { rep->reserve(n); }
//: Raw Memory Access
inline const value_type* data() const { return &(*rep)[0]; }
//: Raw Memory Access
inline value_type* data() { return &(*rep)[0]; }
//: Insert x at the indicated position in the vector
//!wheredef: Container
inline iterator
insert (iterator position, const value_type& x = value_type()) {
return iterator(rep->insert(position.base(), x), position.index()+1);
/* JGS, not sure about what to do with the index here */
}
inline IndexArrayRef nz_struct() const { return IndexArrayRef(*this); }
inline self& adjust_index(size_type i) {
first += i;
return *this;
}
protected:
rep_ptr rep;
size_type first;
};
//: External 1-D Container
//!category: containers
//!component: type
//
// This is similar to dense1D, except that the memory is provided
// by the user. This allows for interoperability with other array
// packages and even with Fortran.
//
//!definition: linalg_vec.h
//!tparam: T - The element type.
//!tparam: NN - The static size of the Vector, 0 if dynamic size
//!tparam: SizeT - The size type to use - size_t
//!tparam: DiffT - ptrdiff_t
//!models: Vector
//!example: dot_prod.cc, apply_givens.cc, euclid_norm.cc, max_index.cc
template <class T, int NN = 0, class SizeType=unsigned int>
class external_vec {
typedef external_vec self;
public:
enum { N = NN };
typedef external_vec<int> IndexArray; /* JGS */
/* Type Definitions */
//: The vector is dense
typedef dense_tag sparsity;
//: Scaled type for the vector
typedef scaled1D< self > scaled_type;
typedef SizeType size_type;
typedef int difference_type;
//: The element type
typedef T value_type;
//: The reference to the value type
typedef T& reference;
//: The pointer ot the value type
typedef T* pointer;
//: The const reference type
typedef const T& const_reference;
//: The const pointer to the value type
typedef const T* const_pointer;
#if defined( _MSVCPP_ )
typedef dense_iterator<T, 0, 0, size_type> iterator;
typedef dense_iterator<T, 1, 0, size_type> const_iterator;
/*
#elif defined ( _MSVCPP7_ )
/// used std::_Ptrit in order to support iterator_traits for
/// pointers masquerading as iterators as per std::vector and std::basic_string - BEL
//
typedef std::_Ptrit<value_type, difference_type, pointer, reference, pointer, reference> ptr_iterator;
typedef std::_Ptrit<value_type, difference_type, const_pointer, const_reference, pointer, reference> ptr_const_iterator;
typedef dense_iterator<ptr_iterator,0,size_type> iterator;
typedef dense_iterator<ptr_const_iterator,0,size_type> const_iterator;
*/
#else
typedef dense_iterator<T*,0,size_type> iterator;
typedef dense_iterator<const T*,0,size_type> const_iterator;
#endif
//: The reverse iterator type
typedef reverse_iter<iterator> reverse_iterator;
//: The const reverse iterator type
typedef reverse_iter<const_iterator> const_reverse_iterator;
//:
//!wheredef: Vector
typedef self IndexArrayRef;
//: The type for the subrange vector
//!wheredef: Vector
typedef self subrange_type;
typedef std::pair<size_type, size_type> range;
//: This is a 1D container
typedef oned_tag dimension;
/* Constructors */
//: Default Constructor
inline external_vec() : rep(0), size_(0), first(0) { }
//: External Data Contructor
inline external_vec(T* data)
: rep(data), size_(N), first(0) { }
//: Preallocated Memory Constructor with optional non-zero starting index
inline external_vec(T* data, size_type n, size_type start = 0)
: rep(data), size_(n), first(start) { }
//: Copy Constructor
inline external_vec(const self& x)
: rep(x.rep), size_(x.size_), first(x.first) { }
//: Assignment
inline self& operator=(const self& x) {
rep = x.rep; size_ = x.size_; first = x.first; return *this;
}
//: Destructor
inline ~external_vec() { }
/* Access Methods */
/* Iterator Access Methods */
//: Return an iterator pointing to the beginning of the vector
//!wheredef: Container
inline iterator begin() { return iterator(rep, 0, first); }
//: Return an iterator pointing past the end of the vector
//!wheredef: Container
inline iterator end() { return iterator(rep, size(), first); }
//: Return a const iterator pointing to the begining of the vector
//!wheredef: Container
inline const_iterator begin() const {
return const_iterator(rep, 0, first);
}
//: Return a const iterator pointing past the end of the vector
//!wheredef: Container
inline const_iterator end() const{
return const_iterator(rep, size(), first);
}
//: Return a reverse iterator pointing to the last element of the vector
//!wheredef: Reversible Container
inline reverse_iterator rbegin() {
return reverse_iterator(end());
}
//: Return a reverse iterator pointing past the end of the vector
//!wheredef: Reversible Container
inline reverse_iterator rend() { return reverse_iterator(begin()); }
//: Return a const reverse iterator pointing to the last element of the vector
//!wheredef: Reversible Container
inline const_reverse_iterator rbegin() const {
return const_reverse_iterator(end());
}
//: Return a const reverse iterator pointing past the end of the vector
//!wheredef: Reversible Container
inline const_reverse_iterator rend() const{
return const_reverse_iterator(begin());
}
/* Element Access Methods */
//: Return a reference to the element with the ith index
//!wheredef: Vector
inline reference operator[](size_type i) { return rep[i - first]; }
//: Return a const reference to the element with the ith index
//!wheredef: Vector
inline const_reference operator[](size_type i) const { return rep[i - first]; }
//: Return a subrange vector with start at s and finish at f
//!wheredef: Vector
inline subrange_type operator()(size_type s, size_type f) const {
return subrange_type(rep + s - first, f - s, 0);
}
inline subrange_type operator()(range r) MTL_THROW_ASSERTION {
return subrange_type(data() + r.first, r.second - r.first, 0);
}
/* Size Methods */
//: The size of the vector
//!wheredef: Container
inline size_type size() const { return N ? N : size_; }
//: The number of non-zeroes in the vector
//!wheredef: Vector
inline size_type nnz() const { return size(); }
//: Resize the vector to size n
#if 0
inline void resize(size_type n) {
if (rep) delete [] rep;
size_ = n;
rep = new T[size_];
}
#else
inline void resize(size_type n) { size_ = n; }
inline void clear() { size_ = 0; }
#endif
//: Raw Memory Access
inline value_type* data() const { return rep; }
inline self& adjust_index(size_type i) {
first += i;
return *this;
}
//: Push x onto the end of the vector, increasing the size
// This function does not allocation memory.
// Better hope enough memory is already there!
void push_back(const T& x) {
rep[size_] = x;
++size_;
}
protected:
T* rep;
size_type size_;
size_type first;
};
//: blah
//!noindex:
template <int N>
struct __make_external {
template <class T>
inline external_vec<T,N> operator()(T* x) {
return external_vec<T,N>(x);
}
};
/* For converting static arrays into MTL vectors */
#define array_to_vec(x) mtl::__make_external<sizeof(x)/sizeof(*x)>()(x)
template <class Container>
inline linalg_vec<Container>
vec(const Container& x)
{
return linalg_vec<Container>(x);
}
} /* namespace mtl */
#endif
dense_iterator.h
(application/octet-stream, 15.5 KB)
//
// Copyright 1997, 1998, 1999 University of Notre Dame.
// Authors: Andrew Lumsdaine, Jeremy G. Siek, Lie-Quan Lee
//
// This file is part of the Matrix Template Library
//
// You should have received a copy of the License Agreement for the
// Matrix Template Library along with the software; see the
// file LICENSE. If not, contact Office of Research, University of Notre
// Dame, Notre Dame, IN 46556.
//
// Permission to modify the code and to distribute modified code is
// granted, provided the text of this NOTICE is retained, a notice that
// the code was modified is included with the above COPYRIGHT NOTICE and
// with the COPYRIGHT NOTICE in the LICENSE file, and that the LICENSE
// file is distributed with the modified code.
//
// LICENSOR MAKES NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED.
// By way of example, but not limitation, Licensor MAKES NO
// REPRESENTATIONS OR WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY
// PARTICULAR PURPOSE OR THAT THE USE OF THE LICENSED SOFTWARE COMPONENTS
// OR DOCUMENTATION WILL NOT INFRINGE ANY PATENTS, COPYRIGHTS, TRADEMARKS
// OR OTHER RIGHTS.
//
//
//
//===========================================================================
#ifndef MTL_DENSE_ITERATOR_H
#define MTL_DENSE_ITERATOR_H
#include "mtl/mtl_iterator.h"
#include "mtl/meta_if.h"
#include "mtl/mtl_config.h"
namespace mtl {
//: dense iterator
//
// An iterator for dense contiguous container that keeps track of the index.
//
//!category: iterators, adaptors
//!component: type
//!definition:dense_iterator.h
//!tparam: RandomAccessIterator - the base iterator
//!models: RandomAccessIterator? (with index())
#if defined ( _MSVCPP_ )
struct _bogus { };
/* The inheritance from Ranit is a VC++ workaround for not having
a working iterator traits */
template <class T, int isConst, int IND_OFFSET=0, class SizeType=int>
class dense_iterator : public std::_Ranit<T,SizeType> {
typedef dense_iterator self;
typedef typename IF<isConst, const T*, T*>::RET RandomAccessIterator;
public:
//: The value type
typedef T value_type;
//: This is a random access iterator
typedef std::random_access_iterator_tag iterator_category;
//: The type for differences between iterators
typedef int difference_type;
typedef int distance_type;
//: The type for pointers to the value type
typedef IF<isConst, const T*, T*>::RET pointer;
//: The type for references to the value type
typedef IF<isConst, const T&, T&>::RET reference;
typedef difference_type Distance;
typedef SizeType size_type;
/*
protected:
*/
RandomAccessIterator start;
size_type pos;
size_type start_index;
public:
//: Return the index of the current element
//!wheredef: IndexedIterator
inline size_type index() const {
return pos + start_index + IND_OFFSET;
}
//: Default Constructor
inline dense_iterator() : pos(0), start_index(0) {}
//: Constructor from underlying iterator
inline dense_iterator(RandomAccessIterator s, size_type i, size_type first_index = 0)
: start(s), pos(i), start_index(first_index) { }
//: Copy Constructor
inline dense_iterator (const self& x)
: start(x.start), pos(x.pos), start_index(x.start_index) {}
typedef typename IF<isConst, dense_iterator<T,0,IND_OFFSET,SizeType>,
_bogus >::RET NonConst;
inline dense_iterator(const NonConst& x)
: start(x.start), pos(x.pos), start_index(x.start_index) {}
//: Assignment operator
inline self& operator=(const self& x) {
start = x.start;
pos = x.pos;
start_index = x.start_index;
return *this;
}
//: Destructor
inline ~dense_iterator () { }
//: Access the underlying iterator
inline RandomAccessIterator base() const { return start + pos; }
inline operator RandomAccessIterator() const { return start + pos; }
//: Dereference operator
inline reference operator*() const { return *(start + pos); }
//: Member access operator
inline pointer operator-> () const { return start + pos; }
//: Pre-increment operator
inline self& operator++ () { ++pos; return *this; }
//: Post-increment operator
inline self operator++ (int) { self tmp = *this; ++pos; return tmp; }
//: Pre-decrement operator
inline self& operator-- () { --pos; return *this; }
//: Post-decrement operator
inline self operator-- (int) { self tmp = *this; --pos; return tmp; }
//: Add iterator and distance n
inline self operator+ (Distance n) const { return self(start, pos + n); }
//: Add distance n to this iterator
inline self& operator+= (Distance n) { pos += n; return *this; }
//: Subtract iterator and distance n
inline self operator- (Distance n) const { return self(start, pos - n); }
//: Return the difference between two iterators
inline difference_type operator- (const self& x) const {
return base() - x.base(); }
//: Subtract distance n from this iterator
inline self& operator-= (Distance n) { pos -= n; return *this; }
//: Return whether this iterator is not equal to iterator x
inline bool operator!= (const self& x) const { return pos != x.pos; }
//: Return whether this iterator is less than iterator x
inline bool operator < (const self& x) const { return pos < x.pos; }
//: Return whether this iterator is greater than iterator x
inline bool operator > (const self& x) const { return pos > x.pos; }
//: Return whether this iterator is equal to iterator x
inline bool operator== (const self& x) const { return pos == x.pos; }
//: Return whether this iterator is less than or equal to iterator x
inline bool operator<= (const self& x) const { return pos <= x.pos; }
//: Return whether this iterator is greater than or equal to iterator x
inline bool operator>= (const self& x) const { return pos >= x.pos; }
//: Equivalent to *(i + n)
#if 0 /* caused ambiguity with + op for VC++ */
inline reference operator[] (Distance n) const {
return *(start + pos + n);
}
#endif
};
template <class T, int isConst,int OS, class ST>
inline
dense_iterator<T,isConst,OS,ST>
operator+ (typename dense_iterator<T,isConst,OS,ST>::size_type n,
const dense_iterator<T,isConst,OS,ST> &x)
{
return dense_iterator<T,isConst,OS,ST>(x.base(), n);
}
/*
#elif defined( _MSVCPP7_ )
// while the iterator_traits behaviour has been fixed with MSVC 7, it still has
// a problem deducing template arguments in some instances. It works if dense_iterator
// is derived from rather than containing RandomAccessIterator. - BEL
template <class RandomAccessIterator, int IND_OFFSET=0, class SizeType=int>
class dense_iterator : public RandomAccessIterator {
typedef dense_iterator self;
public:
//: The value type
typedef typename std::iterator_traits<RandomAccessIterator>::value_type value_type;
//: This is a random access iterator
typedef typename std::iterator_traits<RandomAccessIterator>::iterator_category iterator_category;
//: The type for differences between iterators
typedef typename std::iterator_traits<RandomAccessIterator>::difference_type difference_type;
//: The type for pointers to the value type
typedef typename std::iterator_traits<RandomAccessIterator>::pointer pointer;
//: The type for references to the value type
typedef typename std::iterator_traits<RandomAccessIterator>::reference reference;
typedef difference_type Distance;
typedef SizeType size_type;
//protected:
// if we derive from RandomAccessIterator we shouldn't embed another
// I implemented start() to make this easier to compare to the non-MSVC version - BEL
//
inline RandomAccessIterator& start() { return *this; }
inline const RandomAccessIterator& start() const { return *this; }
size_type pos;
size_type start_index;
public:
//: Return the index of the current element
//!wheredef: IndexedIterator
inline size_type index() const {
return pos + start_index + IND_OFFSET;
}
//: Default Constructor
inline dense_iterator() : pos(0), start_index(0) {}
//: Constructor from underlying iterator
inline dense_iterator(RandomAccessIterator s,
size_type i, size_type first_index = 0)
: RandomAccessIterator(s), pos(i), start_index(first_index) { }
//: Copy Constructor
inline dense_iterator (const self& x)
: RandomAccessIterator(x), pos(x.pos), start_index(x.start_index) {}
template <class SELF>
inline dense_iterator (const SELF& x)
: RandomAccessIterator(x), pos(x.pos), start_index(x.start_index) {}
//: Assignment operator
inline self& operator=(const self& x) {
RandomAccessIterator::operator=( x );
pos = x.pos;
start_index = x.start_index;
return *this;
}
//: Destructor
inline ~dense_iterator () { }
//: Access the underlying iterator
inline RandomAccessIterator base() const { return start() + pos; }
inline operator RandomAccessIterator() const { return start() + pos; }
//: Dereference operator
inline reference operator*() const { return *(start() + pos); }
//: Member access operator
inline pointer operator-> () const { return start() + pos; }
//: Pre-increment operator
inline self& operator++ () { ++pos; return *this; }
//: Post-increment operator
inline self operator++ (int) { self tmp = *this; ++pos; return tmp; }
//: Pre-decrement operator
inline self& operator-- () { --pos; return *this; }
//: Post-decrement operator
inline self operator-- (int) { self tmp = *this; --pos; return tmp; }
//: Add iterator and distance n
inline self operator+ (Distance n) const { return self(start(), pos + n); }
//: Add distance n to this iterator
inline self& operator+= (Distance n) { pos += n; return *this; }
//: Subtract iterator and distance n
inline self operator- (Distance n) const { return self(start(), pos - n); }
//: Return the difference between two iterators
inline difference_type operator- (const self& x) const {
return base() - x.base(); }
//: Subtract distance n from this iterator
inline self& operator-= (Distance n) { pos -= n; return *this; }
//: Return whether this iterator is not equal to iterator x
inline bool operator!= (const self& x) const { return pos != x.pos; }
//: Return whether this iterator is less than iterator x
inline bool operator < (const self& x) const { return pos < x.pos; }
//: Return whether this iterator is greater than iterator x
inline bool operator > (const self& x) const { return pos > x.pos; }
//: Return whether this iterator is equal to iterator x
inline bool operator== (const self& x) const { return pos == x.pos; }
//: Return whether this iterator is less than or equal to iterator x
inline bool operator<= (const self& x) const { return pos <= x.pos; }
//: Return whether this iterator is greater than or equal to iterator x
inline bool operator>= (const self& x) const { return pos >= x.pos; }
//: Equivalent to *(i + n)
inline reference operator[] (Distance n) const {
return *(start() + pos + n);
}
};
template <class T, int OS, class ST>
inline
dense_iterator<T>
operator+ (typename dense_iterator<T,OS,ST>::size_type n,
const dense_iterator<T,OS,ST> &x)
{
return dense_iterator<T,OS,ST>(x.base(), n);
}
*/
#else // other compilers
template <class RandomAccessIterator, int IND_OFFSET=0, class SizeType=int>
class dense_iterator {
typedef dense_iterator self;
public:
//: The value type
typedef typename std::iterator_traits<RandomAccessIterator>::value_type value_type;
//: This is a random access iterator
typedef typename std::iterator_traits<RandomAccessIterator>::iterator_category iterator_category;
//: The type for differences between iterators
typedef typename std::iterator_traits<RandomAccessIterator>::difference_type difference_type;
//: The type for pointers to the value type
typedef typename std::iterator_traits<RandomAccessIterator>::pointer pointer;
//: The type for references to the value type
typedef typename std::iterator_traits<RandomAccessIterator>::reference reference;
typedef difference_type Distance;
typedef SizeType size_type;
/*
protected:
*/
RandomAccessIterator start;
size_type pos;
size_type start_index;
public:
//: Return the index of the current element
//!wheredef: IndexedIterator
inline size_type index() const {
return pos + start_index + IND_OFFSET;
}
//: Default Constructor
inline dense_iterator() : pos(0), start_index(0) {}
//: Constructor from underlying iterator
inline dense_iterator(RandomAccessIterator s,
size_type i, size_type first_index = 0)
: start(s), pos(i), start_index(first_index) { }
//: Copy Constructor
inline dense_iterator (const self& x)
: start(x.start), pos(x.pos), start_index(x.start_index) {}
template <class SELF>
inline dense_iterator (const SELF& x)
: start(x.start), pos(x.pos), start_index(x.start_index) {}
//: Assignment operator
inline self& operator=(const self& x) {
start = x.start;
pos = x.pos;
start_index = x.start_index;
return *this;
}
//: Destructor
inline ~dense_iterator () { }
//: Access the underlying iterator
inline RandomAccessIterator base() const { return start + pos; }
inline operator RandomAccessIterator() const { return start + pos; }
//: Dereference operator
inline reference operator*() const { return *(start + pos); }
//: Member access operator
inline pointer operator-> () const { return start + pos; }
//: Pre-increment operator
inline self& operator++ () { ++pos; return *this; }
//: Post-increment operator
inline self operator++ (int) { self tmp = *this; ++pos; return tmp; }
//: Pre-decrement operator
inline self& operator-- () { --pos; return *this; }
//: Post-decrement operator
inline self operator-- (int) { self tmp = *this; --pos; return tmp; }
//: Add iterator and distance n
inline self operator+ (Distance n) const { return self(start, pos + n); }
//: Add distance n to this iterator
inline self& operator+= (Distance n) { pos += n; return *this; }
//: Subtract iterator and distance n
inline self operator- (Distance n) const { return self(start, pos - n); }
//: Return the difference between two iterators
inline difference_type operator- (const self& x) const {
return base() - x.base(); }
//: Subtract distance n from this iterator
inline self& operator-= (Distance n) { pos -= n; return *this; }
//: Return whether this iterator is not equal to iterator x
inline bool operator!= (const self& x) const { return pos != x.pos; }
//: Return whether this iterator is less than iterator x
inline bool operator < (const self& x) const { return pos < x.pos; }
//: Return whether this iterator is greater than iterator x
inline bool operator > (const self& x) const { return pos > x.pos; }
//: Return whether this iterator is equal to iterator x
inline bool operator== (const self& x) const { return pos == x.pos; }
//: Return whether this iterator is less than or equal to iterator x
inline bool operator<= (const self& x) const { return pos <= x.pos; }
//: Return whether this iterator is greater than or equal to iterator x
inline bool operator>= (const self& x) const { return pos >= x.pos; }
//: Equivalent to *(i + n)
inline reference operator[] (Distance n) const {
return *(start + pos + n);
}
};
template <class T, int OS, class ST>
inline
dense_iterator<T>
operator+ (typename dense_iterator<T,OS,ST>::size_type n,
const dense_iterator<T,OS,ST> &x)
{
return dense_iterator<T,OS,ST>(x.base(), n);
}
#endif
} /* namespace mtl */
#endif