Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- // solve in parallel Ax = b, there A is a coefficientsMatrix and b is a freeMembersVector
- std::vector<double> solveParallelSeidel(matrix &a,
- std::vector<double> &b,
- double eps,
- int approximationsCount) {
- int size;
- int rank;
- MPI_Comm_size(MPI_COMM_WORLD, &size);
- MPI_Comm_rank(MPI_COMM_WORLD, &rank);
- MPI_Status status;
- int root = 0;
- int tag = 0;
- int procNumber = 0;
- int indexInProc = 0;
- int elementsCounter = 0;
- int n = static_cast<int>(b.size());
- int elementsPerBlock = n / size;
- int remainElements = n % size + elementsPerBlock;
- int elementsCount = (rank == size - 1 ? remainElements : elementsPerBlock);
- double currSum = 0.;
- double globalSum = 0.;
- double oldNorm = 0.;
- double newNorm = 0.;
- double currEps = 0.;
- char endOfProgram = 0;
- std::vector<double> resultX;
- std::vector<double> x(elementsCount);
- matrix currA;
- // Init data on all Proc
- if (rank == root) {
- resultX.resize(n);
- for (int i = 1; i < size; i++) {
- elementsCounter = (i == size - 1 ? remainElements : elementsPerBlock);
- MPI_Send(a[i * elementsPerBlock], elementsCounter * n, MPI_DOUBLE, i, tag, MPI_COMM_WORLD);
- }
- } else {
- currA.init(n, elementsCount);
- MPI_Recv (currA[0], elementsCount * n, MPI_DOUBLE, root, tag, MPI_COMM_WORLD, &status);
- }
- //main process
- do {
- newNorm = 0.;
- for (int i = 0; i < n; i++) {
- currSum = 0;
- procNumber = i / elementsPerBlock;
- if (procNumber == size) procNumber--;
- indexInProc = i - procNumber * elementsPerBlock;
- for (int j = 0; j < elementsCount; j++) {
- currSum -= currA[i][j] * x[j];
- }
- //MPI_Reduce(&currSum, &globalSum, 1, MPI_DOUBLE, MPI_SUM, root, MPI_COMM_WORLD);
- if (rank == root) {
- newNorm += globalSum * globalSum;
- resultX[i] = (globalSum + b[i]) / a[i][i];
- globalSum = 0.;
- if (procNumber != root) {
- MPI_Send (&resultX[i], 1, MPI_DOUBLE, procNumber, tag, MPI_COMM_WORLD);
- } else {
- x[indexInProc] = resultX[i];
- }
- } else if (rank == procNumber) {
- MPI_Recv (&x[indexInProc], 1, MPI_DOUBLE, root, tag, MPI_COMM_WORLD, &status);
- }
- }
- if (rank == root) {
- std::cerr << std::endl;
- currEps = abs(sqrt(newNorm) - sqrt(oldNorm));
- approximationsCount--;
- oldNorm = newNorm;
- if (currEps < eps || !approximationsCount) endOfProgram = 1;
- MPI_Bcast (&endOfProgram, 1, MPI_BYTE, root, MPI_COMM_WORLD);
- }
- } while (!endOfProgram);
- return resultX;
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement