Skip to content
Snippets Groups Projects
Commit 8ec5fe46 authored by Oliver Sander's avatar Oliver Sander Committed by sander@PCPOOL.MI.FU-BERLIN.DE
Browse files

after some cleanup of UG the linearsolver.cc hack is not necessary anymore

[[Imported from SVN: r1649]]
parent 2b9b2e9d
No related branches found
No related tags found
No related merge requests found
...@@ -16,7 +16,7 @@ staticrod_SOURCES = staticrod.cc ...@@ -16,7 +16,7 @@ staticrod_SOURCES = staticrod.cc
staticrod2_SOURCES = staticrod2.cc staticrod2_SOURCES = staticrod2.cc
rod3d_SOURCES = rod3d.cc rod3d_SOURCES = rod3d.cc
dirneucoupling_SOURCES = dirneucoupling.cc linearsolver.cc dirneucoupling_SOURCES = dirneucoupling.cc
dirneucoupling_CXXFLAGS = $(UG_CPPFLAGS) $(AMIRAMESH_CPPFLAGS) $(MPI_CPPFLAGS) -I$(IPOPT_DIR)/IPOPT/include \ dirneucoupling_CXXFLAGS = $(UG_CPPFLAGS) $(AMIRAMESH_CPPFLAGS) $(MPI_CPPFLAGS) -I$(IPOPT_DIR)/IPOPT/include \
$(LAPACKPP_CPPFLAGS) $(LAPACKPP_CPPFLAGS)
dirneucoupling_LDADD = $(UG_LDFLAGS) $(AMIRAMESH_LDFLAGS) $(UG_LIBS) $(AMIRAMESH_LIBS) $(MPI_LDFLAGS) \ dirneucoupling_LDADD = $(UG_LDFLAGS) $(AMIRAMESH_LDFLAGS) $(UG_LIBS) $(AMIRAMESH_LIBS) $(MPI_LDFLAGS) \
......
// The following line switches of the f2c.h in UG. Boy, is this disgusting!
#define F2C_INCLUDE
#include "linearsolver.hh"
#include "lapackpp.h"
using namespace Dune;
// Solve a small linear system using lapack++
void linearSolver(const Dune::Matrix<Dune::FieldMatrix<double,1,1> >& A,
Dune::BlockVector<Dune::FieldVector<double,1> >& x,
const FieldVector<double,6>& b)
{
assert(A.N()==6);
assert(A.M()%3 == 0);
LaGenMatDouble matrix(A.N(),A.M());
for (int i=0; i<A.N(); i++)
for (int j=0; j<A.M(); j++)
matrix(i,j) = A[i][j];
LaVectorDouble X(A.M());
for (int i=0; i<A.M(); i++)
X(i) = x[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<A.M(); i++)
x[i] = X(i);
}
// Compute the svd using lapack++
void lapackSVD(const FieldMatrix<double,3,3>& A,
FieldMatrix<double,3,3>& U,
FieldVector<double,3>& sigma,
FieldMatrix<double,3,3>& VT)
{
LaGenMatDouble lpA(3,3), lpU(3,3), lpVT(3,3);
LaVectorDouble lpSigma(3);
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
lpA(i,j) = A[i][j];
LaSVD_IP(lpA, lpSigma, lpU, lpVT);
for (int i=0; i<3; i++) {
sigma[i] = lpSigma(i);
for (int j=0; j<3; j++) {
U[i][j] = lpU(i,j);
VT[i][j] = lpVT(i,j);
}
}
}
#ifndef LINEAR_SOLVER_HH
#define LINEAR_SOLVER_HH
#include <dune/common/fmatrix.hh>
#include <dune/istl/matrix.hh>
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,
Dune::FieldMatrix<double,3,3>& U,
Dune::FieldVector<double,3>& sigma,
Dune::FieldMatrix<double,3,3>& VT);
#endif
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment