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