diff --git a/ex3_benchmarks/geom.h b/ex3_benchmarks/geom.h new file mode 100644 index 0000000..14328fa --- /dev/null +++ b/ex3_benchmarks/geom.h @@ -0,0 +1,381 @@ +#ifndef GEOM_FILE +#define GEOM_FILE +#include +#include // function; C++11 +#include +#include + +/** + * Basis class for finite element meshes. + */ +class Mesh +{ +public: + /** + * Constructor initializing the members with default values. + * + * @param[in] ndim space dimensions (dimension for coordinates) + * @param[in] nvert_e number of vertices per element (dimension for connectivity) + * @param[in] ndof_e degrees of freedom per element (= @p nvert_e for linear elements) + */ + explicit Mesh(int ndim, int nvert_e = 0, int ndof_e = 0); + + /** + * Destructor. + * + * See clang warning on + * weak-vtables. + */ + virtual ~Mesh(); + + /** + * Number of finite elements in (sub)domain. + * @return number of elements. + */ + int Nelems() const + { + return _nelem; + } + + /** + * Global number of vertices for each finite element. + * @return number of vertices per element. + */ + int NverticesElements() const + { + return _nvert_e; + } + + /** + * Global number of degrees of freedom (dof) for each finite element. + * @return degrees of freedom per element. + */ + int NdofsElement() const + { + return _ndof_e; + } + + /** + * Number of vertices in mesh. + * @return number of vertices. + */ + int Nnodes() const + { + return _nnode; + } + + /** + * Space dimension. + * @return number of dimensions. + */ + int Ndims() const + { + return _ndim; + } + + /** + * (Re-)Allocates memory for the element connectivity and redefines the appropriate dimensions. + * + * @param[in] nelem number of elements + * @param[in] nvert_e number of vertices per element + */ + void Resize_Connectivity(int nelem, int nvert_e) + { + SetNelem(nelem); // number of elements + SetNverticesElement(nvert_e); // vertices per element + _ia.resize(nelem * nvert_e); + } + + /** + * Read connectivity information (g1,g2,g3)_i. + * @return convectivity vector [nelems*ndofs]. + */ + const std::vector &GetConnectivity() const + { + return _ia; + } + + /** + * Access/Change connectivity information (g1,g2,g3)_i. + * @return convectivity vector [nelems*ndofs]. + */ + std::vector &GetConnectivity() + { + return _ia; + } + + /** + * (Re-)Allocates memory for the element connectivity and redefines the appropriate dimensions. + * + * @param[in] nnodes number of nodes + * @param[in] ndim space dimension + */ + void Resize_Coords(int nnodes, int ndim) + { + SetNnode(nnodes); // number of nodes + SetNdim(ndim); // space dimension + _xc.resize(nnodes * ndim); + } + + /** + * Read coordinates of vertices (x,y)_i. + * @return coordinates vector [nnodes*2]. + */ + const std::vector &GetCoords() const + { + return _xc; + } + + /** + * Access/Change coordinates of vertices (x,y)_i. + * @return coordinates vector [nnodes*2]. + */ + std::vector &GetCoords() + { + return _xc; + } + + /** + * Calculate values in vector @p v via function @p func(x,y) + * @param[in] v vector + * @param[in] func function of (x,y) returning a double value. + */ + void SetValues(std::vector &v, const std::function &func) const; + + /** + * Prints the information for a finite element mesh + */ + void Debug() const; + + /** + * Determines the indices of those vertices with Dirichlet boundary conditions + * @return index vector. + */ + virtual std::vector Index_DirichletNodes() const = 0; + + /** + * Write vector @p v toghether with its mesh information to an ASCii file @p fname. + * + * The data are written in C-style. + * + * @param[in] fname file name + * @param[in] v vector + */ + void Write_ascii_matlab(std::string const &fname, std::vector const &v) const; + + /** + * Visualize @p v together with its mesh information via matlab or octave. + * + * Comment/uncomment those code lines in method Mesh:Visualize (geom.cpp) + * that are supported on your system. + * + * @param[in] v vector + * + * @warning matlab files ascii_read_meshvector.m visualize_results.m + * must be in the executing directory. + */ + void Visualize(std::vector const &v) const; + + +protected: + void SetNelem(int nelem) + { + _nelem = nelem; + } + + void SetNverticesElement(int nvert) + { + _nvert_e = nvert; + } + + void SetNdofsElement(int ndof) + { + _ndof_e = ndof; + } + + void SetNnode(int nnode) + { + _nnode = nnode; + } + + void SetNdim(int ndim) + { + _ndim = ndim; + } + +private: + int _nelem; //!< number elements + int _nvert_e; //!< number of vertices per element + int _ndof_e; //!< degrees of freedom (d.o.f.) per element + int _nnode; //!< number nodes/vertices + int _ndim; //!< space dimension of the problem (1, 2, or 3) + std::vector _ia; //!< element connectivity + std::vector _xc; //!< coordinates +}; + +/** + * 2D finite element mesh of the square consiting of linear triangular elements. + */ +class Mesh_2d_3_square: public Mesh +{ +public: + /** + * Generates the f.e. mesh for the unit square. + * + * @param[in] nx number of discretization intervals in x-direction + * @param[in] ny number of discretization intervals in y-direction + * @param[in] myid my MPI-rank / subdomain + * @param[in] procx number of ranks/subdomains in x-direction + * @param[in] procy number of processes in y-direction + */ + Mesh_2d_3_square(int nx, int ny, int myid = 0, int procx = 1, int procy = 1); + + /** + * Destructor + */ + ~Mesh_2d_3_square() override + {} + + /** + * Set solution vector based on a tensor product grid in the rectangle. + * @param[in] u solution vector + */ + void SetU(std::vector &u) const; + + /** + * Set right hand side (rhs) vector on a tensor product grid in the rectangle. + * @param[in] f rhs vector + */ + void SetF(std::vector &f) const; + + /** + * Determines the indices of those vertices with Dirichlet boundary conditions + * @return index vector. + */ + std::vector Index_DirichletNodes() const override; + + /** + * Stores the values of vector @p u of (sub)domain into a file @p name for further processing in gnuplot. + * The file stores rowise the x- and y- coordinates together with the value from @p u . + * The domain [@p xl, @p xr] x [@p yb, @p yt] is discretized into @p nx x @p ny intervals. + * + * @param[in] name basename of file name (file name will be extended by the rank number) + * @param[in] u local vector + * + * @warning Assumes tensor product grid in unit square; rowise numbered + * (as generated in class constructor). + * The output is provided for tensor product grid visualization + * ( similar to Matlab-surf() ). + * + * @see Mesh_2d_3_square + */ + void SaveVectorP(std::string const &name, std::vector const &u) const; + + // here will still need to implement in the class + // GetBound(), AddBound() + // or better a generalized way with indices and their appropriate ranks for MPI communication + +private: + /** + * Determines the coordinates of the dicretization nodes of the domain [@p xl, @p xr] x [@p yb, @p yt] + * which is discretized into @p nx x @p ny intervals. + * + * @param[in] ny number of discretization intervals in y-direction + * @param[in] xl x-coordinate of left boundary + * @param[in] xr x-coordinate of right boundary + * @param[in] yb y-coordinate of lower boundary + * @param[in] yt y-coordinate of upper boundary + * @param[out] xc coordinate vector of length 2n with x(2*k,2*k+1) as coodinates of node k + */ + + void GetCoordsInRectangle(int nx, int ny, double xl, double xr, double yb, double yt, + double xc[]); + /** + * Determines the element connectivity of linear triangular elements of a FEM discretization + * of a rectangle using @p nx x @p ny equidistant intervals for discretization. + * @param[in] nx number of discretization intervals in x-direction + * @param[in] ny number of discretization intervals in y-direction + * @param[out] ia element connectivity matrix with ia(3*s,3*s+1,3*s+2) as node numbers od element s + */ + void GetConnectivityInRectangle(int nx, int ny, int ia[]); + +private: + int _myid; //!< my MPI rank + int _procx; //!< number of MPI ranks in x-direction + int _procy; //!< number of MPI ranks in y-direction + std::array _neigh; //!< MPI ranks of neighbors (negative: no neighbor but b.c.) + int _color; //!< red/black coloring (checker board) of subdomains + + double _xl; //!< x coordinate of lower left corner of square + double _xr; //!< x coordinate of lower right corner of square + double _yb; //!< y coordinate or lower left corner of square + double _yt; //!< y coordinate of upper right corner of square + int _nx; //!< number of intervals in x-direction + int _ny; //!< number of intervals in y-direction +}; + +// #################### still some old code (--> MPI) ############################ +/** + * Copies the values of @p w corresponding to boundary @p ib + * onto vector s. South (ib==1), East (ib==2), North (ib==3), West (ib==4). + * The vector @p s has to be long enough!! + * @param[in] ib my local boundary + * @param[in] nx number of discretization intervals in x-direction + * @param[in] ny number of discretization intervals in y-direction + * @param[in] w vector for all nodes of local discretization + * @param[out] s short vector with values on boundary @p ib +*/ +// GH_NOTE: Absicherung bei s !! +void GetBound(int ib, int nx, int ny, double const w[], double s[]); + + +/** + * Computes @p w := @p w + @p s at the interface/boundary nodes on the + * boundary @p ib . South (ib==1), East (ib==2), North (ib==3), West (ib==4) + * @param[in] ib my local boundary + * @param[in] nx number of discretization intervals in x-direction + * @param[in] ny number of discretization intervals in y-direction + * @param[in,out] w vector for all nodes of local discretization + * @param[in] s short vector with values on boundary @p ib +*/ +void AddBound(int ib, int nx, int ny, double w[], double const s[]); + +// #################### Mesh from Matlab ############################ +/** + * 2D finite element mesh of the square consiting of linear triangular elements. + */ +class Mesh_2d_3_matlab: public Mesh +{ +public: + /** + * Reads mesh data from a binary file. + * + * File format, see ascii_write_mesh.m + * + * @param[in] fname file name + */ + explicit Mesh_2d_3_matlab(std::string const &fname); + + /** + * Determines the indices of those vertices with Dirichlet boundary conditions. + * @return index vector. + * + * @warning All boundary nodes are considered as Dirchlet nodes. + */ + std::vector Index_DirichletNodes() const override; + +private: + /** + * Determines the indices of those vertices with Dirichlet boundary conditions + * @return index vector. + */ + int Nnbedges() const + { + return static_cast(bedges.size()); + } + + std::vector bedges; //!< boundary edges [nbedges][2] storing start/end vertex + +}; + +#endif diff --git a/ex3_benchmarks/getmatrix.cpp b/ex3_benchmarks/getmatrix.cpp new file mode 100644 index 0000000..479caa5 --- /dev/null +++ b/ex3_benchmarks/getmatrix.cpp @@ -0,0 +1,347 @@ +#include "getmatrix.h" +#include "userset.h" + +#include +#include +#include +#include +#include +#include +#include +using namespace std; + + +// general routine for lin. triangular elements + +void CalcElem(int const ial[3], double const xc[], double ske[3][3], double fe[3]) +//void CalcElem(const int* __restrict__ ial, const double* __restrict__ xc, double* __restrict__ ske[3], double* __restrict__ fe) +{ + const int i1 = 2 * ial[0], i2 = 2 * ial[1], i3 = 2 * ial[2]; + const double x13 = xc[i3 + 0] - xc[i1 + 0], y13 = xc[i3 + 1] - xc[i1 + 1], + x21 = xc[i1 + 0] - xc[i2 + 0], y21 = xc[i1 + 1] - xc[i2 + 1], + x32 = xc[i2 + 0] - xc[i3 + 0], y32 = xc[i2 + 1] - xc[i3 + 1]; + const double jac = fabs(x21 * y13 - x13 * y21); + + ske[0][0] = 0.5 / jac * (y32 * y32 + x32 * x32); + ske[0][1] = 0.5 / jac * (y13 * y32 + x13 * x32); + ske[0][2] = 0.5 / jac * (y21 * y32 + x21 * x32); + ske[1][0] = ske[0][1]; + ske[1][1] = 0.5 / jac * (y13 * y13 + x13 * x13); + ske[1][2] = 0.5 / jac * (y21 * y13 + x21 * x13); + ske[2][0] = ske[0][2]; + ske[2][1] = ske[1][2]; + ske[2][2] = 0.5 / jac * (y21 * y21 + x21 * x21); + + const double xm = (xc[i1 + 0] + xc[i2 + 0] + xc[i3 + 0]) / 3.0, + ym = (xc[i1 + 1] + xc[i2 + 1] + xc[i3 + 1]) / 3.0; + //fe[0] = fe[1] = fe[2] = 0.5 * jac * FunctF(xm, ym) / 3.0; + fe[0] = fe[1] = fe[2] = 0.5 * jac * fNice(xm, ym) / 3.0; +} + + +// general routine for lin. triangular elements, +// non-symm. matrix +// node numbering in element: a s c e n d i n g indices !! +// GH: deprecated +void AddElem(int const ial[3], double const ske[3][3], double const fe[3], + int const id[], int const ik[], double sk[], double f[]) +{ + for (int i = 0; i < 3; ++i) + { + const int ii = ial[i], // row ii (global index) + id1 = id[ii], // start and + id2 = id[ii + 1]; // end of row ii in matrix + int ip = id1; + for (int j = 0; j < 3; ++j) // no symmetry assumed + { + const int jj = ial[j]; + bool not_found = true; + do // find entry jj (global index) in row ii + { + not_found = (ik[ip] != jj); + ++ip; + } + while (not_found && ip < id2); + +#ifndef NDEBUG // compiler option -DNDEBUG switches off the check + if (not_found) // no entry found !! + { + cout << "Error in AddElem: (" << ii << "," << jj << ") [" + << ial[0] << "," << ial[1] << "," << ial[2] << "]\n"; + assert(!not_found); + } +#endif + sk[ip - 1] += ske[i][j]; + } + f[ii] += fe[i]; + } +} + + +// ---------------------------------------------------------------------------- + + + + +// #################################################################### + +CRS_Matrix::CRS_Matrix(Mesh const &mesh) + : _mesh(mesh), _nrows(0), _nnz(0), _id(0), _ik(0), _sk(0) +{ + Derive_Matrix_Pattern(); + return; +} + +void CRS_Matrix::Derive_Matrix_Pattern() +{ + int const nelem(_mesh.Nelems()); + int const ndof_e(_mesh.NdofsElement()); + auto const &ia(_mesh.GetConnectivity()); +// Determine the number of matrix rows + _nrows = *max_element(ia.cbegin(), ia.cbegin() + ndof_e * nelem); + ++_nrows; // node numberng: 0 ... nnode-1 + assert(*min_element(ia.cbegin(), ia.cbegin() + ndof_e * nelem) == 0); // numbering starts with 0 ? + +// Collect for each node those nodes it is connected to (multiple entries) +// Detect the neighboring nodes + vector< list > cc(_nrows); // cc[i] is the list of nodes a node i is connected to + for (int i = 0; i < nelem; ++i) + { + int const idx = ndof_e * i; + for (int k = 0; k < ndof_e; ++k) + { + list &cck = cc.at(ia[idx + k]); + cck.insert( cck.end(), ia.cbegin() + idx, ia.cbegin() + idx + ndof_e ); + } + } +// Delete the multiple entries + _nnz = 0; + for (auto &it : cc) + { + it.sort(); + it.unique(); + _nnz += static_cast(it.size()); + // cout << it.size() << " :: "; copy(it->begin(),it->end(), ostream_iterator(cout," ")); cout << endl; + } + +// CSR data allocation + _id.resize(_nrows + 1); // Allocate memory for CSR row pointer + _ik.resize(_nnz); // Allocate memory for CSR column index vector + +// copy CSR data + _id[0] = 0; // begin of first row + for (size_t i = 0; i < cc.size(); ++i) + { + //cout << i << " " << nid.at(i) << endl;; + const list &ci = cc.at(i); + const auto nci = static_cast(ci.size()); + _id[i + 1] = _id[i] + nci; // begin of next line + copy(ci.begin(), ci.end(), _ik.begin() + _id[i] ); + } + + assert(_nnz == _id[_nrows]); + _sk.resize(_nnz); // Allocate memory for CSR column index vector + return; +} + + +void CRS_Matrix::Debug() const +{ +// ID points to first entry of row +// no symmetry assumed + cout << "\nMatrix (nnz = " << _id[_nrows] << ")\n"; + + for (int row = 0; row < _nrows; ++row) + { + cout << "Row " << row << " : "; + int const id1 = _id[row]; + int const id2 = _id[row + 1]; + for (int j = id1; j < id2; ++j) + { + cout.setf(ios::right, ios::adjustfield); + cout << "[" << setw(2) << _ik[j] << "] " << setw(4) << _sk[j] << " "; + } + cout << endl; + } + return; +} + +void CRS_Matrix::CalculateLaplace(vector &f) +{ + assert(_mesh.NdofsElement() == 3); // only for triangular, linear elements + //cout << _nnz << " vs. " << _id[_nrows] << " " << _nrows<< endl; + assert(_nnz == _id[_nrows]); + + for (int k = 0; k < _nrows; ++k) + { + _sk[k] = 0.0; + } + for (int k = 0; k < _nrows; ++k) + { + f[k] = 0.0; + } + + double ske[3][3], fe[3]; + // Loop over all elements + auto const nelem = _mesh.Nelems(); + auto const &ia = _mesh.GetConnectivity(); + auto const &xc = _mesh.GetCoords(); + + for (int i = 0; i < nelem; ++i) + { + CalcElem(ia.data() + 3 * i, xc.data(), ske, fe); + //AddElem(ia.data()+3 * i, ske, fe, _id.data(), _ik.data(), _sk.data(), f.data()); // GH: deprecated + AddElem_3(ia.data() + 3 * i, ske, fe, f); + } + + //Debug(); + + return; +} + +void CRS_Matrix::ApplyDirichletBC(std::vector const &u, std::vector &f) +{ + double const PENALTY = 1e6; + auto const idx = _mesh.Index_DirichletNodes(); + int const nidx = static_cast(idx.size()); + + for (int row = 0; row < nidx; ++row) + { + int const k = idx[row]; + int const id1 = fetch(k, k); // Find diagonal entry of row + assert(id1 >= 0); + _sk[id1] += PENALTY; // matrix weighted scaling feasible + f[k] += PENALTY * u[k]; + } + + return; +} + +void CRS_Matrix::GetDiag(vector &d) const +{ + assert( _nrows == static_cast(d.size()) ); + + for (int row = 0; row < _nrows; ++row) + { + const int ia = fetch(row, row); // Find diagonal entry of row + assert(ia >= 0); + d[row] = _sk[ia]; + } + return; +} + +bool CRS_Matrix::Compare2Old(int nnode, int const id[], int const ik[], double const sk[]) const +{ + bool bn = (nnode == _nrows); // number of rows + if (!bn) + { + cout << "######### Error: " << "number of rows" << endl; + } + + bool bz = (id[nnode] == _nnz); // number of non zero elements + if (!bz) + { + cout << "######### Error: " << "number of non zero elements" << endl; + } + + bool bd = equal(id, id + nnode + 1, _id.cbegin()); // row starts + if (!bd) + { + cout << "######### Error: " << "row starts" << endl; + } + + bool bk = equal(ik, ik + id[nnode], _ik.cbegin()); // column indices + if (!bk) + { + cout << "######### Error: " << "column indices" << endl; + } + + bool bv = equal(sk, sk + id[nnode], _sk.cbegin()); // values + if (!bv) + { + cout << "######### Error: " << "values" << endl; + } + + return bn && bz && bd && bk && bv; +} + + +void CRS_Matrix::Mult(vector &w, vector const &u) const +{ + assert( _nrows == static_cast(w.size()) ); + assert( w.size() == u.size() ); + + for (int row = 0; row < _nrows; ++row) + { + double wi = 0.0; + for (int ij = _id[row]; ij < _id[row + 1]; ++ij) + { + wi += _sk[ij] * u[ _ik[ij] ]; + } + w[row] = wi; + } + return; +} + +void CRS_Matrix::Defect(vector &w, + vector const &f, vector const &u) const +{ + assert( _nrows == static_cast(w.size()) ); + assert( w.size() == u.size() && u.size() == f.size() ); + + for (int row = 0; row < _nrows; ++row) + { + double wi = f[row]; + for (int ij = _id[row]; ij < _id[row + 1]; ++ij) + { + wi -= _sk[ij] * u[ _ik[ij] ]; + } + w[row] = wi; + } + return; +} + +int CRS_Matrix::fetch(int const row, int const col) const +{ + int const id2 = _id[row + 1]; // end and + int ip = _id[row]; // start of recent row (global index) + + while (ip < id2 && _ik[ip] != col) // find index col (global index) + { + ++ip; + } + if (ip >= id2) + { + ip = -1; +#ifndef NDEBUG // compiler option -DNDEBUG switches off the check + cout << "No column " << col << " in row " << row << endl; + assert(ip >= id2); +#endif + } + return ip; +} + + +void CRS_Matrix::AddElem_3(int const ial[3], double const ske[3][3], double const fe[3], vector &f) +{ + for (int i = 0; i < 3; ++i) + { + const int ii = ial[i]; // row ii (global index) + for (int j = 0; j < 3; ++j) // no symmetry assumed + { + const int jj = ial[j]; // column jj (global index) + int ip = fetch(ii, jj); // find column entry jj in row ii +#ifndef NDEBUG // compiler option -DNDEBUG switches off the check + if (ip < 0) // no entry found !! + { + cout << "Error in AddElem: (" << ii << "," << jj << ") [" + << ial[0] << "," << ial[1] << "," << ial[2] << "]\n"; + assert(ip >= 0); + } +#endif + _sk[ip] += ske[i][j]; + } + f[ii] += fe[i]; + } +} + diff --git a/ex3_benchmarks/getmatrix.h b/ex3_benchmarks/getmatrix.h new file mode 100644 index 0000000..922606c --- /dev/null +++ b/ex3_benchmarks/getmatrix.h @@ -0,0 +1,178 @@ +#ifndef GETMATRIX_FILE +#define GETMATRIX_FILE + +#include "geom.h" +#include + +/** + * Calculates the element stiffness matrix @p ske and the element load vector @p fe + * of one triangular element with linear shape functions. + * @param[in] ial node indices of the three element vertices + * @param[in] xc vector of node coordinates with x(2*k,2*k+1) as coodinates of node k + * @param[out] ske element stiffness matrix + * @param[out] fe element load vector + */ +void CalcElem(int const ial[3], double const xc[], double ske[3][3], double fe[3]); + +/** + * Adds the element stiffness matrix @p ske and the element load vector @p fe + * of one triangular element with linear shape functions to the appropriate positions in + * the symmetric stiffness matrix, stored as CSR matrix K(@p sk,@p id, @p ik) + * + * @param[in] ial node indices of the three element vertices + * @param[in] ske element stiffness matrix + * @param[in] fe element load vector + * @param[out] sk vector non-zero entries of CSR matrix + * @param[in] id index vector containing the first entry in a CSR row + * @param[in] ik column index vector of CSR matrix + * @param[out] f distributed local vector storing the right hand side + * + * @warning Algorithm requires indices in connectivity @p ial in ascending order. + * Currently deprecated. +*/ +void AddElem(int const ial[3], double const ske[3][3], double const fe[3], + int const id[], int const ik[], double sk[], double f[]); + + +// ##################################################################### +/** + * Square matrix in CRS format (compressed row storage; also named CSR), + * see an introduction. + */ +class CRS_Matrix +{ + public: + /** + * Intializes the CRS matrix structure from the given discetization in @p mesh. + * + * The sparse matrix pattern is generated but the values are 0. + * + * @param[in] mesh given discretization + * + * @warning A reference to the discretization @p mesh is stored inside this class. + * Therefore, changing @p mesh outside requires also + * to call method @p Derive_Matrix_Pattern explicitely. + * + * @see Derive_Matrix_Pattern + */ + explicit CRS_Matrix(Mesh const & mesh); + + /** + * Destructor. + */ + ~CRS_Matrix() + {} + + /** + * Generates the sparse matrix pattern and overwrites the existing pattern. + * + * The sparse matrix pattern is generated but the values are 0. + */ + void Derive_Matrix_Pattern(); + + /** + * Calculates the entries of f.e. stiffness matrix and load/rhs vector @p f for the Laplace operator in 2D. + * No memory is allocated. + * + * @param[in,out] f (preallocated) rhs/load vector + */ + void CalculateLaplace(std::vector &f); + + /** + * Applies Dirichlet boundary conditions to stiffness matrix and to load vector @p f. + * The penalty method + * is used for incorporating the given values @p u. + * + * @param[in] u (global) vector with Dirichlet data + * @param[in,out] f load vector + */ + void ApplyDirichletBC(std::vector const &u, std::vector &f); + + /** + * Extracts the diagonal elemenst of the sparse matrix. + * + * @param[in,out] d (prellocated) vector of diagonal elements + */ + void GetDiag(std::vector &d) const; + + /** + * Performs the matrix-vector product w := K*u. + * + * @param[in,out] w resulting vector (preallocated) + * @param[in] u vector + */ + void Mult(std::vector &w, std::vector const &u) const; + + /** + * Calculates the defect/residuum w := f - K*u. + * + * @param[in,out] w resulting vector (preallocated) + * @param[in] f load vector + * @param[in] u vector + */ + void Defect(std::vector &w, + std::vector const &f, std::vector const &u) const; + + /** + * Number rows in matrix. + * @return number of rows. + */ + int Nrows() const + {return _nrows;} + + /** + * Show the matrix entries. + */ + void Debug() const; + + /** + * Finds in a CRS matrix the access index for an entry at row @p row and column @p col. + * + * @param[in] row row index + * @param[in] col column index + * @return index for element (@p row, @p col). If no appropriate entry exists then -1 will be returned. + * + * @warning assert() stops the function in case that matrix element (@p row, @p col) doesn't exist. + */ + int fetch(int row, int col) const; + + /** + * Adds the element stiffness matrix @p ske and the element load vector @p fe + * of one triangular element with linear shape functions to the appropriate positions in + * the stiffness matrix, stored as CSR matrix K(@p sk,@p id, @p ik). + * + * @param[in] ial node indices of the three element vertices + * @param[in] ske element stiffness matrix + * @param[in] fe element load vector + * @param[in,out] f distributed local vector storing the right hand side + * + * @warning Algorithm assumes linear triangular elements (ndof_e==3). + */ + void AddElem_3(int const ial[3], double const ske[3][3], double const fe[3], std::vector &f); + + /** + * Compare @p this CRS matrix with an external CRS matrix stored in C-Style. + * + * The method prints statements on differences found. + * + * @param[in] nnode row number of external matrix + * @param[in] id start indices of matrix rows of external matrix + * @param[in] ik column indices of external matrix + * @param[in] sk non-zero values of external matrix + * + * @return true iff all data are identical. + */ + bool Compare2Old(int nnode, int const id[], int const ik[], double const sk[]) const; + + private: + Mesh const & _mesh; //!< reference to discretization + int _nrows; //!< number of rows in matrix + int _nnz; //!< number of non-zero entries + std::vector _id; //!< start indices of matrix rows + std::vector _ik; //!< column indices + std::vector _sk; //!< non-zero values + +}; + + +#endif diff --git a/ex3_benchmarks/jacsolve.cpp b/ex3_benchmarks/jacsolve.cpp new file mode 100644 index 0000000..30766d2 --- /dev/null +++ b/ex3_benchmarks/jacsolve.cpp @@ -0,0 +1,61 @@ +#include "vdop.h" +#include "getmatrix.h" +#include "jacsolve.h" + +#include +#include +#include +#include +using namespace std; + +// ##################################################################### +// const int neigh[], const int color, +// const MPI::Intracomm& icomm, +void JacobiSolve(CRS_Matrix const &SK, vector const &f, vector &u) +{ + const double omega = 1.0; + const int maxiter = 1000; + const double tol = 1e-5, // tolerance + tol2 = tol * tol; // tolerance^2 + + int nrows = SK.Nrows(); // number of rows == number of columns + assert( nrows == static_cast(f.size()) && f.size() == u.size() ); + + cout << endl << " Start Jacobi solver for " << nrows << " d.o.f.s" << endl; + // Choose initial guess + for (int k = 0; k < nrows; ++k) + { + u[k] = 0.0; // u := 0 + } + + vector dd(nrows); // matrix diagonal + vector r(nrows); // residual + vector w(nrows); // correction + + SK.GetDiag(dd); // dd := diag(K) + ////DebugVector(dd);{int ijk; cin >> ijk;} + + // Initial sweep + SK.Defect(r, f, u); // r := f - K*u + + vddiv(w, r, dd); // w := D^{-1}*r + double sigma0 = dscapr(w, r); // s0 := + + // Iteration sweeps + int iter = 0; + double sigma = sigma0; + while ( sigma > tol2 * sigma0 && maxiter > iter) + { + ++iter; + vdaxpy(u, u, omega, w ); // u := u + om*w + SK.Defect(r, f, u); // r := f - K*u + vddiv(w, r, dd); // w := D^{-1}*r + sigma = dscapr(w, r); // s0 := +// cout << "Iteration " << iter << " : " << sqrt(sigma/sigma0) << endl; + } + cout << "aver. Jacobi rate : " << exp(log(sqrt(sigma / sigma0)) / iter) << " (" << iter << " iter)" << endl; + cout << "final error: " << sqrt(sigma / sigma0) << " (rel) " << sqrt(sigma) << " (abs)\n"; + + return; +} + diff --git a/ex3_benchmarks/jacsolve.h b/ex3_benchmarks/jacsolve.h new file mode 100644 index 0000000..dfa3802 --- /dev/null +++ b/ex3_benchmarks/jacsolve.h @@ -0,0 +1,18 @@ +#ifndef JACSOLVE_FILE +#define JACSOLVE_FILE +#include "getmatrix.h" +#include + + +/** + * Solves linear system of equations K @p u = @p f via the Jacobi iteration. + * We use a distributed symmetric CSR matrix @p SK and initial guess of the + * solution is set to 0. + * @param[in] SK CSR matrix + * @param[in] f distributed local vector storing the right hand side + * @param[out] u accumulated local vector storing the solution. +*/ +void JacobiSolve(CRS_Matrix const &SK, std::vector const &f, std::vector &u); + + +#endif