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

Linux Cross Reference
Tina4/src/math/numeric/pentadiag.c

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 ] ~

This page was automatically generated by the LXR engine.
Visit the LXR main site for more information.