home *** CD-ROM | disk | FTP | other *** search
/ DP Tool Club 26 / CD_ASCQ_26_1295.iso / vrac / laipe.zip / CBANDED.C next >
C/C++ Source or Header  |  1995-07-31  |  7KB  |  266 lines
  1. /*
  2. ;
  3. ;  Program to demo parallel performance of the following two LAIPE functions
  4. ;  (1) Decompose_CSP_S that decomposes a constantly-banded, symmetric,
  5. ;      and positive definite matrix into the product of [L][L]^T.
  6. ;  (2) Substitute_CSP_S that implements forward and backward substitutions.
  7. ;
  8. ;  If a multiple-processor computer is available, these two functions will
  9. ;  be implemented in a loop. The loop bound begins with one to the number
  10. ;  of system processors, and each loop sets the loop counter as the number
  11. ;  of employed processors. Such that, each loop executes on different
  12. ;  number of processors. Keep computer STANDING ALONE while running this
  13. ;  program, and elapsed time with respect to different number of processors
  14. ;  will be reported.
  15. ;
  16. ;  Running this program, user has to provide the matrix order, i.e., 10000,
  17. ;  and the lower bandwidth, i.e., 600. Large problem may have a better
  18. ;  performance. It is inappropriate to try this program for a small problem
  19. ;  where no sufficient operations may be distributed onto multiple
  20. ;  processors. For a small bandwidth problem, try a multiple-entry solver.
  21. ;
  22. */
  23.  
  24. /*  initial values  */
  25. #define LIMIT 6100000
  26.  
  27. /*  head file  */
  28. #include <stdlib.h>
  29. #include <stdio.h>
  30.  
  31. /*  function protype  */
  32. float crandom( int *seed);
  33. void GenerateData( float *A,
  34.                    float *X,
  35.                    float *Solution,
  36.                    int *N,
  37.                    int *M
  38.                  );
  39.  
  40.  
  41. /*  main function  */
  42. main()
  43. {
  44.    static float A[LIMIT];
  45.    float Zero, Lower, Upper, R4TEMP;
  46.    int N, M, IndexA, IndexX, IndexSolution;
  47.    int processors, CPUs, KK, JJ;
  48.    int NoGood, I;
  49.    float second;
  50.    Zero = 0.0001;
  51.  
  52. /*  enter order of the system  */
  53.  
  54.    printf("Enter order (i.e., 10000 or something else): ");
  55.    scanf("%d", &N);
  56.  
  57. /*  enter half bandwidth  */
  58.  
  59.    printf("Enter half bandwidth (i.e., 600 or something else): ");
  60.    scanf("%d", &M);
  61.  
  62. /*  memory distribution  */
  63.  
  64.    IndexA = 0;
  65.    IndexX = IndexA + (N - 1) * M + N;
  66.    IndexSolution = IndexX + N;
  67.  
  68. /*  check memory space  */
  69.  
  70.    if( IndexSolution + N >= LIMIT )
  71.    {
  72.        printf("Memory overflow\n");
  73.        exit;
  74.    }
  75.  
  76. /*  get number of system processors  */
  77.  
  78.    GetCPUs( &CPUs );
  79.  
  80. /*  set number of processor to execute the following statements  */
  81.  
  82.    for ( processors = 1; processors <= CPUs; processors ++ )
  83.    {
  84.       SetEmployedProcessors( &processors );
  85.  
  86. /*
  87. ;  generate left side matrix and right side vector
  88. */
  89.       printf("\n");
  90.       printf("Generate data in a pseudo-random procedure...");
  91.       GenerateData( &A[IndexA],
  92.                     &A[IndexX],
  93.                     &A[IndexSolution],
  94.                     &N,
  95.                     &M);
  96. /*  output number of employed processors  */
  97.  
  98.       printf("\nNumber of processors employed: %d\n", processors);
  99.  
  100. /*  start collecting elapsed time  */
  101.  
  102.       CollectElapsedTime();
  103.       
  104. /*  decompose in parallel  */
  105.  
  106.       Decompose_CSP_S( &A[IndexA],
  107.                        &N,
  108.                        &M,
  109.                        &Zero,
  110.                        &NoGood
  111.                      );
  112.  
  113.  
  114. /*  output elapsed time  */
  115.  
  116.       GetElapsedTime(&second);
  117.       printf("Decompose>>> Elapsed Time (seconds): %8.2f\n", second);
  118.  
  119. /*  stop if the system is not positive definite  */
  120.  
  121.       if( NoGood ==  1)
  122.       {
  123.          printf("The system is not positive definite.\n");
  124.          exit;
  125.       }
  126.  
  127. /*  start collecting elapsed time  */
  128.  
  129.       CollectElapsedTime();
  130.  
  131. /*  substitute in parallel  */
  132.  
  133.       Substitute_CSP_S( &A[IndexA],
  134.                         &N,
  135.                         &M,
  136.                         &A[IndexX]
  137.                       );
  138.  
  139. /*  output elapsed time  */
  140.  
  141.       GetElapsedTime(&second);
  142.       printf("Substitute>>> Elapsed Time (seconds): %8.2f\n", second);
  143.  
  144. /*  calculate upper and lower bounds of relative error  */
  145.  
  146.       KK = IndexX;
  147.       JJ = IndexSolution;
  148.       R4TEMP = ( A[KK] - A[JJ] ) / A[JJ];
  149.       Lower = R4TEMP;
  150.       Upper = R4TEMP;
  151.       for ( I = 2; I <= N; I ++)
  152.       {
  153.          JJ ++;
  154.          KK ++;
  155.          R4TEMP = ( A[KK]-A[JJ] ) / A[JJ];
  156.          if( R4TEMP < Lower)
  157.          {
  158.             Lower = R4TEMP;
  159.          }
  160.          else if (R4TEMP > Upper )
  161.          {
  162.             Upper = R4TEMP;
  163.          }
  164.       }
  165.  
  166. /*  output lower bound and upper bound of the relative error  */
  167.  
  168.       printf("Lower and upper bound of the relative error:\n");
  169.       printf("%6.12f    %6.12f\n", Lower, Upper);
  170.    }
  171. }
  172.  
  173. void GenerateData( float *A,
  174.                    float *X,
  175.                    float *Solution,
  176.                    int *N,
  177.                    int *M
  178.                  )
  179. /*
  180. ;
  181. ;  routine to generate the left side matrix, and the right side vector
  182. ;  Parameters:
  183. ;     1.A: <R4> left side matrix, dimension((N-1)*M+N)
  184. ;     2.X: <R4> right side vector, dimension(N)
  185. ;     3.Solution: <R4> return the solution for checking accuracy,
  186. ;                      dimension(N)
  187. ;     4.N: <I4> order
  188. ;     5.M: <I4> half band width
  189. ;
  190. */
  191. {
  192.    int SEED, K, J, I, II, IJ;
  193.    SEED = 456789;
  194.  
  195. /* left side matrix  */
  196.  
  197.    for ( I = 0; I < (( *N - 1 ) * *M + *N ); I ++ )
  198.        A[I] = crandom( &SEED ) * 100.0;
  199.    for ( I = 1; I <= *N; I ++)
  200.    {
  201.       II = (I - 1) * ( *M + 1);
  202.       K = I - *M;
  203.       if( K < 1 ) K = 1;
  204.       for ( J = K; J < I; J ++ )
  205.       {
  206.          IJ = ( J - 1 ) * *M + I - 1;
  207.          A[II] = A[II] + A[IJ];
  208.       }
  209.       K = I + *M;
  210.       if( K > *N ) K = *N;
  211.       IJ = (I - 1) * *M + I;
  212.       for ( J = I + 1; J <= K; J ++ )
  213.       {
  214.          A[II] = A[II] + A[IJ];
  215.          IJ ++;
  216.       }
  217.       A[II] = A[II] * 100.0;
  218.    }
  219.  
  220. /*  solution  */
  221.    
  222.    for ( I = 0; I < *N; I ++ )
  223.    {
  224.       Solution[I] = crandom( &SEED ) * 500.0;
  225.    }
  226.  
  227. /*  right side vector  */
  228.  
  229.    II = 0;
  230.    for ( I = 1; I <= *N; I ++ )
  231.    {
  232.       II = (I - 1) * ( *M + 1);
  233.       X[I-1] = A[II] * Solution[I-1];
  234.       K = I - *M;
  235.       if( K < 1 ) K = 1;
  236.       for ( J = K; J < I; J ++)
  237.       {
  238.           IJ = ( J - 1) * *M + I - 1;
  239.           X[I-1] = X[I-1] + A[IJ] * Solution[J-1];
  240.       }
  241.       K = I + *M;
  242.       if( K > *N ) K = *N;
  243.       IJ= (I - 1) * *M + I;
  244.       for ( J = I + 1; J <= K; J ++)
  245.       {
  246.           X[I-1] = X[I-1] + A[IJ] * Solution[J-1];
  247.           IJ ++;
  248.       }
  249.    }
  250.    return;
  251. }
  252. float crandom(int *seed)
  253. {
  254. #define m 2147483647
  255. #define a 16807
  256.    double rtempa, rtempb;
  257.    float rtemp;
  258.    int itemp;
  259.    rtempa = (double) *seed * a;
  260.    itemp = (int) (rtempa / (double) m);
  261.    rtempb = (double) m * (double) itemp;
  262.    *seed = (int) (rtempa - rtempb);
  263.    rtemp = (double) *seed / (double) m;
  264.    return rtemp;
  265. }
  266.