Task 1 (Conjugate gradient with MPI) Assignment

Scaffold Head

//MPI CG SOLVER #include <iostream> #include <cstdint> #include <vector> #include <string> #include <algorithm> #include <limits> #include <mpi.h> #include <numeric> #define EPSILON (0.001) template <class vector_T> void operator+=(vector_T& vec, const vector_T& vec2){ if(vec.size() != vec2.size()){ std::cout << "in operator '+=': dimensions are not compatible!" << std::endl; exit(1); } for(int i=0;i<vec.size();++i) vec[i]+=vec2[i]; } template <class vector_T> void operator-=(vector_T& vec, const vector_T& vec2){ if(vec.size() != vec2.size()){ std::cout << "in operator '-=': dimensions are not compatible!" << std::endl; exit(1); } for(int i=0;i<vec.size();++i) vec[i]-=vec2[i]; } template <class vector_T> vector_T operator+(const vector_T& vec, const vector_T& vec2){ if(vec.size() != vec2.size()){ std::cout << "in operator '+': dimensions are not compatible!" << std::endl; exit(1); } vector_T tmp(vec.size()); for(int i=0;i<vec.size();++i) tmp[i]=vec[i]+vec2[i]; return tmp; } template <class vector_T> vector_T operator-(const vector_T& vec, const vector_T& vec2){ if(vec.size() != vec2.size()){ std::cout << "in operator '-': dimensions are not compatible!" << std::endl; exit(1); } vector_T tmp(vec.size()); for(int i=0;i<vec.size();++i) tmp[i]=vec[i]-vec2[i]; return tmp; } template <class vector_T> vector_T operator/(const vector_T& vec, const vector_T& vec2){ if(vec.size() != vec2.size() || vec.size() != 1){ std::cout << "in operator '/': dimensions are not compatible or different from 1!" << std::endl; exit(1); } if(vec.size() == 1 && vec2.size() == 1) return vector_T{vec[0]/vec2[0]}; } // This function is able to compute matrix * vector, vector * vector, vector * scalar and matrix * scalar products. // Obviously you should make sure that your dimension are actually compatible for the desired operation template <class vector_T> vector_T operator*(const vector_T& mat, const vector_T& vec){ size_t n = vec.size(); if(n==0){ std::cout << "You are trying to multiply with an empty vector! This should not happen!" << std::endl; exit(1); } size_t m = mat.size()/n; if(mat.size()%n != 0){ std::cout << "in operator '*': dimensions are not compatible!" << std::endl; exit(1); } // we assume that mat has m*n entries where m will be the output dimension of this operation. // if this is not the case something bad will probably happen! vector_T tmp(m); for(int i=0;i<m;++i){ for(int j=0;j<n;++j){ tmp[i] += mat[i*n+j]*vec[j]; } } return tmp; } int main(int argc, char * argv[]){ MPI::Init(argc, argv); // Matrix of size 0.75GB (Matrix kept relatively small due to timeout reasons) // Note: We need a square matrix for the CG-method to succeed const size_t global_m = 10000; const size_t num_ranks = MPI::COMM_WORLD.Get_size(); const size_t rank = MPI::COMM_WORLD.Get_rank(); const bool is_root = rank == 0; // Why do we multiplicate by is_root? // Setup the complete equation system std::vector<double> global_A(global_m*global_m*is_root); std::vector<double> global_x(global_m*is_root); std::vector<double> global_b(global_m*is_root); // Fill data into the linear equation system if(is_root){ // We create a diagonal matrix, // That way we can ensure that the equation system is solvable for(int i = 0; i< global_m; i++){ global_A[i*global_m+i] = i%100+1; } // Fill b with values std::iota(global_b.begin(), global_b.end(),1); }

Scaffold Foot

if(is_root) printf("#Final error: %f \n", err_v/global_m); if(is_root) printf("#Performed %lu Iterations\n", k); if (is_root) { std::cout << "# Final error: " << err_v/global_m << std::endl; std::cout << "\nParallel programming is " << ((err_v/global_m < EPSILON && k < 50) ? "fun!" : "error-prone!") << std::endl; } MPI::Finalize(); }

Start time:

Do 14 Nov 2019 12:00:00

End time:

Mi 20 Nov 2019 11:00:00

General test timeout:

10.0 seconds

Tests

Parallel programming is fun!