1 subroutine i3LU (Diag, r, code) 2c 3c---------------------------------------------------------------------- 4c 5c This routine performs a LU factorization/solve of a set of matrices, 6c used for block diagonal preconditioning in the iterative driver. 7 8c input: 9c Diag (nshg,nflow,nflow) : block diagonal (symmetric storage) 10c r (nshg,nflow) : residual 11c code : operation code 12c .eq. 'LU_Fact ', Cholesky Factor 13c .eq. 'forward ', forward reduction 14c .eq. 'backward', backward substitution 15c .eq. 'product ', product Diag.r 16c 17c output: 18c Diag (nshg,nflow,nflow) : LU decomp. of block diagonal 19c r (nshg,nflow) : reduced residual 20c 21c 22c Note: the formulation used here is taken from Golub's "Matrix 23c Computations" Book (1984), pages 82-85 algorithm P5.1-3. 24c 25c L and U overwrite Diag 26c 27c The diagonal terms (i,i) of U are stored in inverted form. 28c 29c 30c---------------------------------------------------------------------- 31c 32 include "common.h" 33c 34 dimension Diag(nshg,nflow,nflow), r(nshg,nflow) 35c 36 character*8 code 37 38c 39c.... perform LU decomposition with the Diagonal terms inverted 40c 41 if (code .eq. 'LU_Fact ') then 42c 43 Diag(:,1,1) = one / Diag(:,1,1) 44c 45 Diag(:,2,1) = Diag(:,1,1) * Diag(:,2,1) 46 Diag(:,3,1) = Diag(:,1,1) * Diag(:,3,1) 47 Diag(:,4,1) = Diag(:,1,1) * Diag(:,4,1) 48 Diag(:,5,1) = Diag(:,1,1) * Diag(:,5,1) 49c 50 Diag(:,2,2) = Diag(:,2,2) - Diag(:,2,1) * Diag(:,1,2) 51 Diag(:,2,3) = Diag(:,2,3) - Diag(:,2,1) * Diag(:,1,3) 52 Diag(:,2,4) = Diag(:,2,4) - Diag(:,2,1) * Diag(:,1,4) 53 Diag(:,2,5) = Diag(:,2,5) - Diag(:,2,1) * Diag(:,1,5) 54 Diag(:,2,2) = one / Diag(:,2,2) 55c 56 Diag(:,3,2) = Diag(:,2,2)*(Diag(:,3,2) 57 & - Diag(:,3,1) * Diag(:,1,2)) 58 Diag(:,4,2) = Diag(:,2,2)*(Diag(:,4,2) 59 & - Diag(:,4,1) * Diag(:,1,2)) 60 Diag(:,5,2) = Diag(:,2,2)*(Diag(:,5,2) 61 & - Diag(:,5,1) * Diag(:,1,2)) 62c 63 64 Diag(:,3,3) = Diag(:,3,3) - Diag(:,3,1) * Diag(:,1,3) 65 & - Diag(:,3,2) * Diag(:,2,3) 66 Diag(:,3,4) = Diag(:,3,4) - Diag(:,3,1) * Diag(:,1,4) 67 & - Diag(:,3,2) * Diag(:,2,4) 68 Diag(:,3,5) = Diag(:,3,5) - Diag(:,3,1) * Diag(:,1,5) 69 & - Diag(:,3,2) * Diag(:,2,5) 70 Diag(:,3,3) = one / Diag(:,3,3) 71c 72 Diag(:,4,3) = Diag(:,3,3) *(Diag(:,4,3) 73 & - Diag(:,4,1) * Diag(:,1,3) 74 & - Diag(:,4,2) * Diag(:,2,3)) 75 Diag(:,5,3) = Diag(:,3,3) *(Diag(:,5,3) 76 & - Diag(:,5,1) * Diag(:,1,3) 77 & - Diag(:,5,2) * Diag(:,2,3)) 78c 79 Diag(:,4,4) = Diag(:,4,4) - Diag(:,4,1) * Diag(:,1,4) 80 & - Diag(:,4,2) * Diag(:,2,4) 81 & - Diag(:,4,3) * Diag(:,3,4) 82 Diag(:,4,4) = one / Diag(:,4,4) 83c 84 Diag(:,5,4) = Diag(:,4,4) *(Diag(:,5,4) 85 & - Diag(:,5,1) * Diag(:,1,4) 86 & - Diag(:,5,2) * Diag(:,2,4) 87 & - Diag(:,5,3) * Diag(:,3,4)) 88c 89 Diag(:,5,5) = Diag(:,5,5) - Diag(:,5,1) * Diag(:,1,5) 90 & - Diag(:,5,2) * Diag(:,2,5) 91 & - Diag(:,5,3) * Diag(:,3,5) 92 & - Diag(:,5,4) * Diag(:,4,5) 93 Diag(:,5,5) = one / Diag(:,5,5) 94c 95c INACCURATE NOW ! flops = flops + 110*nshg 96c 97 return 98 endif 99c 100c.... perform forward reduction 101c 102 if (code .eq. 'forward ') then 103c 104c r(:,1) = r(:,1) !no-op 105 r(:,2) = r(:,2) - Diag(:,2,1) * r(:,1) 106 r(:,3) = r(:,3) - Diag(:,3,1) * r(:,1) 107 & - Diag(:,3,2) * r(:,2) 108 r(:,4) = r(:,4) - Diag(:,4,1) * r(:,1) 109 & - Diag(:,4,2) * r(:,2) 110 & - Diag(:,4,3) * r(:,3) 111 r(:,5) = r(:,5) - Diag(:,5,1) * r(:,1) 112 & - Diag(:,5,2) * r(:,2) 113 & - Diag(:,5,3) * r(:,3) 114 & - Diag(:,5,4) * r(:,4) 115c 116c.... flop count 117c 118c INACCURATE ! flops = flops + 25*nshg 119c 120 return 121 endif 122c 123c.... perform backward substitution 124c 125 if (code .eq. 'backward') then 126c 127 r(:,5) = Diag(:,5,5) * r(:,5) 128 r(:,4) = Diag(:,4,4) * ( r(:,4) 129 & - r(:,5) * Diag(:,4,5) ) 130 r(:,3) = Diag(:,3,3) * ( r(:,3) 131 & - r(:,5) * Diag(:,3,5) 132 & - r(:,4) * Diag(:,3,4) ) 133 r(:,2) = Diag(:,2,2) * ( r(:,2) 134 & - r(:,5) * Diag(:,2,5) 135 & - r(:,4) * Diag(:,2,4) 136 & - r(:,3) * Diag(:,2,3) ) 137 r(:,1) = Diag(:,1,1) * ( r(:,1) 138 & - r(:,5) * Diag(:,1,5) 139 & - r(:,4) * Diag(:,1,4) 140 & - r(:,3) * Diag(:,1,3) 141 & - r(:,2) * Diag(:,1,2) ) 142c 143c.... flop count 144c 145! flops = flops + 25*nshg 146 147 return 148 endif 149c 150c.... perform product U.r 151c 152 if (code .eq. 'product ') then 153c 154 r(:,1) = r(:,1) / Diag(:,1,1) + r(:,2) * Diag(:,1,2) + 155 & r(:,3) * Diag(:,1,3) + r(:,4) * Diag(:,1,4) + 156 & r(:,5) * Diag(:,1,5) 157 r(:,2) = r(:,2) / Diag(:,2,2) + 158 & r(:,3) * Diag(:,2,3) + r(:,4) * Diag(:,2,4) + 159 & r(:,5) * Diag(:,2,5) 160 r(:,3) = r(:,3) / Diag(:,3,3) + 161 & r(:,4) * Diag(:,3,4) + 162 & r(:,5) * Diag(:,3,5) 163 r(:,4) = r(:,4) / Diag(:,4,4) + 164 & r(:,5) * Diag(:,4,5) 165 r(:,5) = r(:,5) / Diag(:,5,5) 166c 167c.... flop count 168c 169! flops = flops + 40*nshg 170 171 return 172 endif 173c 174 call error ('i3LU ', code, 0) 175c 176c.... return 177c 178c 179111 format(5(e14.7,2x)) 180 return 181 end 182