White R.E. Computational Mathematics: Models, Methods, and Analysis with MATLAB and MPI

Подождите немного. Документ загружается.

372 CHAPTER 9. KRYLOV METHODS FOR AX = D

85. mag(k) = rho;

86. end

87. % End of gmres loop.

88. % h(1:k,1:k) is upper triangular matrix in QR.

89. y = h(1:k,1:k)\g(1:k);

90. % Form linear combination.

91. for i=1:k

92. x(2:nx,2:ny) = x(2:nx,2:ny) + v(2:nx,2:ny,i)*y(i);

93. end

94. k

95. semilogy(mag)

96. x((nx+1)/2,(nx+1)/2)

97. % mesh(x)

98. % eig(h(1:k,1:k))

With the SSOR preconditioner convergence of the above code is attained in

25 iterations, and 127 iterations are required with no preconditioner. Larger

numbers of iterations require more storage for the increasing number of basis

vectors. One alternative is to restart the iteration and to use the last iterate as

an initial guess for the restarted GMRES. This is examined in the next section.

9.5.1 Exercises

1. Experiment with the parameters nx, errtol and w in the code pcgmres.m.

2. Experiment with the parameters a1, a2 and a3 in the code pcgmres.m.

3. Verify the calculations with and without the SSOR preconditioner. Com-

pare the SSOR preconditioner with others such as block diagonal or ADI pre-

conditioning.

9.6 GMRES(m) and MPI

In order to avoid storage of the basis vectors that are constructed in the GMRES

method, after doing a number of iterates one can restart the GMRES iteration

using the last GMRES iterate as the initial iterate of the new GMRES iteration.

GMRES(m) Method.

let

{

be an initial guess for the solution

for i = 1, imax

for k = 1, m

ﬁnd

{

via GMRES

test for convergence

endloop

{

= {

endloop.

9.6. GMRES(M) AND MPI 373

The following is a partial listing of an MPI implementation of GMRES(m).

It solves the same partial di

erential equation as in the previous section where

the M

ATLA B code pcgmres.m used GMRES. Lines 1-66 are the initialization of

the code. The outer loop of GMRES(m) is executed in the while loop in lines

66-256. The inner loop is expected in lines 135-230, and here the restart m is

given by kmax. The new initial guess is deﬁned in lines 112-114 where the new

initial residual is computed. The GMRES implementation is similar to that used

in the M

AT LAB code pcgmres.m. The additive Schwarz SSOR preconditioner

is also used, but here it changes with the number of processors. Concurrent

calculations used to do the matrix products, dot products and vector updates

are similar to the MPI code cgssormpi.f.

MPI/Fortran Code gmresmmpi.f

1. program gmres

2.! This code approximates the solution of

3.! -u_xx - u_yy + a1 u_x + a2 u_y + a3 u = f

4.! GMRES(m) is used with a SSOR verson of the

5.! Schwarz additive preconditioner.

6.! The sparse matrix product, dot products and updates

7.! are also done in parallel.

8. implicit none

9. include ’mpif.h’

10. real, dimension(0:1025,0:1025,1:51):: v

11. real, dimension(0:1025,0:1025):: r,b,x,rhat

12. real, dimension(0:1025):: xx,yy

13. real, dimension(1:51,1:51):: h

14. real, dimension(1:51):: g,c,s,y,mag

15. real:: errtol,rho,hik,hipk,nu,gk,gkp,w,t0,timef,tend

16. real :: loc_rho,loc_ap,loc_error,temp

17. real :: hh,a1,a2,a3,ac,ae,aw,an,as,rac

18. integer :: nx,ny,n,kmax,k,i,j,mmax,m,sbn

19. integer :: my_rank,proc,source,dest,tag,ierr,loc_n

20. integer :: status(mpi_status_size),bn,en

Lines 21-56 initialize arrays and are not listed

57. call mpi_init(ierr)

58. call mpi_comm_rank(mpi_comm_world,my_rank,ierr)

59. call mpi_comm_size(mpi_comm_world,proc,ierr)

60. loc_n = (n-1)/proc

61. bn = 1+(my_rank)*loc_n

62. en = bn + loc_n -1

63. call mpi_barrier(mpi_comm_world,ierr)

64. if (my_rank.eq.0) then

65. t0 = timef()

66. end if

67.! Begin restart loop.

374 CHAPTER 9. KRYLOV METHODS FOR AX = D

68. do while ((rho

Aerrtol).and.(m?mmax))

69. m = m+1

70. h = 0.0

71. v= 0.0

72. c= 0.0

73. s= 0.0

74. g = 0.0

75. y = 0.0

76.! Matrix vector product for the initial residual.

77.! First, exchange information between processors.

Lines 78-111 are not listed.

112. r(1:nx-1,bn:en) = b(1:nx-1,bn:en)+aw*x(0:nx-2,bn:en)+&

113. ae*x(2:nx,bn:en)+as*x(1:nx-1,bn-1:en-1)+&

114. an*x(1:nx-1,bn+1:en+1)-ac*x(1:nx-1,bn:en)

115.! This preconditioner changes with the number of processors!

Lines 116-126 are not listed.

127. r(1:n-1,bn:en) = rhat(1:n-1,bn:en)

128. loc_rho = (sum(r(1:nx-1,bn:en)*r(1:nx-1,bn:en)))

129. call mpi_allreduce(loc_rho,rho,1,mpi_real,mpi_sum,&

130. mpi_comm_world,ierr)

131. rho = sqrt(rho)

132. g(1) =rho

133. v(1:nx-1,bn:en,1)=r(1:nx-1,bn:en)/rho

134. k=0

135.! Begin gmres loop.

136. do while((rho

A errtol).and.(k ? kmax))

137. k=k+1

138.! Matrix vector product.

139.! First, exchange information between processors.

Lines 140-173 are not listed.

174. v(1:nx-1,bn:en,k+1 = -aw*v(0:nx-2,bn:en,k)&

175. -ae*v(2:nx,bn:en,k)-as*v(1:nx-1,bn-1:en-1,k)&

176. -an*v(1:nx-1,bn+1:en+1,k)+ac*v(1:nx-1,bn:en,k)

177.! This preconditioner changes with the number of processors!

Lines 178-188 are not listed.

189. v(1:n-1,bn:en,k+1) = rhat(1:n-1,bn:en)

190.! Begin modiﬁed GS. May need to reorthogonalize.

191. do j=1,k

192. temp = sum(v(1:nx-1,bn:en,j)*v(1:nx-1,bn:en,k+1))

193. call mpi_allreduce(temp,h(j,k),1,mpi_real,&

194. mpi_sum,mpi_comm_world,ierr)

195. v(1:nx-1,bn:en,k+1) = v(1:nx-1,bn:en,k+1)-&

196. h(j,k)*v(1:nx-1,bn:en,j)

197. end do

198. temp = (sum(v(1:nx-1,bn:en,k+1)*v(1:nx-1,bn:en,k+1)))

9.6. GMRES(M) AND MPI 375

199. call mpi_allreduce(temp,h(k+1,k),1,mpi_real,&

200. mpi_sum,mpi_comm_world,ierr)

201. h(k+1,k) = sqrt(h(k+1,k))

202. if (h(k+1,k).gt.0.0.or.h(k+1,k).lt.0.0) then

203. v(1:nx-1,bn:en,k+1)=v(1:nx-1,bn:en,k+1)/h(k+1,k)

204. end if

205. if (k

A1) then

206.! Apply old Givens rotations to h(1:k,k).

207. do i=1,k-1

208. hik = c(i)*h(i,k)-s(i)*h(i+1,k)

209. hipk = s(i)*h(i,k)+c(i)*h(i+1,k)

210. h(i,k) = hik

211. h(i+1,k) = hipk

212. end do

213. end if

214. nu = sqrt(h(k,k)**2 + h(k+1,k)**2)

215.! May need better Givens implementation.

216.! Deﬁne and Apply new Givens rotations to h(k:k+1,k).

217. if (nu.gt.0.0) then

218. c(k) =h(k,k)/nu

219. s(k) =-h(k+1,k)/nu

220. h(k,k) =c(k)*h(k,k)-s(k)*h(k+1,k)

221. h(k+1,k) =0

222. gk =c(k)*g(k) -s(k)*g(k+1)

223. gkp =s(k)*g(k) +c(k)*g(k+1)

224. g(k) = gk

225. g(k+1) = gkp

226. end if

227. rho = abs(g(k+1))

228. mag(k) = rho

229.! End of gmres loop.

230. end do

231.! h(1:k,1:k) is upper triangular matrix in QR.

232. y(k) = g(k)/h(k,k)

233. do i = k-1,1,-1

234. y(i) = g(i)

235. do j = i+1,k

236. y(i) = y(i) -h(i,j)*y(j)

237. end do

238. y(i) = y(i)/h(i,i)

239. end do

240.! Form linear combination.

241. do i = 1,k

242. x(1:nx-1,bn:en) = x(1:nx-1,bn:en) + v(1:nx-1,bn:en,i)*y(i)

243. end do

376 CHAPTER 9. KRYLOV METHODS FOR AX = D

Table 9.6.1: MPI Times for gmresmmpi.f

p time iteration

2 358.7 10,9

4 141.6 9,8

8 096.6 10,42

16 052.3 10,41

32 049.0 12,16

244.! Send the local solutions to processor zero.

245. if (my_rank.eq.0) then

246. do source = 1,proc-1

247. sbn = 1+(source)*loc_n

248. call mpi_recv(x(0,sbn),(n+1)*loc_n,mpi_real,&

249. source,50,mpi_comm_world,status,ierr)

250. end do

251. else

252. call mpi_send(x(0,bn),(n+1)*loc_n,mpi_real,0,50,&

253. mpi_comm_world,ierr)

254. end if

255. ! End restart loop.

256. end do

257. if (my_rank.eq.0) then

258. tend = timef()

259. print*, m, mag(k)

260. print*, m,k,x(513,513)

261. print*, ’time =’, tend

262. end if

263. call mpi_ﬁnalize(ierr)

264. end program

The Table 9.6.1 contains computations for

q = 1025 using z = 1=8= The

computation times are in seconds, and note the number of iterations changes

with the number of processors. The restarts are after 50 inner iterations, and

the iterations in the third column are (o uter, inner) so that the total is outer *

50 + inner.

9.6.1 Exercises

1. Examine the full code gmresmmpi.f and identify the concurrent computa-

tions. Also study the communications that are required to do the matrix-vector

2. Verify the computations in Table 9.6.1. Experiment with di

erent num-

ber of iterations used before restarting GMRES.

product, which are similar to those used in Section 6.6 and illustrated in Figures

6.6.1 and 6.6.2.

9.6. GMRES(M) AND MPI 377

3. Experiment with variations on the SSOR preconditioner and include dif-

ferent

q and $=

4. Experiment with variations of the SSOR preconditioner to include the

use of a coarse mesh in the additive Schwarz preconditioner.

5. Use an ADI preconditioner in place of the SSOR preconditioner.

Bibliography

[1] E. Anderson, Z. Bai, C. Bischof, J. Demmel, J. Dongarra, J. Du Croz,

A. Greenbaum, S. Hammarling, A. McKenney, S. Ostrouchov and D.

Sorensen, LAPACK Users’ Guide, SIAM, 2nd ed., 1995.

[2] Edward Beltrami, Mathematical Models for Society and Biology, Academic

Press, 2002.

[3] M. Bertero and P. Boccacci, Introduction to Inverse Problems in Imaging,

IOP Publishing, Bristol, UK, 1998.

[4] Richard J. Burden and Douglas J. Faires, Numerical Analysis, Brooks Cole,

7th ed., 2000.

[5] Edmond Chow and Yousef Saad, Approximate inverse techniques for block-

partitioned matrices, SIAM J. Sci. Comp., vol. 18, no. 6, pp. 1657-1675,

Nov. 1997.

[6] Jack J. Dongarra, Iain S. Du

, Danny C. Sorensen and Henk A. van der

Vorst, Numerical Linear Algebra for High-Performance Computers, SIAM,

1998.

[7] Loyd D. Fosdick, Elizabeth J. Jessup, Carolyn J. C. Schauble and Gitta

Domik, Introduction to High-Performance Scientiﬁc Computing, MIT

Press, 1996.

[8] William Gropp, Ewing Lusk, Anthony Skjellum and Rajeev Thahur, Using

MPI 2nd Edition: Portable Parallel Programming with Message Passing

Interface, MIT Press, 2nd ed., 1999.

[9] Marcus J. Grote and Thomas Huckle, Parallel preconditioning with sparse

approximate inverses, SIAM J. Sci. Comp., vol. 18, no. 3, pp. 838-853,

May 1997.

[10] Michael T. Heath, Scientiﬁc Computing, Second Edition, McGraw-Hill,

2001.

[11] C. T. Kelley, Iterative Methods for Linear and Nonlinear Equations, SIAM,

1995.

379

380 BIBLIOGRAPHY

[12] David E. Keyes, Yousef Saad and Donald G. Truhlar (editors), Domain-

Based Parallelism and Problem Decomposition Methods in Computional

Science and Engineering, SIAM, 1995.

[13] Rubin H. Landau and M. J. Paez, Computational Physics, Problem Solving

with Computers, John Wiley, 1997.

[15] Nathan Mattor, Timothy J. Williams and Dennis W. Hewett, Algorithm

for solving tridiagonal matrix problems in parallel, Parallel Computing, vol.

21, pp. 1769-1782, 1995.

[16] Carl D. Meyer, Matrix Analysis and Applied Linear Algebra, SIAM, 2000.

[17] NPACI (National Parternship for Advanced Computational Infrastruc-

[18] NCSC (North Carolina Computing Center), NCSC User Guide,

[19] Akira Okubo and Simon A. Levin, Di

usion and Ecological Problems:

Modern Perspectives, 2nd ed., Springer-Verlag, 2001.

[20] James J. Ortega, Introduction to Parallel and Vector Solution of Linear

Systems, Plenum Press, 1988.

[21] P. S. Pacheco, Parall el Programming with MPI, Morgan-Kaufmann, 1996.

[23] G. D. Smith, Numerical Solution of Partial Di

erential Equations, Oxford,

3rd ed., 1985.

[24] J. Stoer and R. Bulirsch, Introduction to Numerical Analysis, Springer-

Verlag, 1992.

[26] C. R. Vogel, Computation Methods for Inverse Problems, SIAM, 2002.

[27] R. E. White, An Introduction to the Finite Element Method with Applica-

tions to Nonlinear Problems, Wiley, 1985.

[28] Paul Wilmott, Sam Howison and Je

 Dewynne, The Mathematics of Fi-

nancial Derivatives, Cambridge, 1995.

[29] Joseph L. Zachary, Introduction to Scientiﬁc Programming: Computational

Problem Solving Using Maple and C, Springer-Verlag, 1996.

[14] The MathWorks Inc., http://www.mathworks.com.

ture), http://www.npaci.edu.

http://www.ncsc.org/usersupport/USERGUIDE/toc.html.

[22] Shodor Education Foundation, Inc., http://www.shodor.org.

[25] UMAP Journal, http://www.comap.com.