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

Linux Cross Reference
Tina6/tina-libs/tina/math/mathNum_pentadiag.c

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

  1 /**********
  2  *
  3  * Copyright (c) 2003, Division of Imaging Science and Biomedical Engineering,
  4  * University of Manchester, UK.  All rights reserved.
  5  * 
  6  * Redistribution and use in source and binary forms, with or without modification, 
  7  * are permitted provided that the following conditions are met:
  8  * 
  9  *   . Redistributions of source code must retain the above copyright notice, 
 10  *     this list of conditions and the following disclaimer.
 11  *    
 12  *   . Redistributions in binary form must reproduce the above copyright notice,
 13  *     this list of conditions and the following disclaimer in the documentation 
 14  *     and/or other materials provided with the distribution.
 15  * 
 16  *   . Neither the name of the University of Manchester nor the names of its
 17  *     contributors may be used to endorse or promote products derived from this 
 18  *     software without specific prior written permission.
 19  * 
 20  * 
 21  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
 22  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
 23  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
 24  * ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE 
 25  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
 26  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
 27  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 28  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
 29  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 
 30  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 
 31  * POSSIBILITY OF SUCH DAMAGE.
 32  *
 33  **********
 34  *
 35  * Program :    TINA
 36  * File    :  $Source: /home/tina/cvs/tina-libs/tina/math/mathNum_pentadiag.c,v $
 37  * Date    :  $Date: 2005/01/23 19:10:21 $
 38  * Version :  $Revision: 1.5 $
 39  * CVS Id  :  $Id: mathNum_pentadiag.c,v 1.5 2005/01/23 19:10:21 paul Exp $
 40  *
 41  * Author  :  Legacy TINA
 42  *
 43  * Notes : Solve cyclic pentadiagonal system [abcde]x = s
 44  *
 45  *********
 46 */
 47 /** 
 48  *  @file
 49  *  @brief  Function for solving cyclic pentadiagonal system [abcde]x = s.      
 50  *
 51  *
 52  * 
 53 */
 54 
 55 
 56 #include "mathNum_pentadiag.h"
 57 
 58 #if HAVE_CONFIG_H
 59 #include <config.h>
 60 #endif
 61 
 62 #include <tina/sys/sysDef.h>
 63 #include <tina/sys/sysPro.h>
 64 
 65 /**
 66 solves cyclic pentadiagonal system [abcde]x = s;
 67  - puts answer in s
 68 Originally written in Pascal, so index shift of 1.
 69 **/
 70 
 71 void            pentadiag(int n, double *a0, double *b0, double *c0, double *d0, double *e0, double *s0)
 72 {
 73         int             i, i1, i2;
 74         int             n1 = n - 1, n2 = n - 2, n3 = n - 3, n4 = n - 4;
 75         double          v, g, h;
 76         double         *a = a0 - 1, *b = b0 - 1, *c = c0 - 1, *d = d0 - 1;
 77         double         *e = e0 - 1, *s = s0 - 1;
 78 
 79         g = a[2];
 80         a[2] = 0.0;
 81         h = 0.0;
 82         c[1] = 1.0 / c[1];
 83 
 84         for (i = 2; i <= n - 4; i++)
 85         {
 86                 i1 = i - 1;
 87                 i2 = i + 1;
 88                 v = b[i] * c[i1];
 89                 c[i] -= v * d[i1];
 90                 d[i] -= v * e[i1];
 91                 b[i] = g - v * b[i1];
 92                 a[i] -= v * a[i1];
 93                 s[i] -= v * s[i1];
 94                 v = a[i2] * c[i1];
 95                 b[i2] -= v * d[i1];
 96                 c[i2] -= v * e[i1];
 97                 a[i2] = -v * a[i1];
 98                 g = -v * b[i1];
 99                 s[i2] -= v * s[i1];
100                 v = e[n1] * c[i1];
101                 c[n1] -= v * a[i1];
102                 d[n1] -= v * b[i1];
103                 e[n1] = h - v * d[i1];
104                 h = -v * e[i1];
105                 s[n1] -= v * s[i1];
106                 v = d[n] * c[i1];
107                 d[n] = e[n] - v * d[i1];
108                 e[n] = -v * e[i1];
109                 b[n] -= v * a[i1];
110                 c[n] -= v * b[i1];
111                 s[n] -= v * s[i1];
112                 c[i] = 1.0 / c[i];
113         }
114 
115         e[n3] += a[n3];
116         a[n3] = g;
117         a[n1] = a[n1] + h;
118         i = n - 3;
119         v = b[n3] * c[n4];
120         c[n3] -= v * d[n4];
121         s[n3] -= v * s[n4];
122         d[n3] -= v * e[n4];
123         e[n3] -= v * a[n4];
124         a[n3] -= v * b[n4];
125         v = a[n2] * c[n4];
126         b[n2] -= v * d[n4];
127         c[n2] -= v * e[n4];
128         s[n2] -= v * s[n4];
129         d[n2] -= v * a[n4];
130         e[n2] -= v * b[n4];
131         v = e[n1] * c[n4];
132         a[n1] -= v * d[n4];
133         b[n1] -= v * e[n4];
134         s[n1] -= v * s[n4];
135         c[n1] -= v * a[n4];
136         d[n1] -= v * b[n4];
137         v = d[n] * c[n4];
138         e[n] -= v * d[n4];
139         a[n] -= v * e[n4];
140         b[n] -= v * a[n4];
141         c[n] -= v * b[n4];
142         s[n] -= v * s[n4];
143         c[n3] = 1.0 / c[n3];
144         i = n - 2;
145         v = b[n2] * c[n3];
146         c[n2] -= v * d[n3];
147         s[n2] -= v * s[n3];
148         d[n2] -= v * e[n3];
149         e[n2] -= v * a[n3];
150         v = a[n1] * c[n3];
151         b[n1] -= v * d[n3];
152         c[n1] -= v * e[n3];
153         s[n1] -= v * s[n3];
154         d[n1] -= v * a[n3];
155         v = e[n] * c[n3];
156         a[n] -= v * d[n3];
157         b[n] -= v * e[n3];
158         s[n] -= v * s[n3];
159         c[n] -= v * a[n3];
160         c[n2] = 1.0 / c[n2];
161         i = n - 1;
162         v = b[n1] * c[n2];
163         c[n1] -= v * d[n2];
164         s[n1] -= v * s[n2];
165         d[n1] -= v * e[n2];
166         v = a[n] * c[n2];
167         b[n] -= v * d[n2];
168         c[n] -= v * e[n2];
169         s[n] -= v * s[n2];
170         c[n1] = 1.0 / c[n1];
171         i = n;
172         v = b[n] * c[n1];
173         c[n] -= v * d[n1];
174         s[n] -= v * s[n1];
175         c[n] = 1.0 / c[n];
176         s[n] *= c[n];
177         s[n1] = (s[n1] - d[n1] * s[n]) * c[n1];
178         s[n2] = (s[n2] - d[n2] * s[n1] - e[n2] * s[n]) * c[n2];
179         s[n3] = (s[n3] - d[n3] * s[n2] - e[n3] * s[n1] - a[n3] * s[n]) * c[n3];
180 
181         for (i = n - 4; i >= 1; i--)
182                 s[i] = (s[i] - d[i] * s[i + 1] - e[i] * s[i + 2] - a[i] * s[n1] - b[i] * s[n]) * c[i];
183 }
184 
185 /**
186 Forms product of cyclic pentadiagonal matrix and vector:
187 [abcde]x = y
188 **/
189 
190 void            pentaprod(int n, double *a, double *b, double *c, double *d, double *e, double *x, double *y)
191 {
192         int             i, n1 = n - 1, n2 = n - 2, n3 = n - 3, n4 = n - 4;
193 
194         y[0] = c[0] * x[0] + d[0] * x[1] + e[0] * x[2] + a[0] * x[n2] + b[0] * x[n1];
195         y[1] = b[1] * x[0] + c[1] * x[1] + d[1] * x[2] + e[1] * x[3] + a[1] * x[n1];
196         for (i = 2; i < n2; i++)
197                 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];
198         y[n2] = e[n2] * x[0] + a[n2] * x[n4] + b[n2] * x[n3] + c[n2] * x[n2] + d[n2] * x[n1];
199         y[n1] = d[n1] * x[0] + e[n1] * x[1] + a[n1] * x[n3] + b[n1] * x[n2] + c[n1] * x[n1];
200 }
201 
202 /** temporary solution of cyclic tridiagonal **/
203 
204 void            triadiag(int n, double *b0, double *c0, double *d0, double *s0)
205 {
206         double         *a0 = dvector_alloc(0, n);
207         double         *e0 = dvector_alloc(0, n);
208 
209         pentadiag(n, a0, b0, c0, d0, e0, s0);
210         dvector_free(a0, 0);
211         dvector_free(e0, 0);
212 }
213 

~ [ 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.