1 subroutine Au2MFG (ypre, y, ac, x, 2 & rmes, res, Dy, 3 & uBrg, BDiag, iBC, BC, 4 & iper, ilwork, 5 & shp, shgl, 6 & shpb, shglb) 7! 8!---------------------------------------------------------------------- 9! 10! This routine performs a matrix-vector product for the Matrix-Free 11! Implicit/Iterative solver using a two-sided scheme and subtracts 12! it from the residual. 13! 14! input: 15! y (nshg,ndof) : Y-variables 16! ypre (nshg,ndof) : preconditioned Y-variables 17! x (numnp,nsd) : node coordinates 18! rmes (nshg,nflow) : modified residual 19! res (nshg,nflow) : residual 20! Dy (nshg,nflow) : solution step 21! BDiag (nshg,nflow,nflow) : block-diagonal preconditioner 22! iBC (nshg) : BC codes 23! BC (nshg,ndofBC) : BC constraint parameters 24! 25! output: 26! uBrg (nshg,nflow) : Krylov vector ( R - A x ) 27! 28! 29! Zdenek Johan, Winter 1991. (Fortran 90) 30!---------------------------------------------------------------------- 31! 32 include "common.h" 33! 34 dimension y(nshg,ndof), ypre(nshg,nflow), 35 & x(numnp,nsd), ac(nshg,ndof), 36 & rmes(nshg,nflow), ytmp(nshg,nflow), 37 & res(nshg,nflow), Dy(nshg,nflow), 38 & uBrg(nshg,nflow), BDiag(nshg,nflow,nflow), 39 & iBC(nshg), BC(nshg,ndofBC), 40 & iper(nshg), Dy2(nshg,nflow) 41! 42 dimension uBtmp1(nshg,nflow), uBtmp2(nshg,nflow), 43 & tmpBC(nshg), ilwork(nlwork) 44! 45 dimension shp(MAXTOP,maxsh,MAXQPT), 46 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 47 & shpb(MAXTOP,maxsh,MAXQPT), 48 & shglb(MAXTOP,nsd,maxsh,MAXQPT) 49 50!$$$ dimension shp(nshl,ngauss), shgl(nsd,nshl,ngauss), 51!$$$ & shpb(nshl,ngaussb), shglb(nsd,nshl,ngaussb) 52! 53!.... compute the finite difference interval 54! 55 Dy2 = Dy**2 56 call sumgat (Dy2, nflow, summed, ilwork) 57 eps = epsM**pt66 / sqrt(summed) 58! 59!.... calculate R(y + eps x) 60! 61 uBtmp1 = zero 62! 63! call yshuffle(ypre, 'new2old ') 64! 65 uBrg = ypre + eps * Dy 66! 67 call i3LU (BDiag, uBrg, 'backward') 68! 69 call yshuffle(uBrg, 'old2new ') 70! 71 call itrBC (uBrg, uBrg, iBC, BC, iper, ilwork) 72! 73 call itrRes (uBrg, y, 74 & x, shp, 75 & shgl, iBC, 76 & BC, shpb, 77 & shglb, uBtmp1, 78 & iper, ilwork, ac) 79!Added ac to the end if itrRes, but not tested - Nicholas 80! 81 call i3LU (BDiag, uBtmp1, 'forward ') 82! 83!.... calculate R(y - eps x) 84! 85 uBtmp2 = zero 86! 87! call yshuffle(ypre, 'new2old ') 88! 89 uBrg = ypre - eps * Dy 90! 91 call i3LU (BDiag, uBrg, 'backward') 92! 93 call yshuffle(uBrg, 'old2new ') 94 95 call itrBC (uBrg, uBrg, iBC, BC, iper, ilwork) 96! 97 call itrRes (uBrg, y, 98 & x, shp, 99 & shgl, iBC, 100 & BC, shpb, 101 & shglb, uBtmp2, 102 & iper, ilwork, ac) 103! 104 call i3LU (BDiag, uBtmp2, 'forward ') 105! 106!.... compute R(y) - ( R(y + eps x) - R(y - eps x) ) / 2 eps 107! 108 uBrg = res - ( uBtmp1 - uBtmp2 ) / (two * eps) 109! 110!.... flop count 111! 112 ! flops = flops + 8*nflow*nshg+nshg 113! 114!.... end 115! 116 return 117 end 118