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

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

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

  1 /**********
  2  *
  3  * This file is part of the TINA Open Source Image Analysis Environment
  4  * henceforth known as TINA
  5  *
  6  * TINA is free software; you can redistribute it and/or modify
  7  * it under the terms of the GNU General Public License as
  8  * published by the Free Software Foundation.
  9  *
 10  * TINA is distributed in the hope that it will be useful,
 11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
 12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 13  * GNU General Public License for more details.
 14  *
 15  * You should have received a copy of the GNU General Public License
 16  * along with TINA; if not, write to the Free Software Foundation, Inc.,
 17  * 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 18  *
 19  * ANY users of TINA who require exemption from the existing licence must
 20  * negotiate a new licence with Dr. Neil.A.Thacker, the sole agent for
 21  * the University of Manchester.
 22  *
 23  **********
 24  *
 25  * Program :    TINA
 26  * File    :  $Source: /home/tina/cvs/tina-libs/tina/math/mathNum_solve.c,v $
 27  * Date    :  $Date: 2005/01/23 19:10:21 $
 28  * Version :  $Revision: 1.6 $
 29  * CVS Id  :  $Id: mathNum_solve.c,v 1.6 2005/01/23 19:10:21 paul Exp $
 30  *
 31  * Author  :  Legacy TINA
 32  *
 33  * Notes :
 34  *
 35  *********
 36 */
 37 /** 
 38  *  @file
 39  *  @brief  LU decomposition and backsubstitution (for a Matrix).       
 40  *
 41  *  See Numerical Recipes (3rd Ed) Ch 2.3 pp 43-50 for more details regarding the method.
 42  * 
 43 */
 44 
 45 #include "mathNum_solve.h"
 46 
 47 #if HAVE_CONFIG_H
 48 #include <config.h>
 49 #endif
 50 
 51 #include <math.h>
 52 #include <tina/sys/sysDef.h>
 53 #include <tina/sys/sysPro.h>
 54 #include <tina/math/math_MatrDef.h>
 55 #include <tina/math/math_VecDef.h>
 56 #include <tina/math/math_VecPro.h>
 57 
 58 #define TINY 1.0e-20;
 59 
 60 void            ludcmp(double **a, int n, int *indx, double *d)
 61 {
 62         int             i, imax=0, j, k;
 63         double          big, dum, sum, temp;
 64         double         *vv;
 65 
 66         vv = dvector_alloc(0, n);
 67         *d = 1.0;
 68         for (i = 0; i < n; i++)
 69         {
 70                 big = 0.0;
 71                 for (j = 0; j < n; j++)
 72                         if ((temp = fabs(a[i][j])) > big)
 73                                 big = temp;
 74                 if (big == 0.0)
 75                         big = TINY;
 76                 vv[i] = 1.0 / big;
 77         }
 78         for (j = 0; j < n; j++)
 79         {
 80                 for (i = 0; i < j; i++)
 81                 {
 82                         sum = a[i][j];
 83                         for (k = 0; k < i; k++)
 84                                 sum -= a[i][k] * a[k][j];
 85                         a[i][j] = sum;
 86                 }
 87                 big = 0.0;
 88                 for (i = j; i < n; i++)
 89                 {
 90                         sum = a[i][j];
 91                         for (k = 0; k < j; k++)
 92                                 sum -= a[i][k] * a[k][j];
 93                         a[i][j] = sum;
 94                         if ((dum = vv[i] * fabs(sum)) >= big)
 95                         {
 96                                 big = dum;
 97                                 imax = i;
 98                         }
 99                 }
100                 if (j != imax)
101                 {
102                         for (k = 0; k < n; k++)
103                         {
104                                 dum = a[imax][k];
105                                 a[imax][k] = a[j][k];
106                                 a[j][k] = dum;
107                         }
108                         *d = -(*d);
109                         vv[imax] = vv[j];
110                 }
111                 indx[j] = imax;
112                 if (a[j][j] == 0.0)
113                         a[j][j] = TINY;
114                 if (j != n - 1)
115                 {
116                         dum = 1.0 / (a[j][j]);
117                         for (i = j + 1; i < n; i++)
118                                 a[i][j] *= dum;
119                 }
120         }
121         dvector_free(vv, 0);
122 }
123 
124 #undef TINY
125 
126 
127 void            lubksb(double **a, int n, int *indx, double *b)
128 {
129         int             i, ii = -1, ip, j;
130         double          sum;
131 
132         for (i = 0; i < n; i++)
133         {
134                 ip = indx[i];
135                 sum = b[ip];
136                 b[ip] = b[i];
137                 if (ii >= 0)
138                         for (j = ii; j <= i - 1; j++)
139                                 sum -= a[i][j] * b[j];
140                 else if (sum)
141                         ii = i;
142                 b[i] = sum;
143         }
144         for (i = n - 1; i >= 0; i--)
145         {
146                 sum = b[i];
147                 for (j = i + 1; j < n; j++)
148                         sum -= a[i][j] * b[j];
149                 b[i] = sum / a[i][i];
150         }
151 }
152 
153 Vector         *matrix_solve(Matrix * mat, Vector * y)
154 {
155         Vector         *x;
156         double        **a, *b, d;
157         int            *index, n;
158 
159         if (mat == NULL || y == NULL)
160                 return (NULL);
161         if (mat->m != y->n || mat->n != y->n)
162                 return (NULL);
163 
164         n = mat->m;
165         a = darray_dcopy(mat->el.double_v, 0, 0, n, n);
166         x = vector_copy(y);
167         b = x->data;
168         index = ivector_alloc(0, n);
169 
170         ludcmp(a, n, index, &d);
171         lubksb(a, n, index, b);
172 
173         darray_free(a, 0, 0, n, n);
174         ivector_free(index, 0);
175 
176         return (x);
177 }
178 

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