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