From 68bde4ab6ee2994ac87bee9d6ab35141d1e3f08c Mon Sep 17 00:00:00 2001 From: Usman Jamshed Date: Sat, 29 Apr 2023 19:52:23 -0700 Subject: [PATCH] Using csr with prints --- .gitignore | 1 + distributed_pcg.cpp | 346 +++++++++++++++++++++++++++++--------------- 2 files changed, 227 insertions(+), 120 deletions(-) create mode 100644 .gitignore diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..d163863 --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +build/ \ No newline at end of file diff --git a/distributed_pcg.cpp b/distributed_pcg.cpp index 8239917..d21ee56 100644 --- a/distributed_pcg.cpp +++ b/distributed_pcg.cpp @@ -7,61 +7,64 @@ #include -typedef Eigen::SparseMatrix SpMat; // declares a column-major sparse matrix type of double +typedef Eigen::SparseMatrix SpMat; // declares a row-major sparse matrix type of double typedef Eigen::Triplet T; -class MapMatrix{ +class MapMatrix +{ public: - typedef std::pair N2; + typedef std::pair N2; - std::map data; + std::map data; int nbrow; int nbcol; public: - MapMatrix(const int& nr, const int& nc): - nbrow(nr), nbcol(nc) {}; - - MapMatrix(const MapMatrix& m): - nbrow(m.nbrow), nbcol(m.nbcol), data(m.data) {}; - - MapMatrix& operator=(const MapMatrix& m){ - if(this!=&m){ - nbrow=m.nbrow; - nbcol=m.nbcol; - data=m.data; - } - return *this; + MapMatrix(const int &nr, const int &nc) : nbrow(nr), nbcol(nc){}; + + MapMatrix(const MapMatrix &m) : nbrow(m.nbrow), nbcol(m.nbcol), data(m.data){}; + + MapMatrix &operator=(const MapMatrix &m) + { + if (this != &m) + { + nbrow = m.nbrow; + nbcol = m.nbcol; + data = m.data; + } + return *this; } - int NbRow() const {return nbrow;} - int NbCol() const {return nbcol;} + int NbRow() const { return nbrow; } + int NbCol() const { return nbcol; } - double operator()(const int& j, const int& k) const { - auto search = data.find(std::make_pair(j,k)); - if(search!=data.end()) return search->second; + double operator()(const int &j, const int &k) const + { + auto search = data.find(std::make_pair(j, k)); + if (search != data.end()) + return search->second; return 0; } - double& Assign(const int& j, const int& k) { - return data[std::make_pair(j,k)]; + double &Assign(const int &j, const int &k) + { + return data[std::make_pair(j, k)]; } // parallel matrix-vector product with distributed vector xi - std::vector operator*(const std::vector& xi) const { - - + std::vector operator*(const std::vector &xi) const + { std::vector x(NbCol()); - std::copy(xi.begin(),xi.end(),x.begin()); - + std::copy(xi.begin(), xi.end(), x.begin()); - std::vector b(NbRow(),0.); - for(auto it=data.begin(); it!=data.end(); ++it){ + std::vector b(NbRow(), 0.); + for (auto it = data.begin(); it != data.end(); ++it) + { int j = (it->first).first; - int k = (it->first).second; + int k = (it->first).second; double Mjk = it->second; - b[j] += Mjk*x[k]; + b[j] += Mjk * x[k]; } return b; @@ -71,59 +74,80 @@ class MapMatrix{ #include // parallel scalar product (u,v) (u and v are distributed) -double operator,(const std::vector& u, const std::vector& v){ - assert(u.size()==v.size()); - double sp=0.; - for(int j=0; j &u, const std::vector &v) +{ + assert(u.size() == v.size()); + double sp = 0.; + for (int j = 0; j < u.size(); j++) + { + sp += u[j] * v[j]; + } - return sp; + return sp; } // norm of a vector u -double Norm(const std::vector& u) { - return sqrt((u,u)); +double Norm(const std::vector &u) +{ + return sqrt((u, u)); } // addition of two vectors u+v -std::vector operator+(const std::vector& u, const std::vector& v){ - assert(u.size()==v.size()); - std::vector w=u; - for(int j=0; j operator+(const std::vector &u, const std::vector &v) +{ + assert(u.size() == v.size()); + std::vector w = u; + for (int j = 0; j < u.size(); j++) + { + w[j] += v[j]; + } return w; } // multiplication of a vector by a scalar a*u -std::vector operator*(const double& a, const std::vector& u){ +std::vector operator*(const double &a, const std::vector &u) +{ std::vector w(u.size()); - for(int j=0; j& u, const std::vector& v){ - assert(u.size()==v.size()); - for(int j=0; j &u, const std::vector &v) +{ + assert(u.size() == v.size()); + for (int j = 0; j < u.size(); j++) + { + u[j] += v[j]; + } } /* block Jacobi preconditioner: perform forward and backward substitution using the Cholesky factorization of the local diagonal block computed by Eigen */ -std::vector prec(const Eigen::SimplicialCholesky>& P, const std::vector& u){ +std::vector prec(const Eigen::SimplicialCholesky> &P, const std::vector &u) +{ Eigen::VectorXd b(u.size()); - for (int i=0; i x(u.size()); - for (int i=0; i& b, - std::vector& x, - double tol=1e-6) { +void CG(const MapMatrix &A, + const std::vector &b, + std::vector &x, + double tol = 1e-6) +{ assert(b.size() == A.NbRow()); - x.assign(b.size(),0.); + x.assign(b.size(), 0.); int rank; MPI_Comm_rank(MPI_COMM_WORLD, &rank); // Get the rank of the process @@ -132,123 +156,205 @@ void CG(const MapMatrix& A, // get the local diagonal block of A std::vector> coefficients; - for(auto it=A.data.begin(); it!=A.data.end(); ++it){ + for (auto it = A.data.begin(); it != A.data.end(); ++it) + { int j = (it->first).first; int k = (it->first).second; - if (k >= 0 && k < n) coefficients.push_back(Eigen::Triplet(j,k,it->second)); + if (k >= 0 && k < n) + coefficients.push_back(Eigen::Triplet(j, k, it->second)); } // compute the Cholesky factorization of the diagonal block for the preconditioner - Eigen::SparseMatrix B(n,n); + Eigen::SparseMatrix B(n, n); B.setFromTriplets(coefficients.begin(), coefficients.end()); Eigen::SimplicialCholesky> P(B); - std::vector r=b, z=prec(P,r), p=z, Ap=A*p; - double np2=(p,Ap), alpha=0.,beta=0.; - double nr = sqrt((z,r)); - double epsilon = tol*nr; + std::vector r = b, z = prec(P, r), p = z, Ap = A * p; + double np2 = (p, Ap), alpha = 0., beta = 0.; + double nr = sqrt((z, r)); + double epsilon = tol * nr; - std::vector res = A*x; - res += (-1)*b; - - double rres = sqrt((res,res)); + std::vector res = A * x; + res += (-1) * b; + + double rres = sqrt((res, res)); int num_it = 0; - while(rres>1e-5) { - alpha = (nr*nr)/(np2); - x += (+alpha)*p; - r += (-alpha)*Ap; - z = prec(P,r); - nr = sqrt((z,r)); - beta = (nr*nr)/(alpha*np2); - p = z+beta*p; - Ap=A*p; - np2=(p,Ap); - - rres = sqrt((r,r)); + while (rres > 1e-5) + { + alpha = (nr * nr) / (np2); + x += (+alpha) * p; + r += (-alpha) * Ap; + z = prec(P, r); + nr = sqrt((z, r)); + beta = (nr * nr) / (alpha * np2); + p = z + beta * p; + Ap = A * p; + np2 = (p, Ap); + + rres = sqrt((r, r)); num_it++; - if(rank == 0 && !(num_it%1)) { + if (rank == 0 && !(num_it % 1)) + { std::cout << "iteration: " << num_it << "\t"; - std::cout << "residual: " << rres << "\n"; + std::cout << "residual: " << rres << "\n"; } } } // Command Line Option Processing -int find_arg_idx(int argc, char** argv, const char* option) { - for (int i = 1; i < argc; ++i) { - if (strcmp(argv[i], option) == 0) { - return i; - } +int find_arg_idx(int argc, char **argv, const char *option) +{ + for (int i = 1; i < argc; ++i) + { + if (strcmp(argv[i], option) == 0) + { + return i; } - return -1; + } + return -1; } -int find_int_arg(int argc, char** argv, const char* option, int default_value) { - int iplace = find_arg_idx(argc, argv, option); +int find_int_arg(int argc, char **argv, const char *option, int default_value) +{ + int iplace = find_arg_idx(argc, argv, option); - if (iplace >= 0 && iplace < argc - 1) { - return std::stoi(argv[iplace + 1]); - } + if (iplace >= 0 && iplace < argc - 1) + { + return std::stoi(argv[iplace + 1]); + } - return default_value; + return default_value; } -int main(int argc, char* argv[]) { +int main(int argc, char *argv[]) +{ MPI_Init(&argc, &argv); // Initialize the MPI environment - + int size; MPI_Comm_size(MPI_COMM_WORLD, &size); // Get the number of processes - + int rank; MPI_Comm_rank(MPI_COMM_WORLD, &rank); // Get the rank of the process - if (find_arg_idx(argc, argv, "-h") >= 0) { - std::cout << "-N : side length of the sparse matrix" << std::endl; - return 0; - } - - + if (rank == 0) + std::cout << "Number of procs " << size << std::endl; + if (find_arg_idx(argc, argv, "-h") >= 0) + { + std::cout << "-N : side length of the sparse matrix" << std::endl; + return 0; + } int N = find_int_arg(argc, argv, "-N", 100000); // global size - assert(N%size == 0); - int n = N/size; // number of local rows + assert(N % size == 0); + int n = N / size; // number of local rows // row-distributed matrix - MapMatrix A(n,N); + double map_time = MPI_Wtime(); + + MapMatrix A(n, N); + + int offset = n * rank; + + // NEW MATRIX + // Making sparse matrix using eigen instead of map. + std::vector tripletList; + for (int i = 0; i < n; i++) + { + int j = i; + int v_ij = 2.0; + tripletList.push_back(T(i, j, v_ij)); + // A.Assign(i, i) = 2.0; + + if (offset + i - 1 >= 0) + { + j = i - 1; + v_ij = -1; + tripletList.push_back(T(i, j, v_ij)); + // A.Assign(i, i - 1) = -1; + } + + if (offset + i + 1 < N) + { + j = i + 1; + v_ij = -1; + tripletList.push_back(T(i, j, v_ij)); + // A.Assign(i, i + 1) = -1; + } + + if (offset + i + N < N) + { + j = i + N; + v_ij = -1; + tripletList.push_back(T(i, j, v_ij)); + // A.Assign(i, i + N) = -1; + } - int offset = n*rank; + if (offset + i - N >= 0) + { + j = i - N; + v_ij = -1; + tripletList.push_back(T(i, j, v_ij)); + // A.Assign(i, i - N) = -1; + } + } + SpMat M(n, N); + M.setFromTriplets(tripletList.begin(), tripletList.end()); + + if (rank == 0) + { + std::cout << "Print matrix " << std::endl; + std::cout << M << std::endl; + } + // ORIGINAL IMPLEMENTATION // local rows of the 1D Laplacian matrix; local column indices start at -1 for rank > 0 - for (int i=0; i= 0) A.Assign(i,i - 1) = -1; - if (offset + i + 1 < N) A.Assign(i,i + 1) = -1; - if (offset + i + N < N) A.Assign(i, i + N) = -1; - if (offset + i - N >= 0) A.Assign(i, i - N) = -1; + for (int i = 0; i < n; i++) + { + A.Assign(i, i) = 2.0; + if (offset + i - 1 >= 0) + A.Assign(i, i - 1) = -1; + if (offset + i + 1 < N) + A.Assign(i, i + 1) = -1; + if (offset + i + N < N) + A.Assign(i, i + N) = -1; + if (offset + i - N >= 0) + A.Assign(i, i - N) = -1; } + // prints map + for (const auto &elem : A.data) + { + std::cout << elem.first.first << " " << elem.first.second << " " << elem.second << "\n"; + } + + MPI_Barrier(MPI_COMM_WORLD); + if (rank == 0) + std::cout << "wall time for Map: " << MPI_Wtime() - map_time << std::endl; + // initial guess - std::vector x(n,0); + std::vector x(n, 0); // right-hand side - std::vector b(n,1); + std::vector b(n, 1); MPI_Barrier(MPI_COMM_WORLD); double time = MPI_Wtime(); - CG(A,b,x); + CG(A, b, x); MPI_Barrier(MPI_COMM_WORLD); - if (rank == 0) std::cout << "wall time for CG: " << MPI_Wtime()-time << std::endl; + if (rank == 0) + std::cout << "wall time for CG: " << MPI_Wtime() - time << std::endl; - std::vector r = A*x + (-1)*b; + std::vector r = A * x + (-1) * b; - double err = Norm(r)/Norm(b); - if (rank == 0) std::cout << "|Ax-b|/|b| = " << err << std::endl; + double err = Norm(r) / Norm(b); + if (rank == 0) + std::cout << "|Ax-b|/|b| = " << err << std::endl; MPI_Finalize(); // Finalize the MPI environment