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