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

Linux Cross Reference
Tina6/tina-libs/tina/math/mathMatv_lu.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/mathMatv_lu.c,v $
 27  * Date    :  $Date: 2005/01/09 17:49:25 $
 28  * Version :  $Revision: 1.6 $
 29  * CVS Id  :  $Id: mathMatv_lu.c,v 1.6 2005/01/09 17:49:25 paul Exp $
 30  *
 31  * Author  :  Legacy TINA
 32  *
 33  * Notes : Routines using LU decoposition of square matrices.
 34  * modifications of numerical recipies equivalents.
 35  *
 36  *********
 37 */
 38 /** 
 39  *  @file
 40  *  @brief Routines for/using LU decomposition of square matrices.      
 41  *
 42  *  Versions of Numerical Recipes equivalents (see Chapter 2.3 of 3rd ed).
 43  * 
 44 */
 45 
 46 #include "mathMatv_lu.h"
 47 
 48 #if HAVE_CONFIG_H
 49 #include <config.h>
 50 #endif
 51 
 52 #include <math.h>
 53 #include <float.h>
 54 #include <tina/sys/sysDef.h>
 55 #include <tina/sys/sysPro.h>
 56 #include <tina/math/math_MatvDef.h>
 57 #include <tina/math/mathMatv_mat.h>
 58 
 59 static void     lu_decomp(int n, double **a, int *indx, double *d)
 60 {
 61         int             i, imax=0, j, k;
 62         double          big, dum, sum, temp;
 63         double         *v = tvector_alloc(0, n, double);
 64         double          tiny = 1.0 / FLT_MAX;
 65 
 66         *d = 1.0;
 67         for (i = 0; i < n; i++)
 68         {
 69                 big = 0.0;
 70                 for (j = 0; j < n; j++)
 71                         if ((temp = fabs(a[i][j])) > big)
 72                                 big = temp;
 73                 if (big == 0.0)
 74                         big = tiny;
 75                 v[i] = 1.0 / big;
 76         }
 77         for (j = 0; j < n; j++)
 78         {
 79                 for (i = 0; i < j; i++)
 80                 {
 81                         sum = a[i][j];
 82                         for (k = 0; k < i; k++)
 83                                 sum -= a[i][k] * a[k][j];
 84                         a[i][j] = sum;
 85                 }
 86                 big = 0.0;
 87                 for (i = j; i < n; i++)
 88                 {
 89                         sum = a[i][j];
 90                         for (k = 0; k < j; k++)
 91                                 sum -= a[i][k] * a[k][j];
 92                         a[i][j] = sum;
 93                         if ((dum = v[i] * fabs(sum)) >= big)
 94                         {
 95                                 big = dum;
 96                                 imax = i;
 97                         }
 98                 }
 99                 if (j != imax)
100                 {
101                         for (k = 0; k < n; k++)
102                         {
103                                 dum = a[imax][k];
104                                 a[imax][k] = a[j][k];
105                                 a[j][k] = dum;
106                         }
107                         *d = -(*d);
108                         v[imax] = v[j];
109                 }
110                 indx[j] = imax;
111                 if (a[j][j] == 0.0)
112                         a[j][j] = tiny;
113                 if (j != n - 1)
114                 {
115                         dum = 1.0 / (a[j][j]);
116                         for (i = j + 1; i < n; i++)
117                                 a[i][j] *= dum;
118                 }
119         }
120         tvector_free(v, 0, double);
121 }
122 
123 static void     lu_backsub(int n, double **a, int *indx, double *b)
124 {
125         int             i, k = -1, ip, j;
126         double          sum;
127 
128         for (i = 0; i < n; i++)
129         {
130                 ip = indx[i];
131                 sum = b[ip];
132                 b[ip] = b[i];
133                 if (k >= 0)
134                         for (j = k; j <= i - 1; j++)
135                                 sum -= a[i][j] * b[j];
136                 else if (sum)
137                         k = i;
138                 b[i] = sum;
139         }
140         for (i = n - 1; i >= 0; i--)
141         {
142                 sum = b[i];
143                 for (j = i + 1; j < n; j++)
144                         sum -= a[i][j] * b[j];
145                 b[i] = sum / a[i][i];
146         }
147 }
148 
149 void            mat_solve(Mat * a, Vec * y)
150 {
151         double          d;
152         int             n = a->n;
153         int            *index = tvector_alloc(0, n, int);
154         if (a == NULL || y == NULL)
155                 return;
156         lu_decomp(n, a->el, index, &d);
157         lu_backsub(n, a->el, index, y->el);
158         tvector_free(index, 0, int);
159 }
160 
161 double          mat_det(Mat * a)
162 {
163         double          det, d;
164         int             i, n = a->n;
165         int            *index = tvector_alloc(0, n, int);
166         Mat            *a1 = mat_make(a->m, a->n);
167 
168         mat_copy(a1, a);
169         lu_decomp(n, a1->el, index, &d);
170         /* determinant set to d NAT 9/4/99 */
171         for (det = d, i = 0; i < n; i++)
172                 det *= a1->el[i][i];
173         tvector_free(index, 0, int);
174         mat_free(a1);
175         return (det);
176 }
177 
178 double          darray_det(double **a, int n)
179 {
180         double          det, d;
181         int             i;
182         int            *index = tvector_alloc(0, n, int);
183         double        **a1 = darray_copy((char **) a, 0, 0, n, n);
184 
185         lu_decomp(n, a1, index, &d);
186         /* determinant set to d NAT 9/4/99 */
187         for (det = d, i = 0; i < n; i++)
188                 det *= a1[i][i];
189         tvector_free(index, 0, int);
190         darray_free(a1, 0, 0, n, n);
191         return (det);
192 }
193 

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