Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Related Pages  

more/num/matrixalgo.h

Go to the documentation of this file.
00001 //  Copyright (C) 1998--2002  Petter Urkedal
00002 //
00003 //  This file is free software; you can redistribute it and/or modify
00004 //  it under the terms of the GNU General Public License as published by
00005 //  the Free Software Foundation; either version 2 of the License, or
00006 //  (at your option) any later version.
00007 //
00008 //  This file is distributed in the hope that it will be useful,
00009 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00011 //  GNU General Public License for more details.
00012 //
00013 //  You should have received a copy of the GNU General Public License
00014 //  along with this program; if not, write to the Free Software
00015 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00016 //
00017 //  As a special exception, you may use this file as part of a free
00018 //  software library without restriction.  Specifically, if other files
00019 //  instantiate templates or use macros or inline functions from this
00020 //  file, or you compile this file and link it with other files to
00021 //  produce an executable, this file does not by itself cause the
00022 //  resulting executable to be covered by the GNU General Public
00023 //  License.  This exception does not however invalidate any other
00024 //  reasons why the executable file might be covered by the GNU General
00025 //  Public License.
00026 
00027 //  $Id: matrixalgo.h,v 1.2 2002/09/01 12:02:10 petter_urkedal Exp $
00028 
00029 
00030 #ifndef MORE_NUM_MATRIXALGO_H
00031 #define MORE_NUM_MATRIXALGO_H
00032 
00033 ///\if unmaintained
00034 #include <iterator>
00035 #include <algorithm>
00036 #include <numeric>
00037 #include <stdexcept>
00038 #include <functional>
00039 #include <more/math/math.h>
00040 #include <more/bits/optype.h>
00041 #include <more/num/matrix.h>
00042 
00043 
00044 namespace more {
00045 namespace num {
00046 
00047     // --- matrix ---
00048 
00049     template<typename T, typename Sig>
00050     inline matrix<T, Sig>
00051     operator+(const matrix<T, Sig>& A, const matrix<T, Sig>& B) {
00052         assert(A.size(0) == B.size(0) && A.size(1) == B.size(1));
00053         matrix<T, Sig> C(A.size(0), A.size(1));
00054         std::transform(A.begin(), A.end(), B.begin(), C.begin(),
00055                        std::plus<T>());
00056         return C;
00057     }
00058     template<typename T, typename Sig>
00059     inline matrix<T, Sig>
00060     operator-(const matrix<T, Sig>& A, const matrix<T, Sig>& B) {
00061         assert(A.size(0) == B.size(0) && A.size(1) == B.size(1));
00062         matrix<T, Sig> C(A.size(0), A.size(1));
00063         std::transform(A.begin(), A.end(), B.begin(), C.begin(),
00064                        std::minus<T>());
00065         return C;
00066     }
00067     template<typename T> inline
00068     typename scalar_type_<T>::eval trace(matrix<T> const& A) {
00069         typename scalar_type_<T>::eval tr = 0;
00070         assert(A.size(0) == A.size(1));
00071         for(int i = 0; i < A.size(0); ++i) tr += trace(A(i,i));
00072         return tr;
00073     }
00074 
00075 
00076     // --- mvector ---
00077 
00078     template<typename T>
00079     inline typename scalar_type_<T>::eval
00080     dot(mvector<T> const& A, mvector<T> const& B) {
00081         typedef typename mvector<T>::const_iterator iterator;
00082         if(A.size() != B.size())
00083             throw std::out_of_range(
00084                 "dot product applied to vectors with different sizes");
00085         typename scalar_type_<T>::eval sum = zero;
00086         for(iterator u = A.begin(), v = B.begin(), uend = A.end();
00087             u != uend; ++u, ++v) sum += dot(*u, *v);
00088         return sum;
00089     }
00090 //     template<typename T> struct norm_type_< mvector<T> >
00091 //      { typedef norm_type_<T>::eval eval; };
00092 
00093     template<typename T>
00094     inline typename norm_type_<T>::eval norm(const mvector<T>& A) {
00095         typedef typename mvector<T>::const_iterator iterator;
00096         typename norm_type_<T>::eval sum = zero;
00097         for(iterator u = A.begin(), uend = A.end();
00098             u != uend; ++u) sum += norm(*u);
00099         return sum;
00100     }
00101 
00102     template<typename T>
00103     inline typename norm_type_<T>::eval
00104     abs(const mvector<T>& A) { return more::sqrt(norm(A)); }
00105 
00106     template<typename T> inline T sum(const mvector<T>& x)
00107         { return std::accumulate(x.begin(), x.end(), (T)0); }
00108 
00109     template<typename T> inline
00110     void normalize(mvector<T>& f) {
00111         norm_type_<T>::eval c = 0;
00112         for(int s = 0; s < f.size(); ++s) c += norm(f(s));
00113         c = 1/sqrt(c);
00114         for(int s = 0; s < f.size(); ++s) f(s) *= c;
00115     }
00116 
00117     template<typename T> inline mvector<T>
00118     operator-(const mvector<T>& A, const mvector<T>& B) {
00119         assert(A.size(0) == B.size(0));
00120         mvector<T> C(A.size(0));
00121         std::transform(A.begin(), A.end(), B.begin(), C.begin(),
00122                        std::minus<T>());
00123         return C;
00124     }
00125 
00126     template<typename T> inline mvector<T>
00127     operator+(const mvector<T>& A, const mvector<T>& B) {
00128         assert(A.size(0) == B.size(0));
00129         mvector<T> C(A.size(0));
00130         std::transform(A.begin(), A.end(), B.begin(), C.begin(),
00131                        std::plus<T>());
00132         return C;
00133     }
00134 
00135     template<typename T, typename U>
00136     inline mvector< typename times_type_<T, U>::eval >
00137     operator*(const T& a, const mvector<U>& B) {
00138         typedef mvector< typename times_type_<T, U>::eval > vec_type;
00139         vec_type C(B.size(0));
00140         typename mvector<U>::const_iterator v = B.begin();
00141         for(typename vec_type::iterator u = C.begin(); u != C.end(); ++u) {
00142             *u = a * (*v);
00143             ++v;
00144         }
00145         return C;
00146     }
00147 }} // namespace more::num
00148 ///\endif
00149 
00150 #endif

Generated on Sat Sep 7 19:11:19 2002 for more with Doxygen 1.2.13.1. Doxygen 1.2.13.1 is written and copyright 1997-2002 by Dimitri van Heesch.