Метод Холецкого решения СЛАУ (C++, MPI)

  • В этой теме 0 ответов, 1 участник, последнее обновление 1 месяц, 3 недели назад сделано Васильев Владимир Сергеевич.
Просмотр 0 веток ответов
  • Автор
    Сообщения
    • #6797
      @admin
      StudLance.ru

      // Примеры запуска программы:
      // ./prog -g 100
      // Генерируется матрица порядка 100, печатается вся обратная матрица
      // ./prog -f test.txt
      // Программа читает матрицу из файла test.txt, печатается вся обратная матрица
      // ./prog -g 100 5 или
      // ./prog -f test.txt 5
      // Печатается блок размера 5
      // Сначала выводится время по каждому процессу на исполнение алгоритма
      // Далее обратная матрица, норма невязки и 
      // время по каждому процессу на вычисление нормы невязки
      #include <stdio.h>
      #include <math.h>
      #include <stdlib.h>
      #include <time.h>
      #include <errno.h>
      #include <string.h>
      #include <mpi.h>
      
      #define EPS 1e-16
      
      
      
      int processInput(int argc, char *argv[], int *nPtr,
                       int *wayPtr, int *blockSizePtr);
      int processFile1(char *argv[], FILE **filePtr, int *nPtr);
      void processFile(char **argv, FILE **filePtr, int n, int rank, int size,
                      int nCols, double *A, double *buff);
      
      void generateMatrix(int n, int rank, int size, int nCols, double *A);
      void setId(int n, int rank, int size, int nCols, double *Ans);
      void transposeRL(int n, int rank, int size, double *A, double *buff);
      void transposeLR(int n, int rank, int size, double *A, double *buff);
      
      int choleskyDecomposition(int n, int rank, int size, int nCols, double *A,
                                double *d, double *buff);
      
      void findInverse(int n, int rank, int size, int nCols,
                        double *A, double *d, double *Ans, double *buff);
      struct timespec diff(struct timespec start, struct timespec end);
      
      void calcResidue(int n, int rank, int size, int nCols,
                       double *A, double *Ans, double *buff, double *d, double *res);
      void printAnswer(int n, int rank, int size,
                       double *Ans, double *buff, int blockSize);
      void printTime(int rank, int size, double t);
      
      int main(int argc, char **argv)
      {
          FILE *file;
          int n;
          int blockSize = -1;
          int checker, way;
          double *A, *Ans;
          double *d;
          
          double *buff;
      //    double res;
          int rank, size;
          int nCols;
      
          double t;
          
          MPI_Init(&argc, &argv);
          MPI_Comm_rank(MPI_COMM_WORLD, &rank);
          MPI_Comm_size(MPI_COMM_WORLD, &size);
          if (rank == 0) {
              checker = processInput(argc, argv, &n, &way, &blockSize);
              if (checker != 0) {
                  MPI_Abort(MPI_COMM_WORLD, checker);
              }
              if (way == 2) {
                  checker = processFile1(argv, &file, &n);
                  if (checker != 0)
                      MPI_Abort(MPI_COMM_WORLD, checker);
              }
              if (blockSize == -1 || blockSize >= n - 1)
                  blockSize = n;
          }
          MPI_Bcast(&n, 1, MPI_INT, 0, MPI_COMM_WORLD);
          MPI_Bcast(&way, 1, MPI_INT, 0, MPI_COMM_WORLD);
          MPI_Bcast(&blockSize, 1, MPI_INT, 0, MPI_COMM_WORLD);
          if (rank < n % size)
              nCols = n/size + 1;
          else
              nCols = n/size;
          A = (double*)malloc(n*nCols*sizeof(double));
          d = (double*)malloc(n*sizeof(double));
          Ans = (double*)malloc(n*nCols*sizeof(double));
          buff = (double*)malloc(n*sizeof(double));
          if (A == NULL || d == NULL || Ans == NULL || buff == NULL) {
              if (A) free(A);
              if (d) free(d);
              if (Ans) free(Ans);
              if (buff) free(buff);
              MPI_Abort(MPI_COMM_WORLD, -1);
          }
          if (way == 1)
              generateMatrix(n, rank, size, nCols, A);
          else
              processFile(argv, &file, n, rank, size, nCols, A, buff);
          setId(n, rank, size, nCols, Ans);
          
          MPI_Barrier(MPI_COMM_WORLD);
          t = MPI_Wtime();
      
          choleskyDecomposition(n, rank, size, nCols, A, d, buff);
          transposeRL(n, rank, size, A, buff);
          findInverse(n, rank, size, nCols, A, d, Ans, buff);
          transposeLR(n, rank, size, Ans, buff);    
      
          t = MPI_Wtime() - t;
          if (rank == 0)
              printf("%lf ", t);
      
          printAnswer(n, rank, size, Ans, buff, blockSize);
      
          if (way == 1)
              generateMatrix(n, rank, size, nCols, A);
          else
              processFile(argv, &file, n, rank, size, nCols, A, buff);
          MPI_Barrier(MPI_COMM_WORLD);
      //    t = MPI_Wtime();
      //    calcResidue(n, rank, size, nCols, A, Ans, buff, d, &res);
      //   t = MPI_Wtime() - t;
      //    if (rank == 0)
      //        printf("%e\n", res);
      //    printTime(rank, size, t);
      
          free(A);
          free(d);
          free(Ans);
          free(buff);
          MPI_Finalize();
          return 0;
      }
      
      
      
      int processInput(int argc, char *argv[], int *nPtr,
                       int *wayPtr, int *blockSizePtr)
      {
          char *end;
          if (argc != 3 && argc != 4) {
              return -1; // Incorrect arguments
          }
          if (strlen(argv[1]) != 2)
              return -1;
          if (argv[1][0] != '-' || (argv[1][1] != 'g' && argv[1][1] != 'f'))
              return -1;
          if (argv[1][1] == 'g') {
              *wayPtr = 1;
              errno = 0;
              *nPtr = strtol(argv[2], &end, 10);
              if (errno != 0 || *end != '\0' || *nPtr <= 0) {
                  return -1;
              }
          } else if (argv[1][1] == 'f')
              *wayPtr = 2;
          if (argc == 4) {
              errno = 0;
              *blockSizePtr = strtol(argv[3], &end, 10);
              if (errno != 0 || *end != '\0' || *blockSizePtr < 0)
                  return -1;
          }
          return 0;
      }
      
      int processFile1(char *argv[], FILE **filePtr, int *nPtr)
      {
          int checker = 1;
          *filePtr = fopen(argv[2], "r");
          if (*filePtr == NULL)
              return -2;
          checker = fscanf(*filePtr, "%d", nPtr);
          if (checker == EOF) {
              fclose(*filePtr);
              return -3;
          } else if (checker == 0) {
              fclose(*filePtr);
              return -4;
          }
          if (*nPtr <= 0) {
              fclose(*filePtr);
              return -5;
          }
          fclose(*filePtr);
          return 0;
      }
      
      
      void processFile(char **argv, FILE **filePtr, int n, int rank, int size,
                      int nCols, double *A, double *buff)
      {
          int i, j;
          MPI_Status status;
          if (rank == 0) {
              *filePtr = fopen(argv[2], "r");
              if (*filePtr == NULL)
                  MPI_Abort(MPI_COMM_WORLD, -2);
              fscanf(*filePtr, "%d", &i);
              for (i = 0; i < n; i++) {
                  if (i % size == rank) {
                      for (j = 0; j < n; j++)
                          if (fscanf(*filePtr, "%lf", A+((i/size)*n+j)) != 1) {
                              fclose(*filePtr);
                              MPI_Abort(MPI_COMM_WORLD, -4);
                          }
                  } else {
                      for (j = 0; j < n; j++)
                          if (fscanf(*filePtr, "%lf", buff+j) != 1) {
                              fclose(*filePtr);
                              MPI_Abort(MPI_COMM_WORLD, -4);
                          }
                      MPI_Ssend(buff, n, MPI_DOUBLE, i % size, 0, MPI_COMM_WORLD);
                  }
              }
          } else {
              for (i = 0; i < nCols; i++)
                  MPI_Recv(A+(i*n), n, MPI_DOUBLE, 0, MPI_ANY_TAG,
                           MPI_COMM_WORLD, &status);
          }
          MPI_Barrier(MPI_COMM_WORLD);
          for (i = 0; i < n; i++) {
              if (rank == i % size) {
                  MPI_Bcast(A+((i/size)*n), i+1, MPI_DOUBLE, i%size, MPI_COMM_WORLD);
                  for (j = 0; j < i; j++)
                      if (j % size == rank &&
                          fabs(A[(j/size)*n+i]- A[(i/size)*n+j]) > EPS)
                          MPI_Abort(MPI_COMM_WORLD, -1);
              } else {
                  MPI_Bcast(buff, i+1, MPI_DOUBLE, i%size, MPI_COMM_WORLD);
                  for (j = 0; j < i; j++)
                      if (j % size == rank &&
                          fabs(A[(j/size)*n+i] - buff[j]) > EPS)
                          MPI_Abort(MPI_COMM_WORLD, -1);
              }
          }
         
      }
      
      void generateMatrix(int n, int rank, int size, int nCols, double *A)
      {
          int i, j;
          if ((nCols-1)*size+rank == n - 1) {
              for (i = 0; i < nCols-1; i++)
                  for (j = 0; j < n - 1; j++)
                      if (i*size+rank == j)
                          A[i*n+j] = 1;
                      else
                          A[i*n+j] = 0;
              for (i = 0; i < n; i++)
                  A[(nCols-1)*n+i] = i + 1;
          } else {
              for (i = 0; i < nCols; i++)
                  for (j = 0; j < n - 1; j++)
                      if (i*size+rank == j)
                          A[i*n+j] = 1;
                      else
                          A[i*n+j] = 0;
          }
          for (i = 0; i < nCols; i++)
              A[i*n+n-1] = i*size + rank + 1;
      }
      
      
      
      int choleskyDecomposition(int n, int rank, int size, int nCols, double *A,
                                double *d, double *buff)
      {
          double temp;
          int i, j, k;
          int start;
          for (i = 0; i < n; i++) {
              if (i % size == rank) {
                  temp = 0;
                  for (k = 0; k < i; k++)
                      temp += A[(i/size)*n+k] * A[(i/size)*n+k] * d[k];
                  d[i] = A[(i/size)*n+i] - temp;
                  A[(i/size)*n+i] = 1;
                  if (fabs(d[i]) < EPS)
                      MPI_Abort(MPI_COMM_WORLD, -8);
                  MPI_Bcast(A+(n*(i/size)), i+1, MPI_DOUBLE, i%size, MPI_COMM_WORLD);
                  MPI_Bcast(d+i, 1, MPI_DOUBLE, i%size, MPI_COMM_WORLD);
                  for (j = i/size+1; j < nCols; j++) {
                      temp = 0;
                      for (k = 0; k < i; k++)
                          temp += A[(i/size)*n+k] * d[k] * A[j*n+k];
                      A[j*n+i] = (A[j*n+i] - temp) / d[i];
                  }
              } else {
                  MPI_Bcast(buff, i+1, MPI_DOUBLE, i%size, MPI_COMM_WORLD);
                  MPI_Bcast(d+i, 1, MPI_DOUBLE, i%size, MPI_COMM_WORLD);
                  if (rank < i % size)
                      start = i / size + 1;
                  else
                      start = i / size;
                  for (j = start; j < nCols; j++) {
                      temp = 0;
                      for (k = 0; k < i; k++)
                          temp += buff[k] * d[k] * A[j*n+k];
                      A[j*n+i] = (A[j*n+i] - temp) / d[i];
                  }
              }
          }
          return 0;
      }
      
      
      
      void findInverse(int n, int rank, int size, int nCols,
                        double *A, double *d, double *Ans, double *buff)
      {
          int i, j, k;
          double temp;
          for (i = 0; i < n; i++) {
              if (i % size == rank) {
                  MPI_Bcast(A+((i/size)*n+i), n-i, MPI_DOUBLE,
                            rank, MPI_COMM_WORLD);
                  for (j = 0; j < nCols && j*size+rank <= i; j++) {
                      temp = Ans[j*n+i];
                      for (k = i + 1; k < n; k++)
                          Ans[j*n+k] -= temp * A[(i/size)*n+k];
                  }
              } else {
                  MPI_Bcast(buff+i, n-i, MPI_DOUBLE, i % size, MPI_COMM_WORLD);
                  for (j = 0; j < nCols && j*size+rank <= i; j++) {
                      temp = Ans[j*n+i];
                      for (k = i + 1; k < n; k++)
                          Ans[j*n+k] -= temp * buff[k];
                  }
              }
          }
          for (i = 0; i < nCols; i++)
              for (j = i*size+rank; j < n; j++)
                  Ans[i*n+j] /= d[j];
          for (i = n - 1; i >= 0; i--) {
              if (i % size == rank) {
                  MPI_Bcast(A+((i/size)*n), i+1, MPI_DOUBLE, rank, MPI_COMM_WORLD);
                  for (j = 0; j < nCols && j * size + rank < i; j++) {
                      temp = Ans[j*n+i];
                      for (k = j * size + rank; k < i; k++)
                          Ans[j*n+k] -= temp * A[(i/size)*n+k];
                  }
              } else {
                  MPI_Bcast(buff, i+1, MPI_DOUBLE, i % size, MPI_COMM_WORLD);
                  for (j = 0; j < nCols && j * size + rank < i; j++) {
                      temp = Ans[j*n+i];
                      for (k = j * size + rank; k < i; k++)
                          Ans[j*n+k] -= temp * buff[k];
                  }
              }
          }
      }
      
      
      void calcResidue(int n, int rank, int size, int nCols,
                       double *A, double *Ans, double *buff, double *d, double *res)
      {
          int i, j, k;
          double temp = 0, s = 0;
          for (i = 0; i < n; i++) {
              s = 0;
              if (i % size == rank) {
                  MPI_Bcast(A+((i/size)*n), n, MPI_DOUBLE, rank, MPI_COMM_WORLD);
                  for (j = 0; j < nCols; j++) {
                      temp = 0;
                      for (k = 0; k < n; k++)
                          temp += A[(i/size)*n+k] * Ans[j*n+k];
                      if (j * size + rank == i)
                          temp -= 1;
                      s += fabs(temp);
                  }
              } else {
                  MPI_Bcast(buff, n, MPI_DOUBLE, i%size, MPI_COMM_WORLD);
                  for (j = 0; j < nCols; j++) {
                      temp = 0;
                      for (k = 0; k < n; k++)
                          temp += buff[k] * Ans[j*n+k];
                      if (j * size + rank == i)
                          temp -= 1;
                      s += fabs(temp);
                  }
              }
              MPI_Allreduce(&s, d+i, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
          }
          s = 0;
          for (i = 0; i < n; i++)
              s += d[i];
          *res = s;
      }
      
      
      void printAnswer(int n, int rank, int size,
                       double *Ans, double *buff, int blockSize)
      {
          int i, j;
          MPI_Status status;
          if (blockSize == 0)
              return;
          if (rank == 0) {
              if (blockSize <= n - 2) {
                  for (i = 0; i < blockSize; i++) {
                      if (i % size == rank) {
                          for (j = 0; j < blockSize; j++)
                              printf("%lf ", Ans[(i/size)*n+j]);
                          printf("... %lf\n", Ans[(i/size)*n+n-1]);
                      } else {
                          MPI_Recv(buff, n, MPI_DOUBLE, i%size,
                                   MPI_ANY_TAG, MPI_COMM_WORLD, &status);
                          for (j = 0; j < blockSize; j++)
                              printf("%lf ", buff[j]);
                          printf("... %lf\n", buff[n-1]);
                      }
                  }
                  for (j = 0; j < blockSize + 2; j++)
                      printf("... ");
                  printf("\n");
                  if ((n-1) % size == rank) {
                      for (j = 0; j < blockSize; j++)
                          printf("%lf ", Ans[((n-1)/size)*n+j]);
                      printf("... %lf\n", Ans[((n-1)/size)*n+n-1]);
                  } else {
                      MPI_Recv(buff, n, MPI_DOUBLE, (n-1)%size,
                               MPI_ANY_TAG, MPI_COMM_WORLD, &status);
                      for (j = 0; j < blockSize; j++)
                          printf("%lf ", buff[j]);
                      printf("... %lf\n", buff[n-1]);
                      
                  }
              } else {
                  for (i = 0; i < n; i++) {
                      if (i % size == rank) {
                          for (j = 0; j < n; j++)
                              printf("%lf ", Ans[(i/size)*n+j]);
                          printf("\n");
                      } else {
                          MPI_Recv(buff, n, MPI_DOUBLE, i%size,
                                   MPI_ANY_TAG, MPI_COMM_WORLD, &status);
                          for (j = 0; j < n; j++)
                              printf("%lf ", buff[j]);
                          printf("\n");
                      }
                  }
              }
          } else {
              if (blockSize <= n - 2) {
                  for (i = 0; i < blockSize; i++)
                      if (i % size == rank)
                          MPI_Ssend(Ans+((i/size)*n), n, MPI_DOUBLE, 0, 0,
                                    MPI_COMM_WORLD);
                  if ((n-1) % size == rank)
                      MPI_Ssend(Ans+(((n-1)/size)*n), n, MPI_DOUBLE, 0, 0,
                                MPI_COMM_WORLD);
              } else {
                  for (i = 0; i < n; i++)
                      if (i % size == rank)
                          MPI_Ssend(Ans+((i/size)*n), n, MPI_DOUBLE, 0, 0,
                                    MPI_COMM_WORLD);
              }
          }
      }
      
      
      void printTime(int rank, int size, double t)
      {
          int i;
          double temp;
          MPI_Status status;
          if (rank == 0) {
              printf("%d %lf\n", rank, t);
              for (i = 1; i < size; i++) {
                  MPI_Recv(&temp, 1, MPI_DOUBLE, i, MPI_ANY_TAG,
                           MPI_COMM_WORLD, &status);
                  printf("%d %lf\n", i, temp);
              }
          } else {
              MPI_Ssend(&t, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
          }
      }
      
      void setId(int n, int rank, int size, int nCols, double *Ans)
      {
          int i, j;
          for (i = 0; i < nCols; i++)
              for (j = 0; j < n; j++)
                  if (j != i*size+rank)
                      Ans[i*n+j] = 0;
                  else
                      Ans[i*n+j] = 1;
      }
      
      void transposeRL(int n, int rank, int size, double *A, double *buff)
      {
          int i, j;
          for (i = 0; i < n; i++) {
              if (rank == i % size) {
                  MPI_Bcast(A+((i/size)*n), i+1, MPI_DOUBLE, i%size, MPI_COMM_WORLD);
                  for (j = 0; j < i; j++)
                      if (j % size == rank)
                          A[(j/size)*n+i] = A[(i/size)*n+j];
              } else {
                  MPI_Bcast(buff, i+1, MPI_DOUBLE, i%size, MPI_COMM_WORLD);
                  for (j = 0; j < i; j++)
                      if (j % size == rank)
                          A[(j/size)*n+i] = buff[j];
              }
          }
      }
      
      void transposeLR(int n, int rank, int size, double *A, double *buff)
      {
          int i, j;
          for (i = 0; i < n; i++) {
              if (rank == i % size) {
                  MPI_Bcast(A+((i/size)*n+i), n-i,
                            MPI_DOUBLE, i%size, MPI_COMM_WORLD);
                  for (j = i+1; j < n; j++)
                      if (j % size == rank)
                          A[(j/size)*n+i] = A[(i/size)*n+j];
              } else {
                  MPI_Bcast(buff+i, n-i, MPI_DOUBLE, i%size, MPI_COMM_WORLD);
                  for (j = i+1; j < n; j++)
                      if (j % size == rank)
                          A[(j/size)*n+i] = buff[j];
              }
          }
      }

      StudLance.ru

Просмотр 0 веток ответов
  • Для ответа в этой теме необходимо авторизоваться.
×