cholersky.c
上传用户:hhyinxing
上传日期:2013-09-10
资源大小:833k
文件大小:3k
源码类别:

并行计算

开发平台:

Unix_Linux

  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  4. #include "mpi.h"
  5. #define  MAX_PROCESSOR_NUM  12   /*set the max number of the Processor*/
  6. #define  MAX_ARRAY_SIZE     32   /*set the max size of the array*/
  7. int main(int argc, char *argv[])
  8. {
  9.     double a[MAX_ARRAY_SIZE][MAX_ARRAY_SIZE], g[MAX_ARRAY_SIZE][MAX_ARRAY_SIZE];
  10.     int n;
  11.     double transTime = 0,tempCurrentTime, beginTime;
  12.     MPI_Status status;
  13.     int rank, size;
  14.     FILE *fin;
  15.     int i, j, k;
  16.     MPI_Init(&argc, &argv);
  17.     MPI_Comm_rank(MPI_COMM_WORLD,&rank);
  18.     MPI_Comm_size(MPI_COMM_WORLD,&size);
  19.     if(rank == 0)
  20.     {
  21. fin = fopen("dataIn.txt","r");
  22.         if (fin == NULL)     /* have no input source file */
  23.         {
  24.             puts("Not find input data file");     
  25.         puts("Please create a file "dataIn.txt"");
  26.             puts("<example for dataIn.txt> ");
  27.             puts("4");
  28.             puts("1  3  4  5");
  29.             puts("3  5  6  1");
  30.             puts("4  6  2  7");
  31.             puts("5  1  7  8");
  32.             puts("nArter's default data are running for youn");
  33.         }
  34.         else
  35.         {
  36.             fscanf(fin, "%d", &n);
  37.     if ((n < 1)||(n > MAX_ARRAY_SIZE)) /* the matrix input error */
  38.             {
  39.                 puts("Input the Matrix's size is out of range!");
  40.                 exit(-1);
  41.             }
  42.             for(i = 0; i < n; i ++) /* get the data of the matric in */ 
  43.             {
  44.                 for(j = 0; j < n; j ++) fscanf(fin,"%lf", &a[i][j]);
  45.             }
  46.         }
  47.         /* put out the matrix */
  48.         puts("Cholersky Decomposion");
  49.         puts("Input Matrix A from dataIn.txt");
  50.         for(i = 0; i < n; i ++)
  51.         {
  52.             for(j = 0; j < n; j ++) printf("%9.5f  ", a[i][j]);
  53.             printf("n");
  54.         }
  55.         printf("n");
  56.     }
  57.     MPI_Bcast(&n, 1, MPI_INT, 0, MPI_COMM_WORLD);
  58.     for(k = 0; k < n; k ++)
  59.     {
  60.         /* gathering the result ,and then broacasting to each processor */
  61.         MPI_Bcast(a[k], (n-k)*MAX_ARRAY_SIZE, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  62.                                                                                                             
  63.         for(i = k+rank; i < n; i += size)
  64.         {
  65.             for(j = 0; j < k; j ++)
  66.             {
  67.                 g[i][j] = a[i][j];
  68.             }
  69.             if (i == k)
  70.             {
  71.                 for(j = k; j < n; j ++) g[i][j] = a[i][j]/sqrt(a[k][k]);
  72.             }
  73.             else
  74.             {
  75.                 g[i][k] = a[i][k]/sqrt(a[k][k]);
  76.                 for(j = k+1; j < n; j ++) g[i][j] = a[i][j] - a[i][k]*a[k][j]/a[k][k];
  77.             }
  78.         }
  79.         /* use the Cholersky Algorithm */
  80.         for(i = k +rank; i < n; i ++)
  81.         {
  82.             MPI_Send(g[i], n, MPI_DOUBLE, 0, k*1000+i, MPI_COMM_WORLD);
  83.         }
  84.         if(rank == 0)
  85.         {
  86.             for(j = 0; j < size; j ++)
  87.             {
  88.                 for(i = k + j; i < n; i += size)
  89.                 {
  90.                     MPI_Recv(a[i], n, MPI_DOUBLE, j, k*1000+i, MPI_COMM_WORLD, &status);
  91.                 }
  92.             }
  93.         }
  94.     }
  95.     if (rank == 0)
  96.     {
  97.         puts("After Cholersky Discomposion");
  98.         puts("Output Matrix G");
  99.         for(i =0; i < n; i ++)
  100.         {
  101.             for(j = 0; j < i; j ++) printf("           ");
  102.             for(j = i; j < n; j ++) printf("%9.5f  ", a[i][j]);
  103.             printf("n");
  104.         } /* output the result */
  105.     }
  106.     MPI_Finalize();/* end of the program */
  107. }