dune-localfunctions
2.2.0
|
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