dune-localfunctions  2.2.1
lfematrix.hh
Go to the documentation of this file.
1 #ifndef DUNE_ALGLIB_MATRIX_HH
2 #define DUNE_ALGLIB_MATRIX_HH
3 
4 #include <cassert>
5 #include <vector>
6 
7 #include "field.hh"
8 
9 #if HAVE_ALGLIB
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)
13 #endif
14 
15 namespace Dune
16 {
17 
18  template< class F, bool aligned = false >
19  class LFEMatrix;
20 
21  template< class F, bool aligned >
22  class LFEMatrix
23  {
25  typedef std::vector< F > Row;
26  typedef std::vector<Row> RealMatrix;
27 
28  public:
29  typedef F Field;
30 
31  operator const RealMatrix & () const
32  {
33  return matrix_;
34  }
35 
36  operator RealMatrix & ()
37  {
38  return matrix_;
39  }
40 
41  template <class Vector>
42  void row( const unsigned int row, Vector &vec ) const
43  {
44  assert(row<rows());
45  for (int i=0;i<cols();++i)
46  field_cast(matrix_[row][i], vec[i]);
47  }
48 
49  const Field &operator() ( const unsigned int row, const unsigned int col ) const
50  {
51  assert(row<rows());
52  assert(col<cols());
53  return matrix_[ row ][ col ];
54  }
55 
56  Field &operator() ( const unsigned int row, const unsigned int col )
57  {
58  assert(row<rows());
59  assert(col<cols());
60  return matrix_[ row ][ col ];
61  }
62 
63  unsigned int rows () const
64  {
65  return rows_;
66  }
67 
68  unsigned int cols () const
69  {
70  return cols_;
71  }
72 
73  const Field *rowPtr ( const unsigned int row ) const
74  {
75  assert(row<rows());
76  return &(matrix_[row][0]);
77  }
78 
79  Field *rowPtr ( const unsigned int row )
80  {
81  assert(row<rows());
82  return &(matrix_[row][0]);
83  }
84 
85  void resize ( const unsigned int rows, const unsigned int cols )
86  {
87  matrix_.resize(rows);
88  for (unsigned int i=0;i<rows;++i)
89  matrix_[i].resize(cols);
90  rows_ = rows;
91  cols_ = cols;
92  }
93 
94  bool invert ()
95  {
96  assert( rows() == cols() );
97  std::vector<unsigned int> p(rows());
98  for (unsigned int j=0;j<rows();++j)
99  p[j] = j;
100  for (unsigned int j=0;j<rows();++j)
101  {
102  // pivot search
103  unsigned int r = j;
104  Field max = std::abs( (*this)(j,j) );
105  for (unsigned int i=j+1;i<rows();++i)
106  {
107  if ( std::abs( (*this)(i,j) ) > max )
108  {
109  max = std::abs( (*this)(i,j) );
110  r = i;
111  }
112  }
113  if (max == Zero<Field>())
114  return false;
115  // row swap
116  if (r > j)
117  {
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] );
121  }
122  // transformation
123  Field hr = Unity<Field>()/(*this)(j,j);
124  for (unsigned int i=0;i<rows();++i)
125  (*this)(i,j) *= hr;
126  (*this)(j,j) = hr;
127  for (unsigned int k=0;k<cols();++k)
128  {
129  if (k==j) continue;
130  for (unsigned int i=0;i<rows();++i)
131  {
132  if (i==j) continue;
133  (*this)(i,k) -= (*this)(i,j)*(*this)(j,k);
134  }
135  (*this)(j,k) *= -hr;
136  }
137  }
138  // column exchange
139  Row hv(rows());
140  for (unsigned int i=0;i<rows();++i)
141  {
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];
146  }
147  return true;
148  }
149 
150  private:
151  RealMatrix matrix_;
152  unsigned int cols_,rows_;
153  };
154 
155 #if HAVE_ALGLIB
156  template< unsigned int precision, bool aligned >
157  class LFEMatrix< amp::ampf< precision >, aligned >
158  {
159  public:
160  typedef amp::ampf< precision > Field;
161  private:
162  typedef LFEMatrix< amp::ampf< precision >, aligned > This;
163  typedef ap::template_2d_array< Field, aligned > RealMatrix;
164 
165  public:
166  operator const RealMatrix & () const
167  {
168  return matrix_;
169  }
170 
171  operator RealMatrix & ()
172  {
173  return matrix_;
174  }
175 
176  template <class Vector>
177  void row( const unsigned int row, Vector &vec ) const
178  {
179  assert(row<rows());
180  for (unsigned int i=0;i<cols();++i)
181  field_cast(matrix_(row,i), vec[i]);
182  }
183 
184  const Field &operator() ( const unsigned int row, const unsigned int col ) const
185  {
186  assert(row<rows());
187  assert(col<cols());
188  return matrix_( row, col );
189  }
190 
191  Field &operator() ( const unsigned int row, const unsigned int col )
192  {
193  assert(row<rows());
194  assert(col<cols());
195  return matrix_( row, col );
196  }
197 
198  unsigned int rows () const
199  {
200  return matrix_.gethighbound( 1 )+1;
201  }
202 
203  unsigned int cols () const
204  {
205  return matrix_.gethighbound( 2 )+1;
206  }
207 
208  const Field *rowPtr ( const unsigned int row ) const
209  {
210  assert(row<rows());
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();
215  }
216 
217  Field *rowPtr ( const unsigned int row )
218  {
219  assert(row<rows());
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();
224  }
225 
226  void resize ( const unsigned int rows, const unsigned int cols )
227  {
228  matrix_.setbounds( 0, rows-1, 0, cols-1 );
229  }
230 
231  bool invert ()
232  {
233  assert( rows() == cols() );
234  int info;
235  matinv::matinvreport< precision > report;
236  matinv::rmatrixinverse< precision >( matrix_, rows(), info, report );
237  return (info >= 0);
238  }
239 
240  private:
241  RealMatrix matrix_;
242  };
243 #endif
244 
245  template< class Field, bool aligned >
246  inline std::ostream &operator<<(std::ostream &out, const LFEMatrix<Field,aligned> &mat)
247  {
248  for (unsigned int r=0;r<mat.rows();++r)
249  {
250  out << field_cast<double>(mat(r,0));
251  for (unsigned int c=1;c<mat.cols();++c)
252  {
253  out << " , " << field_cast<double>(mat(r,c));
254  }
255  out << std::endl;
256  }
257  return out;
258  }
259 }
260 
261 #endif // #ifndef DUNE_ALGLIB_MATRIX_HH