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