1#include "applu_share.h" 2#include "applu_macros.h" 3 4void blts(int k) { 5 6//c compute the regular-sparse, block lower triangular solution: 7 8 9//c local variables 10 11 int i, j, m; 12 int iex; 13 double tmp, tmp1; 14 double tmat[5*5]; 15 16//c receive data from north and west 17 18 iex = 0; 19 exchange_1(rsd, k, iex); 20 21 for (j=jst; j<=jend; j++) { 22 for (i=ist; i<=iend; i++) { 23 for (m=1; m<=5; m++) { 24 25 rsd( m, i, j, k ) = rsd( m, i, j, k ) 26 - omega * ( a( m, 1, i, j ) * rsd( 1, i, j, k-1 ) 27 + a( m, 2, i, j ) * rsd( 2, i, j, k-1 ) 28 + a( m, 3, i, j ) * rsd( 3, i, j, k-1 ) 29 + a( m, 4, i, j ) * rsd( 4, i, j, k-1 ) 30 + a( m, 5, i, j ) * rsd( 5, i, j, k-1 ) ); 31 } 32 } 33 } 34 35 for (j=jst; j<=jend; j++) { 36 for (i=ist; i<=iend; i++) { 37 for (m=1; m<=5; m++) { 38 39 rsd( m, i, j, k ) = rsd( m, i, j, k ) 40 - omega * ( b( m, 1, i, j ) * rsd( 1, i, j-1, k ) 41 + c( m, 1, i, j ) * rsd( 1, i-1, j, k ) 42 + b( m, 2, i, j ) * rsd( 2, i, j-1, k ) 43 + c( m, 2, i, j ) * rsd( 2, i-1, j, k ) 44 + b( m, 3, i, j ) * rsd( 3, i, j-1, k ) 45 + c( m, 3, i, j ) * rsd( 3, i-1, j, k ) 46 + b( m, 4, i, j ) * rsd( 4, i, j-1, k ) 47 + c( m, 4, i, j ) * rsd( 4, i-1, j, k ) 48 + b( m, 5, i, j ) * rsd( 5, i, j-1, k ) 49 + c( m, 5, i, j ) * rsd( 5, i-1, j, k ) ); 50 } 51 52//c diagonal block inversion 53//c 54//c forward elimination 55 56 for (m=1; m<=5; m++) { 57 tmat( m, 1 ) = d( m, 1, i, j ); 58 tmat( m, 2 ) = d( m, 2, i, j ); 59 tmat( m, 3 ) = d( m, 3, i, j ); 60 tmat( m, 4 ) = d( m, 4, i, j ); 61 tmat( m, 5 ) = d( m, 5, i, j ); 62 } 63 64 tmp1 = 1.0 / tmat( 1, 1 ); 65 tmp = tmp1 * tmat( 2, 1 ); 66 tmat( 2, 2 ) = tmat( 2, 2 ) 67 - tmp * tmat( 1, 2 ); 68 tmat( 2, 3 ) = tmat( 2, 3 ) 69 - tmp * tmat( 1, 3 ); 70 tmat( 2, 4 ) = tmat( 2, 4 ) 71 - tmp * tmat( 1, 4 ); 72 tmat( 2, 5 ) = tmat( 2, 5 ) 73 - tmp * tmat( 1, 5 ); 74 rsd( 2, i, j, k ) = rsd( 2, i, j, k ) 75 - rsd( 1, i, j, k ) * tmp; 76 77 tmp = tmp1 * tmat( 3, 1 ); 78 tmat( 3, 2 ) = tmat( 3, 2 ) 79 - tmp * tmat( 1, 2 ); 80 tmat( 3, 3 ) = tmat( 3, 3 ) 81 - tmp * tmat( 1, 3 ); 82 tmat( 3, 4 ) = tmat( 3, 4 ) 83 - tmp * tmat( 1, 4 ); 84 tmat( 3, 5 ) = tmat( 3, 5 ) 85 - tmp * tmat( 1, 5 ); 86 rsd( 3, i, j, k ) = rsd( 3, i, j, k ) 87 - rsd( 1, i, j, k ) * tmp; 88 89 tmp = tmp1 * tmat( 4, 1 ); 90 tmat( 4, 2 ) = tmat( 4, 2 ) 91 - tmp * tmat( 1, 2 ); 92 tmat( 4, 3 ) = tmat( 4, 3 ) 93 - tmp * tmat( 1, 3 ); 94 tmat( 4, 4 ) = tmat( 4, 4 ) 95 - tmp * tmat( 1, 4 ); 96 tmat( 4, 5 ) = tmat( 4, 5 ) 97 - tmp * tmat( 1, 5 ); 98 rsd( 4, i, j, k ) = rsd( 4, i, j, k ) 99 - rsd( 1, i, j, k ) * tmp; 100 101 tmp = tmp1 * tmat( 5, 1 ); 102 tmat( 5, 2 ) = tmat( 5, 2 ) 103 - tmp * tmat( 1, 2 ); 104 tmat( 5, 3 ) = tmat( 5, 3 ) 105 - tmp * tmat( 1, 3 ); 106 tmat( 5, 4 ) = tmat( 5, 4 ) 107 - tmp * tmat( 1, 4 ); 108 tmat( 5, 5 ) = tmat( 5, 5 ) 109 - tmp * tmat( 1, 5 ); 110 rsd( 5, i, j, k ) = rsd( 5, i, j, k ) 111 - rsd( 1, i, j, k ) * tmp; 112 113 tmp1 = 1.0 / tmat( 2, 2 ); 114 tmp = tmp1 * tmat( 3, 2 ); 115 tmat( 3, 3 ) = tmat( 3, 3 ) 116 - tmp * tmat( 2, 3 ); 117 tmat( 3, 4 ) = tmat( 3, 4 ) 118 - tmp * tmat( 2, 4 ); 119 tmat( 3, 5 ) = tmat( 3, 5 ) 120 - tmp * tmat( 2, 5 ); 121 rsd( 3, i, j, k ) = rsd( 3, i, j, k ) 122 - rsd( 2, i, j, k ) * tmp; 123 124 tmp = tmp1 * tmat( 4, 2 ); 125 tmat( 4, 3 ) = tmat( 4, 3 ) 126 - tmp * tmat( 2, 3 ); 127 tmat( 4, 4 ) = tmat( 4, 4 ) 128 - tmp * tmat( 2, 4 ); 129 tmat( 4, 5 ) = tmat( 4, 5 ) 130 - tmp * tmat( 2, 5 ); 131 rsd( 4, i, j, k ) = rsd( 4, i, j, k ) 132 - rsd( 2, i, j, k ) * tmp; 133 134 tmp = tmp1 * tmat( 5, 2 ); 135 tmat( 5, 3 ) = tmat( 5, 3 ) 136 - tmp * tmat( 2, 3 ); 137 tmat( 5, 4 ) = tmat( 5, 4 ) 138 - tmp * tmat( 2, 4 ); 139 tmat( 5, 5 ) = tmat( 5, 5 ) 140 - tmp * tmat( 2, 5 ); 141 rsd( 5, i, j, k ) = rsd( 5, i, j, k ) 142 - rsd( 2, i, j, k ) * tmp; 143 144 tmp1 = 1.0 / tmat( 3, 3 ); 145 tmp = tmp1 * tmat( 4, 3 ); 146 tmat( 4, 4 ) = tmat( 4, 4 ) 147 - tmp * tmat( 3, 4 ); 148 tmat( 4, 5 ) = tmat( 4, 5 ) 149 - tmp * tmat( 3, 5 ); 150 rsd( 4, i, j, k ) = rsd( 4, i, j, k ) 151 - rsd( 3, i, j, k ) * tmp; 152 153 tmp = tmp1 * tmat( 5, 3 ); 154 tmat( 5, 4 ) = tmat( 5, 4 ) 155 - tmp * tmat( 3, 4 ); 156 tmat( 5, 5 ) = tmat( 5, 5 ) 157 - tmp * tmat( 3, 5 ); 158 rsd( 5, i, j, k ) = rsd( 5, i, j, k ) 159 - rsd( 3, i, j, k ) * tmp; 160 161 tmp1 = 1.0 / tmat( 4, 4 ); 162 tmp = tmp1 * tmat( 5, 4 ); 163 tmat( 5, 5 ) = tmat( 5, 5 ) 164 - tmp * tmat( 4, 5 ); 165 rsd( 5, i, j, k ) = rsd( 5, i, j, k ) 166 - rsd( 4, i, j, k ) * tmp; 167 168//c back substitution 169 170 rsd( 5, i, j, k ) = rsd( 5, i, j, k ) 171 / tmat( 5, 5 ); 172 173 rsd( 4, i, j, k ) = rsd( 4, i, j, k ) 174 - tmat( 4, 5 ) * rsd( 5, i, j, k ); 175 rsd( 4, i, j, k ) = rsd( 4, i, j, k ) 176 / tmat( 4, 4 ); 177 178 rsd( 3, i, j, k ) = rsd( 3, i, j, k ) 179 - tmat( 3, 4 ) * rsd( 4, i, j, k ) 180 - tmat( 3, 5 ) * rsd( 5, i, j, k ); 181 rsd( 3, i, j, k ) = rsd( 3, i, j, k ) 182 / tmat( 3, 3 ); 183 184 rsd( 2, i, j, k ) = rsd( 2, i, j, k ) 185 - tmat( 2, 3 ) * rsd( 3, i, j, k ) 186 - tmat( 2, 4 ) * rsd( 4, i, j, k ) 187 - tmat( 2, 5 ) * rsd( 5, i, j, k ); 188 rsd( 2, i, j, k ) = rsd( 2, i, j, k ) 189 / tmat( 2, 2 ); 190 191 rsd( 1, i, j, k ) = rsd( 1, i, j, k ) 192 - tmat( 1, 2 ) * rsd( 2, i, j, k ) 193 - tmat( 1, 3 ) * rsd( 3, i, j, k ) 194 - tmat( 1, 4 ) * rsd( 4, i, j, k ) 195 - tmat( 1, 5 ) * rsd( 5, i, j, k ); 196 rsd( 1, i, j, k ) = rsd( 1, i, j, k ) 197 / tmat( 1, 1 ); 198 199 } 200 } 201 202//c send data to east and south 203 204 iex = 2; 205 exchange_1(rsd, k, iex); 206 207 return; 208} 209 210 211 212