diff --git a/linearsolver.cc b/linearsolver.cc index d8dac993b5ed4cdbcf9e467f550b1b9533267f9f..2c84d9ab321c4c07f6e3dadd18d30280d4a7095c 100644 --- a/linearsolver.cc +++ b/linearsolver.cc @@ -7,30 +7,30 @@ using namespace Dune; // Solve a small linear system using lapack++ -void linearSolver(const FieldMatrix<double,6,12>& A, - FieldVector<double,12>& x, +void linearSolver(const Dune::Matrix<Dune::FieldMatrix<double,1,1> >& A, + Dune::BlockVector<Dune::FieldVector<double,1> >& x, const FieldVector<double,6>& b) { - int N = 6; - int M = 12; + assert(A.N()==6); + assert(A.M()%3 == 0); - LaGenMatDouble matrix(6,12); + LaGenMatDouble matrix(A.N(),A.M()); - for (int i=0; i<N; i++) - for (int j=0; j<M; j++) + for (int i=0; i<A.N(); i++) + for (int j=0; j<A.M(); j++) matrix(i,j) = A[i][j]; - LaVectorDouble X(M); - for (int i=0; i<M; i++) + LaVectorDouble X(A.M()); + for (int i=0; i<A.M(); i++) X(i) = x[i]; - LaVectorDouble B(N); - for (int i=0; i<N; i++) + LaVectorDouble B(A.N()); + for (int i=0; i<A.N(); i++) B(i) = b[i]; LaLinearSolve(matrix, X, B); - for (int i=0; i<M; i++) + for (int i=0; i<A.M(); i++) x[i] = X(i); } diff --git a/linearsolver.hh b/linearsolver.hh index b66376e1105d70e35dd9583246a2914282feb05c..bf6ecea6d9fb391bda424335687835dd82771ee7 100644 --- a/linearsolver.hh +++ b/linearsolver.hh @@ -2,9 +2,10 @@ #define LINEAR_SOLVER_HH #include <dune/common/fmatrix.hh> +#include <dune/istl/matrix.hh> -void linearSolver(const Dune::FieldMatrix<double,6,12>& A, - Dune::FieldVector<double,12>& x, +void linearSolver(const Dune::Matrix<Dune::FieldMatrix<double,1,1> >& A, + Dune::BlockVector<Dune::FieldVector<double,1> >& x, const Dune::FieldVector<double,6>& b); void lapackSVD(const Dune::FieldMatrix<double,3,3>& A,