1 subroutine i3pre (BDtmp, Binv, EGmass, ilwork ) 2c 3c--------------------------------------------------------------------- 4c This is the initialization routine for the EBE-GMRES solver. 5c It pre-preconditions the LHS mass matrix and sets up the 6c EBE preconditioners. (pre-preconditioning is block diagonal 7c scaling). 8c 9c input: 10c BDtmp (nshg,nflow,nflow) : block diagonal scaling matrix 11c which is already LU factored. 12c EGmass (numel,nedof,nedof) : element mass matrix 13c 14c output: 15c EGmass (numel,nedof,nedof) : pre-preconditioned (scaled) mass matrix 16c Binv (numel,nedof,nedof) : EBE preconditioner for each element 17c 18c--------------------------------------------------------------------- 19c 20 use pointer_data 21 22 include "common.h" 23 24 dimension BDtmp(nshg,nflow,nflow), 25 & EGmass(numel,nedof,nedof), 26ctoomuchmemory & Binv(numel,nedof,nedof), 27 & ilwork(nlwork) 28c 29 dimension BDiagl(numel,nshape,nflow,nflow), 30 & BDiag(nshg,nflow,nflow) 31 32 BDiag = BDtmp 33 34 if (numpe > 1) then 35 call commu (BDiag , ilwork, nflow*nflow , 'out') 36 endif 37c 38c.... block-diagonal pre-precondition LHS 39c 40c 41c.... loop over element blocks 42c 43 do iblk = 1, nelblk 44 iel = lcblk(1,iblk) 45 nenl = lcblk(5,iblk) 46 npro = lcblk(1,iblk+1) - iel 47 n = iel - 1 48 inum = iel + npro - 1 49 nshl = lcblk(10,iblk) 50c 51c.... localize block diagonal scaling matrices 52c 53 call local (BDiag, BDiagl(iel:inum,:,:,:), abs(mien(iblk)%p), 54 & nflow*nflow, 'gather ' ) 55c 56c.... loop through local nodes and reduce all columns and rows 57c 58 do inode = 1, nshl 59 i = (inode - 1) * nflow ! EGmass dof offset 60c 61c.... reduce by columns, (left block diagonal preconditioning) 62c 63c EGmass <-- inverse(L^tilde) EGmass 64c 65 do j = 1, nedof 66 do iv = 1, npro 67 EGmass(n+iv,i+1,j) = EGmass(n+iv,i+1,j) 68c 69 EGmass(n+iv,i+2,j) = (EGmass(n+iv,i+2,j) 70 & - BDiagl(n+iv,inode,2,1) * EGmass(n+iv,i+1,j)) 71c 72 EGmass(n+iv,i+3,j) = (EGmass(n+iv,i+3,j) 73 & - BDiagl(n+iv,inode,3,1) * EGmass(n+iv,i+1,j) 74 & - BDiagl(n+iv,inode,3,2) * EGmass(n+iv,i+2,j)) 75c 76 EGmass(n+iv,i+4,j) = (EGmass(n+iv,i+4,j) 77 & - BDiagl(n+iv,inode,4,1) * EGmass(n+iv,i+1,j) 78 & - BDiagl(n+iv,inode,4,2) * EGmass(n+iv,i+2,j) 79 & - BDiagl(n+iv,inode,4,3) * EGmass(n+iv,i+3,j)) 80c 81 EGmass(n+iv,i+5,j) = (EGmass(n+iv,i+5,j) 82 & - BDiagl(n+iv,inode,5,1) * EGmass(n+iv,i+1,j) 83 & - BDiagl(n+iv,inode,5,2) * EGmass(n+iv,i+2,j) 84 & - BDiagl(n+iv,inode,5,3) * EGmass(n+iv,i+3,j) 85 & - BDiagl(n+iv,inode,5,4) * EGmass(n+iv,i+4,j)) 86 enddo 87 enddo 88 enddo 89 90 do inode = 1, nshl 91 i = (inode - 1) * nflow ! EGmass dof offset 92 93c 94c.... reduce by rows, (right block diagonal preconditioning) 95c 96c EGmass <-- EGmass inverse(U^tilde) 97c 98 do j = 1, nedof 99 do iv = 1, npro 100 EGmass(n+iv,j,i+1) = BDiagl(n+iv,inode,1,1) * 101 & (EGmass(n+iv,j,i+1)) 102c 103 EGmass(n+iv,j,i+2) = BDiagl(n+iv,inode,2,2) * ( 104 & EGmass(n+iv,j,i+2) 105 & - BDiagl(n+iv,inode,1,2) * EGmass(n+iv,j,i+1)) 106c 107 EGmass(n+iv,j,i+3) = BDiagl(n+iv,inode,3,3) * ( 108 & EGmass(n+iv,j,i+3) 109 & - BDiagl(n+iv,inode,1,3) * EGmass(n+iv,j,i+1) 110 & - BDiagl(n+iv,inode,2,3) * EGmass(n+iv,j,i+2)) 111c 112 EGmass(n+iv,j,i+4) = BDiagl(n+iv,inode,4,4) * ( 113 & EGmass(n+iv,j,i+4) 114 & - BDiagl(n+iv,inode,1,4) * EGmass(n+iv,j,i+1) 115 & - BDiagl(n+iv,inode,2,4) * EGmass(n+iv,j,i+2) 116 & - BDiagl(n+iv,inode,3,4) * EGmass(n+iv,j,i+3)) 117c 118 EGmass(n+iv,j,i+5) = BDiagl(n+iv,inode,5,5) * ( 119 & EGmass(n+iv,j,i+5) 120 & - BDiagl(n+iv,inode,1,5) * EGmass(n+iv,j,i+1) 121 & - BDiagl(n+iv,inode,2,5) * EGmass(n+iv,j,i+2) 122 & - BDiagl(n+iv,inode,3,5) * EGmass(n+iv,j,i+3) 123 & - BDiagl(n+iv,inode,4,5) * EGmass(n+iv,j,i+4)) 124 enddo 125 enddo 126c 127c.... end loops over row and column nodes 128c 129 enddo 130c 131c.... end of element blocks loop 132c 133 enddo 134c 135c.... calculate non-symmetric Cholesky EBE preconditioners 136c 137ctoomuchmemory Binv = EGmass 138c$$$ if (iPcond .eq. 2) then 139c$$$ call itrPr2 (ieneg, lcblk, Binv, ubBgl, ubBgl, 'LDU_Fact') 140c$$$ endif 141c 142c.... end of (pre process), return 143c 144 return 145 146 end 147c 148c 149c 150 subroutine i3preSclr (Diag, Dinv, EGmassT, ilwork ) 151c 152c--------------------------------------------------------------------- 153c This is the initialization routine for the EBE-GMRES solver. 154c It pre-preconditions the LHS mass matrix and sets up the 155c EBE preconditioners. (pre-preconditioning is block diagonal 156c scaling). 157c 158c input: 159c Diag (numnp,nflow,nflow) : diagonal scaling matrix 160c EGmass (numel,nedof,nedof) : element mass matrix 161c 162c output: 163c EGmass (numel,nedof,nedof) : pre-preconditioned (scaled) mass matrix 164c Dinv (numel,nedof,nedof) : EBE preconditioner for each element 165c 166c--------------------------------------------------------------------- 167c 168 use pointer_data 169 170 include "common.h" 171 172 dimension EGmassT(numel,nshape,nshape), 173 & Dinv(nshg), ilwork(nlwork) 174c 175 dimension Dinvl(numel,nshape), Diag(nshg) 176c 177 Dinv = one/Diag 178c 179 if (numpe > 1) then 180 call commu (Dinv , ilwork, 1 , 'out') 181 endif 182c 183c.... loop over element blocks 184c 185 do iblk = 1, nelblk 186 iel = lcblk(1,iblk) 187 nenl = lcblk(5,iblk) 188 npro = lcblk(1,iblk+1) - iel 189 n = iel - 1 190 inum = iel + npro - 1 191 nshl = lcblk(10,iblk) 192c 193c.... localize diagonal scaling matrices 194c 195 call local (Dinv, Dinvl(iel:inum,:), mien(iblk)%p, 196 & 1, 'gather ' ) 197c 198c.... loop through and reduce all columns 199c 200 do icol = 1, nshl 201 do irow = 1, nshl 202 do iv = 1, npro 203 EGmassT(n+iv,irow,icol) = EGmassT(n+iv,irow,icol) 204 & *Dinvl(n+iv,icol) 205 enddo 206 enddo 207 enddo 208c 209c.... end of element blocks loop 210c 211 enddo 212c 213c.... end of (pre process), return 214c 215 return 216 217 end 218 219 220 221 222