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