home *** CD-ROM | disk | FTP | other *** search
/ DP Tool Club 26 / CD_ASCQ_26_1295.iso / vrac / laipe.zip / FBANDED.FOR < prev    next >
Text File  |  1995-07-31  |  7KB  |  258 lines
  1. C
  2. C
  3. C  Program to demo parallel performance of the following LAIPE subroutines
  4. C  (1) Decompose_CSP_S that decomposes a constantly-banded, symmetric, 
  5. C      positive definitive matrix into the produce of [L][L]^T.
  6. C  (2) Substitute_CSP_S that implements forward and backward substitutions.
  7. C
  8. C  If a multiple-processor computer is available, these two subroutines
  9. C  will be implemented in a loop. The loop bound begins with one to the
  10. C  number of system processors, and each loop sets the loop counter
  11. C  as the number of employed processors. Such that, each loop executes on
  12. C  different number of processors. Keep computer STANDING ALONE while
  13. C  running this program, and elapsed time with respect to different number
  14. C  of processors will be reported. And, the speedup can be seen.
  15. C
  16. C  Running this program, user has to provide the matrix order, i.e., 10000,
  17. C  and the lower bandwidth, i.e., 600. Larger problem may have a better
  18. C  performance. It is inappropriate to try this program for a small problem
  19. C  where no sufficient operations may be distributed onto multiple
  20. C  processors. For a small bandwidth problem, try a multiple-entry solver.
  21. C
  22. C  define global variable
  23. C
  24.       PARAMETER (LIMIT=6100000)
  25.       REAL*4 A(LIMIT)
  26. C
  27. C  variables for number of processors
  28. C  where Processors : number of employed processors
  29. C        CPUs: number of system processors
  30. C
  31.       INTEGER*4 Processors,CPUs
  32. C
  33. C  working variables
  34. C
  35.       INTEGER*4 I,JJ,KK
  36. C
  37. C  variables for calculating relative error
  38. C  where Lower : lower bound of relative error
  39. C        Upper : upper bound of relative error
  40. C        R4TEMP : working variable
  41. C
  42.       REAL*4 Lower,Upper,R4TEMP
  43. C
  44. C  variable for collecting elapsed time
  45. C
  46.       REAL*4 Second
  47. C
  48. C  variable for problem size
  49. C  where N : matrix order
  50. C        M : lower bandwidth
  51. C
  52.       INTEGER*4 N,M
  53. C
  54. C  variables for memory distribution
  55. C  where IndexA : index to the left side matrix [A]
  56. C        IndexX : index to the right side vector {X} and
  57. C                 the solution computed
  58. C        IndexSolution : index to the exact solution for checking accuracy
  59. C
  60.       INTEGER*4 IndexA,IndexX,IndexSolution
  61. C
  62. C  execution flag
  63. C
  64.       LOGICAL*4 NoGood
  65. C
  66. C  numerical zero
  67. C
  68.       REAL*4 ZERO
  69.       DATA ZERO/0.0001/
  70. C
  71. C  enter order of the system
  72. C
  73.       WRITE(*,'('' Enter order (i.e. 10000 or something else): '',$)')
  74.       READ(*,*) N
  75. C
  76. C  enter half bandwidth
  77. C
  78.       WRITE(*,
  79.      &      '('' Enter half bandwidth (i.e. 600 or something else): '',
  80.      &      $)')
  81.       READ(*,*) M
  82. C
  83. C  memory distribution
  84. C
  85.       IndexA=1
  86.       IndexX=IndexA+(N-1)*M+N
  87.       IndexSolution=IndexX+N
  88. C
  89. C  check memory space
  90. C
  91.       IF(IndexSolution+N.GT.LIMIT) THEN
  92.          WRITE(*,'('' Memory overflow'')')
  93.          STOP
  94.       END IF
  95. C
  96. C  get number of CPUs on board
  97. C
  98.       CALL GetCPUs(CPUs)
  99. C
  100. C  set number of processor to execute the following statements
  101. C
  102.       DO Processors=1,CPUs
  103.          CALL SetEmployedProcessors(Processors)
  104. C
  105. C  generate left side matrix and right side vector
  106. C
  107.          Write(*,*)
  108.          Write(*,'('' Generate data in a pseudo-random procedure...'',
  109.      &             $)')
  110.          CALL GenerateData(A(IndexA),A(IndexX),A(IndexSolution),N,M)
  111. C
  112. C  output number of employed processors
  113. C
  114.          WRITE(*,'(/,'' Number of processors employed: '',I2)')
  115.      &                                               Processors
  116. C
  117. C  start collecting timing result spent in decomposition
  118. C
  119.          CALL CollectElapsedTime
  120. C      
  121. C  decompose in parallel
  122. C
  123.          CALL Decompose_CSP_S
  124.      &                       (
  125.      &                        A(IndexA),
  126.      &                        N,
  127.      &                        M,
  128.      &                        Zero,
  129.      &                        NoGood
  130.      &                       )
  131. C
  132. C  output elapsed time
  133. C
  134.          CALL GetElapsedTime(Second)
  135.          WRITE(*,'('' Decompose>>> Elapsed Time (Seconds): '',
  136.      &             F8.2)') Second
  137. C
  138. C  stop if the system is not positive definite
  139. C
  140.          IF(NoGood) THEN
  141.             WRITE(*,'('' The system is not positive definite.'')')
  142.             STOP
  143.          END IF
  144. C
  145. C  start collecting timing result spent in substitutions
  146. C
  147.          CALL CollectElapsedTime
  148. C
  149. C  substitute in parallel
  150. C
  151.          CALL Substitute_CSP_S(A(IndexA),N,M,A(IndexX))
  152. C
  153. C  output elapsed time
  154. C
  155.          CALL GetElapsedTime(Second)
  156.          WRITE(*,'('' Substitute>>> Elapsed Time (Seconds): '',
  157.      &                                           F8.2)') Second
  158. C
  159. C  calculate relative error
  160. C
  161.          KK=IndexX
  162.          JJ=IndexSolution
  163.          R4TEMP=(A(KK)-A(JJ))/A(JJ)
  164.          Lower=R4TEMP
  165.          Upper=R4TEMP
  166.          DO I=2,N
  167.             JJ=JJ+1
  168.             KK=KK+1
  169.             R4TEMP=(A(KK)-A(JJ))/A(JJ)
  170.             IF(R4TEMP.LT.Lower) THEN
  171.                Lower=R4TEMP
  172.             ELSE IF(R4TEMP.GT.Upper) THEN
  173.                Upper=R4TEMP
  174.             END IF
  175.          END DO
  176. C
  177. C  output lower bound and upper bound of the solution
  178. C
  179.          WRITE(*,*) 'Lower and upper bound of the relative error:'
  180.          WRITE(*,*) Lower,Upper
  181.       END DO
  182.       STOP
  183.       END
  184.       SUBROUTINE GenerateData(A,X,Solution,N,M)
  185. C
  186. C
  187. C  routine to generate the left side matrix, the solution,
  188. C  and the right side vector
  189. C  (A)FORTRAN CALL: CALL GenerateData(A,X,Solution,N,M)
  190. C     1.A: <R4> left side matrix, dimension((N-1)*M+N)
  191. C     2.X: <R4> right side vector, dimension(N)
  192. C     3.Solution: <R4> return the solution for checking accuracy,
  193. C                      dimension(N)
  194. C     4.N: <I4> order
  195. C     5.M: <I4> half band width
  196. C
  197.       INTEGER*4 N,M,SEED
  198.       REAL*4 A(M,1),X(1),Solution(1)
  199. C
  200. C  left side matrix
  201. C
  202.       SEED=456789
  203.       DO I=1,(N-1)*M+N
  204.          A(I,1)=RAN(SEED)*100.0
  205.       END DO
  206.       DO I=1,N
  207.          DO J=MAX0(1,I-M),I-1
  208.             A(I,I)=A(I,I)+A(I,J)
  209.          END DO
  210.          DO J=I+1,MIN0(I+M,N)
  211.             A(I,I)=A(I,I)+A(J,I)
  212.          END DO
  213.          A(I,I)=A(I,I)*100.0
  214.       END DO
  215. C
  216. C  solution
  217. C
  218.       DO I=1,N
  219.          Solution(I)=Ran(Seed)*500.0
  220.       END DO
  221. C
  222. C  right side vector
  223. C
  224.       DO I=1,N
  225.          X(I)=A(I,I)*Solution(I)
  226.          DO J=MAX0(1,I-M),I-1
  227.             X(I)=X(I)+A(I,J)*Solution(J)
  228.          END DO
  229.          DO J=I+1,MIN0(I+M,N)
  230.             X(I)=X(I)+A(J,I)*Solution(J)
  231.          END DO
  232.       END DO
  233. C
  234.       RETURN
  235.       END
  236.       REAL*4 FUNCTION RAN(SEED)
  237. C
  238. C
  239. C  function to generate a pseudo-random number
  240. C  (A)FORTRAN ENTRY: RAN(SEED)
  241. C     1.RAN: <R4> return a random number
  242. C     2.SEED: <I4> a seed (updated)
  243. C
  244.       INTEGER*4 SEED,A
  245.       REAL*8 RTEMPA,RTEMPB
  246.       DATA M/2147483647/,A/16807/
  247. C
  248. C  generate a random number
  249. C
  250.       RTEMPA=DFLOAT(SEED)*A
  251.       ITEMP=INT(RTEMPA/DFLOAT(M))
  252.       RTEMPB=DFLOAT(M)*DFLOAT(ITEMP)
  253.       SEED=INT(RTEMPA-RTEMPB)
  254.       RAN=DFLOAT(SEED)/DFLOAT(M)
  255. C
  256.       RETURN
  257.       END
  258.