dune-localfunctions  2.2.0
lfematrix.hh
Go to the documentation of this file.
00001 #ifndef DUNE_ALGLIB_MATRIX_HH
00002 #define DUNE_ALGLIB_MATRIX_HH
00003 
00004 #include <cassert>
00005 #include <vector>
00006 
00007 #include <dune/localfunctions/utility/field.hh>
00008 // #include <dune/localfunctions/utility/vector.hh>
00009 
00010 #if HAVE_ALGLIB
00011 #include <alglib/amp.h>
00012 #include <alglib/matinv.h>
00013 #warning ALGLIB support is deprecated and will be dropped after DUNE 2.2 (cf. FS#931)
00014 #endif
00015 
00016 namespace Dune
00017 {
00018 
00019   template< class F, bool aligned = false >
00020   class LFEMatrix;
00021 
00022 
00023   template< class F, bool aligned >
00024   class LFEMatrix
00025   {
00026     typedef LFEMatrix< F, aligned > This;
00027     // typedef LFEVector< F > Row;
00028     typedef std::vector< F > Row;
00029     typedef std::vector<Row> RealMatrix;
00030 
00031   public:
00032     typedef F Field;
00033 
00034     operator const RealMatrix & () const
00035     {
00036       return matrix_;
00037     }
00038 
00039     operator RealMatrix & ()
00040     {
00041       return matrix_;
00042     }
00043 
00044     template <class Vector>
00045     void row( const unsigned int row, Vector &vec ) const
00046     {
00047       assert(row<rows());
00048       for (int i=0;i<cols();++i)
00049         field_cast(matrix_[row][i], vec[i]);
00050     }
00051 
00052     const Field &operator() ( const unsigned int row, const unsigned int col ) const
00053     {
00054       assert(row<rows());
00055       assert(col<cols());
00056       return matrix_[ row ][ col ];
00057     }
00058 
00059     Field &operator() ( const unsigned int row, const unsigned int col )
00060     {
00061       assert(row<rows());
00062       assert(col<cols());
00063       return matrix_[ row ][ col ];
00064     }
00065 
00066     unsigned int rows () const
00067     {
00068       return rows_;
00069     }
00070 
00071     unsigned int cols () const
00072     {
00073       return cols_;
00074     }
00075 
00076     const Field *rowPtr ( const unsigned int row ) const
00077     {
00078       assert(row<rows());
00079       return &(matrix_[row][0]);
00080     }
00081 
00082     Field *rowPtr ( const unsigned int row )
00083     {
00084       assert(row<rows());
00085       return &(matrix_[row][0]);
00086     }
00087 
00088     void resize ( const unsigned int rows, const unsigned int cols )
00089     {
00090       matrix_.resize(rows);
00091       for (unsigned int i=0;i<rows;++i)
00092         matrix_[i].resize(cols);
00093       rows_ = rows;
00094       cols_ = cols;
00095     }
00096 
00097     bool invert ()
00098     {
00099       assert( rows() == cols() );
00100       std::vector<unsigned int> p(rows());
00101       for (unsigned int j=0;j<rows();++j)
00102         p[j] = j;
00103       for (unsigned int j=0;j<rows();++j)
00104       {
00105         // pivot search
00106         unsigned int r = j;
00107         Field max = std::abs( (*this)(j,j) ); 
00108         for (unsigned int i=j+1;i<rows();++i)
00109         {
00110           if ( std::abs( (*this)(i,j) ) > max )
00111           {
00112             max = std::abs( (*this)(i,j) ); 
00113             r = i;
00114           }
00115         }
00116         if (max == Zero<Field>())
00117           return false;
00118         // row swap
00119         if (r > j) 
00120         {
00121           for (unsigned int k=0;k<cols();++k)
00122             std::swap( (*this)(j,k), (*this)(r,k) );
00123           std::swap( p[j], p[r] );
00124         }
00125         // transformation
00126         Field hr = Unity<Field>()/(*this)(j,j);
00127         for (unsigned int i=0;i<rows();++i) 
00128           (*this)(i,j) *= hr;
00129         (*this)(j,j) = hr;
00130         for (unsigned int k=0;k<cols();++k) 
00131         {
00132           if (k==j) continue;
00133           for (unsigned int i=0;i<rows();++i)
00134           {
00135             if (i==j) continue;
00136             (*this)(i,k) -= (*this)(i,j)*(*this)(j,k);
00137           }
00138           (*this)(j,k) *= -hr;
00139         }
00140       }
00141       // column exchange
00142       Row hv(rows());
00143       for (unsigned int i=0;i<rows();++i)
00144       {
00145         for (unsigned int k=0;k<rows();++k)
00146           hv[ p[k] ] = (*this)(i,k);
00147         for (unsigned int k=0;k<rows();++k)
00148           (*this)(i,k) = hv[k];
00149       }
00150       return true;
00151     }
00152 
00153   private:
00154     RealMatrix matrix_;
00155     unsigned int cols_,rows_;
00156   };
00157 
00158 #if HAVE_ALGLIB
00159   template< unsigned int precision, bool aligned >
00160   class LFEMatrix< amp::ampf< precision >, aligned >
00161   {
00162   public:
00163     typedef amp::ampf< precision > Field;
00164   private:
00165     typedef LFEMatrix< amp::ampf< precision >, aligned > This;
00166     typedef ap::template_2d_array< Field, aligned > RealMatrix;
00167 
00168   public:
00169     operator const RealMatrix & () const
00170     {
00171       return matrix_;
00172     }
00173 
00174     operator RealMatrix & ()
00175     {
00176       return matrix_;
00177     }
00178 
00179     template <class Vector>
00180     void row( const unsigned int row, Vector &vec ) const
00181     {
00182       assert(row<rows());
00183       for (unsigned int i=0;i<cols();++i)
00184         field_cast(matrix_(row,i), vec[i]);
00185     }
00186 
00187     const Field &operator() ( const unsigned int row, const unsigned int col ) const
00188     {
00189       assert(row<rows());
00190       assert(col<cols());
00191       return matrix_( row, col );
00192     }
00193 
00194     Field &operator() ( const unsigned int row, const unsigned int col )
00195     {
00196       assert(row<rows());
00197       assert(col<cols());
00198       return matrix_( row, col );
00199     }
00200 
00201     unsigned int rows () const
00202     {
00203       return matrix_.gethighbound( 1 )+1;
00204     }
00205 
00206     unsigned int cols () const
00207     {
00208       return matrix_.gethighbound( 2 )+1;
00209     }
00210 
00211     const Field *rowPtr ( const unsigned int row ) const
00212     {
00213       assert(row<rows());
00214       const int lastCol = matrix_.gethighbound( 2 );
00215       ap::const_raw_vector< Field > rowVector = matrix_.getrow( row, 0, lastCol );
00216       assert( (rowVector.GetStep() == 1) && (rowVector.GetLength() == lastCol+1) );
00217       return rowVector.GetData();
00218     }
00219 
00220     Field *rowPtr ( const unsigned int row )
00221     {
00222       assert(row<rows());
00223       const int lastCol = matrix_.gethighbound( 2 );
00224       ap::raw_vector< Field > rowVector = matrix_.getrow( row, 0, lastCol );
00225       assert( (rowVector.GetStep() == 1) && (rowVector.GetLength() == lastCol+1) );
00226       return rowVector.GetData();
00227     }
00228 
00229     void resize ( const unsigned int rows, const unsigned int cols )
00230     {
00231       matrix_.setbounds( 0, rows-1, 0, cols-1 );
00232     }
00233 
00234     bool invert ()
00235     {
00236       assert( rows() == cols() );
00237       int info;
00238       matinv::matinvreport< precision > report;
00239       matinv::rmatrixinverse< precision >( matrix_, rows(), info, report );
00240       return (info >= 0);
00241     }
00242 
00243   private:
00244     RealMatrix matrix_;
00245   };
00246 #endif
00247 
00248   template< class Field, bool aligned >
00249   inline std::ostream &operator<<(std::ostream &out, const LFEMatrix<Field,aligned> &mat) 
00250   {
00251     for (unsigned int r=0;r<mat.rows();++r)
00252     {
00253       out << field_cast<double>(mat(r,0));
00254       for (unsigned int c=1;c<mat.cols();++c) 
00255       {
00256         out << " , " << field_cast<double>(mat(r,c));
00257       }
00258       out << std::endl;
00259     }
00260     return out;
00261   }
00262 }
00263 
00264 #endif // #ifndef DUNE_ALGLIB_MATRIX_HH