SHARE
TWEET

Untitled

a guest Feb 22nd, 2019 65 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  4. #include "Lab3IO.h"
  5.  
  6. #include "timer.h"
  7.  
  8. int main(int argc, char *argv[])
  9. {
  10.     int i, j, k, size;
  11.     double** Au;
  12.     double* X;
  13.     double temp;
  14.     int* index;
  15.     //FILE* fp;
  16.     double start, end;
  17.     int thread_count=0;
  18.  
  19.     int largest;
  20.     int largest_index;
  21.     /*Load the datasize and verify it*/
  22.     Lab3LoadInput(&Au, &size);
  23.  
  24.     thread_count = atoi(argv[1]);
  25.  
  26.     /*Calculate the solution by parallel code*/
  27.     X = CreateVec(size);
  28.     index = malloc(size * sizeof(int));
  29.  
  30.     GET_TIME(start);
  31.     // # pragma omp parallel for num_threads(thread_count)
  32.     for (i = 0; i < size; ++i)
  33.         index[i] = i;
  34.  
  35.     if (size == 1)
  36.         X[0] = Au[0][1] / Au[0][0];
  37.     else
  38.     {
  39.         /*Gaussian elimination*/
  40.         for (k = 0; k < size - 1; ++k){
  41.             /*Pivoting*/
  42.             temp = 0;
  43.             j = 0;
  44.             #pragma omp parallel num_threads(thread_count)
  45.             {
  46.                 #pragma omp for
  47.                 for (i = k; i < size; ++i)
  48.                 {
  49.                     #pragma omp critical
  50.                     if (temp < Au[index[i]][k] * Au[index[i]][k])
  51.                     {
  52.                         temp = Au[index[i]][k] * Au[index[i]][k];
  53.                         j = i;
  54.                     }  
  55.                 }
  56.                 #pragma omp single
  57.                 if (j != k)/*swap*/
  58.                 {
  59.                     i = index[j];
  60.                     index[j] = index[k];
  61.                     index[k] = i;
  62.                 }
  63.  
  64.                 /*calculating*/
  65.                 #pragma omp for private(temp, j)
  66.                 for (i = k + 1; i < size; ++i)
  67.                 {
  68.                     temp = Au[index[i]][k] / Au[index[k]][k];
  69.                     for (j = k; j < size + 1; ++j)
  70.                         Au[index[i]][j] -= Au[index[k]][j] * temp;
  71.                 }  
  72.             }    
  73.         }
  74.  
  75.         /*Jordan elimination*/
  76.         for (k = size - 1; k > 0; --k)
  77.         {
  78.             #pragma omp parallel for num_threads(thread_count) private(temp)
  79.             for (i = k - 1; i >= 0; --i )
  80.             {
  81.                 temp = Au[index[i]][k] / Au[index[k]][k];
  82.                 Au[index[i]][k] -= temp * Au[index[k]][k];
  83.                 Au[index[i]][size] -= temp * Au[index[k]][size];
  84.             }
  85.         }
  86.         /*solution*/
  87.         #pragma omp parallel for num_threads(thread_count)
  88.         for (k=0; k< size; ++k)
  89.         {
  90.             X[k] = Au[index[k]][size] / Au[index[k]][k];
  91.         }
  92.     }    
  93.  
  94.  
  95.     GET_TIME(end);
  96.  
  97.     Lab3SaveOutput(X,size,end-start);
  98.     DestroyVec(X);
  99.     DestroyMat(Au, size);
  100.     free(index);
  101.     return 0;
  102. }
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top