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