1c* ///////////////////////////////////////////////////////////////////////////
2c* @file    cgmgd.f
3c* @author  Michael Holst
4c* @brief   Multigrid-Preconditioned CG.
5c* @version $Id: cgmgd.f,v 1.2 2010/08/12 05:45:31 fetk Exp $
6c* @attention
7c* @verbatim
8c*
9c* PMG -- Parallel algebraic MultiGrid
10c* Copyright (C) 1994-- Michael Holst.
11c*
12c* Michael Holst <mholst@math.ucsd.edu>
13c* University of California, San Diego
14c* Department of Mathematics, 5739 AP&M
15c* 9500 Gilman Drive, Dept. 0112
16c* La Jolla, CA 92093-0112 USA
17c* http://math.ucsd.edu/~mholst
18c*
19c* This file is part of PMG.
20c*
21c* PMG is free software; you can redistribute it and/or modify
22c* it under the terms of the GNU General Public License as published by
23c* the Free Software Foundation; either version 2 of the License, or
24c* (at your option) any later version.
25c*
26c* PMG is distributed in the hope that it will be useful,
27c* but WITHOUT ANY WARRANTY; without even the implied warranty of
28c* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
29c* GNU General Public License for more details.
30c*
31c* You should have received a copy of the GNU General Public License
32c* along with PMG; if not, write to the Free Software
33c* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA
34c*
35c* Linking PMG statically or dynamically with other modules is making a
36c* combined work based on PMG. Thus, the terms and conditions of the GNU
37c* General Public License cover the whole combination.
38c*
39c* SPECIAL GPL EXCEPTION
40c* In addition, as a special exception, the copyright holders of PMG
41c* give you permission to combine the PMG program with free software
42c* programs and libraries that are released under the GNU LGPL or with
43c* code included in releases of ISIM, PMV, PyMOL, SMOL, VMD, and Vision.
44c* Such combined software may be linked with PMG and redistributed together
45c* in original or modified form as mere aggregation without requirement that
46c* the entire work be under the scope of the GNU General Public License.
47c* This special exception permission is also extended to any software listed
48c* in the SPECIAL GPL EXCEPTION clauses by the FEtk and APBS libraries.
49c*
50c* Note that people who make modified versions of PMG are not obligated
51c* to grant this special exception for their modified versions; it is
52c* their choice whether to do so. The GNU General Public License gives
53c* permission to release a modified version without this exception; this
54c* exception also makes it possible to release a modified version which
55c* carries forward this exception.
56c*
57c* @endverbatim
58c* ///////////////////////////////////////////////////////////////////////////
59
60      subroutine fcgmg(nx,ny,nz,x,iz,w0,w1,w2,w3,
61     2   istop,itmax,iters,ierror,nlev,ilev,nlev_real,mgsolv,
62     3   iok,iinfo,epsiln,errtol,omega,nu1,nu2,mgsmoo,
63     4   w4,w5,w6,
64     5   ipc,rpc,pc,ac,cc,fc,tru)
65c* *********************************************************************
66c* purpose:
67c*
68c*    nested iteration for multilevel preconditioned cg
69c*
70c* author:  michael holst
71c* *********************************************************************
72      implicit         none
73c*
74c*    *** other declarations ***
75      integer          ipc(*),iz(50,*),iok,ilev,iinfo,nlev,itmax
76      integer          iters,ierror,level,itmxd,nlevd,iterd,iokd,istop
77      integer          nx,ny,nz,nxf,nyf,nzf,nxc,nyc,nzc,nlev_real,istpd
78      integer          nu1,nu2,mgsmoo,mgsolv
79      double precision epsiln,errd,errtol,omega
80      double precision x(*),w0(*),w1(*),w2(*),w3(*)
81      double precision rpc(*),pc(*),ac(*),cc(*),fc(*),tru(*)
82      double precision w4(*),w5(*),w6(*)
83c*
84c*    *** recover gridsizes ***
85      nxf = nx
86      nyf = ny
87      nzf = nz
88      call mkcors(nlev-1,nxf,nyf,nzf,nxc,nyc,nzc)
89c*
90c*    *** move up grids: interpolate solution to finer, do cgmg ***
91      if (iinfo.ne.0) then
92         write(6,100)'% FCGMG: starting:  ',nxf,nyf,nzf,nxc,nyc,nzc
93 100     format(a,2(2x,' [',i3,',',i3,',',i3,'] '))
94      endif
95      do 10 level = nlev_real, ilev+1, -1
96c*
97c*       *** call mv cycle ***
98         errd  = 1.0e-5
99         itmxd = 1
100         nlevd = nlev_real - level + 1
101         iterd = 0
102         iokd  = 0
103         istpd = 1
104         if (iinfo .ge. 2) iokd = 2
105         call cgmg(nxc,nyc,nzc,x,iz,w0,w1,w2,w3,
106     2      istpd,itmxd,iterd,ierror,nlevd,level,nlev_real,mgsolv,
107     3      iokd,iinfo,epsiln,errtol,omega,nu1,nu2,mgsmoo,
108     4      w4,w5,w6,
109     5      ipc,rpc,pc,ac,cc,fc,tru)
110c*
111c*       *** find new grid size ***
112         call mkfine(1,nxc,nyc,nzc,nxf,nyf,nzf)
113c*
114c*       *** interpolate to next finer grid (use correct bc's) ***
115         call interp(nxc,nyc,nzc,nxf,nyf,nzf,
116     2      x(iz(1,level)),x(iz(1,level-1)),pc(iz(11,level-1)))
117c*
118c*       *** new grid size ***
119         nxc = nxf
120         nyc = nyf
121         nzc = nzf
122 10   continue
123c*
124c*    *** call mv cycle ***
125      level = ilev
126      call cgmg(nxf,nyf,nzf,x,iz,w0,w1,w2,w3,
127     2   istop,itmax,iters,ierror,nlev,level,nlev_real,mgsolv,
128     3   iok,iinfo,epsiln,errtol,omega,nu1,nu2,mgsmoo,
129     4   w4,w5,w6,
130     5   ipc,rpc,pc,ac,cc,fc,tru)
131c*
132c*    *** return and end ***
133      return
134      end
135      subroutine cgmg(nx,ny,nz,x,iz,w0,w1,w2,w3,
136     2   istop,itmax,iters,ierror,nlev,ilev,nlev_real,mgsolv,
137     3   iok,iinfo,epsiln,errtol,omega,nu1,nu2,mgsmoo,
138     4   rr,zz,pp,
139     5   ipc,rpc,pc,ac,cc,fc,tru)
140c* *********************************************************************
141c* purpose:
142c*
143c*    multilevel preconditioned cg
144c*
145c* author:  michael holst
146c* *********************************************************************
147      implicit         none
148c*
149c*    *** other declarations ***
150      integer          ipc(*),iz(50,*),iok,ilev,iinfo,nlev,level,lev
151      integer          itmax,iters,ierror,istop,nu1,nu2,mgsmoo
152      integer          itmax_s,iters_s,ierror_s,iok_s,iinfo_s,istop_s
153      integer          nx,ny,nz,mgsolv,nlev_real
154      double precision errtol,epsiln,rhok1,rhok2,ztmp,omega
155      double precision rsden,rsnrm,orsnrm,xnrm1,xnrm2,xdot
156      double precision x(*),w0(*),w1(*),w2(*),w3(*)
157      double precision rpc(*),pc(*),ac(*),cc(*),fc(*),tru(*)
158      double precision rr(*),zz(*),pp(*)
159c*
160c*    *** recover level information ***
161      level = 1
162      lev   = (ilev-1)+level
163c*
164c*    *** do some i/o if requested ***
165      if (iinfo.ne.0) then
166         write(6,100)'% CGMG: starting:   ',nx,ny,nz
167 100     format(a,(2x,' [',i3,',',i3,',',i3,'] '))
168      endif
169c*
170c*    *** initial wall clock ***
171      if (iok.ne.0) then
172         call prtini(istop)
173         call prtstp(iok,-1,0.0d0,0.0d0,0.0d0)
174      endif
175c*
176c*    **************************************************************
177c*    *** note: if (iok.ne.0) then:  use a stopping test.        ***
178c*    ***       else:  use just the itmax to stop iteration.     ***
179c*    **************************************************************
180c*    *** istop=0 most efficient (whatever it is)                ***
181c*    *** istop=1 relative residual                              ***
182c*    *** istop=2 rms difference of successive iterates          ***
183c*    *** istop=3 relative true error (provided for testing)     ***
184c*    **************************************************************
185c*
186c*    *** compute denominator for stopping criterion ***
187      if (istop .eq. 0) then
188         rsden = 1.0d0
189      elseif (istop .eq. 1) then
190         rsden = xnrm1(nx,ny,nz,fc(iz(1,lev)))
191      elseif (istop .eq. 2) then
192         rsden = dsqrt(dble(nx*ny*nz))
193      elseif (istop .eq. 3) then
194         rsden = xnrm2(nx,ny,nz,tru(iz(1,lev)))
195      elseif (istop .eq. 4) then
196         rsden = xnrm2(nx,ny,nz,tru(iz(1,lev)))
197      elseif (istop .eq. 5) then
198         call matvec(nx,ny,nz,
199     2      ipc(iz(5,lev)),rpc(iz(6,lev)),
200     3      ac(iz(7,lev)),cc(iz(1,lev)),
201     4      tru(iz(1,lev)),w1)
202         rsden = dsqrt(xdot(nx,ny,nz,tru(iz(1,lev)),w1))
203      else
204         print*,'% CGMG: bad istop value... '
205      endif
206      if (rsden.eq.0.0d0) then
207         rsden = 1.0d0
208         print*,'% CGMG: rhs is zero...'
209      endif
210      rsnrm = rsden
211      orsnrm = rsnrm
212      if (iok.ne.0) call prtstp (iok,0,rsnrm,rsden,orsnrm)
213c*
214c* *********************************************************************
215c* *** begin cg iteration
216c* *********************************************************************
217c*
218c*    *** compute the initial residual (always required) ***
219      call mresid(nx,ny,nz,
220     2   ipc(iz(5,lev)),rpc(iz(6,lev)),
221     3   ac(iz(7,lev)),cc(iz(1,lev)),fc(iz(1,lev)),
222     4   x(iz(1,lev)),rr(iz(1,lev)))
223c*
224c* *********************************************************************
225c* ** *** no preconditioning ***
226c* ** call xcopy(nx,ny,nz,rr(iz(1,lev)),zz(iz(1,lev)))
227c* *********************************************************************
228c*    *** multilevel preconditioning ***
229c*
230c*    *** restrict residual to form rhs of coarse grid systems ***
231      call getpre (nx,ny,nz,iz,ilev,nlev_real,rr,pc)
232c*
233c*    *** do a linear multigrid solve of the precond equations ***
234      call azeros(nx,ny,nz,zz(iz(1,lev)))
235      istop_s   = 1
236      itmax_s   = 1
237      iters_s   = 0
238      ierror_s  = 0
239      iinfo_s   = 0
240      iok_s     = 0
241c* ***if ((iinfo .ge. 2) .and. (ilev .eq. 1)) iok_s = 2
242      call mvcs(nx,ny,nz,zz,iz,w0,w1,w2,w3,
243     2   istop_s,itmax_s,iters_s,ierror_s,
244     3   nlev,ilev,nlev_real,mgsolv,
245     4   iok_s,iinfo_s,epsiln,errtol,omega,nu1,nu2,mgsmoo,
246     5   ipc,rpc,pc,ac,cc,rr,tru)
247c* *********************************************************************
248c*
249c*    *** setup ***
250      rhok1 = xdot(nx,ny,nz,zz(iz(1,lev)),rr(iz(1,lev)))
251      if (rhok1 .eq. 0.0d0) goto 99
252c*
253c*    *** do the cg iteration ***
254      iters = 0
255 30   continue
256         iters = iters + 1
257c*
258c*       *** save iterate if stop test will use it on next iter ***
259         if (istop .eq. 2) call xcopy(nx,ny,nz,x(iz(1,lev)),
260     2      tru(iz(1,lev)))
261c*
262c*       *** main computation ***
263         if (iters .eq. 1) then
264            call xcopy(nx,ny,nz,zz(iz(1,lev)),pp)
265         else
266            call xaxpy(nx,ny,nz,(rhok2/rhok1),zz(iz(1,lev)),pp)
267            call xscal(nx,ny,nz,(rhok1/rhok2),pp)
268         endif
269         call matvec(nx,ny,nz,
270     2      ipc(iz(5,lev)),rpc(iz(6,lev)),
271     3      ac(iz(7,lev)),cc(iz(1,lev)),
272     4      pp,w1)
273         ztmp = rhok1 / xdot(nx,ny,nz,pp,w1)
274         call xaxpy(nx,ny,nz,ztmp,pp,x(iz(1,lev)))
275         call xaxpy(nx,ny,nz,(-ztmp),w1,rr(iz(1,lev)))
276c*
277c* *********************************************************************
278c* ***** *** no preconditioning ***
279c* ***** call xcopy(nx,ny,nz,rr(iz(1,lev)),zz(iz(1,lev)))
280c* *********************************************************************
281c*       *** multilevel preconditioning ***
282c*
283c*       *** restrict residual to form rhs of coarse grid systems ***
284         call getpre (nx,ny,nz,iz,ilev,nlev_real,rr,pc)
285c*
286c*       *** do a linear multigrid solve of the precond equations ***
287         call azeros(nx,ny,nz,zz(iz(1,lev)))
288         istop_s   = 1
289         itmax_s   = 1
290         iters_s   = 0
291         ierror_s  = 0
292         iinfo_s   = 0
293         iok_s     = 0
294c* ******if ((iinfo .ge. 2) .and. (ilev .eq. 1)) iok_s = 2
295         call mvcs(nx,ny,nz,zz,iz,w0,w1,w2,w3,
296     2      istop_s,itmax_s,iters_s,ierror_s,
297     3      nlev,ilev,nlev_real,mgsolv,
298     4      iok_s,iinfo_s,epsiln,errtol,omega,nu1,nu2,mgsmoo,
299     5      ipc,rpc,pc,ac,cc,rr,tru)
300c* *********************************************************************
301c*
302c*       *** new residual ***
303         rhok2 = rhok1
304         rhok1 = xdot(nx,ny,nz,zz(iz(1,lev)),rr(iz(1,lev)))
305         if (rhok1 .eq. 0.0d0) goto 99
306c*
307c*       *** compute/check the current stopping test ***
308         orsnrm = rsnrm
309         if (istop .eq. 0) then
310            rsnrm = xnrm1(nx,ny,nz,rr(iz(1,lev)))
311         elseif (istop .eq. 1) then
312            rsnrm = xnrm1(nx,ny,nz,rr(iz(1,lev)))
313         elseif (istop .eq. 2) then
314            call xcopy(nx,ny,nz,tru(iz(1,lev)),w1)
315            call xaxpy(nx,ny,nz,(-1.0d0),x(iz(1,lev)),w1)
316            rsnrm = xnrm1(nx,ny,nz,w1)
317         elseif (istop .eq. 3) then
318            call xcopy(nx,ny,nz,tru(iz(1,lev)),w1)
319            call xaxpy(nx,ny,nz,(-1.0d0),x(iz(1,lev)),w1)
320            rsnrm = xnrm2(nx,ny,nz,w1)
321         elseif (istop .eq. 4) then
322            call xcopy(nx,ny,nz,tru(iz(1,lev)),w1)
323            call xaxpy(nx,ny,nz,(-1.0d0),x(iz(1,lev)),w1)
324            rsnrm = xnrm2(nx,ny,nz,w1)
325         elseif (istop .eq. 5) then
326            call xcopy(nx,ny,nz,tru(iz(1,lev)),w1)
327            call xaxpy(nx,ny,nz,(-1.0d0),x(iz(1,lev)),w1)
328            call matvec(nx,ny,nz,
329     2         ipc(iz(5,lev)),rpc(iz(6,lev)),
330     3         ac(iz(7,lev)),cc(iz(1,lev)),
331     4         w1,w2)
332            rsnrm = dsqrt(xdot(nx,ny,nz,w1,w2))
333         else
334            print*,'% MVCS: bad istop value... '
335         endif
336         if (iok.ne.0) then
337            call prtstp (iok,iters,rsnrm,rsden,orsnrm)
338         endif
339c*
340c*       *** check iteration count and tolerance ***
341         if ((rsnrm/rsden) .le. errtol) goto 99
342         if (iters .ge. itmax) goto 99
343c*
344c*    *** main loop ***
345      goto 30
346c*
347c*    *** return and end ***
348 99   continue
349      return
350      end
351      subroutine getpre(nx,ny,nz,iz,lev,nlev_real,r,pc)
352c* *********************************************************************
353c* purpose:
354c*
355c*    form the residual on all levels to prepare for multilevel prec.
356c*
357c* author:  michael holst
358c* *********************************************************************
359      implicit         none
360      integer          iz(50,*),nx,ny,nz,nlev_real,lev,level
361      integer          nxx,nyy,nzz,nxold,nyold,nzold
362      double precision pc(*),r(*)
363c*
364c*    *** setup ***
365      nxx    = nx
366      nyy    = ny
367      nzz    = nz
368c*
369c*    *** build the (nlev-1) level operators ***
370      do 10 level = lev+1, nlev_real
371         nxold = nxx
372         nyold = nyy
373         nzold = nzz
374         call mkcors(1,nxold,nyold,nzold,nxx,nyy,nzz)
375c*
376c*       *** make the coarse grid rhs functions ***
377         call restrc(nxold,nyold,nzold,nxx,nyy,nzz,
378     2      r(iz(1,level-1)),r(iz(1,level)),pc(iz(11,level-1)))
379 10   continue
380c*
381c*    *** end it ***
382      return
383      end
384