1c* ///////////////////////////////////////////////////////////////////////////
2c* @file    ncgdrvd.f
3c* @author  Michael Holst
4c* @brief   Driver for the nonlinear CG methods.
5c* @version $Id: ncgdrvd.f,v 1.2 2010/08/12 05:45:35 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 ncghsdriv(iparm,rparm,iwork,rwork,u,
61     2   xf,yf,zf,gxcf,gycf,gzcf,a1cf,a2cf,a3cf,ccf,fcf,tcf)
62c* *********************************************************************
63c* purpose:
64c*
65c*    linear/nonlinear conjugate gradient driver (fletcher-reeves).
66c*
67c* author:  michael holst
68c* *********************************************************************
69      implicit         none
70c*
71c*    *** input parameters ***
72      integer          iparm(*),iwork(*)
73      double precision rparm(*),rwork(*),u(*)
74      double precision xf(*),yf(*),zf(*),gxcf(*),gycf(*),gzcf(*)
75      double precision a1cf(*),a2cf(*),a3cf(*),ccf(*),fcf(*),tcf(*)
76c*
77c*    *** variables returned from mgsz ***
78      integer          iretot,iintot
79c*
80c*    *** misc variables ***
81      integer          nrwk,niwk,nx,ny,nz,nlev,ierror,n
82      integer          k_iz,k_w0
83      integer          k_ipc,k_rpc,k_ac,k_cc,k_fc
84      integer          n_iz,n_ipc,n_rpc
85c*
86c*    *** decode some parameters ***
87      nrwk    = iparm(1)
88      niwk    = iparm(2)
89      nx      = iparm(3)
90      ny      = iparm(4)
91      nz      = iparm(5)
92      nlev    = iparm(6)
93      n       = nx * ny * nz
94      n_iz    = 10*(nlev+1)
95      n_ipc   = 100*(nlev+1)
96      n_rpc   = 100*(nlev+1)
97c*
98c*    *** compute required work array sizes ***
99      iintot  = n_iz + n_ipc
100      iretot  = n_rpc + (3*n) + (4*n)
101c*
102c*    *** some more checks on input ***
103      if ((nrwk.lt.iretot) .or. (niwk.lt.iintot)) then
104         print*,'% NCGHSDRIV: real    work space must be: ',iretot
105         print*,'% NCGHSDRIV: integer work space must be: ',iintot
106         ierror = -3
107         iparm(51) = ierror
108         return
109      endif
110c*
111c*    *** iwork offsets ***
112      k_iz  = 1
113      k_ipc = k_iz  + n_iz
114c*
115c*    *** rwork offsets ***
116      k_rpc = 1
117      k_cc  = k_rpc + n_rpc
118      k_fc  = k_cc  + n
119      k_w0  = k_fc  + n
120      k_ac  = k_w0  + n
121c* ***k_ac_after  = k_ac + 4*n
122c*
123c*    *** call the multigrid driver ***
124      call ncghsdriv2(iparm,rparm,nx,ny,nz,u,iwork(k_iz),
125     2   rwork(k_w0),
126     4   iwork(k_ipc),rwork(k_rpc),
127     5   rwork(k_ac),rwork(k_cc),rwork(k_fc),
128     8   xf,yf,zf,gxcf,gycf,gzcf,a1cf,a2cf,a3cf,ccf,fcf,tcf)
129c*
130c*    *** return and end ***
131      return
132      end
133      subroutine ncghsdriv2(iparm,rparm,nx,ny,nz,u,iz,
134     2   w0,
135     3   ipc,rpc,ac,cc,fc,
136     4   xf,yf,zf,gxcf,gycf,gzcf,a1cf,a2cf,a3cf,ccf,fcf,tcf)
137c* *********************************************************************
138c* purpose:
139c*
140c*    linear/nonlinear conjugate gradient driver (fletcher-reeves).
141c*
142c* author:  michael holst
143c* *********************************************************************
144      implicit         none
145c*
146c*    *** input parameters ***
147      integer          iparm(*),ipc(*),iz(*)
148      double precision rparm(*),rpc(*)
149      double precision u(*),w0(*),ac(*),cc(*),fc(*)
150      double precision xf(*),yf(*),zf(*),gxcf(*),gycf(*),gzcf(*)
151      double precision a1cf(*),a2cf(*),a3cf(*),ccf(*),fcf(*),tcf(*)
152c*
153c*    *** misc variables ***
154      integer          mgkey,nlev,itmax,iok,iinfo,istop,ipkey,nu1,nu2
155      integer          nx,ny,nz,ilev,ido,iters,ierror,ibound,mode
156      integer          mgprol,mgcoar,mgsolv,mgdisc
157      double precision epsiln,epsmac,errtol,omegal,omegan
158      double precision pc_dumm
159      double precision bf,oh,tsetupf,tsetupc,tsolve
160c*
161c*    *** decode the iparm array ***
162      nlev   = iparm(6)
163      nu1    = iparm(7)
164      nu2    = iparm(8)
165      mgkey  = iparm(9)
166      itmax  = iparm(10)
167      istop  = iparm(11)
168      iinfo  = iparm(12)
169      ipkey  = iparm(14)
170      mode   = iparm(16)
171      mgprol = iparm(17)
172      mgcoar = iparm(18)
173      mgdisc = iparm(19)
174      mgsolv = iparm(21)
175      errtol = rparm(1)
176      omegal = rparm(9)
177      omegan = rparm(10)
178c*
179c*    *** intitialize the iteration timer ***
180      call prtstp(0,-99,0.0d0,0.0d0,0.0d0)
181c*
182c*    *** build the multigrid data structure in iz ***
183      call buildstr (nx,ny,nz,nlev,iz)
184c*
185c*    *** start timer ***
186      call tstart(bf,oh)
187c*
188c*    *** build op and rhs on fine grid ***
189      ido = 0
190      call buildops (nx,ny,nz,nlev,ipkey,iinfo,ido,iz,
191     2   mgprol,mgcoar,mgsolv,mgdisc,
192     3   ipc,rpc,pc_dumm,ac,cc,fc,
193     4   xf,yf,zf,gxcf,gycf,gzcf,a1cf,a2cf,a3cf,ccf,fcf,tcf)
194c*
195c*    *** stop timer ***
196      call tstop(bf,oh,tsetupf)
197      print*,'% NCGHSDRIV2: fine problem setup time: ',tsetupf
198      tsetupc = 0.0d0
199c*
200c* ******************************************************************
201c* *** this overwrites the rhs array provided by pde specification
202c* ****** compute an algebraically produced rhs for the given tcf ***
203      if ((istop .eq. 4) .or. (istop .eq. 5)) then
204         if ((mode .eq. 1) .or. (mode .eq. 2)) then
205            call nmatvec(nx,ny,nz,ipc,rpc,ac,cc,tcf,fc,w0)
206         else
207            call matvec(nx,ny,nz,ipc,rpc,ac,cc,tcf,fc)
208         endif
209      endif
210c* ******************************************************************
211c*
212c*    *** determine machine epsilon ***
213      epsiln = epsmac(0)
214c*
215c*    *** impose zero dirichlet boundary conditions (now in source fcn) ***
216      call fbound00(nx,ny,nz,u)
217c*
218c*    *** MATLAB ***
219      print*,' cg = [ '
220c*
221c*    *** start timer ***
222      call tstart(bf,oh)
223c*
224c*    *** call specified multigrid method ***
225      if ((mode .eq. 0) .or. (mode .eq. 2)) then
226         print*,'% NCGHSDRIV2: linear mode...'
227         iok  = 1
228         ilev = 1
229         call cghsgo(nx,ny,nz,u,w0,a1cf,a2cf,a3cf,ccf,fcf,
230     2      istop,itmax,iters,ierror,
231     3      iok,iinfo,epsiln,errtol,omegal,
232     4      ipc,rpc,ac,cc,fc,tcf)
233      endif
234      if ((mode .eq. 1) .or. (mode .eq. 2)) then
235         print*,'% NCGHSDRIV2: nonlinear mode...'
236         iok  = 1
237         ilev = 1
238         call ncghsgo(nx,ny,nz,u,w0,a1cf,a2cf,a3cf,ccf,fcf,
239     2      istop,itmax,iters,ierror,
240     3      iok,iinfo,epsiln,errtol,omegan,
241     4      ipc,rpc,ac,cc,fc,tcf)
242      endif
243c*
244c*    *** stop timer ***
245      call tstop(bf,oh,tsolve)
246      print*,'% NCGHSDRIV2: solve time: ',tsolve
247c*
248c*    *** MATLAB ***
249      write(*,100) 'cg_sf',tsetupf,'cg_sc',tsetupc,
250     2   'cg_st',(tsetupf+tsetupc),'cg_so',tsolve
251 100  format(' ];',4(' ',a7,'=',1pe9.3,';'))
252c*
253c*    *** restore boundary conditions ***
254      ibound = 1
255      call fbound(ibound,nx,ny,nz,u,gxcf,gycf,gzcf)
256c*
257c*    *** return and end ***
258      return
259      end
260      subroutine ncghsgo(nx,ny,nz,x,r,p,ap,zk,zkp1,tmp,
261     2   istop,itmax,iters,ierror,
262     3   iok,iinfo,epsiln,errtol,omega,
263     4   ipc,rpc,ac,cc,fc,tru)
264c* *********************************************************************
265c* purpose:
266c*
267c*    nonlinear conjugate gradients (fletcher-reeves).
268c*
269c* author:  michael holst
270c* *********************************************************************
271      implicit         none
272c*
273c*    *** other declarations ***
274      integer          ipc(*),iok,iinfo
275      integer          itmax,iters,ierror
276      integer          istop
277      integer          nx,ny,nz
278      double precision omega,errtol,epsiln
279      double precision rsden,rsnrm,orsnrm,xnrm1,xnrm2,xdot
280      double precision x(*),r(*),p(*),ap(*),zk(*),zkp1(*),tmp(*)
281      double precision rpc(*),ac(*),cc(*),fc(*),tru(*)
282      double precision alpha,rhok2,beta,rhok1
283c*
284cmdir 0 0
285c*
286c*    *** do some i/o if requested ***
287      if (iinfo.ne.0) then
288         write(6,100)'% NCGHSGO: starting:',nx,ny,nz
289 100     format(a,(2x,' [',i3,',',i3,',',i3,'] '))
290      endif
291c*
292c*    *** initial wall clock ***
293      call prtini(istop)
294      call prtstp(iok,-1,0.0d0,0.0d0,0.0d0)
295c*
296c*    **************************************************************
297c*    *** note: if (iok.ne.0) then:  use a stopping test.        ***
298c*    ***       else:  use just the itmax to stop iteration.     ***
299c*    **************************************************************
300c*    *** istop=0 most efficient (whatever it is)                ***
301c*    *** istop=1 relative residual                              ***
302c*    *** istop=2 rms difference of successive iterates          ***
303c*    *** istop=3 relative true error (provided for testing)     ***
304c*    **************************************************************
305c*
306c*    *** compute denominator for stopping criterion ***
307      if (istop .eq. 0) then
308         rsden = 1.0d0
309      elseif (istop .eq. 1) then
310c*       *** compute initial residual with zero initial guess ***
311c*       *** this is analogous to the linear case where one can ***
312c*       *** simply take norm of rhs for a zero initial guess ***
313         call azeros(nx,ny,nz,tmp)
314         call nmresid(nx,ny,nz,ipc,rpc,ac,cc,fc,tmp,r,zk)
315         rsden = xnrm1(nx,ny,nz,r)
316      elseif (istop .eq. 2) then
317         rsden = dsqrt(dble((nx-2)*(ny-2)*(nz-2)))
318      elseif (istop .eq. 3) then
319         rsden = xnrm2(nx,ny,nz,tru)
320      elseif (istop .eq. 4) then
321         rsden = xnrm2(nx,ny,nz,tru)
322      elseif (istop .eq. 5) then
323         call nmatvec(nx,ny,nz,ipc,rpc,ac,cc,tru,r,zk)
324         rsden = dsqrt(xdot(nx,ny,nz,tru,r))
325      else
326         print*,'% NCGHSGO: bad istop value... '
327      endif
328      if (rsden.eq.0.0d0) then
329         rsden = 1.0d0
330         print*,'% NCGHSGO: rhs is zero '
331      endif
332      rsnrm = rsden
333      orsnrm = rsnrm
334      call prtstp (iok,0,rsnrm,rsden,orsnrm)
335c*
336c*    *** now compute residual with the initial guess ***
337      call nmresid(nx,ny,nz,ipc,rpc,ac,cc,fc,x,r,zk)
338c*
339c*    *** setup for the looping ***
340      iters = 0
341 30   continue
342c*
343c*       *** save iterate if stop test will use it on next iter ***
344         if (istop .eq. 2) call xcopy(nx,ny,nz,x,tru)
345c*
346c*       *** form new direction vector from old one and residual ***
347         rhok2 = xdot(nx,ny,nz,r,r)
348         if (iters .eq. 0) then
349            call xcopy(nx,ny,nz,r,p)
350         else
351            beta = rhok2 / rhok1
352            call xaxpy(nx,ny,nz,((1.0d0)/beta),r,p)
353            call xscal(nx,ny,nz,beta,p)
354         endif
355c*
356c*       *** nonlinear case: do a line search ***
357c*       *** (note: "ap,zk,zkp1" passed back from line search as desired) ***
358         call xcopy(nx,ny,nz,r,tmp)
359         call linesearch(nx,ny,nz,alpha,
360     2      ipc,rpc,ac,cc,fc,p,x,tmp,ap,zk,zkp1)
361c*
362c*       *** save rhok2 for next iteration ***
363         rhok1 = rhok2
364c*
365c*       *** update solution in direction p of length alpha ***
366         call xaxpy(nx,ny,nz,alpha,p,x)
367c*
368c*       *** update residual ***
369         call xaxpy(nx,ny,nz,(-alpha),ap,r)
370         call xaxpy(nx,ny,nz,(1.0d0),zk,r)
371         call xaxpy(nx,ny,nz,(-1.0d0),zkp1,r)
372c*
373c* ***** *** switch to descent if necessary ***
374c* ***** itmax_d = 10
375c* ***** iter_d = 0
376c* ***** alpha = -1.0d0
377c* ***** rsnrm_tmp = xnrm1(nx,ny,nz,r)
378c* ***** 18 continue
379c* *****    if ((rsnrm_tmp.lt.rsnrm).or.(iter_d.gt.itmax_d)) then
380c* *****       print*,'% finished with descent: r_o,r_n',rsnrm,rsnrm_tmp
381c* *****       if (iter_d .gt. 0) call xcopy(nx,ny,nz,tmp4,r)
382c* *****       goto 19
383c* *****    endif
384c* *****    print*,'% trying a descent:      r_o,r_n',rsnrm,rsnrm_tmp
385c* *****    call xcopy(nx,ny,nz,tmp2,tmp)
386c* *****    call xaxpy(nx,ny,nz,alpha,tmp3,tmp)
387c* *****    call nmresid(nx,ny,nz,ipc,rpc,ac,cc,fc,tmp,tmp4,ap)
388c* *****    rsnrm_tmp = xnrm1(nx,ny,nz,tmp4)
389c* *****    alpha = alpha / 2.0d0
390c* *****    iter_d = iter_d + 1
391c* ***** goto 18
392c* ***** 19 continue
393c*
394c*       *** some bookkeeping ***
395         iters = iters + 1
396c*
397c*       *** compute/check the current stopping test ***
398         orsnrm = rsnrm
399         if (istop .eq. 0) then
400            rsnrm = xnrm1(nx,ny,nz,r)
401         elseif (istop .eq. 1) then
402            rsnrm = xnrm1(nx,ny,nz,r)
403         elseif (istop .eq. 2) then
404            call xcopy(nx,ny,nz,tru,tmp)
405            call xaxpy(nx,ny,nz,(-1.0d0),x,tmp)
406            rsnrm = xnrm1(nx,ny,nz,tmp)
407         elseif (istop .eq. 3) then
408            call xcopy(nx,ny,nz,tru,tmp)
409            call xaxpy(nx,ny,nz,(-1.0d0),x,tmp)
410            rsnrm = xnrm2(nx,ny,nz,tmp)
411         elseif (istop .eq. 4) then
412            call xcopy(nx,ny,nz,tru,tmp)
413            call xaxpy(nx,ny,nz,(-1.0d0),x,tmp)
414            rsnrm = xnrm2(nx,ny,nz,tmp)
415         elseif (istop .eq. 5) then
416            call xcopy(nx,ny,nz,tru,tmp)
417            call xaxpy(nx,ny,nz,(-1.0d0),x,tmp)
418            call nmatvec(nx,ny,nz,ipc,rpc,ac,cc,tmp,zk,zkp1)
419            rsnrm = dsqrt(xdot(nx,ny,nz,tmp,zk))
420         else
421            print*,'% NCGHSGO: bad istop value... '
422         endif
423         call prtstp (iok,iters,rsnrm,rsden,orsnrm)
424         if ((rsnrm/rsden) .le. errtol) goto 99
425         if (iters .ge. itmax) goto 99
426      goto 30
427c*
428c*    *** return and end ***
429 99   continue
430      return
431      end
432      subroutine cghsgo(nx,ny,nz,x,r,p,ap,zk,zkp1,tmp,
433     2   istop,itmax,iters,ierror,
434     3   iok,iinfo,epsiln,errtol,omega,
435     4   ipc,rpc,ac,cc,fc,tru)
436c* *********************************************************************
437c* purpose:
438c*
439c*    linear conjugate gradients (hestenes-steifel).
440c*
441c* author:  michael holst
442c* *********************************************************************
443      implicit         none
444c*
445c*    *** other declarations ***
446      integer          ipc(*),iok,iinfo
447      integer          itmax,iters,ierror
448      integer          istop
449      integer          nx,ny,nz
450      double precision omega,errtol,epsiln
451      double precision rsden,rsnrm,orsnrm,xnrm1,xnrm2,xdot
452      double precision x(*),r(*),p(*),ap(*),zk(*),zkp1(*),tmp(*)
453      double precision rpc(*),ac(*),cc(*),fc(*),tru(*)
454      double precision alpha,rhok2,beta,rhok1,pAp
455c*
456cmdir 0 0
457c*
458c*    *** do some i/o if requested ***
459      if (iinfo.ne.0) then
460         write(6,100)'% CGHSGO: starting: ',nx,ny,nz
461 100     format(a,(2x,' [',i3,',',i3,',',i3,'] '))
462      endif
463c*
464c*    *** initial wall clock ***
465      call prtini(istop)
466      call prtstp(iok,-1,0.0d0,0.0d0,0.0d0)
467c*
468c*    **************************************************************
469c*    *** note: if (iok.ne.0) then:  use a stopping test.        ***
470c*    ***       else:  use just the itmax to stop iteration.     ***
471c*    **************************************************************
472c*    *** istop=0 most efficient (whatever it is)                ***
473c*    *** istop=1 relative residual                              ***
474c*    *** istop=2 rms difference of successive iterates          ***
475c*    *** istop=3 relative true error (provided for testing)     ***
476c*    **************************************************************
477c*
478c*    *** compute denominator for stopping criterion ***
479      if (istop .eq. 0) then
480         rsden = 1.0d0
481      elseif (istop .eq. 1) then
482         rsden = xnrm1(nx,ny,nz,fc)
483      elseif (istop .eq. 2) then
484         rsden = dsqrt(dble((nx-2)*(ny-2)*(nz-2)))
485      elseif (istop .eq. 3) then
486         rsden = xnrm2(nx,ny,nz,tru)
487      elseif (istop .eq. 4) then
488         rsden = xnrm2(nx,ny,nz,tru)
489      elseif (istop .eq. 5) then
490         call matvec(nx,ny,nz,ipc,rpc,ac,cc,tru,r)
491         rsden = dsqrt(xdot(nx,ny,nz,tru,r))
492      else
493         print*,'% CGHSGO: bad istop value... '
494      endif
495      if (rsden.eq.0.0d0) then
496         rsden = 1.0d0
497         print*,'% CGHSGO: rhs is zero '
498      endif
499      rsnrm = rsden
500      orsnrm = rsnrm
501      call prtstp (iok,0,rsnrm,rsden,orsnrm)
502c*
503c*    *** now compute residual with the initial guess ***
504      call mresid(nx,ny,nz,ipc,rpc,ac,cc,fc,x,r)
505c*
506c*    *** setup for the looping ***
507      iters = 0
508 30   continue
509c*
510c*       *** save iterate if stop test will use it on next iter ***
511         if (istop .eq. 2) call xcopy(nx,ny,nz,x,tru)
512c*
513c*       *** form new direction vector from old one and residual ***
514         rhok2 = xdot(nx,ny,nz,r,r)
515         if (iters .eq. 0) then
516            call xcopy(nx,ny,nz,r,p)
517         else
518            beta = rhok2 / rhok1
519            call xaxpy(nx,ny,nz,((1.0d0)/beta),r,p)
520            call xscal(nx,ny,nz,beta,p)
521         endif
522c*
523c*       *** linear case: alpha which minimizes energy norm of error ***
524         call matvec(nx,ny,nz,ipc,rpc,ac,cc,p,ap)
525         pAp = xdot(nx,ny,nz,p,ap)
526         alpha = rhok2 / pAp
527c*
528c*       *** save rhok2 for next iteration ***
529         rhok1 = rhok2
530c*
531c*       *** update solution in direction p of length alpha ***
532         call xaxpy(nx,ny,nz,alpha,p,x)
533c*
534c*       *** update residual ***
535         call xaxpy(nx,ny,nz,(-alpha),ap,r)
536c*
537c*       *** some bookkeeping ***
538         iters = iters + 1
539c*
540c*       *** compute/check the current stopping test ***
541         orsnrm = rsnrm
542         if (istop .eq. 0) then
543            rsnrm = xnrm1(nx,ny,nz,r)
544         elseif (istop .eq. 1) then
545            rsnrm = xnrm1(nx,ny,nz,r)
546         elseif (istop .eq. 2) then
547            call xcopy(nx,ny,nz,tru,tmp)
548            call xaxpy(nx,ny,nz,(-1.0d0),x,tmp)
549            rsnrm = xnrm1(nx,ny,nz,tmp)
550         elseif (istop .eq. 3) then
551            call xcopy(nx,ny,nz,tru,tmp)
552            call xaxpy(nx,ny,nz,(-1.0d0),x,tmp)
553            rsnrm = xnrm2(nx,ny,nz,tmp)
554         elseif (istop .eq. 4) then
555            call xcopy(nx,ny,nz,tru,tmp)
556            call xaxpy(nx,ny,nz,(-1.0d0),x,tmp)
557            rsnrm = xnrm2(nx,ny,nz,tmp)
558         elseif (istop .eq. 5) then
559            call xcopy(nx,ny,nz,tru,tmp)
560            call xaxpy(nx,ny,nz,(-1.0d0),x,tmp)
561            call matvec(nx,ny,nz,ipc,rpc,ac,cc,tmp,zk)
562            rsnrm = dsqrt(xdot(nx,ny,nz,tmp,zk))
563         else
564            print*,'% CGHSGO: bad istop value... '
565         endif
566         call prtstp (iok,iters,rsnrm,rsden,orsnrm)
567         if ((rsnrm/rsden) .le. errtol) goto 99
568         if (iters .ge. itmax) goto 99
569      goto 30
570c*
571c*    *** return and end ***
572 99   continue
573      return
574      end
575