Commit 710118f0 by Praetorius, Simon

### Remove exercise folder

parent f3fa3658
 # Exercises for the lecture SCPROG in the winter term 2019/20 **Tutor:** Dr. Simon Praetorius (simon.praetorius@tu-dresden.de), user @spraetor There will be an exercise sheet every week with some exercises to be worked on during the tutorial and some to be submitted for review. The submission procedure is briefly described on the lecture repository [README.md](/README.md) page and follows the way of code submission in many (scientific) software projects. Thus, it is a good way of learning and training collaborative coding. Code submission and code storage is organized using the [Git](https://git-scm.com/) version control system and the data is hosted at the MatNat [GitLab](https://gitlab.mn.tu-dresden.de) platform. Access to this platform is established during the first weeks in the semester. Additionally to the submission of your own code, the exercise sheets and lecture notes are provided in this Git repository. \ No newline at end of file
 #include #include "linear_algebra.hh" int main(int argc, char** argv) { using namespace scprog; std::size_t n = 50; // right-hand side vector b DenseVector b(n*n, 1.0); // system matrix DenseMatrix A; laplacian_setup(A,n,n); // solution vector DenseVector x(n*n); // iteration object with max_iter and rtol BasicIteration iter(b, 1000, 1.e-6, 0, 10); // solve the linear system using a conjugate-gradient algorithm int err = cg(A, x, b, iter); if (err) { std::cerr << "ERROR: Linear system could not be solved.\n" << " |b-A*x| = " << iter.resid() << ", |b-A*x|/|b| = " << iter.relresid() << std::endl; std::abort(); } return err; } \ No newline at end of file
 #include #include "linear_algebra.hh" int main(int argc, char** argv) { using namespace scprog; std::size_t n = 50; // right-hand side vector b DenseVector b(n, 1.0); // system matrix DenseMatrix A; laplacian_setup(A,n,n); // solution vector DenseVector x(n); // iteration object with max_iter and rtol BasicIteration iter(b, 0, 1.e-6, 0, 10); // solve the linear system using a conjugate-gradient algorithm int err = cg(A, x, b, iter); if (err) { std::cerr << "ERROR: Linear system could not be solved.\n" << " |b-A*x| = " << iter.resid() << ", |b-A*x|/|b| = " << iter.relresid() << std::endl; std::abort(); } return err; } \ No newline at end of file
 #include #include "linear_algebra.hh" namespace scprog { // set all entries of the vector to value v DenseVector& DenseVector::operator=(value_type v) { for (auto& v_i : data_) v_i = v; return *this; } // perform update-assignment elementwise += DenseVector& DenseVector::operator+=(DenseVector const& that) { assert(size() == that.size()); for (size_type i = 0; i < size(); ++i) data_[i] += that.data_[i]; return *this; } // perform update-assignment elementwise += DenseVector& DenseVector::operator-=(DenseVector const& that) { assert(size() == that.size()); for (size_type i = 0; i < size(); ++i) data_[i] -= that.data_[i]; return *this; } // perform update-assignment elementwise *= with a scalar DenseVector& DenseVector::operator*=(value_type s) { for (size_type i = 0; i < size(); ++i) data_[i] *= s; return *this; } // perform update-assignment elementwise /= with a scalar DenseVector& DenseVector::operator/=(value_type s) { assert(s != value_type(0)); for (size_type i = 0; i < size(); ++i) data_[i] /= s; return *this; } // computes Y = a*X + Y. void DenseVector::axpy(value_type a, DenseVector const& x) { assert(size() == x.size()); for (size_type i = 0; i < data_.size(); ++i) data_[i] += a * x.data_[i]; } // computes Y = a*Y + X. void DenseVector::aypx(value_type a, DenseVector const& x) { assert(size() == x.size()); for (size_type i = 0; i < data_.size(); ++i) data_[i] = a * data_[i] + x.data_[i]; } // return the two-norm ||vector||_2 = sqrt(sum_i v_i^2) typename DenseVector::value_type DenseVector::two_norm() const { using std::sqrt; value_type result = 0; for (auto const& d : data_) result += d * d; return sqrt(result); } // return the infinity-norm ||vector||_inf = max_i(|v_i|) typename DenseVector::value_type DenseVector::inf_norm() const { using std::abs; using std::max; value_type result = 0; for (auto const& d : data_) result = max(result, value_type(abs(d))); return result; } // return v^T*v typename DenseVector::value_type DenseVector::unary_dot() const { value_type result = 0; for (auto const& d : data_) result += d*d; return result; } // return v^T*v2 typename DenseVector::value_type DenseVector::dot(DenseVector const& v2) const { assert(v2.size() == size()); value_type result = 0; for (size_type i = 0; i < size(); ++i) result += data_[i] * v2.data_[i]; return result; } // construct a matrix from initializer lists DenseMatrix::DenseMatrix(std::initializer_list> l) { // 1. determine number of entries size_type columns = 0; size_type rows = l.size(); for (auto const& row : l) { if (columns == 0) columns = row.size(); else assert(columns == row.size()); } // 2. insert entries from initializer lists into matrix data_.reserve(rows*columns); for (auto const& row : l) data_.insert(data_.end(), row.begin(), row.end()); } // perform update-assignment elementwise += DenseMatrix& DenseMatrix::operator+=(DenseMatrix const& that) { assert(rows() == that.rows()); assert(cols() == that.cols()); for (size_type i = 0; i < data_.size(); ++i) data_[i] += that.data_[i]; return *this; } // perform update-assignment elementwise += DenseMatrix& DenseMatrix::operator-=(DenseMatrix const& that) { assert(rows() == that.rows()); assert(cols() == that.cols()); for (size_type i = 0; i < data_.size(); ++i) data_[i] -= that.data_[i]; return *this; } // set all entries to v DenseMatrix& DenseMatrix::operator=(value_type v) { for (auto& A_ij : data_) A_ij = v; return *this; } // matrix-vector product A*x DenseVector operator*(DenseMatrix const& A, DenseVector const& x) { using value_type = typename DenseMatrix::value_type; DenseVector y(A.cols(), value_type(0)); A.mult(x, y); return y; } // computes the matrix-vector product, y = Ax. void DenseMatrix::mult(DenseVector const& x, DenseVector& y) const { assert(x.size() == cols()); assert(y.size() == rows()); for (size_type r = 0; r < rows(); ++r) { y[r] = value_type(0); value_type const* row = (*this)[r]; for (size_type c = 0; c < cols(); ++c) y[r] += row[c]*x[c]; } } // computes v3 = v2 + A * v1. void DenseMatrix::mult_add(DenseVector const& v1, DenseVector const& v2, DenseVector& v3) const { assert(v1.size() == cols()); assert(v2.size() == rows()); assert(v3.size() == rows()); for (size_type r = 0; r < rows(); ++r) { v3[r] = v2[r]; value_type const* row = (*this)[r]; for (size_type c = 0; c < cols(); ++c) v3[r] += row[c]*v1[c]; } } // computes Y = a*X + Y. void DenseMatrix::axpy(value_type a, DenseMatrix const& X) { assert(rows() == X.rows()); assert(cols() == X.cols()); for (size_type i = 0; i < data_.size(); ++i) data_[i] += a * X.data_[i]; } // computes Y = a*Y + X. void DenseMatrix::aypx(value_type a, DenseMatrix const& X) { assert(rows() == X.rows()); assert(cols() == X.cols()); for (size_type i = 0; i < data_.size(); ++i) data_[i] = a * data_[i] + X.data_[i]; } // Setup a matrix according to a Laplacian equation on a 2D-grid using a five-point-stencil. // Results in a matrix A of size (m*n) x (m*n) void laplacian_setup(DenseMatrix& A, std::size_t m, std::size_t n) { A.resize(m*n, m*n); A = 0; for (std::size_t i = 0; i < m; i++) { for (std::size_t j = 0; j < n; j++) { std::size_t row = i * n + j; A(row, row) = 4; if (j < n - 1) A(row, row + 1) = -1; if (i < m - 1) A(row, row + n) = -1; if (j > 0) A(row, row - 1) = -1; if (i > 0) A(row, row - n) = -1; } } } // Iteration finished according to residual value r bool BasicIteration::finished(real_type const& r) { bool result = false; if (converged(r)) result = finished_ = true; if (!result) result = check_max(); print_resid(); return result; } bool BasicIteration::check_max() { if (i_ >= max_iter_) error_ = 1, finished_ = true, err_msg_ = "Too many iterations."; return finished_; } bool BasicIteration::converged() const { if (norm_r0_ == 0) return resid_ <= atol_; // ignore relative tolerance if |r0| is zero return resid_ <= rtol_ * norm_r0_ || resid_ <= atol_; } void BasicIteration::print_resid() { if (!quite_ && i_ % cycle_ == 0) { if (i_ != last_print_) { // Avoid multiple print-outs in same iteration std::cout << "iteration " << i_ << ": resid " << resid() << std::endl; last_print_ = i_; } } } int BasicIteration::error_code() const { using std::pow; if (!suppress_) std::cout << "finished! error code = " << error_ << '\n' << iterations() << " iterations\n" << resid() << " is actual final residual. \n" << relresid() << " is actual relative tolerance achieved. \n" << "Relative tol: " << rtol_ << " Absolute tol: " << atol_ << '\n' << "Convergence: " << pow(relresid(), 1.0 / double(iterations())) << std::endl; return error_; } // Apply the conjugate gradient algorithm to the linear system A*x = b and return the number of iterations int cg(DenseMatrix const& A, DenseVector& x, DenseVector const& b, BasicIteration& iter) { using std::abs; using Vector = DenseVector; using Scalar = typename DenseVector::value_type; using Real = typename BasicIteration::real_type; Scalar rho(0), rho_1(0), alpha(0); Vector p(b), q(b), z(b); Vector r(b - A*x); rho = r.unary_dot(); while (! iter.finished(Real(sqrt(abs(rho))))) { ++iter; if (iter.first()) p = r; else p.aypx(rho / rho_1, r); // p = r + (rho / rho_1) * p; q = A * p; alpha = rho / p.dot(q); x.axpy(alpha, p); // x += alpha * p r.axpy(-alpha, q); // r -= alpha * q rho_1 = rho; rho = r.unary_dot(); // rho = r^T * r } return iter; } } // end namespace scprog \ No newline at end of file
 #include #include #include #include #include #include namespace scprog { /// A contiguous vector with vector-space operations class DenseVector { public: using size_type = std::size_t; using value_type = double; using reference = value_type&; using const_reference = value_type const&; using pointer = value_type*; using const_pointer = value_type const*; // ----- constructors / assignment ------------------------------------------- public: /// default constructor, creates an empty vector of size 0 DenseVector() = default; /// constructor of vector with size s and all entries initialized with value v explicit DenseVector(size_type s, value_type v = value_type{}) : data_(s, v) {} /// constructor with vector entries initialized by initializer_list explicit DenseVector(std::initializer_list l) : data_(l.begin(), l.end()) {} /// set all entries of the vector to value v DenseVector& operator=(value_type v); /// resize vector to size s and fill new entries with value v void resize(size_type s, value_type v = value_type{}) { data_.resize(s, v); } /// return the number of elements in the vector size_type size() const { return data_.size(); } // ----- vector-space operations -------------------------------------------- public: /// perform update-assignment elementwise += DenseVector& operator+=(DenseVector const& that); /// perform update-assignment elementwise += DenseVector& operator-=(DenseVector const& that); /// perform update-assignment elementwise *= with a scalar DenseVector& operator*=(value_type s); /// perform update-assignment elementwise /= with a scalar DenseVector& operator/=(value_type s); // ----- element access functions ------------------------------------------- public: /// return a mutable reference to the vector entry v_i reference operator[](size_type i) { assert(i < data_.size()); return data_[i]; } /// return a const reference to the vector entry v_i const_reference operator[](size_type i) const { assert(i < data_.size()); return data_[i]; } // ----- binary operations --------------------------------------------------- public: /// addition of two vectors friend DenseVector operator+(DenseVector lhs, DenseVector const& rhs) { return lhs += rhs; } /// subtraction of two vectors friend DenseVector operator-(DenseVector lhs, DenseVector const& rhs) { return lhs -= rhs; } /// multiplication of the vector with a scalar from the right, i.e. vec * s friend DenseVector operator*(DenseVector vec, value_type s) { return vec *= s; } /// multiplication of the vector with a scalar from the left, i.e. s * vec friend DenseVector operator*(value_type s, DenseVector vec) { return vec *= s; } /// computes Y = a*X + Y. void axpy(value_type a, DenseVector const& X); /// computes Y = a*Y + X. void aypx(value_type a, DenseVector const& X); // ----- reduction operators ------------------------------------------------ public: /// return the two-norm ||vector||_2 = sqrt(sum_i v_i^2) value_type two_norm() const; /// return the infinity-norm ||vector||_inf = max_i(|v_i|) value_type inf_norm() const; /// return v^T*v value_type unary_dot() const; /// return v^T*v2 value_type dot(DenseVector const& v2) const; // ----- data members ------------------------------------------------------- private: std::vector data_; }; /// A dense matrix with row-wise contiguous storage and matrix-matrix as well as /// matrix-vector operations. class DenseMatrix { public: using size_type = std::size_t; using value_type = double; using reference = value_type&; using const_reference = value_type const&; using pointer = value_type*; using const_pointer = value_type const*; // ----- constructors / assignment ------------------------------------------- public: /// default constructor, creates and empty matrix of size 0x0 DenseMatrix() = default; /// constructor of matrix with rows r, columns c and all entries initialized with value v explicit DenseMatrix(size_type r, size_type c, value_type v = value_type{}) : data_(r*c, v) , rows_(r) , cols_(c) {} /// constructor with matrix entries initialized by initializer_list explicit DenseMatrix(std::initializer_list> l); /// set all entries to v DenseMatrix& operator=(value_type v); /// resize matrix to rows r and columns c and fill new entries with value v void resize(size_type r, size_type c, value_type v = value_type{}) { data_.resize(r*c, v); rows_ = r; cols_ = c; } /// return the number of rows in the matrix size_type rows() const { return rows_; } /// return the number of columns in the matrix