diff --git a/cicecore/cicedynB/dynamics/fgmresD.F90 b/cicecore/cicedynB/dynamics/fgmresD.F90 deleted file mode 100644 index d65e8cba0..000000000 --- a/cicecore/cicedynB/dynamics/fgmresD.F90 +++ /dev/null @@ -1,289 +0,0 @@ - subroutine fgmres_legacy (n,im,rhs,sol,i,vv,w,wk1, wk2, & - gamma,maxits,iout,icode,its,ro) - - use ice_fileunits, only: nu_diag - -!----------------------------------------------------------------------- -! jfl Dec 1st 2006. We modified the routine so that it is double precison. -! Here are the modifications: -! 1) implicit real (a-h,o-z) becomes implicit real*8 (a-h,o-z) -! 2) real bocomes real*8 -! 3) subroutine scopy.f has been changed for dcopy.f -! 4) subroutine saxpy.f has been changed for daxpy.f -! 5) function sdot.f has been changed for ddot.f -! 6) 1e-08 becomes 1d-08 -! -! Be careful with the dcopy, daxpy and ddot code...there is a slight -! difference with the single precision versions (scopy, saxpy and sdot). -! In the single precision versions, the array are declared sightly differently. -! It is written for single precision: -! -! modified 12/3/93, array(1) declarations changed to array(*) -!----------------------------------------------------------------------- - - implicit double precision (a-h,o-z) !jfl modification - integer n, im, maxits, iout, icode - double precision rhs(*), sol(*), vv(n,im+1),w(n,im) - double precision wk1(n), wk2(n), gamma, ro -!----------------------------------------------------------------------- -! flexible GMRES routine. This is a version of GMRES which allows a -! a variable preconditioner. Implemented with a reverse communication -! protocole for flexibility - -! DISTRIBUTED VERSION (USES DISTDOT FOR DDOT) -! explicit (exact) residual norms for restarts -! written by Y. Saad, modified by A. Malevsky, version February 1, 1995 -!----------------------------------------------------------------------- -! This Is A Reverse Communication Implementation. -!------------------------------------------------- -! USAGE: (see also comments for icode below). FGMRES -! should be put in a loop and the loop should be active for as -! long as icode is not equal to 0. On return fgmres will -! 1) either be requesting the new preconditioned vector applied -! to wk1 in case icode.eq.1 (result should be put in wk2) -! 2) or be requesting the product of A applied to the vector wk1 -! in case icode.eq.2 (result should be put in wk2) -! 3) or be terminated in case icode .eq. 0. -! on entry always set icode = 0. So icode should be set back to zero -! upon convergence. -!----------------------------------------------------------------------- -! Here is a typical way of running fgmres: -! -! icode = 0 -! 1 continue -! call fgmres (n,im,rhs,sol,i,vv,w,wk1, wk2,eps,maxits,iout,icode) -! -! if (icode .eq. 1) then -! call precon(n, wk1, wk2) <--- user's variable preconditioning -! goto 1 -! else if (icode .ge. 2) then -! call matvec (n,wk1, wk2) <--- user's matrix vector product. -! goto 1 -! else -! ----- done ---- -! ......... -!----------------------------------------------------------------------- -! list of parameters -!------------------- -! -! n == integer. the dimension of the problem -! im == size of Krylov subspace: should not exceed 50 in this -! version (can be reset in code. looking at comment below) -! rhs == vector of length n containing the right hand side -! sol == initial guess on input, approximate solution on output -! vv == work space of size n x (im+1) -! w == work space of length n x im -! wk1, -! wk2, == two work vectors of length n each used for the reverse -! communication protocole. When on return (icode .ne. 1) -! the user should call fgmres again with wk2 = precon * wk1 -! and icode untouched. When icode.eq.1 then it means that -! convergence has taken place. -! -! eps == tolerance for stopping criterion. process is stopped -! as soon as ( ||.|| is the euclidean norm): -! || current residual||/||initial residual|| <= eps -! -! maxits== maximum number of iterations allowed -! -! iout == output unit number number for printing intermediate results -! if (iout .le. 0) no statistics are printed. -! if (iout .eq. 1) L2norm of 1st ite is printed. -! if (iout .gt. 1) L2norm of all ite are printed. -! -! icode = integer. indicator for the reverse communication protocole. -! ON ENTRY : icode should be set to icode = 0. -! ON RETURN: -! * icode .eq. 1 value means that fgmres has not finished -! and that it is requesting a preconditioned vector before -! continuing. The user must compute M**(-1) wk1, where M is -! the preconditioing matrix (may vary at each call) and wk1 is -! the vector as provided by fgmres upun return, and put the -! result in wk2. Then fgmres must be called again without -! changing any other argument. -! * icode .eq. 2 value means that fgmres has not finished -! and that it is requesting a matrix vector product before -! continuing. The user must compute A * wk1, where A is the -! coefficient matrix and wk1 is the vector provided by -! upon return. The result of the operation is to be put in -! the vector wk2. Then fgmres must be called again without -! changing any other argument. -! * icode .eq. 0 means that fgmres has finished and sol contains -! the approximate solution. -! comment: typically fgmres must be implemented in a loop -! with fgmres being called as long icode is returned with -! a value .ne. 0. -!----------------------------------------------------------------------- -! local variables -- !jfl modif - double precision hh(201,200),c(200),s(200),rs(201),t,ddot,sqrt -! -!------------------------------------------------------------- -! arnoldi size should not exceed 50 in this version.. -! to reset modify sizes of hh, c, s, rs -!------------------------------------------------------------- - - save - data epsmac/1.d-16/ -! -! computed goto -! - goto (100,200,300,11) icode +1 - 100 continue - n1 = n + 1 - its = 0 -!------------------------------------------------------------- -! ** outer loop starts here.. -!--------------compute initial residual vector -------------- -! 10 continue - call dcopy (n, sol, 1, wk1, 1) !jfl modification - icode = 3 - return - 11 continue - do j=1,n - vv(j,1) = rhs(j) - wk2(j) - enddo - 20 ro = ddot(n, vv, 1, vv,1) !jfl modification - ro = sqrt(ro) - if (ro .eq. 0.0d0) goto 999 - t = 1.0d0/ ro - do j=1, n - vv(j,1) = vv(j,1)*t - enddo -! if (its .eq. 0) eps1=eps - if (its .eq. 0) then - r0 = ro - eps1=gamma*ro - endif - - if (iout .gt. 0) write(nu_diag, 199) its, ro!& -! -! initialize 1-st term of rhs of hessenberg system.. -! - rs(1) = ro - i = 0 - 4 i=i+1 - its = its + 1 - i1 = i + 1 - do k=1, n - wk1(k) = vv(k,i) - enddo -! -! return -! - icode = 1 - - return - 200 continue - do k=1, n - w(k,i) = wk2(k) - enddo -! -! call matvec operation -! - icode = 2 - call dcopy(n, wk2, 1, wk1, 1) !jfl modification -! -! return -! - return - 300 continue -! -! first call to ope corresponds to intialization goto back to 11. -! -! if (icode .eq. 3) goto 11 - call dcopy (n, wk2, 1, vv(1,i1), 1) !jfl modification -! -! modified gram - schmidt... -! - do j=1, i - t = ddot(n, vv(1,j), 1, vv(1,i1), 1) !jfl modification - hh(j,i) = t - call daxpy(n, -t, vv(1,j), 1, vv(1,i1), 1) !jfl modification - enddo - t = sqrt(ddot(n, vv(1,i1), 1, vv(1,i1), 1)) !jfl modification - hh(i1,i) = t - if (t .eq. 0.0d0) goto 58 - t = 1.0d0 / t - do k=1,n - vv(k,i1) = vv(k,i1)*t - enddo -! -! done with modified gram schimd and arnoldi step. -! now update factorization of hh -! - 58 if (i .eq. 1) goto 121 -! -! perfrom previous transformations on i-th column of h -! - do k=2,i - k1 = k-1 - t = hh(k1,i) - hh(k1,i) = c(k1)*t + s(k1)*hh(k,i) - hh(k,i) = -s(k1)*t + c(k1)*hh(k,i) - enddo - 121 gam = sqrt(hh(i,i)**2 + hh(i1,i)**2) - if (gam .eq. 0.0d0) gam = epsmac -!-----------#determine next plane rotation #------------------- - c(i) = hh(i,i)/gam - s(i) = hh(i1,i)/gam - rs(i1) = -s(i)*rs(i) - rs(i) = c(i)*rs(i) -! -! determine res. norm. and test for convergence- -! - hh(i,i) = c(i)*hh(i,i) + s(i)*hh(i1,i) - ro = abs(rs(i1)) - if (iout .gt. 1) & - write(nu_diag, 199) its, ro - if (i .lt. im .and. (ro .gt. eps1)) goto 4 -! -! now compute solution. first solve upper triangular system. -! - rs(i) = rs(i)/hh(i,i) - do ii=2,i - k=i-ii+1 - k1 = k+1 - t=rs(k) - do j=k1,i - t = t-hh(k,j)*rs(j) - enddo - rs(k) = t/hh(k,k) - enddo -! -! done with back substitution.. -! now form linear combination to get solution -! - do j=1, i - t = rs(j) - call daxpy(n, t, w(1,j), 1, sol,1) !jfl modification - enddo -! -! test for return -! - if (ro .le. eps1 .or. its .ge. maxits) goto 999 -! -! else compute residual vector and continue.. -! -! goto 10 - - do j=1,i - jj = i1-j+1 - rs(jj-1) = -s(jj-1)*rs(jj) - rs(jj) = c(jj-1)*rs(jj) - enddo - do j=1,i1 - t = rs(j) - if (j .eq. 1) t = t-1.0d0 - call daxpy (n, t, vv(1,j), 1, vv, 1) - enddo -! -! restart outer loop. -! - goto 20 - 999 icode = 0 - - 199 format('monitor_fgmres: iter_fmgres=', i4, ' L2norm=', d26.16) -! - return -!-----end-of-fgmres----------------------------------------------------- -!----------------------------------------------------------------------- - end diff --git a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 index 520f1469d..b0a969114 100644 --- a/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 +++ b/cicecore/cicedynB/dynamics/ice_dyn_vp.F90 @@ -898,23 +898,12 @@ subroutine anderson_solver (icellt, icellu, & Diagu, Diagv, & gamma, im_fgmres, & maxits_fgmres, nbiter, conv) - ! Put FGMRES solution solx,soly in fpfunc vector + ! Put FGMRES solution solx,soly in fpfunc vector (needed for anderson) call arrays_to_vec (nx_block, ny_block, nblocks, & max_blocks, icellu (:), ntot, & indxui (:,:), indxuj(:,:), & solx (:,:,:), soly (:,:,:), & fpfunc(:)) - ! FGMRES linear solver (solution is in fpfunc) - ! fpfunc = sol - ! call fgmres_solver (ntot, bvec, & - ! fpfunc, diagvec, & - ! icellt, icellu, & - ! indxti, indxtj, & - ! indxui, indxuj, & - ! zetaD, & - ! Cb, vrel, & - ! aiu, umassdti) - elseif (fpfunc_andacc == 2) then ! g_2(x) = x - A(x)x + b(x) = x - F(x) endif @@ -1085,199 +1074,6 @@ end subroutine anderson_solver !======================================================================= -! Driver for the FGMRES linear solver - - subroutine fgmres_solver (ntot, bvec, & - sol, diagvec, & - icellt, icellu, & - indxti, indxtj, & - indxui, indxuj, & - zetaD, & - Cb, vrel, & - aiu, umassdti ) - - use ice_blocks, only: nx_block, ny_block - use ice_boundary, only: ice_HaloUpdate - use ice_domain, only: halo_info - use ice_domain_size, only: max_blocks - use ice_flux, only: fm - use ice_grid, only: dxt, dyt, dxhy, dyhx, cxp, cyp, cxm, cym, & - tarear, uarear, tinyarea - use ice_state, only: uvel, vvel - - integer (kind=int_kind), intent(in) :: & - ntot ! size of problem for fgmres (for given cpu) - - real (kind=dbl_kind), dimension (ntot), intent(in) :: & - bvec , & ! RHS vector for FGMRES - diagvec ! diagonal of matrix A for preconditioners - - real (kind=dbl_kind), dimension (ntot), intent(inout) :: & - sol ! solution vector for FGMRES - - integer (kind=int_kind), dimension(max_blocks), intent(in) :: & - icellt , & ! no. of cells where icetmask = 1 - icellu ! no. of cells where iceumask = 1 - - integer (kind=int_kind), dimension (nx_block*ny_block, max_blocks), intent(in) :: & - indxti , & ! compressed index in i-direction - indxtj , & ! compressed index in j-direction - indxui , & ! compressed index in i-direction - indxuj ! compressed index in j-direction - - real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), intent(in) :: & - vrel , & ! coeff for tauw - Cb , & ! seabed stress coeff - aiu , & ! ice fraction on u-grid - umassdti ! mass of U-cell/dte (kg/m^2 s) - - real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks,4), intent(in) :: & - zetaD ! zetaD = 2zeta (viscous coeff) - - ! local variables - - integer (kind=int_kind) :: & - iblk , & ! block index - icode , & ! code for fgmres solver - its , & ! iteration nb for fgmres - fgmres_its , & ! final nb of fgmres iterations - ierr ! code for pgmres preconditioner !phb: needed? - - real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: & - Au , & ! matvec, Fx = bx - Au - Av ! matvec, Fy = by - Av - - real (kind=dbl_kind), allocatable :: & - vv(:,:), ww(:,:) ! work arrays for FGMRES - - real (kind=dbl_kind), allocatable :: & - wk11(:), wk22(:) ! work vectors for FGMRES - - real (kind=dbl_kind) :: & - res_norm ! residual norm for FGMRES - - character(len=*), parameter :: subname = '(fgmres_solver)' - - ! ! Allocate space for FGMRES work arrays - ! allocate(wk11(ntot), wk22(ntot)) - ! allocate(vv(ntot,im_fgmres+1), ww(ntot,im_fgmres)) - ! - ! !----------------------------------------------------------------------- - ! ! prep F G M R E S - ! !----------------------------------------------------------------------- - ! - ! icode = 0 - ! - ! !----------------------------------------------------------------------- - ! ! F G M R E S L O O P - ! !----------------------------------------------------------------------- - ! 1 continue - ! !----------------------------------------------------------------------- - ! - ! - ! call fgmres (ntot,im_fgmres,bvec,sol,its,vv,ww,wk11,wk22, & - ! gamma, maxits_fgmres,monitor_fgmres, & - ! icode, fgmres_its, res_norm) - ! - ! if (icode == 1) then - ! - ! if (precond .eq. 1) then - ! - ! wk22(:)=wk11(:) ! precond=identity - ! - ! elseif (precond .eq. 2) then ! use diagonal of A for precond step - ! - ! call precond_diag (ntot, & - ! diagvec (:), & - ! wk11 (:), wk22 (:) ) - ! - ! elseif (precond .eq. 3) then - ! - ! call pgmres (nx_block, ny_block, nblocks , & - ! max_blocks , icellu (:) , & - ! indxui (:,:) , indxuj (:,:) , & - ! icellt (:) , & - ! indxti (:,:) , indxtj (:,:) , & - ! dxt (:,:,:) , dyt (:,:,:) , & - ! dxhy (:,:,:) , dyhx (:,:,:) , & - ! cxp (:,:,:) , cyp (:,:,:) , & - ! cxm (:,:,:) , cym (:,:,:) , & - ! vrel (:,:,:) , Cb (:,:,:) , & - ! zetaD (:,:,:,:) , & - ! umassdti (:,:,:) , fm (:,:,:) , & - ! uarear (:,:,:) , diagvec(:) , & - ! wk22 (:) , wk11(:) , & - ! ntot , im_pgmres , & - ! epsprecond , maxits_pgmres , & - ! monitor_pgmres , ierr ) - ! endif ! precond - ! - ! goto 1 - ! - ! elseif (icode >= 2) then - ! - ! call vec_to_arrays (nx_block, ny_block, nblocks, & - ! max_blocks, icellu (:), ntot, & - ! indxui (:,:), indxuj(:,:), & - ! wk11 (:), & - ! uvel (:,:,:), vvel (:,:,:)) - ! - ! ! JFL halo update could be in subroutine... - ! !$OMP PARALLEL DO PRIVATE(iblk) - ! do iblk = 1, nblocks - ! fld2(:,:,1,iblk) = uvel(:,:,iblk) - ! fld2(:,:,2,iblk) = vvel(:,:,iblk) - ! enddo - ! !$OMP END PARALLEL DO - ! - ! call ice_HaloUpdate (fld2, halo_info, & - ! field_loc_NEcorner, field_type_vector) - ! - ! !$OMP PARALLEL DO PRIVATE(iblk) - ! do iblk = 1, nblocks - ! uvel(:,:,iblk) = fld2(:,:,1,iblk) - ! vvel(:,:,iblk) = fld2(:,:,2,iblk) - ! enddo - ! !$OMP END PARALLEL DO - ! - ! !$OMP PARALLEL DO PRIVATE(iblk) - ! do iblk = 1, nblocks - ! - ! call matvec (nx_block , ny_block, & - ! icellu (iblk) , icellt (iblk) , & - ! indxui (:,iblk) , indxuj (:,iblk) , & - ! indxti (:,iblk) , indxtj (:,iblk) , & - ! dxt (:,:,iblk) , dyt (:,:,iblk), & - ! dxhy (:,:,iblk) , dyhx (:,:,iblk), & - ! cxp (:,:,iblk) , cyp (:,:,iblk), & - ! cxm (:,:,iblk) , cym (:,:,iblk), & - ! uvel (:,:,iblk) , vvel (:,:,iblk), & - ! vrel (:,:,iblk) , Cb (:,:,iblk), & - ! zetaD (:,:,iblk,:), & - ! umassdti (:,:,iblk) , fm (:,:,iblk), & - ! uarear (:,:,iblk) , & - ! Au (:,:,iblk) , Av (:,:,iblk)) - ! - ! enddo - ! !$OMP END PARALLEL DO - ! - ! ! form wk2 from Au and Av arrays - ! call arrays_to_vec (nx_block, ny_block, nblocks, & - ! max_blocks, icellu (:), ntot, & - ! indxui (:,:), indxuj(:,:), & - ! Au (:,:,:), Av (:,:,:), & - ! wk22(:)) - ! - ! goto 1 - ! - ! endif ! icode - ! - ! deallocate(wk11, wk22, vv, ww) - - end subroutine fgmres_solver - -!======================================================================= - ! Computes the viscous coefficients (in fact zetaD=2*zeta) and dPr/dx. subroutine calc_zeta_Pr (nx_block, ny_block, & diff --git a/cicecore/cicedynB/dynamics/pgmres.F90 b/cicecore/cicedynB/dynamics/pgmres.F90 deleted file mode 100644 index 5799f209d..000000000 --- a/cicecore/cicedynB/dynamics/pgmres.F90 +++ /dev/null @@ -1,325 +0,0 @@ - -!**s/r pgmres - preconditionner for GEM_H : PGmres -! - - subroutine pgmres_legacy(nx_block, ny_block, nblocks, & - max_blocks, icellu, & - indxui, indxuj, & - icellt, & - indxti, indxtj, & - dxt, dyt, & - dxhy, dyhx, & - cxp, cyp, & - cxm, cym, & - vrel, Cb, & - zetaD, & - umassdti, fm, & - uarear, diagvec, & - sol, rhs, & - n, im, & - eps, maxits, & - iout, ierr) - -!----------------------------------------------------------------------- - -! use grid_options -! use prec - use ice_kinds_mod - use ice_dyn_vp, only: matvec, arrays_to_vec, vec_to_arrays, precond_diag - use ice_fileunits, only: nu_diag - - implicit none - -!#include - - integer (kind=int_kind), intent(in) :: & - nx_block, ny_block, & ! block dimensions - nblocks, & ! nb of blocks - max_blocks ! max nb of blocks - - - integer (kind=int_kind), dimension (max_blocks), intent(in) :: & - icellu , & - icellt ! no. of cells where icetmask = 1 - - integer (kind=int_kind), dimension (nx_block*ny_block, max_blocks), & - intent(in) :: & - indxui , & ! compressed index in i-direction - indxuj , & ! compressed index in j-direction - indxti , & ! compressed index in i-direction - indxtj ! compressed index in j-direction - - real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), intent(in) :: & - dxt , & ! width of T-cell through the middle (m) - dyt , & ! height of T-cell through the middle (m) - dxhy , & ! 0.5*(HTE - HTE) - dyhx , & ! 0.5*(HTN - HTN) - cyp , & ! 1.5*HTE - 0.5*HTE - cxp , & ! 1.5*HTN - 0.5*HTN - cym , & ! 0.5*HTE - 1.5*HTE - cxm ! 0.5*HTN - 1.5*HTN - - real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks), & - intent(in) :: & - vrel , & ! coefficient for tauw - Cb , & ! coefficient for basal stress - umassdti, & ! mass of U-cell/dt (kg/m^2 s) - fm , & ! Coriolis param. * mass in U-cell (kg/s) - uarear ! 1/uarea - - real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks,4), & - intent(in) :: & - zetaD ! 2*zeta - - real (kind=dbl_kind), dimension (n), intent(in) :: & - diagvec - - real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: & - utp , & ! x-component of velocity (m/s) - vtp , & ! y-component of velocity (m/s) - Au , & ! matvec, Fx = Au - bx (N/m^2)! jfl - Av ! matvec, Fy = Av - by (N/m^2)! jfl - - integer n, im, maxits, iout, ierr, iblk - real*8 rhs(n), sol(n) ,eps ! wk11, wk22, eps -! Abdessamad Qaddouri - 2018 -! -!revision -! v5.0 - Qaddouri A. - initial version - - real*8 vv(n,im+1), gam,eps1 - real*8 wk(n),r0 -!----------------------------------------------------------------------* - integer kmax,ii,i,j,n1,its,k1,i1,jj,k - parameter (kmax=50) - - real*8 hh(kmax+1,kmax), c(kmax), s(kmax), rs(kmax+1),t - real*8 hhloc(kmax+1,kmax) -!------------------------------------------------------------- -! arnoldi size should not exceed kmax=50 in this version.. -! to reset modify paramter kmax accordingly. -!------------------------------------------------------------- - real*8 epsmac ,ro,ddot - parameter (epsmac=1.d-16) - integer l -! character(len= 9) communicate_S -! communicate_S = "GRID" -! if (Grd_yinyang_L) communicate_S = "MULTIGRID" - - - n1 = n + 1 - its = 0 - sol=0.0 !JFL ...veut-on vraiment mettre sol = 0 ici?????? -!------------------------------------------------------------- -! outer loop starts here.. -!-------------- compute initial residual vector -------------- - do 21 j=1,n - vv(j,1) = rhs(j) - 21 continue - -!------------------------------------------------------------- - 20 continue - ro = ddot(n, vv,1,vv,1) - ro = dsqrt(ro) - - r0=ro - - if (ro .eq. 0.0d0) goto 999 - t = 1.0d0/ ro - do 210 j=1, n - vv(j,1) = vv(j,1)*t - 210 continue - if (its .eq. 0) eps1=eps*ro - if (iout .gt. 0 .and. its .eq. 0)& - write(nu_diag, 199) its, ro ,eps1 -! ** initialize 1-st term of rhs of hessenberg system.. - rs(1) = ro - i = 0 - 4 i=i+1 - its = its + 1 - i1 = i + 1 - - do l=1,n - rhs(l)= 0.0 - wk(l)= vv(l,i) - enddo -! precond - call precond_diag (n, & - diagvec (:), & - wk (:), rhs (:) ) - -! rhs = wk !!! JFL - -! matrix-vector -! call sol_matvec_H JFL - - call vec_to_arrays (nx_block, ny_block, nblocks, & - max_blocks, icellu (:), n, & - indxui (:,:), indxuj(:,:), & - rhs (:), & - utp (:,:,:), vtp (:,:,:)) - - - !$OMP PARALLEL DO PRIVATE(iblk) - do iblk = 1, nblocks - call matvec (nx_block , ny_block, & - icellu (iblk) , icellt (iblk) , & - indxui (:,iblk) , indxuj (:,iblk) , & - indxti (:,iblk) , indxtj (:,iblk) , & - dxt (:,:,iblk) , dyt (:,:,iblk), & - dxhy (:,:,iblk) , dyhx (:,:,iblk), & - cxp (:,:,iblk) , cyp (:,:,iblk), & - cxm (:,:,iblk) , cym (:,:,iblk), & - utp (:,:,iblk) , vtp (:,:,iblk), & - vrel (:,:,iblk) , Cb (:,:,iblk), & - zetaD (:,:,iblk,:), & - umassdti (:,:,iblk) , fm (:,:,iblk), & - uarear (:,:,iblk) , & - Au (:,:,iblk) , Av (:,:,iblk)) - - enddo - !$OMP END PARALLEL DO - - ! form wk2 from Au and Av arrays - call arrays_to_vec (nx_block, ny_block, nblocks, & - max_blocks, icellu (:), n, & - indxui (:,:), indxuj(:,:), & - Au (:,:,:), Av (:,:,:), & - vv(1,i1)) - -! classical gram - schmidt... -! - do 55 j=1, i - hhloc(j,i) = ddot(n, vv(1,j), 1, vv(1,i1), 1) - hh(j,i) = hhloc(j,i) - 55 continue - - do 56 j=1, i - call daxpy(n, -hh(j,i), vv(1,j), 1, vv(1,i1), 1) - 56 continue - t = ddot(n, vv(1,i1), 1, vv(1,i1), 1) -! - t=dsqrt(t) -! - - hh(i1,i) = t - if ( t .eq. 0.0d0) goto 58 - t = 1.0d0/t - do 57 k=1,n - vv(k,i1) = vv(k,i1)*t - 57 continue -! -! done with modified gram schimd and arnoldi step.. -! now update factorization of hh -! - 58 if (i == 1) goto 121 -! -! perfrom previous transformations on i-th column of h -! - do 66 k=2,i - k1 = k-1 - t = hh(k1,i) - hh(k1,i) = c(k1)*t + s(k1)*hh(k,i) - hh(k,i) = -s(k1)*t + c(k1)*hh(k,i) - 66 continue - 121 gam = sqrt(hh(i,i)**2 + hh(i1,i)**2) - -! if gamma is zero then any small value will do... -! will affect only residual estimate -! - if (gam == 0.0d0) gam = epsmac -!-----------#determinenextplane rotation #------------------- - c(i) = hh(i,i)/gam - s(i) = hh(i1,i)/gam - rs(i1) = -s(i)*rs(i) - rs(i) = c(i)*rs(i) - -! -! detrermine residual norm and test for convergence- -! - hh(i,i) = c(i)*hh(i,i) + s(i)*hh(i1,i) - ro = abs(rs(i1)) - if (iout .gt. 0) & - write(nu_diag, 199) its, ro , eps1 - if (i .lt. im .and. (ro .gt. eps1)) goto 4 -! -! now compute solution. first solve upper triangular system. -! - rs(i) = rs(i)/hh(i,i) - do 30 ii=2,i - k=i-ii+1 - k1 = k+1 - t=rs(k) - do 40 j=k1,i - t = t-hh(k,j)*rs(j) - 40 continue - rs(k) = t/hh(k,k) - 30 continue -! -! form linear combination of -!,i)'s to get solution -! - t = rs(1) - do 15 k=1, n - rhs(k) = vv(k,1)*t - 15 continue - do 16 j=2, i - t = rs(j) - do 161 k=1, n - rhs(k) = rhs(k)+t*vv(k,j) - 161 continue - 16 continue -! -! call preconditioner. -! - - do l=1,n - wk(l)= rhs(l) - rhs(l)=0.0 - enddo -! precond - call precond_diag (n, & - diagvec (:), & - wk (:), rhs (:) ) -! rhs = wk !!! JFL - - do 17 k=1, n - sol(k) = sol(k) + rhs(k) - 17 continue -! -! restart outer loop when necessary -! - if (ro .le. eps1) goto 990 - if (its .ge. maxits) goto 991 -! -! else compute residual vector and continue.. -! - do 24 j=1,i - jj = i1-j+1 - rs(jj-1) = -s(jj-1)*rs(jj) - rs(jj) = c(jj-1)*rs(jj) - 24 continue - do 25 j=1,i1 - t = rs(j) - if (j .eq. 1) t = t-1.0d0 - call daxpy (n, t, vv(1,j), 1, vv, 1) - 25 continue - 199 format('monitor_pgmres: iter_pmgres=', i4, ' L2norm=', d26.16, ' epsprecond*initial_L2norm=', d26.6) -! restart outer loop. - goto 20 - 990 ierr = 0 -! write(iout, 198) its, ro/r0 - 198 format(' its =', i4, ' conv =', d20.6) - return - 991 ierr = 1 -! write(iout, 198) its, ro/r0 - - return - 999 continue - ierr = -1 -! write(iout, 198) its, ro/r0 - - return -!-----------------end of pgmres --------------------------------------- -!----------------------------------------------------------------------- - end