---

more/num/matrixalgo.h

Go to the documentation of this file.

//  Copyright (C) 1998--2002  Petter Urkedal
//
//  This file is free software; you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation; either version 2 of the License, or
//  (at your option) any later version.
//
//  This file is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.
//
//  You should have received a copy of the GNU General Public License
//  along with this program; if not, write to the Free Software
//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//
//  As a special exception, you may use this file as part of a free
//  software library without restriction.  Specifically, if other files
//  instantiate templates or use macros or inline functions from this
//  file, or you compile this file and link it with other files to
//  produce an executable, this file does not by itself cause the
//  resulting executable to be covered by the GNU General Public
//  License.  This exception does not however invalidate any other
//  reasons why the executable file might be covered by the GNU General
//  Public License.

//  $Id: matrixalgo.h,v 1.2 2002/09/01 12:02:10 petter_urkedal Exp $


#ifndef MORE_NUM_MATRIXALGO_H
#define MORE_NUM_MATRIXALGO_H

///\if unmaintained
#include <iterator>
#include <algorithm>
#include <numeric>
#include <stdexcept>
#include <functional>
#include <more/math/math.h>
#include <more/bits/optype.h>
#include <more/num/matrix.h>


namespace more {
namespace num {

    // --- matrix ---

    template<typename T, typename Sig>
    inline matrix<T, Sig>
    operator+(const matrix<T, Sig>& A, const matrix<T, Sig>& B) {
        assert(A.size(0) == B.size(0) && A.size(1) == B.size(1));
        matrix<T, Sig> C(A.size(0), A.size(1));
        std::transform(A.begin(), A.end(), B.begin(), C.begin(),
                       std::plus<T>());
        return C;
    }
    template<typename T, typename Sig>
    inline matrix<T, Sig>
    operator-(const matrix<T, Sig>& A, const matrix<T, Sig>& B) {
        assert(A.size(0) == B.size(0) && A.size(1) == B.size(1));
        matrix<T, Sig> C(A.size(0), A.size(1));
        std::transform(A.begin(), A.end(), B.begin(), C.begin(),
                       std::minus<T>());
        return C;
    }
    template<typename T> inline
    typename scalar_type_<T>::eval trace(matrix<T> const& A) {
        typename scalar_type_<T>::eval tr = 0;
        assert(A.size(0) == A.size(1));
        for(int i = 0; i < A.size(0); ++i) tr += trace(A(i,i));
        return tr;
    }


    // --- mvector ---

    template<typename T>
    inline typename scalar_type_<T>::eval
    dot(mvector<T> const& A, mvector<T> const& B) {
        typedef typename mvector<T>::const_iterator iterator;
        if(A.size() != B.size())
            throw std::out_of_range(
                "dot product applied to vectors with different sizes");
        typename scalar_type_<T>::eval sum = zero;
        for(iterator u = A.begin(), v = B.begin(), uend = A.end();
            u != uend; ++u, ++v) sum += dot(*u, *v);
        return sum;
    }
//     template<typename T> struct norm_type_< mvector<T> >
//      { typedef norm_type_<T>::eval eval; };

    template<typename T>
    inline typename norm_type_<T>::eval norm(const mvector<T>& A) {
        typedef typename mvector<T>::const_iterator iterator;
        typename norm_type_<T>::eval sum = zero;
        for(iterator u = A.begin(), uend = A.end();
            u != uend; ++u) sum += norm(*u);
        return sum;
    }

    template<typename T>
    inline typename norm_type_<T>::eval
    abs(const mvector<T>& A) { return more::sqrt(norm(A)); }

    template<typename T> inline T sum(const mvector<T>& x)
        { return std::accumulate(x.begin(), x.end(), (T)0); }

    template<typename T> inline
    void normalize(mvector<T>& f) {
        norm_type_<T>::eval c = 0;
        for(int s = 0; s < f.size(); ++s) c += norm(f(s));
        c = 1/sqrt(c);
        for(int s = 0; s < f.size(); ++s) f(s) *= c;
    }

    template<typename T> inline mvector<T>
    operator-(const mvector<T>& A, const mvector<T>& B) {
        assert(A.size(0) == B.size(0));
        mvector<T> C(A.size(0));
        std::transform(A.begin(), A.end(), B.begin(), C.begin(),
                       std::minus<T>());
        return C;
    }

    template<typename T> inline mvector<T>
    operator+(const mvector<T>& A, const mvector<T>& B) {
        assert(A.size(0) == B.size(0));
        mvector<T> C(A.size(0));
        std::transform(A.begin(), A.end(), B.begin(), C.begin(),
                       std::plus<T>());
        return C;
    }

    template<typename T, typename U>
    inline mvector< typename times_type_<T, U>::eval >
    operator*(const T& a, const mvector<U>& B) {
        typedef mvector< typename times_type_<T, U>::eval > vec_type;
        vec_type C(B.size(0));
        typename mvector<U>::const_iterator v = B.begin();
        for(typename vec_type::iterator u = C.begin(); u != C.end(); ++u) {
            *u = a * (*v);
            ++v;
        }
        return C;
    }
}} // namespace more::num
///\endif

#endif

---