Canoe
Comprehensive Atmosphere N' Ocean Engine
band_back_sub.c
Go to the documentation of this file.
1 /*======================= band_back_sub() ====================================*/
2 
3 /*
4  * Back substitute for a banded matrix using the output from band_decomp().
5  * Adapted from Numerical Recipes in C, 2nd ed., p. 54.
6  * Assumes zero-based indexing.
7  */
8 
9 #undef SWAP
10 #define SWAP(a, b) \
11  { \
12  tmp = (a); \
13  (a) = (b); \
14  (b) = tmp; \
15  }
16 
17 #undef A
18 #define A(i, j) a[(m1 + m2 + 1) * (i) + (j)]
19 #undef AL
20 #define AL(i, j) al[m1 * (i) + (j)]
21 
22 void band_back_sub(int n, int m1, int m2, double *a, double *al, int *index,
23  double *b) {
24  int i, k, l, mm;
25  double tmp;
26  /*
27  * The following are part of DEBUG_MILESTONE(.) statements:
28  */
29  int idbms = 0;
30  static char dbmsname[] = "band_back_sub";
31 
32  mm = m1 + m2 + 1;
33  l = m1;
34  for (k = 0; k < n; k++) {
35  i = index[k];
36  if (i != k) {
37  SWAP(b[k], b[i]);
38  }
39  if (l < n) {
40  l++;
41  }
42  for (i = k + 1; i < l; i++) {
43  b[i] -= AL(k, i - k - 1) * b[k];
44  }
45  }
46  l = 1;
47  for (i = n - 1; i >= 0; i--) {
48  tmp = b[i];
49  for (k = 1; k < l; k++) {
50  tmp -= A(i, k) * b[k + i];
51  }
52  b[i] = tmp / A(i, 0);
53  if (l < mm) {
54  l++;
55  }
56  }
57 
58  return;
59 }
60 
61 /*======================= end of band_back_sub() =============================*/
void band_back_sub(int n, int m1, int m2, double *a, double *al, int *index, double *b)
Definition: band_back_sub.c:22
#define SWAP(a, b)
Definition: band_back_sub.c:10
#define A(i, j)
Definition: band_back_sub.c:18
#define AL(i, j)
Definition: band_back_sub.c:20