1c* /////////////////////////////////////////////////////////////////////////// 2c* @file mgcsd.f 3c* @author Michael Holst 4c* @brief The core linear (correction scheme) multigrid routines. 5c* @version $Id: mgcsd.f,v 1.2 2010/08/12 05:45:33 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 fmvcs(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 ipc,rpc,pc,ac,cc,fc,tru) 64c* ********************************************************************* 65c* purpose: 66c* 67c* nested iteration for a linear multilevel method. 68c* 69c* algorithm: linear multigrid iteration (cs) 70c* 71c* this routine is the full multigrid front-end for a multigrid 72c* v-cycle solver. in other words, at repeatedly calls the v-cycle 73c* multigrid solver on successively finer and finer grids. 74c* 75c* author: michael holst 76c* ********************************************************************* 77 implicit none 78c* 79c* *** other declarations *** 80 integer ipc(*),iz(50,*),iok,ilev,iinfo,nlev,itmax 81 integer iters,ierror,level,itmxd,nlevd,iterd,iokd,istop 82 integer nx,ny,nz,nxf,nyf,nzf,nxc,nyc,nzc,nlev_real,istpd 83 integer nu1,nu2,mgsmoo,iinfod,mgsolv 84 double precision epsiln,errd,errtol,omega 85 double precision x(*),w0(*),w1(*),w2(*),w3(*) 86 double precision rpc(*),pc(*),ac(*),cc(*),fc(*),tru(*) 87c* 88c* *** recover gridsizes *** 89 nxf = nx 90 nyf = ny 91 nzf = nz 92 call mkcors(nlev-1,nxf,nyf,nzf,nxc,nyc,nzc) 93c* 94c* *** move up grids: interpolate solution to finer, do v cycle *** 95 if (iinfo.ne.0) then 96 write(6,100)'% FMVCS: starting: ',nxf,nyf,nzf,nxc,nyc,nzc 97 100 format(a,2(2x,' [',i3,',',i3,',',i3,'] ')) 98 endif 99 do 10 level = nlev_real, ilev+1, -1 100c* 101c* *** call mv cycle *** 102 errd = 1.0e-5 103 itmxd = 1 104 nlevd = nlev_real - level + 1 105 iterd = 0 106 iokd = 2 107 iinfod = iinfo 108 istpd = istop 109 call mvcs(nxc,nyc,nzc,x,iz,w0,w1,w2,w3, 110 2 istpd,itmxd,iterd,ierror,nlevd,level,nlev_real,mgsolv, 111 3 iokd,iinfod,epsiln,errtol,omega,nu1,nu2,mgsmoo, 112 4 ipc,rpc,pc,ac,cc,fc,tru) 113c* 114c* *** find new grid size *** 115 call mkfine(1,nxc,nyc,nzc,nxf,nyf,nzf) 116c* 117c* *** interpolate to next finer grid *** 118 call interp(nxc,nyc,nzc,nxf,nyf,nzf, 119 2 x(iz(1,level)),x(iz(1,level-1)),pc(iz(11,level-1))) 120c* 121c* *** new grid size *** 122 nxc = nxf 123 nyc = nyf 124 nzc = nzf 125 10 continue 126c* 127c* *** call mv cycle *** 128 level = ilev 129 call mvcs(nxf,nyf,nzf,x,iz,w0,w1,w2,w3, 130 2 istop,itmax,iters,ierror,nlev,level,nlev_real,mgsolv, 131 3 iok,iinfo,epsiln,errtol,omega,nu1,nu2,mgsmoo, 132 4 ipc,rpc,pc,ac,cc,fc,tru) 133c* 134c* *** return and end *** 135 return 136 end 137 subroutine mvcs(nx,ny,nz,x,iz,w0,w1,w2,w3, 138 2 istop,itmax,iters,ierror,nlev,ilev,nlev_real,mgsolv, 139 3 iok,iinfo,epsiln,errtol,omega,nu1,nu2,mgsmoo, 140 4 ipc,rpc,pc,ac,cc,fc,tru) 141c* ********************************************************************* 142c* purpose: 143c* 144c* screaming linear multilevel method. 145c* 146c* algorithm: linear multigrid iteration (cs) 147c* 148c* multigrid v-cycle solver. 149c* 150c* input: 151c* (1) fine and coarse grid discrete linear operators: L_h, L_H 152c* (2) fine grid source function: f_h 153c* (3) fine grid approximate solution: u_h 154c* 155c* output: 156c* (1) fine grid improved solution: u_h 157c* 158c* the two-grid algorithm is: 159c* (1) pre-smooth: u1_h = smooth(L_h,f_h,u_h) 160c* (2) restrict defect: d_H = r(L_h(u1_h) - f_h) 161c* (3) solve for correction: c_H = L_H^{-1}(d_H) 162c* (4) prolongate and correct: u2_h = u1_h - p(c_H) 163c* (5) post-smooth: u_h = smooth(L_h,f_h,u2_h) 164c* 165c* (of course, c_H is determined with another two-grid algorithm) 166c* 167c* implementation notes: 168c* (0) "u1_h" must be kept on each level until "c_H" is computed, 169c* and then both are used to compute "u2_h". 170c* (1) "u_h" (and then "u1_h") on all levels is stored in the "x" array. 171c* (2) "d_H" is identically "f_h" for f_h on the next coarser grid. 172c* (3) "c_h" is identically "u_h" for u_h on the next coarser grid. 173c* (4) "d_H" is stored in the "r" array (must be kept for post-smooth). 174c* (5) "f_h" is stored in the "fc" array. 175c* (6) "L_h" on all levels is stored in the "ac" array. 176c* (7) signs may be reveresed; i.e., residual is used in place 177c* of the defect in places, etc. 178c* 179c* author: michael holst 180c* ********************************************************************* 181 implicit none 182c* 183c* *** other declarations *** 184 integer ipc(*),iz(50,*),iok,ilev,iinfo,nlev,level,lev 185 integer itmax,iters,ierror,istop,nu1,nu2,mgsmoo 186 integer itmax_s,iters_s,nuuu,ivariv,mgsmoo_s,iresid 187 integer nx,ny,nz,nxf,nyf,nzf,nxc,nyc,nzc 188 integer lpv,n,m,lda,mgsolv,nlev_real,iadjoint 189 double precision omega,errtol,epsiln,errtol_s 190 double precision rsden,rsnrm,orsnrm,xnrm1,xnrm2,xdot,xnum,xden 191 double precision xdamp 192 double precision x(*),w0(*),w1(*),w2(*),w3(*) 193 double precision rpc(*),pc(*),ac(*),cc(*),fc(*),tru(*) 194c* 195c* *** recover level information *** 196 level = 1 197 lev = (ilev-1)+level 198c* 199c* *** recover gridsizes *** 200 nxf = nx 201 nyf = ny 202 nzf = nz 203 call mkcors(nlev-1,nxf,nyf,nzf,nxc,nyc,nzc) 204c* 205c* *** do some i/o if requested *** 206 if (iinfo.ne.0) then 207 write(6,100)'% MVCS: starting: ',nxf,nyf,nzf,nxc,nyc,nzc 208 100 format(a,2(2x,' [',i3,',',i3,',',i3,'] ')) 209 endif 210c* 211c* *** initial wall clock *** 212 if (iok.ne.0) then 213 call prtini(istop) 214 call prtstp(iok,-1,0.0d0,0.0d0,0.0d0) 215 endif 216c* 217c* ************************************************************** 218c* *** note: if (iok.ne.0) then: use a stopping test. *** 219c* *** else: use just the itmax to stop iteration. *** 220c* ************************************************************** 221c* *** istop=0 most efficient (whatever it is) *** 222c* *** istop=1 relative residual *** 223c* *** istop=2 rms difference of successive iterates *** 224c* *** istop=3 relative true error (provided for testing) *** 225c* ************************************************************** 226c* 227c* *** compute denominator for stopping criterion *** 228 if (iok.ne.0) then 229 if (istop .eq. 0) then 230 rsden = 1.0d0 231 elseif (istop .eq. 1) then 232 rsden = xnrm1(nxf,nyf,nzf,fc(iz(1,lev))) 233 elseif (istop .eq. 2) then 234 rsden = dsqrt(dble(nxf*nyf*nzf)) 235 elseif (istop .eq. 3) then 236 rsden = xnrm2(nxf,nyf,nzf,tru(iz(1,lev))) 237 elseif (istop .eq. 4) then 238 rsden = xnrm2(nxf,nyf,nzf,tru(iz(1,lev))) 239 elseif (istop .eq. 5) then 240 call matvec(nxf,nyf,nzf, 241 2 ipc(iz(5,lev)),rpc(iz(6,lev)), 242 3 ac(iz(7,lev)),cc(iz(1,lev)), 243 4 tru(iz(1,lev)),w1) 244 rsden = dsqrt(xdot(nxf,nyf,nzf,tru(iz(1,lev)),w1)) 245 else 246 print*,'% MVCS: bad istop value... ' 247 endif 248 if (rsden.eq.0.0d0) then 249 rsden = 1.0d0 250 print*,'% MVCS: rhs is zero on finest level ' 251 endif 252 rsnrm = rsden 253 orsnrm = rsnrm 254 call prtstp (iok,0,rsnrm,rsden,orsnrm) 255 endif 256c* 257c* ********************************************************************* 258c* *** solve directly if nlev = 1 259c* ********************************************************************* 260c* 261c* *** solve directly if on the coarse grid *** 262 if (nlev .eq. 1) then 263c* 264c* *** use iterative method? *** 265 if (mgsolv .eq. 0) then 266c* 267c* *** solve on coarsest grid with cghs, mgsmoo_s=4 (no residual) *** 268 iresid = 0 269 iadjoint = 0 270 itmax_s = 100 271 iters_s = 0 272 errtol_s = epsiln 273 mgsmoo_s = 4 274 call azeros(nxf,nyf,nzf,x(iz(1,lev))) 275 call smooth (nxf,nyf,nzf, 276 2 ipc(iz(5,lev)),rpc(iz(6,lev)), 277 3 ac(iz(7,lev)),cc(iz(1,lev)),fc(iz(1,lev)), 278 4 x(iz(1,lev)),w1,w2,w3, 279 5 itmax_s,iters_s,errtol_s,omega, 280 6 iresid,iadjoint,mgsmoo_s) 281c* 282c* *** check for trouble on the coarse grid *** 283 if (iters_s .ge. itmax_s) then 284 print*,'% MVCS: iters on coarse grid: ',iters_s 285 endif 286c* 287c* *** use direct method? *** 288 elseif (mgsolv .eq. 1) then 289c* 290c* *** setup lpv to access the factored/banded operator *** 291 lpv = lev+1 292c* 293c* *** setup for banded format *** 294 n = ipc((iz(5,lpv)-1)+1) 295 m = ipc((iz(5,lpv)-1)+2) 296 lda = ipc((iz(5,lpv)-1)+3) 297c* 298c* *** call dpbsl to solve *** 299 call xcopy_small(nxf,nyf,nzf,fc(iz(1,lev)),w1) 300 call dpbsl(ac(iz(7,lpv)),lda,n,m,w1) 301 call xcopy_large(nxf,nyf,nzf,w1,x(iz(1,lev))) 302 call fbound00(nxf,nyf,nzf,x(iz(1,lev))) 303 else 304 print*,'% MVCS: invalid coarse solver requested...' 305 endif 306c* 307c* *** compute the stopping test *** 308 iters = 1 309 if (iok.ne.0) then 310 orsnrm = rsnrm 311 if (istop .eq. 0) then 312 call mresid(nxf,nyf,nzf, 313 2 ipc(iz(5,lev)),rpc(iz(6,lev)), 314 3 ac(iz(7,lev)),cc(iz(1,lev)),fc(iz(1,lev)), 315 4 x(iz(1,lev)),w1) 316 rsnrm = xnrm1(nxf,nyf,nzf,w1) 317 elseif (istop .eq. 1) then 318 call mresid(nxf,nyf,nzf, 319 2 ipc(iz(5,lev)),rpc(iz(6,lev)), 320 3 ac(iz(7,lev)),cc(iz(1,lev)),fc(iz(1,lev)), 321 4 x(iz(1,lev)),w1) 322 rsnrm = xnrm1(nxf,nyf,nzf,w1) 323 elseif (istop .eq. 2) then 324 call xcopy(nxf,nyf,nzf,tru(iz(1,lev)),w1) 325 call xaxpy(nxf,nyf,nzf,(-1.0d0),x(iz(1,lev)),w1) 326 rsnrm = xnrm1(nxf,nyf,nzf,w1) 327 call xcopy(nxf,nyf,nzf,x(iz(1,lev)),tru(iz(1,lev))) 328 elseif (istop .eq. 3) then 329 call xcopy(nxf,nyf,nzf,tru(iz(1,lev)),w1) 330 call xaxpy(nxf,nyf,nzf,(-1.0d0),x(iz(1,lev)),w1) 331 rsnrm = xnrm2(nxf,nyf,nzf,w1) 332 elseif (istop .eq. 4) then 333 call xcopy(nxf,nyf,nzf,tru(iz(1,lev)),w1) 334 call xaxpy(nxf,nyf,nzf,(-1.0d0),x(iz(1,lev)),w1) 335 rsnrm = xnrm2(nxf,nyf,nzf,w1) 336 elseif (istop .eq. 5) then 337 call xcopy(nxf,nyf,nzf,tru(iz(1,lev)),w1) 338 call xaxpy(nxf,nyf,nzf,(-1.0d0),x(iz(1,lev)),w1) 339 call matvec(nxf,nyf,nzf, 340 2 ipc(iz(5,lev)),rpc(iz(6,lev)), 341 3 ac(iz(7,lev)),cc(iz(1,lev)), 342 4 w1,w2) 343 rsnrm = dsqrt(xdot(nxf,nyf,nzf,w1,w2)) 344 else 345 print*,'% MVCS: bad istop value... ' 346 endif 347 call prtstp (iok,iters,rsnrm,rsden,orsnrm) 348 endif 349c* 350c* *** return now *** 351 goto 99 352 endif 353c* 354c* ********************************************************************* 355c* *** begin mg iteration (note nxf,nyf,nzf changes during loop) 356c* ********************************************************************* 357c* 358c* *** setup for the v-cycle looping *** 359 iters = 0 360 30 continue 361c* 362c* *** finest level initialization *** 363 level = 1 364 lev = (ilev-1)+level 365c* 366c* *** nu1 pre-smoothings on fine grid (with residual) *** 367 iresid = 1 368 iadjoint = 0 369 iters_s = 0 370 errtol_s = 0.0d0 371 nuuu = ivariv (nu1,lev) 372 call smooth(nxf,nyf,nzf, 373 2 ipc(iz(5,lev)),rpc(iz(6,lev)), 374 3 ac(iz(7,lev)),cc(iz(1,lev)),fc(iz(1,lev)), 375 4 x(iz(1,lev)),w2,w3,w1, 376 5 nuuu,iters_s,errtol_s,omega, 377 6 iresid,iadjoint,mgsmoo) 378 call xcopy(nxf,nyf,nzf,w1,w0(iz(1,lev))) 379c* 380c* ********************************************************************* 381c* begin cycling down to coarse grid 382c* ********************************************************************* 383c* 384c* *** go down grids: restrict resid to coarser and smooth *** 385 do 40 level = 2, nlev 386 lev = (ilev-1)+level 387c* 388c* *** find new grid size *** 389 call mkcors(1,nxf,nyf,nzf,nxc,nyc,nzc) 390c* 391c* *** restrict residual to coarser grid *** 392 call restrc(nxf,nyf,nzf,nxc,nyc,nzc, 393 2 w1,w0(iz(1,lev)),pc(iz(11,lev-1))) 394c* 395c* *** new grid size *** 396 nxf = nxc 397 nyf = nyc 398 nzf = nzc 399c* 400c* *** if not on coarsest level yet... *** 401 if (level .ne. nlev) then 402c* 403c* *** nu1 pre-smoothings on this level (with residual) *** 404c* *** (w1 has residual...) *** 405 call azeros(nxf,nyf,nzf,x(iz(1,lev))) 406 iresid = 1 407 iadjoint = 0 408 iters_s = 0 409 errtol_s = 0.0d0 410 nuuu = ivariv (nu1,lev) 411 call smooth(nxf,nyf,nzf, 412 2 ipc(iz(5,lev)),rpc(iz(6,lev)), 413 3 ac(iz(7,lev)),cc(iz(1,lev)),w0(iz(1,lev)), 414 4 x(iz(1,lev)),w2,w3,w1, 415 5 nuuu,iters_s,errtol_s,omega, 416 6 iresid,iadjoint,mgsmoo) 417 endif 418c* 419c* *** end of cycling down to coarse grid loop *** 420 40 continue 421c* 422c* ********************************************************************* 423c* begin coarse grid 424c* ********************************************************************* 425c* 426c* *** coarsest level *** 427 level = nlev 428 lev = (ilev-1)+level 429c* 430c* *** use iterative method? *** 431 if (mgsolv .eq. 0) then 432c* 433c* *** solve on coarsest grid with cghs, mgsmoo_s=4 (no residual) *** 434 iresid = 0 435 iadjoint = 0 436 itmax_s = 100 437 iters_s = 0 438 errtol_s = epsiln 439 mgsmoo_s = 4 440 call azeros(nxf,nyf,nzf,x(iz(1,lev))) 441 call smooth (nxf,nyf,nzf, 442 2 ipc(iz(5,lev)),rpc(iz(6,lev)), 443 3 ac(iz(7,lev)),cc(iz(1,lev)),w0(iz(1,lev)), 444 4 x(iz(1,lev)),w1,w2,w3, 445 5 itmax_s,iters_s,errtol_s,omega, 446 6 iresid,iadjoint,mgsmoo_s) 447c* 448c* *** check for trouble on the coarse grid *** 449 if (iters_s .ge. itmax_s) then 450 print*,'% MVCS: iters on coarse grid: ',iters_s 451 endif 452c* 453c* *** use direct method? *** 454 elseif (mgsolv .eq. 1) then 455c* 456c* *** setup lpv to access the factored/banded operator *** 457 lpv = lev+1 458c* 459c* *** setup for banded format *** 460 n = ipc((iz(5,lpv)-1)+1) 461 m = ipc((iz(5,lpv)-1)+2) 462 lda = ipc((iz(5,lpv)-1)+3) 463c* 464c* *** call dpbsl to solve *** 465 call xcopy_small(nxf,nyf,nzf,w0(iz(1,lev)),w1) 466 call dpbsl(ac(iz(7,lpv)),lda,n,m,w1) 467 call xcopy_large(nxf,nyf,nzf,w1,x(iz(1,lev))) 468 call fbound00(nxf,nyf,nzf,x(iz(1,lev))) 469 else 470 print*,'% MVCS: invalid coarse solver requested...' 471 endif 472c* 473c* ********************************************************************* 474c* begin cycling back to fine grid 475c* ********************************************************************* 476c* 477c* *** move up grids: interpolate resid to finer and smooth *** 478 do 70 level = nlev-1, 1, -1 479 lev = (ilev-1)+level 480c* 481c* *** find new grid size *** 482 call mkfine(1,nxf,nyf,nzf,nxc,nyc,nzc) 483c* 484c* *** interpolate to next finer grid *** 485 call interp(nxf,nyf,nzf,nxc,nyc,nzc, 486 2 x(iz(1,lev+1)),w1,pc(iz(11,lev))) 487c* 488c* *** compute the hackbusch/reusken damping parameter *** 489c* *** which is equivalent to the standard linear cg steplength *** 490 call matvec(nxf,nyf,nzf, 491 2 ipc(iz(5,lev+1)),rpc(iz(6,lev+1)), 492 3 ac(iz(7,lev+1)),cc(iz(1,lev+1)), 493 4 x(iz(1,lev+1)),w2) 494 xnum = xdot(nxf,nyf,nzf,x(iz(1,lev+1)),w0(iz(1,lev+1))) 495 xden = xdot(nxf,nyf,nzf,x(iz(1,lev+1)),w2) 496 xdamp = xnum / xden 497c* 498c* *** new grid size *** 499 nxf = nxc 500 nyf = nyc 501 nzf = nzc 502c* 503c* *** perform the coarse grid correction *** 504CZZZZZ xdamp = 1.0d0 505 call xaxpy(nxf,nyf,nzf,xdamp,w1,x(iz(1,lev))) 506c* 507c* *** nu2 post-smoothings for correction (no residual) *** 508 iresid = 0 509 iadjoint = 1 510 iters_s = 0 511 errtol_s = 0.0d0 512 nuuu = ivariv (nu2,lev) 513 if (level .eq. 1) then 514 call smooth(nxf,nyf,nzf, 515 2 ipc(iz(5,lev)),rpc(iz(6,lev)), 516 3 ac(iz(7,lev)),cc(iz(1,lev)),fc(iz(1,lev)), 517 4 x(iz(1,lev)),w1,w2,w3, 518 5 nuuu,iters_s,errtol_s,omega, 519 6 iresid,iadjoint,mgsmoo) 520 else 521 call smooth(nxf,nyf,nzf, 522 2 ipc(iz(5,lev)),rpc(iz(6,lev)), 523 3 ac(iz(7,lev)),cc(iz(1,lev)),w0(iz(1,lev)), 524 4 x(iz(1,lev)),w1,w2,w3, 525 5 nuuu,iters_s,errtol_s,omega, 526 6 iresid,iadjoint,mgsmoo) 527 endif 528 70 continue 529c* 530c* ********************************************************************* 531c* iteration complete: do some i/o 532c* ********************************************************************* 533c* 534c* *** increment the iteration counter *** 535 iters = iters + 1 536c* 537c* *** compute/check the current stopping test *** 538 if (iok.ne.0) then 539 orsnrm = rsnrm 540 if (istop .eq. 0) then 541 call mresid(nxf,nyf,nzf, 542 2 ipc(iz(5,lev)),rpc(iz(6,lev)), 543 3 ac(iz(7,lev)),cc(iz(1,lev)),fc(iz(1,lev)), 544 4 x(iz(1,lev)),w1) 545 rsnrm = xnrm1(nxf,nyf,nzf,w1) 546 elseif (istop .eq. 1) then 547 call mresid(nxf,nyf,nzf, 548 2 ipc(iz(5,lev)),rpc(iz(6,lev)), 549 3 ac(iz(7,lev)),cc(iz(1,lev)),fc(iz(1,lev)), 550 4 x(iz(1,lev)),w1) 551 rsnrm = xnrm1(nxf,nyf,nzf,w1) 552 elseif (istop .eq. 2) then 553 call xcopy(nxf,nyf,nzf,tru(iz(1,lev)),w1) 554 call xaxpy(nxf,nyf,nzf,(-1.0d0),x(iz(1,lev)),w1) 555 rsnrm = xnrm1(nxf,nyf,nzf,w1) 556 call xcopy(nxf,nyf,nzf,x(iz(1,lev)),tru(iz(1,lev))) 557 elseif (istop .eq. 3) then 558 call xcopy(nxf,nyf,nzf,tru(iz(1,lev)),w1) 559 call xaxpy(nxf,nyf,nzf,(-1.0d0),x(iz(1,lev)),w1) 560 rsnrm = xnrm2(nxf,nyf,nzf,w1) 561 elseif (istop .eq. 4) then 562 call xcopy(nxf,nyf,nzf,tru(iz(1,lev)),w1) 563 call xaxpy(nxf,nyf,nzf,(-1.0d0),x(iz(1,lev)),w1) 564 rsnrm = xnrm2(nxf,nyf,nzf,w1) 565 elseif (istop .eq. 5) then 566 call xcopy(nxf,nyf,nzf,tru(iz(1,lev)),w1) 567 call xaxpy(nxf,nyf,nzf,(-1.0d0),x(iz(1,lev)),w1) 568 call matvec(nxf,nyf,nzf, 569 2 ipc(iz(5,lev)),rpc(iz(6,lev)), 570 3 ac(iz(7,lev)),cc(iz(1,lev)), 571 4 w1,w2) 572 rsnrm = dsqrt(xdot(nxf,nyf,nzf,w1,w2)) 573 else 574 print*,'% MVCS: bad istop value... ' 575 endif 576 call prtstp (iok,iters,rsnrm,rsden,orsnrm) 577 if ((rsnrm/rsden) .le. errtol) goto 99 578 endif 579 if (iters .ge. itmax) goto 91 580 goto 30 581c* 582c* *** problems *** 583 91 continue 584 ierror = 1 585c* 586c* *** return and end *** 587 99 continue 588 return 589 end 590 591