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