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) 155c 156c---------------------------------------------------------------------- 157c 158c This subroutine computes the "optimum" finite difference 159c interval eps for the forward difference scheme 160c 161c Rmod(y + eps u) - Rmod(y) 162c --------------------------- 163c eps 164c 165c where u is the step and Rmod is the modified residual. 166c 167c Note: A good theoretical reference is 'Practical Optimization' by 168c P.E. Gill, W. Murray and M.H. Wright [1981]. 169c 170c input: 171c y (nshg,ndof) : Y-variables 172c ypre (nshg,ndof) : preconditioned Y-variables 173c x (numnp,nsd) : node coordinates 174c rmes (nshg,nflow) : modified residual 175c uBrg (nshg,nflow) : step 176c BDiag (nshg,nflow,nflow) : block-diagonal preconditioner 177c iBC (nshg) : BC codes 178c BC (nshg,ndofBC) : BC constraint parameters 179c engBC (nshg) : energy for BC on density or pressure 180c 181c 182c Zdenek Johan, Winter 1989. 183c Zdenek Johan, Winter 1991. (Fortran 90) 184c---------------------------------------------------------------------- 185c 186 include "common.h" 187c 188 dimension y(nshg,ndof), ypre(nshg,ndof), 189 & x(numnp,nsd), 190 & rmes(nshg,nflow), 191 & uBrg(nshg,nflow), BDiag(nshg,nflow,nflow), 192 & iBC(nshg), BC(nshg,ndofBC), 193 & engBC(nshg), ilwork(nlwork), 194 & iper(nshg) 195c 196 dimension ytmp(nshg,ndof), rtmp(nshg,nflow) 197c 198c.... compute the accuracy (cancellation error) -> epsA 199c 200 rtmp = zero 201c 202 ytmp = ypre 203c 204c call tnanq(ytmp,ndof,"ytmp ") 205 iabres = 1 206c 207 call itrRes (ytmp, y, 208 & x, a(mshp), 209 & a(mshgl), a(mwght), iBC, 210 & BC, engBC, a(mshpb), 211 & a(mshglb), a(mwghtb), rtmp, 212 & iper, ilwork, ac) 213!Added ac to the end if itrRes, but not tested - Nicholas 214c 215 iabres = 0 216c 217 rtmp = rtmp**2 218 call sumgat (rtmp, nflow, summed, ilwork) 219 epsA = (epsM**2) * sqrt(summed) 220c 221c.... compute the norm of the second derivative (truncation error) 222c 223c.... set interval 224c 225 epsSD = sqrt(epsM) 226c 227c.... compute the first residual 228c 229 rtmp = zero 230c 231 ytmp = ypre + epsSD * uBrg 232c 233 call itrRes (ytmp, y, 234 & x, a(mshp), 235 & a(mshgl), a(mwght), iBC, 236 & BC, engBC, a(mshpb), 237 & a(mshglb), a(mwghtb), rtmp, 238 & iper, ilwork, ac) 239!Added ac to the end if itrRes, but not tested - Nicholas 240c 241c.... compute the second residual and add it to the first one 242c 243 ytmp = ypre - epsSD * uBrg 244c 245 call itrRes (ytmp, y, 246 & x, a(mshp), 247 & a(mshgl), a(mwght), iBC, 248 & BC, engBC, a(mshpb), 249 & a(mshglb), a(mwghtb), rtmp, 250 & iper, ilwork, ac) 251!Added ac to the end if itrRes, but not tested - Nicholas 252c 253c.... compute the second derivative and its norm 254c 255 rtmp = (( rtmp - two * rmes ) / epsM)**2 256c 257 call sumgat (rtmp, nflow, summed, ilwork) 258 SDnrm = sqrt(summed) 259c 260c.... compute the 'optimum' interval 261c 262 eGMRES = two * sqrt( epsA / SDnrm ) 263c 264c.... flop count 265c 266 ! flops = flops + 10*nflow*nshg+3*nshg 267c 268c.... end 269c 270 return 271 end 272