1 #ifndef DUNE_ALGLIB_MATRIX_HH
2 #define DUNE_ALGLIB_MATRIX_HH
11 #include <alglib/amp.h>
12 #include <alglib/matinv.h>
13 #warning ALGLIB support is deprecated and will be dropped after DUNE 2.2 (cf. FS#931)
19 template<
class F,
bool aligned = false >
23 template<
class F,
bool aligned >
28 typedef std::vector< F > Row;
29 typedef std::vector<Row> RealMatrix;
34 operator const RealMatrix & ()
const
39 operator RealMatrix & ()
44 template <
class Vector>
45 void row(
const unsigned int row, Vector &vec )
const
48 for (
int i=0;i<
cols();++i)
56 return matrix_[
row ][ col ];
63 return matrix_[
row ][ col ];
79 return &(matrix_[
row][0]);
85 return &(matrix_[
row][0]);
91 for (
unsigned int i=0;i<
rows;++i)
100 std::vector<unsigned int> p(
rows());
101 for (
unsigned int j=0;j<
rows();++j)
103 for (
unsigned int j=0;j<
rows();++j)
107 Field max = std::abs( (*
this)(j,j) );
108 for (
unsigned int i=j+1;i<
rows();++i)
110 if ( std::abs( (*
this)(i,j) ) > max )
112 max = std::abs( (*
this)(i,j) );
121 for (
unsigned int k=0;k<
cols();++k)
122 std::swap( (*
this)(j,k), (*
this)(r,k) );
123 std::swap( p[j], p[r] );
127 for (
unsigned int i=0;i<
rows();++i)
130 for (
unsigned int k=0;k<
cols();++k)
133 for (
unsigned int i=0;i<
rows();++i)
136 (*this)(i,k) -= (*
this)(i,j)*(*
this)(j,k);
143 for (
unsigned int i=0;i<
rows();++i)
145 for (
unsigned int k=0;k<
rows();++k)
146 hv[ p[k] ] = (*
this)(i,k);
147 for (
unsigned int k=0;k<
rows();++k)
148 (*
this)(i,k) = hv[k];
155 unsigned int cols_,rows_;
159 template<
unsigned int precision,
bool aligned >
160 class LFEMatrix< amp::ampf< precision >, aligned >
163 typedef amp::ampf< precision >
Field;
165 typedef LFEMatrix< amp::ampf< precision >, aligned > This;
166 typedef ap::template_2d_array< Field, aligned > RealMatrix;
169 operator const RealMatrix & ()
const
174 operator RealMatrix & ()
179 template <
class Vector>
180 void row(
const unsigned int row, Vector &vec )
const
183 for (
unsigned int i=0;i<
cols();++i)
187 const Field &
operator() (
const unsigned int row,
const unsigned int col )
const
191 return matrix_( row, col );
198 return matrix_( row, col );
201 unsigned int rows ()
const
203 return matrix_.gethighbound( 1 )+1;
206 unsigned int cols ()
const
208 return matrix_.gethighbound( 2 )+1;
211 const Field *
rowPtr (
const unsigned int row )
const
214 const int lastCol = matrix_.gethighbound( 2 );
215 ap::const_raw_vector< Field > rowVector = matrix_.getrow( row, 0, lastCol );
216 assert( (rowVector.GetStep() == 1) && (rowVector.GetLength() == lastCol+1) );
217 return rowVector.GetData();
223 const int lastCol = matrix_.gethighbound( 2 );
224 ap::raw_vector< Field > rowVector = matrix_.getrow( row, 0, lastCol );
225 assert( (rowVector.GetStep() == 1) && (rowVector.GetLength() == lastCol+1) );
226 return rowVector.GetData();
231 matrix_.setbounds( 0, rows-1, 0, cols-1 );
238 matinv::matinvreport< precision > report;
239 matinv::rmatrixinverse< precision >( matrix_,
rows(), info, report );
248 template<
class Field,
bool aligned >
249 inline std::ostream &operator<<(std::ostream &out, const LFEMatrix<Field,aligned> &mat)
251 for (
unsigned int r=0;r<mat.rows();++r)
253 out << field_cast<double>(mat(r,0));
254 for (
unsigned int c=1;c<mat.cols();++c)
264 #endif // #ifndef DUNE_ALGLIB_MATRIX_HH