C -*- Mode: Fortran; -*- c----------------------------------------------------------------- c Ravi Samtaney c Copyright 2004 c Princeton Plasma Physics Laboratory c All Rights Reserved c----------------------------------------------------------------- c $Log$ c----------------------------------------------------------------------- subroutine NewtonSolve c----------------------------------------------------------------- use mesh use mesh_common use properties #ifdef PARALLEL use mpistuff #endif use OptionFlags use profiles c implicit none c c=======Declarations========= c integer:: maxiter, iter double precision:: L2err, LInfErr, tolerance double precision:: lamda, epsilon integer:: ndfdu double precision, allocatable:: dFdU(:,:) double precision, allocatable:: rhs(:) double precision, allocatable:: Wtmp(:) double precision, allocatable:: WtmpI(:) double precision, allocatable:: Vtmp(:,:) integer:: i,j,k,l,m,n,ilo,ihi c c iter=1 maxiter=10 tolerance=1.D-08 lamda=0.1D0 L2err=1.D16 c c epsilon=1.D-05 ilo=1 ihi=2 c ndfdu=2*(nxlocal-1)*(nylocal-1) c ndfdu=2*(nxlocal-11)*nylocal ndfdu=5 allocate(dFdU(ndfdu,ndfdu)) allocate(Vtmp(ndfdu,ndfdu)) allocate(Wtmp(ndfdu)) allocate(rhs(ndfdu)) allocate(WtmpI(ndfdu)) c dfdu(:,:)=0.D0 do i=1,ndfdu,1 dfdu(i,i)=-1.D0 enddo do i=2,ndfdu,1 dfdu(i,i-1)=2.D0 enddo c rhs=0.D0 rhs(1)=1.D0 do i=1,ndfdu,1 write(6,*) 'i',i,dfdu(i,:) enddo call KrylovSolve(rhs,dFdU,ndfdu) c return end c----------------------------------------------------------------------- subroutine KrylovSolve(rhs,dFdU,ndfdu & ) use mesh use mesh_common use properties #ifdef PARALLEL use mpistuff #endif use OptionFlags c implicit none c c=======Declarations========= c double precision:: epsilon, l2norm c integer:: ndfdu double precision:: rhs(ndfdu) double precision:: dfdu(ndfdu,ndfdu) c integer:: nsize double precision,allocatable:: w(:) double precision,allocatable:: Frhs(:) c c GMRES variables integer:: rsize integer::nwork double precision,allocatable::work(:) integer::irc(5) integer:: icntl(8) double precision:: cntl(5) integer:: info(3) double precision::rinfo(2) integer matvec, precondLeft, precondRight, dotProd parameter (matvec=1, precondLeft=2, precondRight=3, dotProd=4) integer revcom, colx, coly, colz, nbscal c integer:: i,j,k,l c c epsilon=1.D-05 c nsize=ndfdu allocate(w(nsize)) allocate(Frhs(nsize)) c rsize=ndfdu c Work space nwork=rsize*rsize+rsize*(nsize+5)+6*nsize+1 c nwork=rsize*rsize+rsize*(nsize+5)+5*nsize+1 write(6,*) 'Unknown space size',nsize, ndfdu write(6,*) 'Work space size',nwork allocate(work(nwork)) work=0.D0 c c Initialize control params to default value call init_DGMRES(icntl,cntl) c control parameters c icntl(1)=6 c icntl(2)=6 icntl(3)=20 c Right Preconditioning c icntl(4)=3 c No preconditioning icntl(4)=0 c icntl(5)=0 c icntl(5)=2 icntl(5)=3 c icntl(6)=1 icntl(7)=rsize icntl(8)=0 c cntl(1)=1.D-09 c cntl(3)=1.0D0 c c c$$$* Error message stream, default is 6 c$$$ icntl(1) = 8 c$$$* Warning message stream, default is 6 c$$$ icntl(2) = 8 c$$$* Convergence history steam, default is 0 (no display) c$$$ icntl(3) = 8 c$$$* Controls the location of preconditioning, no default c$$$* 0 means no preconditioning c$$$ icntl(4) = 0 c$$$* Determines the orthogonalization scheme, default is 0 (MGS) c$$$ icntl(5) = 1 c$$$* 0 means no initial guess is supplied, default is 0 c$$$ icntl(6) = 0 c$$$* Maximum number of iterations, no default c$$$ icntl(7) = nsize*nsize c$$$* Controls residual calculation, default is 0 c$$$ icntl(8) = 0 c$$$ * Convergence tolerance, default is 1.D-5 cntl(1) = 1.D-5 * Normalizing factor alpha, default is 0.D0 cntl(2) = 1.D0 c$$$* Normalizing factor beta, default is 0.D0 c$$$ cntl(3) = 1.D0 c$$$* Normalizing factor alpha**p, default is 0.D0 c$$$ cntl(4) = 0.D0 c$$$* Normalizing factor beta**p, default is 0.D0 c$$$ cntl(5) = 0.D0 l=1 W(:)=0.D0 Frhs(:)=-rhs(:) c work(1:nsize)=w(1:nsize) work(nsize+1:2*nsize)=Frhs(1:nsize) do i=1,nsize,1 c write(6,*) 'RZ=', work(i) enddo c write(6,*) 'MaxFRHS', c & maxval(dabs(Frhs(:))) c 10 call Drive_DGMRES(nsize,nsize,rsize,nwork,work,irc,icntl,cntl, & info,rinfo) c revcom = irc(1) colx = irc(2) coly = irc(3) colz = irc(4) nbscal = irc(5) c c write(6,*) 'REVCOM=', revcom,colx,coly,colz,nbscal * if (revcom.eq.matvec) then * perform the matrix vector product * work(colz) <-- A * work(colx) call dgemv('N',nsize,nsize,1.D0,dFdU,nsize,work(colx),1, & 0.D0,work(colz),1) goto 10 * c else if (revcom.eq.precondLeft) then * perform the left preconditioning * work(colz) <-- M^{-1} * work(colx) c call dcopy(nsize,work(colx),1,work(colz),1) c goto 10 * c else if (revcom.eq.precondRight) then c write(6,*) 'Performing Right Preconditioning' * perform the right preconditioning c call dcopy(nsize,work(colx),1,work(colz),1) c call dgemv('N',nsize,nsize,1.D0,PconInv,nsize,work(colx),1, c & 0.D0,work(colz),1) * c goto 10 else if (revcom.eq.dotProd) then * perform the scalar product * work(colz) <-- work(colx) work(coly) * c write(6,*) 'Performing dot product', colx,coly,colz call dgemv('C',nsize,nbscal,1.D0,work(colx),nsize, & work(coly),1,0.D0,work(colz),1) goto 10 endif w=work(1:nsize) do i=1,nsize,1 write(6,*) 'SOL=', i,w(i), work(i+nsize) enddo c w=w+work(nsize+1:2*nsize) write(6,*) 'info=', info write(6,*) 'rinfo=', rinfo write(6,*) 'irc=', irc c write(6,*) 'Soln=',w w=matmul(dfdu,w)-work(1+nsize:2*nsize) c do i=1,nsize,1 write(6,*) 'RES=', i,w(i) enddo deallocate(work) deallocate(w) deallocate(Frhs) return end c c----------------------------------------------------------------------- subroutine ComputeL2Norm(l2norm,vec,nsize) integer:: nsize double precision:: vec(nsize), l2norm l2norm=dsqrt(sum(vec(:)*vec(:))/nsize) return end c----------------------------------------------------------------------- c----------------------------------------------------------------------- subroutine ComputeL2Norm2(l2norm,vec,rhs,nsize) integer:: nsize double precision:: vec(nsize), rhs(nsize),l2norm l2norm=dsqrt(sum((vec(:)-rhs(:))**2)/nsize) return end c----------------------------------------------------------------------- c