1 subroutine Spsi3pre (BDiag, lhsK, col, row) 2c 3c------------------------------------------------------------------------------ 4c This is the initialization routine for the Sparse-GMRES solver. 5c It pre-preconditions the LHS mass matrix and sets up the 6c sparse preconditioners. (pre-preconditioning is block diagonal 7c scaling). 8c 9c input: 10c BDiag (nshg,nflow,nflow) : block diagonal scaling matrix 11c which is already LU factored. 12c lhsK(nflow*nflow,nnz_tot) : sparse LHS mass matrix 13c 14c output: 15c lhsK(nflow*nflow,nnz_tot) : pre-preconditioned (scaled) mass matrix 16c 17c Nahid Razmara, Spring 2000. (Sparse Matrix) 18c------------------------------------------------------------------------------ 19c 20 use pointer_data 21 22 include "common.h" 23c 24 integer col(nshg+1), row(nnz*nshg) 25 integer sparseloc, c, c1 26 real*8 lhsK(nflow*nflow,nnz_tot) 27 28c 29 30 dimension BDiag(nshg,nflow,nflow) 31 32 33c 34c.... block-diagonal pre-precondition LHS 35c 36c 37c.... reduce by columns, (left block diagonal preconditioning) 38c 39c lhsK <-- inverse(L^tilde) lhsK 40c 41c 42 if(nflow.ne.5) then 43 write(*,*)' spsi3pre.f assumed nflow=5' 44 stop 45 endif 46 do i = 1, nshg 47c 48 do k = col(i), col(i+1)-1 49c 50 lhsK(2,k) = lhsK(2,k) 51 & - BDiag(i,2,1) * lhsK(1,k) 52 lhsK(7,k) = lhsK(7,k) 53 & - BDiag(i,2,1) * lhsK(6,k) 54 lhsK(12,k) = lhsK(12,k) 55 & - BDiag(i,2,1) * lhsK(11,k) 56 lhsK(17,k) = lhsK(17,k) 57 & - BDiag(i,2,1) * lhsK(16,k) 58 lhsK(22,k) = lhsK(22,k) 59 & - BDiag(i,2,1) * lhsK(21,k) 60c 61 lhsK(3,k) = lhsK(3,k) 62 & - BDiag(i,3,1) * lhsK(1,k) 63 & - BDiag(i,3,2) * lhsK(2,k) 64 lhsK(8,k) = lhsK(8,k) 65 & - BDiag(i,3,1) * lhsK(6,k) 66 & - BDiag(i,3,2) * lhsK(7,k) 67 lhsK(13,k) = lhsK(13,k) 68 & - BDiag(i,3,1) * lhsK(11,k) 69 & - BDiag(i,3,2) * lhsK(12,k) 70 lhsK(18,k) = lhsK(18,k) 71 & - BDiag(i,3,1) * lhsK(16,k) 72 & - BDiag(i,3,2) * lhsK(17,k) 73 lhsK(23,k) = lhsK(23,k) 74 & - BDiag(i,3,1) * lhsK(21,k) 75 & - BDiag(i,3,2) * lhsK(22,k) 76c 77 lhsK(4,k) = lhsK(4,k) 78 & - BDiag(i,4,1) * lhsK(1,k) 79 & - BDiag(i,4,2) * lhsK(2,k) 80 & - BDiag(i,4,3) * lhsK(3,k) 81 lhsK(9,k) = lhsK(9,k) 82 & - BDiag(i,4,1) * lhsK(6,k) 83 & - BDiag(i,4,2) * lhsK(7,k) 84 & - BDiag(i,4,3) * lhsK(8,k) 85 lhsK(14,k) = lhsK(14,k) 86 & - BDiag(i,4,1) * lhsK(11,k) 87 & - BDiag(i,4,2) * lhsK(12,k) 88 & - BDiag(i,4,3) * lhsK(13,k) 89 lhsK(19,k) = lhsK(19,k) 90 & - BDiag(i,4,1) * lhsK(16,k) 91 & - BDiag(i,4,2) * lhsK(17,k) 92 & - BDiag(i,4,3) * lhsK(18,k) 93 lhsK(24,k) = lhsK(24,k) 94 & - BDiag(i,4,1) * lhsK(21,k) 95 & - BDiag(i,4,2) * lhsK(22,k) 96 & - BDiag(i,4,3) * lhsK(23,k) 97c 98 lhsK(5,k) = lhsK(5,k) 99 & - BDiag(i,5,1) * lhsK(1,k) 100 & - BDiag(i,5,2) * lhsK(2,k) 101 & - BDiag(i,5,3) * lhsK(3,k) 102 & - BDiag(i,5,4) * lhsK(4,k) 103 lhsK(10,k) = lhsK(10,k) 104 & - BDiag(i,5,1) * lhsK(6,k) 105 & - BDiag(i,5,2) * lhsK(7,k) 106 & - BDiag(i,5,3) * lhsK(8,k) 107 & - BDiag(i,5,4) * lhsK(9,k) 108 lhsK(15,k) = lhsK(15,k) 109 & - BDiag(i,5,1) * lhsK(11,k) 110 & - BDiag(i,5,2) * lhsK(12,k) 111 & - BDiag(i,5,3) * lhsK(13,k) 112 & - BDiag(i,5,4) * lhsK(14,k) 113 lhsK(20,k) = lhsK(20,k) 114 & - BDiag(i,5,1) * lhsK(16,k) 115 & - BDiag(i,5,2) * lhsK(17,k) 116 & - BDiag(i,5,3) * lhsK(18,k) 117 & - BDiag(i,5,4) * lhsK(19,k) 118 lhsK(25,k) = lhsK(25,k) 119 & - BDiag(i,5,1) * lhsK(21,k) 120 & - BDiag(i,5,2) * lhsK(22,k) 121 & - BDiag(i,5,3) * lhsK(23,k) 122 & - BDiag(i,5,4) * lhsK(24,k) 123 enddo 124 enddo 125 126 do i = 1, nshg 127c 128 do k = col(i), col(i+1)-1 129 j = row(k) 130c 131c.... reduce by rows, (right block diagonal preconditioning) 132c 133c lhsK <-- lhsK inverse(U^tilde) 134c 135 136 137 lhsK(1,k) = BDiag(j,1,1)*lhsK(1,k) 138 lhsK(2,k) = BDiag(j,1,1)*lhsK(2,k) 139 lhsK(3,k) = BDiag(j,1,1)*lhsK(3,k) 140 lhsK(4,k) = BDiag(j,1,1)*lhsK(4,k) 141 lhsK(5,k) = BDiag(j,1,1)*lhsK(5,k) 142 143 lhsK(6,k) = BDiag(j,2,2)*( lhsK(6,k) 144 & - BDiag(j,1,2) * lhsK(1,k) ) 145 lhsK(7,k) = BDiag(j,2,2)*( lhsK(7,k) 146 & - BDiag(j,1,2) * lhsK(2,k) ) 147 lhsK(8,k) = BDiag(j,2,2)*( lhsK(8,k) 148 & - BDiag(j,1,2) * lhsK(3,k) ) 149 lhsK(9,k) = BDiag(j,2,2)*( lhsK(9,k) 150 & - BDiag(j,1,2) * lhsK(4,k) ) 151 lhsK(10,k) = BDiag(j,2,2)*( lhsK(10,k) 152 & - BDiag(j,1,2) * lhsK(5,k) ) 153c 154 lhsK(11,k) = BDiag(j,3,3)*( lhsK(11,k) 155 & - BDiag(j,1,3) * lhsK(1,k) 156 & - BDiag(j,2,3) * lhsK(6,k) ) 157 lhsK(12,k) = BDiag(j,3,3)*( lhsK(12,k) 158 & - BDiag(j,1,3) * lhsK(2,k) 159 & - BDiag(j,2,3) * lhsK(7,k) ) 160 lhsK(13,k) = BDiag(j,3,3)*( lhsK(13,k) 161 & - BDiag(j,1,3) * lhsK(3,k) 162 & - BDiag(j,2,3) * lhsK(8,k) ) 163 lhsK(14,k) = BDiag(j,3,3)*( lhsK(14,k) 164 & - BDiag(j,1,3) * lhsK(4,k) 165 & - BDiag(j,2,3) * lhsK(9,k) ) 166 lhsK(15,k) = BDiag(j,3,3)*( lhsK(15,k) 167 & - BDiag(j,1,3) * lhsK(5,k) 168 & - BDiag(j,2,3) * lhsK(10,k) ) 169c 170 lhsK(16,k) = BDiag(j,4,4)*( lhsK(16,k) 171 & - BDiag(j,1,4) * lhsK(1,k) 172 & - BDiag(j,2,4) * lhsK(6,k) 173 & - BDiag(j,3,4) * lhsK(11,k) ) 174 lhsK(17,k) = BDiag(j,4,4)*( lhsK(17,k) 175 & - BDiag(j,1,4) * lhsK(2,k) 176 & - BDiag(j,2,4) * lhsK(7,k) 177 & - BDiag(j,3,4) * lhsK(12,k) ) 178 lhsK(18,k) = BDiag(j,4,4)*( lhsK(18,k) 179 & - BDiag(j,1,4) * lhsK(3,k) 180 & - BDiag(j,2,4) * lhsK(8,k) 181 & - BDiag(j,3,4) * lhsK(13,k) ) 182 lhsK(19,k) = BDiag(j,4,4)*( lhsK(19,k) 183 & - BDiag(j,1,4) * lhsK(4,k) 184 & - BDiag(j,2,4) * lhsK(9,k) 185 & - BDiag(j,3,4) * lhsK(14,k) ) 186 lhsK(20,k) = BDiag(j,4,4)*( lhsK(20,k) 187 & - BDiag(j,1,4) * lhsK(5,k) 188 & - BDiag(j,2,4) * lhsK(10,k) 189 & - BDiag(j,3,4) * lhsK(15,k) ) 190c 191 lhsK(21,k) = BDiag(j,5,5)*( lhsK(21,k) 192 & - BDiag(j,1,5) * lhsK(1,k) 193 & - BDiag(j,2,5) * lhsK(6,k) 194 & - BDiag(j,3,5) * lhsK(11,k) 195 & - BDiag(j,4,5) * lhsK(16,k) ) 196 lhsK(22,k) = BDiag(j,5,5)*( lhsK(22,k) 197 & - BDiag(j,1,5) * lhsK(2,k) 198 & - BDiag(j,2,5) * lhsK(7,k) 199 & - BDiag(j,3,5) * lhsK(12,k) 200 & - BDiag(j,4,5) * lhsK(17,k) ) 201 lhsK(23,k) = BDiag(j,5,5)*( lhsK(23,k) 202 & - BDiag(j,1,5) * lhsK(3,k) 203 & - BDiag(j,2,5) * lhsK(8,k) 204 & - BDiag(j,3,5) * lhsK(13,k) 205 & - BDiag(j,4,5) * lhsK(18,k) ) 206 lhsK(24,k) = BDiag(j,5,5)*( lhsK(24,k) 207 & - BDiag(j,1,5) * lhsK(4,k) 208 & - BDiag(j,2,5) * lhsK(9,k) 209 & - BDiag(j,3,5) * lhsK(14,k) 210 & - BDiag(j,4,5) * lhsK(19,k) ) 211 lhsK(25,k) = BDiag(j,5,5)*( lhsK(25,k) 212 & - BDiag(j,1,5) * lhsK(5,k) 213 & - BDiag(j,2,5) * lhsK(10,k) 214 & - BDiag(j,3,5) * lhsK(15,k) 215 & - BDiag(j,4,5) * lhsK(20,k) ) 216 enddo 217 enddo 218c 219 return 220 221 end 222 223 224