1 subroutine itrFDI (ypre, y, ac, x, 2 & rmes, uBrg, BDiag, 3 & iBC, BC, iper, 4 & ilwork, shp, shgl, 5 & shpb, shglb) 6c 7c---------------------------------------------------------------------- 8c 9c This subroutine computes the "optimum" finite difference 10c interval eps for the forward difference scheme 11c 12c Rmod(y + eps u) - Rmod(y) 13c --------------------------- 14c eps 15c 16c where u is the step and Rmod is the modified residual. 17c 18c Note: A good theoretical reference is 'Practical Optimization' by 19c P.E. Gill, W. Murray and M.H. Wright [1981]. 20c 21c input: 22c y (nshg,ndof) : Y-variables 23c ypre (nshg,ndof) : preconditioned Y-variables 24c x (numnp,nsd) : node coordinates 25c rmes (nshg,nflow) : modified residual 26c uBrg (nshg,nflow) : step 27c BDiag (nshg,nflow,nflow) : block-diagonal preconditioner 28c iBC (nshg) : BC codes 29c BC (nshg,ndofBC) : BC constraint parameters 30c 31c 32c 33c Zdenek Johan, Winter 1989. 34c Zdenek Johan, Winter 1991. (Fortran 90) 35c---------------------------------------------------------------------- 36c 37 include "common.h" 38c 39 dimension y(nshg,ndof), ypre(nshg,nflow), 40 & x(numnp,nsd), ac(nshg,ndof), 41 & rmes(nshg,nflow), 42 & uBrg(nshg,nflow), BDiag(nshg,nflow,nflow), 43 & iBC(nshg), BC(nshg,ndofBC), 44 & ilwork(nlwork), 45 & iper(nshg) 46c 47 dimension ytmp(nshg,nflow), rtmp(nshg,nflow) 48 dimension tmpy(nshg,ndof) 49c 50 dimension shp(MAXTOP,maxsh,MAXQPT), 51 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 52 & shpb(MAXTOP,maxsh,MAXQPT), 53 & shglb(MAXTOP,nsd,maxsh,MAXQPT) 54c 55c.... compute the accuracy (cancellation error) -> epsA 56c 57 rtmp = zero 58c 59 ytmp = ypre 60c 61c call yshuffle(ytmp, 'new2old ') 62c 63 call i3LU (BDiag, ytmp, 'backward') 64c 65 call yshuffle(ytmp, 'old2new ') 66c 67 iabres = 1 68c 69 call itrRes (ytmp, y, 70 & x, shp, 71 & shgl, iBC, 72 & BC, shpb, 73 & shglb, rtmp, 74 & iper, ilwork, 75 & ac) 76c 77 iabres = 0 78c 79 call i3LU (BDiag, rtmp, 'forward ') 80c 81 rtmp = rtmp**2 82 call sumgat (rtmp, nflow, summed, ilwork) 83 epsA = (epsM**2) * sqrt(summed) 84c 85c.... compute the norm of the second derivative (truncation error) 86c 87c.... set interval 88c 89 epsSD = sqrt(epsM) 90c 91c.... compute the first residual 92c 93 rtmp = zero 94c 95c call yshuffle(ypre, 'new2old ') 96c 97 ytmp = ypre + epsSD * uBrg 98c 99 call i3LU (BDiag, ytmp, 'backward') 100c 101 call yshuffle(ytmp, 'old2new ') 102c 103 call itrRes (ytmp, y, 104 & x, shp, 105 & shgl, iBC, 106 & BC, shpb, 107 & shglb, rtmp, 108 & iper, ilwork, ac) 109!Added ac to the end if itrRes, but not tested - Nicholas 110c 111c.... compute the second residual and add it to the first one 112c 113 ytmp = ypre - epsSD * uBrg 114c 115c call yshuffle(ypre, 'old2new ') 116c 117 call i3LU (BDiag, ytmp, 'backward') 118c 119 call yshuffle(ytmp, 'old2new ') 120 call itrRes (ytmp, y, 121 & x, shp, 122 & shgl, iBC, 123 & BC, shpb, 124 & shglb, rtmp, 125 & iper, ilwork, ac) 126!Added ac to the end if itrRes, but not tested - Nicholas 127c 128 call i3LU (BDiag, rtmp, 'forward ') 129c 130c.... compute the second derivative and its norm 131c 132 rtmp = (( rtmp - two * rmes ) / epsM)**2 133c 134 call sumgat (rtmp, nflow, summed, ilwork) 135 SDnrm = sqrt(summed) 136c 137c.... compute the 'optimum' interval 138c 139 eGMRES = two * sqrt( epsA / SDnrm ) 140c 141c.... flop count 142c 143 ! flops = flops + 10*nflow*nshg+3*nshg 144c 145c.... end 146c 147 return 148 end 149 150 151 subroutine itrFDISclr (y, ypre, x, 152 & rmes, uBrg, BDiag, 153 & iBC, BC, engBC, iper, 154 & ilwork, shp, shgl, wght, 155 & shpb, shglb, wghtb) 156c 157c---------------------------------------------------------------------- 158c 159c This subroutine computes the "optimum" finite difference 160c interval eps for the forward difference scheme 161c 162c Rmod(y + eps u) - Rmod(y) 163c --------------------------- 164c eps 165c 166c where u is the step and Rmod is the modified residual. 167c 168c Note: A good theoretical reference is 'Practical Optimization' by 169c P.E. Gill, W. Murray and M.H. Wright [1981]. 170c 171c input: 172c y (nshg,ndof) : Y-variables 173c ypre (nshg,ndof) : preconditioned Y-variables 174c x (numnp,nsd) : node coordinates 175c rmes (nshg,nflow) : modified residual 176c uBrg (nshg,nflow) : step 177c BDiag (nshg,nflow,nflow) : block-diagonal preconditioner 178c iBC (nshg) : BC codes 179c BC (nshg,ndofBC) : BC constraint parameters 180c engBC (nshg) : energy for BC on density or pressure 181c 182c 183c Zdenek Johan, Winter 1989. 184c Zdenek Johan, Winter 1991. (Fortran 90) 185c---------------------------------------------------------------------- 186c 187 include "common.h" 188c 189 dimension y(nshg,ndof), ypre(nshg,ndof), 190 & x(numnp,nsd), 191 & rmes(nshg,nflow), 192 & uBrg(nshg,nflow), BDiag(nshg,nflow,nflow), 193 & iBC(nshg), BC(nshg,ndofBC), 194 & engBC(nshg), ilwork(nlwork), 195 & iper(nshg) 196c 197 dimension ytmp(nshg,ndof), rtmp(nshg,nflow) 198c 199c.... compute the accuracy (cancellation error) -> epsA 200c 201 rtmp = zero 202c 203 ytmp = ypre 204c 205c call tnanq(ytmp,ndof,"ytmp ") 206 iabres = 1 207c 208 call itrRes (ytmp, y, 209 & x, shp, 210 & shgl, wght, iBC, 211 & BC, engBC, shpb, 212 & shglb, wghtb, rtmp, 213 & iper, ilwork, ac) 214!Added ac to the end if itrRes, but not tested - Nicholas 215c 216 iabres = 0 217c 218 rtmp = rtmp**2 219 call sumgat (rtmp, nflow, summed, ilwork) 220 epsA = (epsM**2) * sqrt(summed) 221c 222c.... compute the norm of the second derivative (truncation error) 223c 224c.... set interval 225c 226 epsSD = sqrt(epsM) 227c 228c.... compute the first residual 229c 230 rtmp = zero 231c 232 ytmp = ypre + epsSD * uBrg 233c 234 call itrRes (ytmp, y, 235 & x, shp, 236 & shgl, wght, iBC, 237 & BC, engBC, shpb, 238 & shglb, wghtb, rtmp, 239 & iper, ilwork, ac) 240!Added ac to the end if itrRes, but not tested - Nicholas 241c 242c.... compute the second residual and add it to the first one 243c 244 ytmp = ypre - epsSD * uBrg 245c 246 call itrRes (ytmp, y, 247 & x, shp, 248 & shgl, wght, iBC, 249 & BC, engBC, shpb, 250 & shglb, wghtb, rtmp, 251 & iper, ilwork, ac) 252!Added ac to the end if itrRes, but not tested - Nicholas 253c 254c.... compute the second derivative and its norm 255c 256 rtmp = (( rtmp - two * rmes ) / epsM)**2 257c 258 call sumgat (rtmp, nflow, summed, ilwork) 259 SDnrm = sqrt(summed) 260c 261c.... compute the 'optimum' interval 262c 263 eGMRES = two * sqrt( epsA / SDnrm ) 264c 265c.... flop count 266c 267 ! flops = flops + 10*nflow*nshg+3*nshg 268c 269c.... end 270c 271 return 272 end 273