COCONUT API: linalg.h Source File

00001 // Linear Algebra Wrapper Header -*- C++ -*-
00002 
00003 // $Id: linalg.h 875 2009-02-13 13:10:41Z dferi $
00004 // Copyright (C) 2001-2003 Hermann Schichl
00005 //
00006 // This file is part of the COCONUT API.  This library
00007 // is free software; you can redistribute it and/or modify it under the
00008 // terms of the Library GNU General Public License as published by the
00009 // Free Software Foundation; either version 2, or (at your option)
00010 // any later version.
00011 
00012 // This library is distributed in the hope that it will be useful,
00013 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00014 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015 // Library GNU General Public License for more details.
00016 
00017 // As a special exception, you may use this file as part of a free software
00018 // library without restriction.  Specifically, if other files instantiate
00019 // templates or use macros or inline functions from this file, or you compile
00020 // this file and link it with other files to produce an executable, this
00021 // file does not by itself cause the resulting executable to be covered by
00022 // the Library GNU General Public License.  This exception does not however
00023 // invalidate any other reasons why the executable file might be covered by
00024 // the Library GNU General Public License.
00025 
00028 #ifndef __LINALG_H_
00029 #define __LINALG_H_
00030 
00031 #if 0
00032 #include <Eigen2/Core>
00033 
00034 namespace linalg {
00035 
00037 
00044 template <class _Tp>
00045 class matrix : public Eigen::DynamicSparseMatrix<_Tp, Eigen::RowMajor>
00046 {
00047   typedef matrix<_Tp> _Self;
00048 public:
00050   typedef typename Eigen::DynamicSparseMatrix<_Tp, Eigen::RowMajor> _Base;
00052   typedef typename Eigen::SparseVector Row;
00053 
00054 #if defined(_MSC_VER)
00055 public:
00056 #else
00057 protected:
00058 #endif
00059 
00062   linalg_output_format output_format;
00063 
00064 public:
00066   matrix() : _Base(), output_format(linalg_output_sparse) {}
00068   matrix(const _Self& __m) : _Base((_Base)__m), output_format(__m.output_format) {}
00069 
00071   matrix(const size_t& n, const size_t& m)
00072                 : _Base(n,m), output_format(linalg_output_sparse) {}
00073 
00075   void setOutputFormat(const linalg_output_format& f) { output_format=f; }
00076 
00078   const linalg_output_format& getOutputFormat() const { return output_format; }
00079 
00080 public:
00082   _Self& operator= (const _Self& a)
00083   {
00084     (*(_Base*)this).operator=((_Base)a);
00085         this->output_format = a.output_format;
00086         return *this;
00087   }
00088 
00089   //arithmetic
00091   _Self& operator+= (const _Self& a)
00092   {
00093     *this += a;
00094     return *this;
00095   }
00096 
00098   _Self& operator*= (_Tp a)
00099   {
00100     *this *= a;
00101     return *this;
00102   }
00103 
00105   bool operator== (const _Self& a) const;
00106 
00108   bool operator!= (const _Self& a) const;
00109 
00110   //end of modifications
00111 
00112 #if !defined(_MSC_VER)
00113   template <typename _Tq>
00114   friend std::ostream& std::operator<<(std::ostream &s, const matrix<_Tq> &A);
00115 #endif
00116 };
00117 
00118 //arithemtic ops
00120 
00122 template <class _Tp>
00123 matrix<_Tp> operator+ (const matrix<_Tp>& A, const matrix<_Tp>& B)
00124 {
00125   unsigned int rmax=std::max(A.nrows(),B.nrows());
00126   unsigned int cmax=std::max(A.ncols(),B.ncols());
00127   matrix<_Tp> out(rmax,cmax);
00128   typename matrix<_Tp>::const_iterator i;
00129   typename matrix<_Tp>::OneD::const_iterator j;
00130   for (i = A.begin(); i < A.end(); ++i){
00131      for (j = (*i).begin(); j < (*i).end(); ++j)
00132        out(j.row(),j.column())=*j;
00133   }
00134   for (i = B.begin(); i < B.end(); ++i){
00135      for (j = (*i).begin(); j < (*i).end(); ++j)
00136        out(j.row(),j.column())+=*j;
00137   }
00138   return out;
00139 }
00140 
00142 
00145 template <class _Tp>
00146 matrix<_Tp> operator* (_Tp a, const matrix<_Tp>& B)
00147 {
00148   matrix<_Tp> out=B.new_instance();
00149   scale(out,a);
00150   return out;
00151 }
00152 
00153 template <class _Tp>
00154 bool matrix<_Tp>::operator==(const matrix<_Tp> &A) const
00155 {
00156   if(this == &A)
00157     return true;
00158   else
00159   {
00160     if(this->nrows() != A.nrows() || this->ncols() != A.ncols())
00161       return false;
00162     typename matrix<_Tp>::const_iterator       i, Ai;
00163     typename matrix<_Tp>::OneD::const_iterator j, Aj;
00164     typename matrix<_Tp>::const_iterator       i_max(this->end()), Ai_max(A.end());
00165     typename matrix<_Tp>::OneD::const_iterator j_max, Aj_max;
00166 
00167     if(this->nrows() > 0 && this->ncols() > 0)
00168     {
00169       for(i = this->begin(), Ai = A.begin(); i != i_max && Ai != Ai_max; ++i, ++Ai)
00170       {
00171         j_max = ((*i).end());
00172         Aj_max = ((*Ai).end());
00173         for(j = (*i).begin(), Aj = (*Ai).begin(); j != j_max && Aj != Aj_max;
00174             ++j, ++Aj)
00175         {
00176           if(j.row() != Aj.row() || j.column() != Aj.column() || *j != *Aj)
00177             return false;
00178         }
00179         if((j == j_max) != (Aj == Aj_max))
00180           return false;
00181       }
00182       if((i == i_max) != (Ai == Ai_max))
00183         return false;
00184     }
00185     return true;
00186   }
00187   return true; // to make the compiler happy
00188 }
00189 
00190 template <class _Tp>
00191 bool matrix<_Tp>::operator!=(const matrix<_Tp> &A) const
00192 {
00193   return !(this->operator==(A));
00194 }
00195 
00196 1} // namespace linalg
00197 
00198 
00199 #else
00200 #include <mtl/mtl.h>
00201 #include <mtl/matrix.h>
00202 #include <mtl/scaled1D.h>
00203 #include <mtl/utils.h>
00204 #include <mtl/blais.h>
00205 #include <typeinfo>
00206 
00208 
00211 namespace linalg {
00212 
00215 #define __LINALG_RANGE_CHECK    0
00216 
00218 
00223 typedef enum {
00227   linalg_output_sparse,
00231   linalg_output_full
00232 } linalg_output_format;
00233 
00234 /********/
00235 /*MATRIX*/
00236 /********/
00237 
00239 
00246 template <class _Tp>
00247 class matrix : public mtl::matrix<_Tp, mtl::rectangle<>, mtl::compressed<>, mtl::row_major>::type
00248 {
00249   typedef matrix<_Tp> _Self;
00250 public:
00252   typedef typename mtl::matrix<_Tp, mtl::rectangle<>, mtl::compressed<>,
00253                    mtl::row_major>::type _Base;
00255   typedef typename _Base::OneD Row;
00256 
00257 #if defined(_MSC_VER)
00258 public:
00259 #else
00260 protected:
00261 #endif
00262 
00265   linalg_output_format output_format;
00266 
00267 public:
00269   matrix() : _Base(), output_format(linalg_output_sparse) {}
00271   matrix(const _Self& __m) : output_format(__m.output_format) {
00272         (*(_Base*)this).operator=(__m.new_instance());
00273   }
00275   matrix(const size_t& n, const size_t& m)
00276                 : _Base(n,m), output_format(linalg_output_sparse) {}
00277 
00278   //**modified by Ferenc Domes
00279   //output
00281   void setOutputFormat(const linalg_output_format& f) { output_format=f; }
00282 
00284   const linalg_output_format& getOutputFormat() const { return output_format; }
00285 
00286 private:
00291   _Base new_instance() const
00292   {
00293     _Base out(this->nrows(),this->ncols());
00294     typename _Self::const_iterator i;
00295     typename _Self::OneD::const_iterator j;
00296     for (i = this->begin(); i < this->end(); ++i){
00297        for (j = (*i).begin(); j < (*i).end(); ++j)
00298          out(j.row(),j.column())=*j;
00299     }
00300     return out;
00301   }
00302 
00303 public:
00304   //OPERATORS
00305   //dereference
00307   _Self& operator= (const _Self& a)
00308   {
00309     (*(_Base*)this).operator=(a.new_instance());
00310         this->output_format = a.output_format;
00311         return *this;
00312   }
00313 
00314   //arithmetic
00316   _Self& operator+= (const _Self& a)
00317   {
00318     _Self t = *this + a;
00319     *this = t;
00320     return *this;
00321   }
00322 
00324   _Self& operator*= (_Tp a)
00325   {
00326     _Self t = a * *this;
00327     *this = t;
00328     return *this;
00329   }
00330 
00332   bool operator== (const _Self& a) const;
00333 
00335   bool operator!= (const _Self& a) const;
00336 
00337   //end of modifications
00338 
00339 #if !defined(_MSC_VER)
00340   template <typename _Tq>
00341   friend std::ostream& std::operator<<(std::ostream &s, const matrix<_Tq> &A);
00342 #endif
00343 };
00344 //**modified by Ferenc Domes!!!
00345 //MATRIX OPERATORS
00346 
00347 //arithemtic ops
00349 
00351 template <class _Tp>
00352 matrix<_Tp> operator+ (const matrix<_Tp>& A, const matrix<_Tp>& B)
00353 {
00354   unsigned int rmax=std::max(A.nrows(),B.nrows());
00355   unsigned int cmax=std::max(A.ncols(),B.ncols());
00356   matrix<_Tp> out(rmax,cmax);
00357   typename matrix<_Tp>::const_iterator i;
00358   typename matrix<_Tp>::OneD::const_iterator j;
00359   for (i = A.begin(); i < A.end(); ++i){
00360      for (j = (*i).begin(); j < (*i).end(); ++j)
00361        out(j.row(),j.column())=*j;
00362   }
00363   for (i = B.begin(); i < B.end(); ++i){
00364      for (j = (*i).begin(); j < (*i).end(); ++j)
00365        out(j.row(),j.column())+=*j;
00366   }
00367   return out;
00368 }
00369 
00371 
00374 template <class _Tp>
00375 matrix<_Tp> operator* (_Tp a, const matrix<_Tp>& B)
00376 {
00377   matrix<_Tp> out=B.new_instance();
00378   scale(out,a);
00379   return out;
00380 }
00381 
00382 template <class _Tp>
00383 bool matrix<_Tp>::operator==(const matrix<_Tp> &A) const
00384 {
00385   if(this == &A)
00386     return true;
00387   else
00388   {
00389     if(this->nrows() != A.nrows() || this->ncols() != A.ncols())
00390       return false;
00391     typename matrix<_Tp>::const_iterator       i, Ai;
00392     typename matrix<_Tp>::OneD::const_iterator j, Aj;
00393     typename matrix<_Tp>::const_iterator       i_max(this->end()), Ai_max(A.end());
00394     typename matrix<_Tp>::OneD::const_iterator j_max, Aj_max;
00395 
00396     if(this->nrows() > 0 && this->ncols() > 0)
00397     {
00398       for(i = this->begin(), Ai = A.begin(); i != i_max && Ai != Ai_max; ++i, ++Ai)
00399       {
00400         j_max = ((*i).end());
00401         Aj_max = ((*Ai).end());
00402         for(j = (*i).begin(), Aj = (*Ai).begin(); j != j_max && Aj != Aj_max;
00403             ++j, ++Aj)
00404         {
00405           if(j.row() != Aj.row() || j.column() != Aj.column() || *j != *Aj)
00406             return false;
00407         }
00408         if((j == j_max) != (Aj == Aj_max))
00409           return false;
00410       }
00411       if((i == i_max) != (Ai == Ai_max))
00412         return false;
00413     }
00414     return true;
00415   }
00416   return true; // to make the compiler happy
00417 }
00418 
00419 template <class _Tp>
00420 bool matrix<_Tp>::operator!=(const matrix<_Tp> &A) const
00421 {
00422   return !(this->operator==(A));
00423 }
00424 
00425 //**end modifications
00426 
00427 /**********/
00428 /*C_MATRIX*/
00429 /**********/
00431 
00438 template <class _Tp>
00439 class c_matrix : public mtl::matrix<_Tp, mtl::rectangle<>, mtl::compressed<>, mtl::column_major>::type
00440 {
00441   typedef c_matrix<_Tp> _Self;
00442 public:
00444   typedef typename mtl::matrix<_Tp, mtl::rectangle<>, mtl::compressed<>, mtl::column_major>::type _Base;
00446   typedef typename _Base::OneD Col;
00447 
00448 protected:
00452   linalg_output_format output_format;
00453 
00454 public:
00456   c_matrix() : _Base(), output_format(linalg_output_sparse) {}
00457   c_matrix(const _Self& __m) : output_format(__m.output_format) {
00458     *this = __m.new_instance();
00459   }
00460   c_matrix(const size_t& n, const size_t& m)
00461         : _Base(n,m), output_format(linalg_output_sparse)  {}
00462 
00464   void setOutputFormat(const linalg_output_format& f) { output_format=f; }
00465 
00467   const linalg_output_format& getOutputFormat() const { return output_format; }
00468 
00469 private:
00474   _Self new_instance() const
00475   {
00476     _Self out(this->nrows(),this->ncols());
00477     typename _Self::const_iterator i;
00478     typename _Self::OneD::const_iterator j;
00479     for (i = this->begin(); i < this->end(); ++i){
00480        for (j = (*i).begin(); j < (*i).end(); ++j)
00481          out(j.row(),j.column())=*j;
00482     }
00483     return out;
00484   }
00485 
00486 public:
00487   //OPERATORS
00488   //dereference
00489 
00490   //arithmetic
00492   _Self& operator+= (const _Self& a)
00493   {
00494     _Self t = *this + a;
00495     *this = t;
00496     return *this;
00497   }
00498 
00500   _Self& operator*= (_Tp a)
00501   {
00502     _Self t = a * *this;
00503     *this = t;
00504     return *this;
00505   }
00506   
00508   bool operator== (const _Self& a) const;
00509 
00511   bool operator!= (const _Self& a) const;
00512 
00513   //end of modifications
00514 
00515 #if !defined(_MSC_VER)
00516   template <typename _Tq>
00517   friend std::ostream& std::operator<<(std::ostream &s, const c_matrix<_Tq> &A);
00518 #endif
00519 };
00520 
00521 //arithemtic ops
00523 
00525 template <class _Tp>
00526 c_matrix<_Tp> operator+ (const c_matrix<_Tp>& A, const c_matrix<_Tp>& B)
00527 {
00528   unsigned int rmax=std::max(A.nrows(),B.nrows());
00529   unsigned int cmax=std::max(A.ncols(),B.ncols());
00530   c_matrix<_Tp> out(rmax,cmax);
00531   typename c_matrix<_Tp>::const_iterator i;
00532   typename c_matrix<_Tp>::OneD::const_iterator j;
00533   for (i = A.begin(); i < A.end(); ++i){
00534      for (j = (*i).begin(); j < (*i).end(); ++j)
00535        out(j.row(),j.column())=*j;
00536   }
00537   for (i = B.begin(); i < B.end(); ++i){
00538      for (j = (*i).begin(); j < (*i).end(); ++j)
00539        out(j.row(),j.column())+=*j;
00540   }
00541   return out;
00542 }
00543 
00545 
00548 template <class _Tp>
00549 c_matrix<_Tp> operator* (_Tp a, const c_matrix<_Tp>& B)
00550 {
00551   c_matrix<_Tp> out=B.new_instance();
00552   scale(out,a);
00553   return out;
00554 }
00555 
00556 template <class _Tp>
00557 bool c_matrix<_Tp>::operator==(const c_matrix<_Tp> &A) const
00558 {
00559   if(this == &A)
00560     return true;
00561   else
00562   {
00563     if(this->nrows() != A.nrows() || this->ncols() != A.ncols())
00564       return false;
00565     typename c_matrix<_Tp>::const_iterator       i, Ai;
00566     typename c_matrix<_Tp>::OneD::const_iterator j, Aj;
00567     typename c_matrix<_Tp>::const_iterator     i_max(this->end()), Ai_max(A.end());
00568     typename c_matrix<_Tp>::OneD::const_iterator j_max, Aj_max;
00569 
00570     if(this->nrows() > 0 && this->ncols() > 0)
00571     {
00572       for(i = this->begin(), Ai = A.begin(); i != i_max && Ai != Ai_max; ++i, ++Ai)
00573       {
00574         j_max = ((*i).end());
00575         Aj_max = ((*Ai).end());
00576         for(j = (*i).begin(), Aj = (*Ai).begin(); j != j_max && Aj != Aj_max;
00577             ++j, ++Aj)
00578         {
00579           if(j.row() != Aj.row() || j.column() != Aj.column() || *j != *Aj)
00580             return false;
00581         }
00582         if((j == j_max) != (Aj == Aj_max))
00583           return false;
00584       }
00585       if((i == i_max) != (Ai == Ai_max))
00586         return false;
00587     }
00588     return true;
00589   }
00590   return true; // to make the compiler happy
00591 }
00592 
00593 template <class _Tp>
00594 bool c_matrix<_Tp>::operator!=(const c_matrix<_Tp> &A) const
00595 {
00596   return !(this->operator==(A));
00597 }
00598 
00599 /***************/
00600 /*SPARSE_VECTOR*/
00601 /***************/
00603 
00610 template <class _Tp>
00611 class sparse_vector : public mtl::compressed1D<_Tp>
00612 {
00613   typedef sparse_vector<_Tp> _Self;
00615   typedef mtl::compressed1D<_Tp> _Base;
00616 
00617   typedef typename _Base::const_reference const_reference;
00618   typedef typename _Base::reference reference;
00619   typedef typename _Base::size_type size_type;
00620   typedef typename std::vector<_Tp> values_vec;
00621   typedef typename std::vector<size_type> indices_vec;
00622 
00623 //**modified by Ferenc Domes
00624 #if defined(_MSC_VER)
00625 public:
00626 #else
00627 protected:
00628 #endif
00629 
00631   linalg_output_format output_format;
00632 
00633 public:
00635   sparse_vector() : _Base(), output_format(linalg_output_sparse) {}
00637   sparse_vector(const _Self& __m): output_format(__m.output_format)
00638   {
00639     this->size_ = __m.size_;
00640     this->values = new values_vec(*__m.values);
00641     this->indices = new indices_vec(*__m.indices);
00642   }
00644   sparse_vector(const size_t& n) :_Base(n), output_format(linalg_output_sparse) {}
00645 
00647   void setOutputFormat(const linalg_output_format& f) { output_format=f; }
00648 
00650   const linalg_output_format& getOutputFormat() const { return output_format; }
00651   
00652   //OPERATORS
00655   _Self& operator=(const _Self& x)
00656   {
00657     this->size_ = x.size_;
00658     this->values = new values_vec(*x.values);
00659     this->indices = new indices_vec(*x.indices);
00660     this->output_format = x.output_format;
00661     return *this;
00662   }
00663 
00664   //dereference
00666   inline const_reference operator[](size_type i) const MTL_THROW_ASSERTION
00667   {
00668     MTL_ASSERT(i < this->size(), "compressed1D::operator[]");
00669     return const_reference(*this, i);
00670   }
00672   inline reference operator[](size_type i)
00673   {
00674     if (i>=this->size())
00675       this->resize(i+1);
00676     return reference(*this,i);
00677   }
00678 
00679   //arithmetic
00681   _Self &operator+= (const _Self& a)
00682   {
00683     _Self t = *this+a;
00684     *this = t;
00685     return *this;
00686   }  
00688   _Self &operator*= (const _Tp& a)
00689   {
00690     _Self t = a * *this;
00691     *this = t;
00692     return *this;
00693   }
00694   
00696   _Tp &operator*= (const _Self& a)
00697   {
00698     _Self t = a * *this;
00699     *this = t;
00700     return *this;
00701   }
00702 //end of modifications
00703 
00704 #if ! defined(_MSC_VER)
00705   template <typename _Tq>
00706   friend std::ostream& std::operator<<(std::ostream &s,
00707                                        const sparse_vector<_Tq> &A);
00708 #endif
00709 };
00710 
00711 //VECTOR OPERATORS
00712 //comparsion
00714 
00716 template <class _Tp>
00717 bool operator==(const sparse_vector<_Tp>& a, const sparse_vector<_Tp>& b)
00718 {
00719   typename sparse_vector<_Tp>::const_iterator _b1, _b2, _e1, _e2;
00720 
00721   _b1 = a.begin();
00722   _b2 = b.begin();
00723   _e1 = a.end();
00724   _e2 = b.end();
00725   while(_b1 != _e1 || _b2 != _e2)
00726   {
00727     if(_b1 == _e1)
00728     {
00729       if(*_b2 == _Tp())
00730       {
00731         ++_b2;
00732         continue;
00733       }
00734       else
00735         return false;
00736     }
00737     else if(_b2 == _e2)
00738     {
00739       if(*_b1 == _Tp())
00740       {
00741         ++_b1;
00742         continue;
00743       }
00744       else
00745         return false;
00746     }
00747     else if(_b1.index() < _b2.index())
00748     {
00749       if(*_b1 == _Tp())
00750       {
00751         ++_b1;
00752         continue;
00753       }
00754       else
00755         return false;
00756     }
00757     else if(_b2.index() < _b1.index())
00758     {
00759       if(*_b2 == _Tp())
00760       {
00761         ++_b2;
00762         continue;
00763       }
00764       else
00765         return false;
00766     }
00767     else if(*_b1 != *_b2)
00768       return false;
00769     ++_b1;
00770     ++_b2;
00771   }
00772   return true;
00773 }
00774 
00775 //**modified by Domes Ferenc NOT FINISHED!!!
00776 //arithmetic
00778 
00780 template <class _Tp> 
00781 sparse_vector<_Tp> operator+ (const sparse_vector<_Tp>& a,
00782                               const sparse_vector<_Tp>& b)
00783 {
00784   unsigned int a_size=a.size(),b_size=b.size(),i=0;
00785   unsigned int maxs=std::max(a_size,b_size), mins=std::min(a_size,b_size);
00786   sparse_vector<_Tp> out(maxs);
00787   for (;i<mins;i++){
00788     if((a[i] != _Tp()) || (b[i]!=_Tp()))
00789       out[i]=a[i]+b[i];
00790   }
00791   if (a_size>b_size){
00792     for(;i<maxs;i++)
00793       if(a[i] != _Tp())
00794         out[i]=a[i];
00795   }
00796   else if (a_size<b_size){
00797     for(;i<maxs;i++)
00798       if(b[i] != _Tp())
00799         out[i]=b[i];
00800   }
00801   return out;
00802 }
00803 
00805 
00807 template <class _Tp> 
00808 _Tp operator* (const sparse_vector<_Tp>& a, const sparse_vector<_Tp>& b)
00809 {
00810   //DOT-PRODUCT
00811   typename sparse_vector<_Tp>::const_iterator _i,_j;
00812   _Tp out;
00813   for (_i=a.begin(), _j=b.begin(); _i != a.end() && _j != b.end();){
00814     unsigned int ii=_i.index(),ji=_j.index();
00815     if (ii==ji){
00816       out+=(*_i)*(*_j);
00817       ++_i;++_j;
00818     }
00819     else if (ii<ji)
00820       ++_i;
00821     else
00822       ++_j;
00823   }
00824   return out;
00825 }
00826 
00828 
00830 template <class _Tp>
00831 sparse_vector<_Tp> operator* (_Tp a, const sparse_vector<_Tp>& b)
00832 {
00833   typename sparse_vector<_Tp>::const_iterator _i;
00834   sparse_vector<_Tp> out(b.size());
00835   for (_i=b.begin();_i< b.end();_i++)
00836       out[_i.index()]=a*(*_i);
00837   return out;
00838 }
00839 //**end modifications
00840 
00842 
00846 template <class _TS>
00847 class __linalg_cvec : public std::vector<_TS>
00848 {
00849 private:
00850   typedef std::vector<_TS> _Base;
00851   typedef __linalg_cvec<_TS> _Self;
00852   const std::vector<_TS>& vec;
00853 
00854 public:
00856 
00857   enum { N = 0 };
00858 
00859   typedef typename std::vector<_TS>::value_type value_type;
00860   typedef typename std::vector<_TS>::reference reference;
00861   typedef typename std::vector<_TS>::pointer pointer;
00862   typedef typename std::vector<_TS>::const_reference const_reference;
00863   typedef typename std::vector<_TS>::const_pointer const_pointer;
00864   typedef typename std::vector<_TS>::size_type size_type;
00865   typedef typename std::vector<_TS>::difference_type difference_type;
00866   typedef typename std::vector<_TS>::allocator_type allocator_type;
00867   typedef typename std::vector<_TS>::iterator iterator;
00868   typedef typename std::vector<_TS>::const_iterator const_iterator;
00869   typedef typename std::vector<_TS>::reverse_iterator reverse_iterator;
00870   typedef typename std::vector<_TS>::const_reverse_iterator const_reverse_iterator;
00871 
00872   typedef mtl::dense_tag sparsity;
00873   typedef mtl::oned_tag dimension;
00874 
00875   typedef mtl::scaled1D<_Self> scaled_type;
00876 
00877   typedef mtl::light1D<_TS> subrange_type;
00878   typedef _Self IndexArray;
00879   typedef _Self IndexArrayRef;
00881 
00883   __linalg_cvec(const std::vector<_TS>& s) : vec(s) {}
00885   ~__linalg_cvec() {}
00886 
00888 
00890   iterator begin() { return vec.begin(); }
00891   iterator end() { return vec.end(); }
00892 
00893   const_iterator begin() const { return vec.begin(); }
00894   const_iterator end() const { return vec.end(); }
00896 
00899   unsigned int size() const { return vec.size(); }
00900 };
00901 
00903 
00907 template <class _TS>
00908 class __linalg_vec : public std::vector<_TS>
00909 {
00910 private:
00911   typedef std::vector<_TS> _Base;
00912   typedef __linalg_vec<_TS> _Self;
00913   std::vector<_TS>& vec;
00914 
00915 public:
00917 
00918   enum { N = 0 };
00919 
00920   typedef typename std::vector<_TS>::value_type value_type;
00921   typedef typename std::vector<_TS>::reference reference;
00922   typedef typename std::vector<_TS>::pointer pointer;
00923   typedef typename std::vector<_TS>::const_reference const_reference;
00924   typedef typename std::vector<_TS>::const_pointer const_pointer;
00925   typedef typename std::vector<_TS>::size_type size_type;
00926   typedef typename std::vector<_TS>::difference_type difference_type;
00927   typedef typename std::vector<_TS>::allocator_type allocator_type;
00928   typedef typename std::vector<_TS>::iterator iterator;
00929   typedef typename std::vector<_TS>::const_iterator const_iterator;
00930   typedef typename std::vector<_TS>::reverse_iterator reverse_iterator;
00931   typedef typename std::vector<_TS>::const_reverse_iterator const_reverse_iterator;
00932 
00933   typedef mtl::dense_tag sparsity;
00934   typedef mtl::oned_tag dimension;
00935   
00936   typedef mtl::scaled1D<_Self> scaled_type;
00937 
00938   typedef mtl::light1D<_TS> subrange_type;
00939   typedef _Self IndexArray;
00940   typedef _Self IndexArrayRef;
00942 
00944   __linalg_vec(std::vector<_TS>& s) : vec(s) {}
00946   ~__linalg_vec() {}
00947 
00949 
00951   iterator begin() { return vec.begin(); }
00952   iterator end() { return vec.end(); }
00953 
00954   const_iterator begin() const { return vec.begin(); }
00955   const_iterator end() const { return vec.end(); }
00957 
00960   unsigned int size() const { return vec.size(); }
00961 };
00962 
00964 
00966 template <class _TS>
00967 inline _TS linalg_dot(const std::vector<_TS>& a, const std::vector<_TS>& b, _TS s)
00968 {
00969   return mtl::dot(__linalg_cvec<_TS>(a), __linalg_cvec<_TS>(b), s);
00970 }
00971 
00972 template <class _Va, class _TS>
00973 inline _TS linalg_dot(const _Va& a, const std::vector<_TS>& b, _TS s)
00974 {
00975   __linalg_cvec<_TS> _b(b);
00976   return mtl::dot(a, _b, s);
00977 }
00978 
00979 template <class _Vb, class _TS>
00980 inline _TS linalg_dot(const std::vector<_TS>& a, const _Vb& b, _TS s)
00981 {
00982   return mtl::dot(__linalg_cvec<_TS>(a), b, s);
00983 }
00984 
00985 template <class _Va, class _Vb, class _TS>
00986 inline _TS linalg_dot(const _Va& a, const _Vb& b, _TS s)
00987 {
00988   return mtl::dot(a, b, s);
00989 }
00991 
00993 
00995 template <class _TS>
00996 inline void linalg_add(const std::vector<_TS>& a, const std::vector<_TS>& b, std::vector<_TS>& c)
00997 {
00998   mtl::add(__linalg_cvec<_TS>(a), __linalg_cvec<_TS>(b), __linalg_vec<_TS>(c));
00999 }
01000 
01001 template <class _Va, class _TS>
01002 inline void linalg_add(const _Va& a, const std::vector<_TS>& b, std::vector<_TS>& c)
01003 {
01004   mtl::add((a), __linalg_cvec<_TS>(b), __linalg_vec<_TS>(c));
01005 }
01006 
01007 template <class _Vb, class _TS>
01008 inline void linalg_add(const std::vector<_TS>& a, const _Vb& b, std::vector<_TS>& c)
01009 {
01010   mtl::add(__linalg_cvec<_TS>(a), (b), __linalg_vec<_TS>(c));
01011 }
01012 
01013 template <class _TS, class _Vc>
01014 inline void linalg_add(const std::vector<_TS>& a, const std::vector<_TS>& b, _Vc& c)
01015 {
01016   mtl::add(__linalg_cvec<_TS>(a), __linalg_cvec<_TS>(b), (c));
01017 }
01018 
01019 template <class _Va, class _Vb, class _TS>
01020 inline void linalg_add(const _Va& a, const _Vb& b, std::vector<_TS>& c)
01021 {
01022   mtl::add((a), (b), __linalg_vec<_TS>(c));
01023 }
01024 
01025 template <class _Va, class _TS, class _Vc>
01026 inline void linalg_add(const _Va& a, const std::vector<_TS>& b, _Vc& c)
01027 {
01028   mtl::add((a), __linalg_cvec<_TS>(b), (c));
01029 }
01030 
01031 template <class _TS, class _Vb, class _Vc>
01032 inline void linalg_add(const std::vector<_TS>& a, const _Vb& b, _Vc& c)
01033 {
01034   mtl::add(__linalg_cvec<_TS>(a), (b), (c));
01035 }
01036 
01037 template <class _Va, class _Vb, class _Vc>
01038 inline void linalg_add(const _Va& a, const _Vb& b, _Vc& c)
01039 {
01040   mtl::add((a), (b), (c));
01041 }
01043 
01044 // matrix vector multiplication
01046 
01048 template <class _Matrix, class _TS>
01049 inline void linalg_matvec(const _Matrix& A, const std::vector<_TS>& a,
01050                           std::vector<_TS>& c)
01051 {
01052   mtl::mult((A), __linalg_cvec<_TS>(a), __linalg_vec<_TS>(c));
01053 }
01054 
01055 template <class _Matrix, class _TS, class _Va>
01056 inline void linalg_matvec(const _Matrix& A, const _Va& a,
01057                           std::vector<_TS>& c)
01058 {
01059   mtl::mult((A), (a), __linalg_vec<_TS>(c));
01060 }
01061 
01062 template <class _Matrix, class _TS, class _Vc>
01063 inline void linalg_matvec(const _Matrix& A, const std::vector<_TS>& a,
01064                           _Vc& c)
01065 {
01066   mtl::mult((A), __linalg_cvec<_TS>(a), (c));
01067 }
01068 
01069 template <class _Matrix, class _Va, class _Vc>
01070 inline void linalg_matvec(const _Matrix& A, const _Va& a, _Vc& c)
01071 {
01072   mtl::mult((A), (a), (c));
01073 }
01075 
01077 
01079 template <class _Matrix, class _TS>
01080 inline void linalg_matvecv(const _Matrix& A, const std::vector<_TS>& a,
01081                           const std::vector<_TS>& b, std::vector<_TS>& c)
01082 {
01083   mtl::mult((A), __linalg_cvec<_TS>(a), __linalg_cvec<_TS>(b),
01084             __linalg_vec<_TS>(c));
01085 }
01086 
01087 template <class _Matrix, class _TS, class _Va>
01088 inline void linalg_matvecv(const _Matrix& A, const _Va& a,
01089                           const std::vector<_TS>& b, std::vector<_TS>& c)
01090 {
01091   mtl::mult((A), (a), __linalg_cvec<_TS>(b), __linalg_vec<_TS>(c));
01092 }
01093 
01094 template <class _Matrix, class _TS, class _Vb>
01095 inline void linalg_matvecv(const _Matrix& A, const std::vector<_TS>& a,
01096                           const _Vb& b, std::vector<_TS>& c)
01097 {
01098   mtl::mult((A), __linalg_cvec<_TS>(a), (b), __linalg_vec<_TS>(c));
01099 }
01100 
01101 template <class _Matrix, class _TS, class _Vc>
01102 inline void linalg_matvecv(const _Matrix& A, const std::vector<_TS>& a,
01103                           const std::vector<_TS>& b, _Vc& c)
01104 {
01105   mtl::mult((A), __linalg_cvec<_TS>(a), __linalg_cvec<_TS>(b), (c));
01106 }
01107 
01108 template <class _Matrix, class _Va, class _Vb, class _TS>
01109 inline void linalg_matvecv(const _Matrix& A, _Va& a,
01110                           const _Vb& b, const std::vector<_TS>& c)
01111 {
01112   mtl::mult((A), (a), (b), __linalg_vec<_TS>(c));
01113 }
01114 
01115 template <class _Matrix, class _Va, class _TS, class _Vc>
01116 inline void linalg_matvecv(const _Matrix& A, const _Va& a,
01117                           const std::vector<_TS>& b, _Vc& c)
01118 {
01119   mtl::mult((A), (a), __linalg_cvec<_TS>(b), (c));
01120 }
01121 
01122 template <class _Matrix, class _TS, class _Vb, class _Vc>
01123 inline void linalg_matvecv(const _Matrix& A, const std::vector<_TS>& a,
01124                           const _Vb& b, _Vc& c)
01125 {
01126   mtl::mult((A), __linalg_cvec<_TS>(a), (b), (c));
01127 }
01128 
01129 template <class _Matrix, class _Va, class _Vb, class _Vc>
01130 inline void linalg_matvecv(const _Matrix& A, const _Va& a, const _Vb& b,
01131                           _Vc& c)
01132 {
01133   mtl::mult((A), (a), (b), (c));
01134 }
01136 
01139 #define linalg_trans(A) mtl::trans((A))
01140 
01142 #define linalg_scale(X, Y) mtl::scaled(X, Y)
01143 
01144 // other functions
01147 template <class _Ve>
01148 inline std::vector<_Ve>& linalg_ssum(std::vector<_Ve>& a, _Ve b, const std::vector<_Ve>& c)
01149 {
01150 #if __LINALG_RANGE_CHECK
01151   if(a.size() != b.size())
01152     throw "linalg: a and c have different size in linalg_ssum";
01153 #endif
01154   typename std::vector<_Ve>::iterator _a(a.begin()), e_a(a.end());
01155   typename std::vector<_Ve>::const_iterator _c(c.begin());
01156 
01157   while(_a != e_a)
01158     *_a++ += b * *_c++;
01159   return a;
01160 }
01161 
01164 template <class _Ve>
01165 inline std::vector<_Ve>& linalg_smult(std::vector<_Ve>& a, _Ve b)
01166 {
01167   typename std::vector<_Ve>::iterator _a(a.begin()), e_a(a.end());
01168 
01169   while(_a != e_a)
01170     *_a++ *= b;
01171   return a;
01172 }
01173 
01174 } // namespace linalg
01175 
01176 #endif
01177 
01178 
01179 #endif /* __LINALG_H_ */