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