1              PROGRAM EXAMPLE

2              INCLUDE 'mpif.h'

3              PARAMETER (NDIM=3001)

4              INTEGER my_id,np,ierr, comm

5              REAL*8 A(NDIM,NDIM),C(NDIM),COND      

6              REAL*8 SUM

7              CALL MPI_Init(ierr)

8              comm = MPI_COMM_WORLD

9              CALL MPI_Comm_rank(comm , my_id, ierr)

10             CALL MPI_Comm_size(comm , np, ierr)

11             Write (*,*) 'Work ',np,' is ',my_id

12             DO III=1000,2000,500

13               N=III

14               DO J=1,N

15                 DO I=J,N

16                   a(i,j)=0.d0

17                   a(j,i)=0.d0

18                 END DO

19               END DO           

20             DO I=1+my_id,N,np

21               DO J=1,N

22                 a(i,j)=1.d0

23               END DO

24               a(I,I)=1.001d0

25             END DO       

26             SUM=0.d0

27             DO I=1,N

28               SUM=SUM+dreal(I)

29             END DO

30             DO I=1,N

31               c(I)=SUM+i*0.001d0   

32             END DO

33             CALL SolveMatrix(NDIM,N,A,C,COND,np,my_id,comm)

34             write(*,*) 'Number of condision ',cond

35             END DO

36             CALL MPI_FINALIZE(ierr)  

37             STOP

38             END

39        

40             SUBROUTINE SolveMatrix(NDIM,N,A,C,COND,np,my_id,comm)

41             INCLUDE 'mpif.h'

42             INTEGER  NDIM,N,ierr

43             INTEGER  IPVT(N)

44             INTEGER  np,my_id,J

45                 INTEGER  comm

46             REAL*8   COND

47             REAL*8   A(NDIM,N),C(N),WORK(N)

48             DOUBLE PRECISION t1,t2,t3

49             CALL DEC(NDIM,N,A,COND,IPVT,WORK,np,my_id,comm)

50             CALL SOLV(NDIM,N,A,C,IPVT,np,my_id,comm)

51             RETURN

52             END

53        

54             SUBROUTINE DEC(NDIM,N,A,COND,IPVT,WORK,np,my_id,comm)

55             INCLUDE 'mpif.h'

56             INTEGER NDIM,N

57             INTEGER IPVT(N),M,ierr

58             INTEGER np,my_id,WorkProc,SendProc

59             INTEGER status

60             INTEGER comm

61             REAL*8  A(NDIM,N),WORK(N)

62             REAL*8  COND

63             REAL*8  NWork,ANORM,T,Mold

64             IPVT(N)=1

65             IF(N.EQ.1) GO TO 80

66             NM1=N-1

67             ANORM=0.D0

68             DO 10 J=1+my_id,N,np

69               T=0.D0

70               DO 5 I=1,N

71                 T=T+ABS(A(I,J))

72       5       CONTINUE

73               IF(T.GT.ANORM) ANORM=T

74       10    CONTINUE

75             CALL MPI_ALLREDUCE(ANORM,T,1,MPI_REAL8,MPI_MAX,

76           *      comm,ierr)

77             ANORM=T

78             DO 35 K=1,NM1

79               KP1=K+1

80               M=K

81               CALL Number(K,np,WorkProc)

82               IF (my_id.eq.WorkProc) Then

83                 DO 15 I=KP1,N

84                   IF(ABS(A(I,K)).GT.ABS(A(M,K))) M=I

85       15        END DO

86               END IF

87               Mold=DReal(M)

88               CALL MPI_Bcast(Mold,1,MPI_REAL8,WorkProc,comm,ierr)

89               M=Int(Mold)

90               IPVT(K)=M

91               IF(M.NE.K) IPVT(N)=-IPVT(N)

92               IF (my_id.eq.WorkProc) Then

93                 T=A (M,K)

94                 A(M,K)=A(K,K)

95                 A(K,K)=T

96               END IF   

97               IF(T.EQ.0.D0) GO TO 37

98               IF (my_id.eq.WorkProc) Then

99                 DO 20 I=KP1,N

100              A(I,K)=-A(I,K)/T

101  20        CONTINUE

102          END IF

103  37      CALL MPI_Bcast(A(1,K),N,MPI_REAL8,WorkProc, comm,ierr)

104          IF(A(K,K).EQ.0.D0) GO TO 35

105          DO 30 J=KP1,N

106            T=0.d0

107            CALL Number(J,np,SendProc)

108            IF (my_id.eq.SendProc) Then

109              T=A(M,J)

110              A(M,J)=A(K,J)

111              A(K,J)=T

1            END IF

113            IF(T.EQ.0.D0) GO TO 30

114            IF (my_id.eq.SendProc) Then

115              DO 25 I=KP1,N

116                A(I,J)=A(I,J)+A(I,K)*T

117  25          CONTINUE

118            END IF

119  30      CONTINUE

120  35    CONTINUE

121        CALL Number(N,np,SendProc)

122        CALL MPI_Bcast(A(1,N),N,MPI_REAL8,SendProc, comm,ierr)

123        DO 50 K=1,N

124          T=0.D0

125          IF(K.EQ.1) GO TO 45

126          KM1=K-1

127          DO 40 I=1,KM1

128            T=T+A(I,K)*WORK(I)

129  40      CONTINUE

130  45      EK=  1.0D0

131          IF(T.LT.0.0) EK=-1.0D0

132          IF(A(K,K).EQ.0.0) GO TO 90

133          WORK(K)=-(EK+T)/A(K,K)

134  50    CONTINUE

135        DO 60 KB=1,NM1

136          K=N-KB

137          T=0.0D0

138          KP1=K+1

139            DO 55 I=KP1,N

140              T=T+A(I,K)*WORK(K)

141  55        CONTINUE

142            WORK(K)=T

143          M=IPVT(K)

144          IF(M.EQ.K) GO TO 60

145          T=WORK(M)

146          WORK(M)=WORK(K)

147          WORK(K)=T

148  60    CONTINUE

149        YNORM=0.0D0

150        DO 65 I=1,N

151          YNORM=YNORM+ABS(WORK(I))

152  65    CONTINUE

153        CALL SOLV(NDIM,N,A,WORK,IPVT,np,my_id,comm)

154        ZNORM=0.0D0

155        DO 70 I=1,N

156          ZNORM=ZNORM+ABS(WORK(I))

157  70    CONTINUE

158        COND=ANORM*ZNORM/YNORM

159        IF(COND.LT.1.D0) COND=1.0D0

160        RETURN

161  80    COND=1.D0

162        IF(A(1,1).NE.0.D0) RETURN

163  90    COND=1.0D+32

164        RETURN

165        END

166   

167        SUBROUTINE SOLV(NDIM,N,A,B,IPVT,np,my_id,comm)

168        INCLUDE 'mpif.h'

169        INTEGER NDIM,N

170        INTEGER IPVT(N),np,my_id,ierr

171        INTEGER comm

172        REAL*8  A(NDIM,N),B(N)

173        IF(N.EQ.1) GO TO 50

174        NM1=N-1

175        DO 20 K=1,NM1

176          KP1=K+1

177          M=IPVT(K)

178          T=B(M)

179          B(M)=B(K)

180          B(K)=T

181          DO 10 I=KP1,N

182            B(I)=B(I)+A(I,K)*T

183  10      CONTINUE

184  20    CONTINUE

185        DO 40 KB=1,NM1

186          KM1=N-KB

187          K=KM1+1

188          B(K)=B(K)/A(K,K)

189          T=-B(K)

190          DO 30 I=1,KM1

191            B(I)=B(I)+A(I,K)*T

192  30      CONTINUE

193  40    CONTINUE

194  50    B(1)=B(1)/A(1,1)

195        RETURN

196        END

197   

198        SUBROUTINE Number(p,s,aa)

199        INTEGER n, k, i, Works, p, s

200        n=-1

201        k=0

202        DO i=1,p+1

203          n=n+1

204          k=k+1

205          if (n.eq.s) then

206            n=0

207          end if

208          if (k1.eq.p) then

209            Works=n

210          end if

211        END DO

210        RETURN

213        END

Рис. 58. Листинг программы, реализующей параллельный алгоритм решения СЛАУ методом Гаусса