~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

```  1 /**@(#)Solve cyclic pentadiagonal system [abcde]x = s;
2  */
3
4 #include <tina/sys.h>
5 #include <tina/sysfuncs.h>
6
7 /**
8 solves cyclic pentadiagonal system [abcde]x = s;
9  - puts answer in s
10 Originally written in Pascal, so index shift of 1.
11 **/
12
13 void    pentadiag(int n, double *a0, double *b0, double *c0, double *d0, double *e0, double *s0)
14 {
15     int     i, i1, i2;
16     int     n1 = n - 1, n2 = n - 2, n3 = n - 3, n4 = n - 4;
17     double  v, g, h;
18     double *a = a0 - 1, *b = b0 - 1, *c = c0 - 1, *d = d0 - 1;
19     double *e = e0 - 1, *s = s0 - 1;
20
21     g = a[2];
22     a[2] = 0.0;
23     h = 0.0;
24     c[1] = 1.0 / c[1];
25
26     for (i = 2; i <= n - 4; i++)
27     {
28         i1 = i - 1;
29         i2 = i + 1;
30         v = b[i] * c[i1];
31         c[i] -= v * d[i1];
32         d[i] -= v * e[i1];
33         b[i] = g - v * b[i1];
34         a[i] -= v * a[i1];
35         s[i] -= v * s[i1];
36         v = a[i2] * c[i1];
37         b[i2] -= v * d[i1];
38         c[i2] -= v * e[i1];
39         a[i2] = -v * a[i1];
40         g = -v * b[i1];
41         s[i2] -= v * s[i1];
42         v = e[n1] * c[i1];
43         c[n1] -= v * a[i1];
44         d[n1] -= v * b[i1];
45         e[n1] = h - v * d[i1];
46         h = -v * e[i1];
47         s[n1] -= v * s[i1];
48         v = d[n] * c[i1];
49         d[n] = e[n] - v * d[i1];
50         e[n] = -v * e[i1];
51         b[n] -= v * a[i1];
52         c[n] -= v * b[i1];
53         s[n] -= v * s[i1];
54         c[i] = 1.0 / c[i];
55     }
56
57     e[n3] += a[n3];
58     a[n3] = g;
59     a[n1] = a[n1] + h;
60     i = n - 3;
61     v = b[n3] * c[n4];
62     c[n3] -= v * d[n4];
63     s[n3] -= v * s[n4];
64     d[n3] -= v * e[n4];
65     e[n3] -= v * a[n4];
66     a[n3] -= v * b[n4];
67     v = a[n2] * c[n4];
68     b[n2] -= v * d[n4];
69     c[n2] -= v * e[n4];
70     s[n2] -= v * s[n4];
71     d[n2] -= v * a[n4];
72     e[n2] -= v * b[n4];
73     v = e[n1] * c[n4];
74     a[n1] -= v * d[n4];
75     b[n1] -= v * e[n4];
76     s[n1] -= v * s[n4];
77     c[n1] -= v * a[n4];
78     d[n1] -= v * b[n4];
79     v = d[n] * c[n4];
80     e[n] -= v * d[n4];
81     a[n] -= v * e[n4];
82     b[n] -= v * a[n4];
83     c[n] -= v * b[n4];
84     s[n] -= v * s[n4];
85     c[n3] = 1.0 / c[n3];
86     i = n - 2;
87     v = b[n2] * c[n3];
88     c[n2] -= v * d[n3];
89     s[n2] -= v * s[n3];
90     d[n2] -= v * e[n3];
91     e[n2] -= v * a[n3];
92     v = a[n1] * c[n3];
93     b[n1] -= v * d[n3];
94     c[n1] -= v * e[n3];
95     s[n1] -= v * s[n3];
96     d[n1] -= v * a[n3];
97     v = e[n] * c[n3];
98     a[n] -= v * d[n3];
99     b[n] -= v * e[n3];
100     s[n] -= v * s[n3];
101     c[n] -= v * a[n3];
102     c[n2] = 1.0 / c[n2];
103     i = n - 1;
104     v = b[n1] * c[n2];
105     c[n1] -= v * d[n2];
106     s[n1] -= v * s[n2];
107     d[n1] -= v * e[n2];
108     v = a[n] * c[n2];
109     b[n] -= v * d[n2];
110     c[n] -= v * e[n2];
111     s[n] -= v * s[n2];
112     c[n1] = 1.0 / c[n1];
113     i = n;
114     v = b[n] * c[n1];
115     c[n] -= v * d[n1];
116     s[n] -= v * s[n1];
117     c[n] = 1.0 / c[n];
118     s[n] *= c[n];
119     s[n1] = (s[n1] - d[n1] * s[n]) * c[n1];
120     s[n2] = (s[n2] - d[n2] * s[n1] - e[n2] * s[n]) * c[n2];
121     s[n3] = (s[n3] - d[n3] * s[n2] - e[n3] * s[n1] - a[n3] * s[n]) * c[n3];
122
123     for (i = n - 4; i >= 1; i--)
124         s[i] = (s[i] - d[i] * s[i + 1] - e[i] * s[i + 2] - a[i] * s[n1] - b[i] * s[n]) * c[i];
125 }
126
127 /**
128 Forms product of cyclic pentadiagonal matrix and vector:
129 [abcde]x = y
130 **/
131
132 void    pentaprod(int n, double *a, double *b, double *c, double *d, double *e, double *x, double *y)
133 {
134     int     i, n1 = n - 1, n2 = n - 2, n3 = n - 3, n4 = n - 4;
135
136     y[0] = c[0] * x[0] + d[0] * x[1] + e[0] * x[2] + a[0] * x[n2] + b[0] * x[n1];
137     y[1] = b[1] * x[0] + c[1] * x[1] + d[1] * x[2] + e[1] * x[3] + a[1] * x[n1];
138     for (i = 2; i < n2; i++)
139         y[i] = a[i] * x[i - 2] + b[i] * x[i - 1] + c[i] * x[i] + d[i] * x[i + 1] + e[i] * x[i + 2];
140     y[n2] = e[n2] * x[0] + a[n2] * x[n4] + b[n2] * x[n3] + c[n2] * x[n2] + d[n2] * x[n1];
141     y[n1] = d[n1] * x[0] + e[n1] * x[1] + a[n1] * x[n3] + b[n1] * x[n2] + c[n1] * x[n1];
142 }
143
144 /** temporary solution of cyclic tridiagonal **/
145
146 void    triadiag(int n, double *b0, double *c0, double *d0, double *s0)
147 {
148     double *a0 = dvector_alloc(0, n);
149     double *e0 = dvector_alloc(0, n);
150
151     pentadiag(n, a0, b0, c0, d0, e0, s0);
152     dvector_free((void *) a0, 0);
153     dvector_free((void *) e0, 0);
154 }
155 ```

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~