| 245 | } |
| 246 | |
| 247 | void gbsl(real **a, integer n, integer smu, integer ml, integer *p, |
| 248 | real *b) |
| 249 | { |
| 250 | integer k, l, i, first_row_k, last_row_k; |
| 251 | real mult, *diag_k; |
| 252 | |
| 253 | /* Solve Ly = Pb, store solution y in b */ |
| 254 | |
| 255 | for (k=0; k < n-1; k++) { |
| 256 | l = p[k]; |
| 257 | mult = b[l]; |
| 258 | if (l != k) { |
| 259 | b[l] = b[k]; |
| 260 | b[k] = mult; |
| 261 | } |
| 262 | diag_k = a[k]+smu; |
| 263 | last_row_k = MIN(n-1,k+ml); |
| 264 | for (i=k+1; i <= last_row_k; i++) |
| 265 | b[i] += mult * diag_k[i-k]; |
| 266 | } |
| 267 | |
| 268 | /* Solve Ux = y, store solution x in b */ |
| 269 | |
| 270 | for (k=n-1; k >= 0; k--) { |
| 271 | diag_k = a[k]+smu; |
| 272 | first_row_k = MAX(0,k-smu); |
| 273 | b[k] /= (*diag_k); |
| 274 | mult = -b[k]; |
| 275 | for (i=first_row_k; i <= k-1; i++) |
| 276 | b[i] += mult*diag_k[i-k]; |
| 277 | } |
| 278 | } |
| 279 | |
| 280 | void bandzero(real **a, integer n, integer mu, integer ml, integer smu) |
| 281 | { |