1 #ifndef DUNE_ALGLIB_MATRIX_HH
2 #define DUNE_ALGLIB_MATRIX_HH
10 #include <alglib/amp.h>
11 #include <alglib/matinv.h>
12 #warning ALGLIB support is deprecated and will be dropped after DUNE 2.2 (cf. FS#931)
18 template<
class F,
bool aligned = false >
21 template<
class F,
bool aligned >
25 typedef std::vector< F > Row;
26 typedef std::vector<Row> RealMatrix;
31 operator const RealMatrix & ()
const
36 operator RealMatrix & ()
41 template <
class Vector>
42 void row(
const unsigned int row, Vector &vec )
const
45 for (
int i=0;i<
cols();++i)
53 return matrix_[
row ][ col ];
60 return matrix_[
row ][ col ];
76 return &(matrix_[
row][0]);
82 return &(matrix_[
row][0]);
88 for (
unsigned int i=0;i<
rows;++i)
97 std::vector<unsigned int> p(
rows());
98 for (
unsigned int j=0;j<
rows();++j)
100 for (
unsigned int j=0;j<
rows();++j)
104 Field max = std::abs( (*
this)(j,j) );
105 for (
unsigned int i=j+1;i<
rows();++i)
107 if ( std::abs( (*
this)(i,j) ) > max )
109 max = std::abs( (*
this)(i,j) );
118 for (
unsigned int k=0;k<
cols();++k)
119 std::swap( (*
this)(j,k), (*
this)(r,k) );
120 std::swap( p[j], p[r] );
124 for (
unsigned int i=0;i<
rows();++i)
127 for (
unsigned int k=0;k<
cols();++k)
130 for (
unsigned int i=0;i<
rows();++i)
133 (*this)(i,k) -= (*
this)(i,j)*(*
this)(j,k);
140 for (
unsigned int i=0;i<
rows();++i)
142 for (
unsigned int k=0;k<
rows();++k)
143 hv[ p[k] ] = (*
this)(i,k);
144 for (
unsigned int k=0;k<
rows();++k)
145 (*
this)(i,k) = hv[k];
152 unsigned int cols_,rows_;
156 template<
unsigned int precision,
bool aligned >
157 class LFEMatrix< amp::ampf< precision >, aligned >
160 typedef amp::ampf< precision >
Field;
162 typedef LFEMatrix< amp::ampf< precision >, aligned > This;
163 typedef ap::template_2d_array< Field, aligned > RealMatrix;
166 operator const RealMatrix & ()
const
171 operator RealMatrix & ()
176 template <
class Vector>
177 void row(
const unsigned int row, Vector &vec )
const
180 for (
unsigned int i=0;i<
cols();++i)
184 const Field &
operator() (
const unsigned int row,
const unsigned int col )
const
188 return matrix_( row, col );
195 return matrix_( row, col );
198 unsigned int rows ()
const
200 return matrix_.gethighbound( 1 )+1;
203 unsigned int cols ()
const
205 return matrix_.gethighbound( 2 )+1;
208 const Field *
rowPtr (
const unsigned int row )
const
211 const int lastCol = matrix_.gethighbound( 2 );
212 ap::const_raw_vector< Field > rowVector = matrix_.getrow( row, 0, lastCol );
213 assert( (rowVector.GetStep() == 1) && (rowVector.GetLength() == lastCol+1) );
214 return rowVector.GetData();
220 const int lastCol = matrix_.gethighbound( 2 );
221 ap::raw_vector< Field > rowVector = matrix_.getrow( row, 0, lastCol );
222 assert( (rowVector.GetStep() == 1) && (rowVector.GetLength() == lastCol+1) );
223 return rowVector.GetData();
228 matrix_.setbounds( 0, rows-1, 0, cols-1 );
235 matinv::matinvreport< precision > report;
236 matinv::rmatrixinverse< precision >( matrix_,
rows(), info, report );
245 template<
class Field,
bool aligned >
246 inline std::ostream &operator<<(std::ostream &out, const LFEMatrix<Field,aligned> &mat)
248 for (
unsigned int r=0;r<mat.rows();++r)
250 out << field_cast<double>(mat(r,0));
251 for (
unsigned int c=1;c<mat.cols();++c)
261 #endif // #ifndef DUNE_ALGLIB_MATRIX_HH