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

Linux Cross Reference
Tina6/tina-libs/tina/math/mathNum_cholesky.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_cholesky.c,v $
 27  * Date    :  $Date: 2003/10/06 12:29:47 $
 28  * Version :  $Revision: 1.5 $
 29  * CVS Id  :  $Id: mathNum_cholesky.c,v 1.5 2003/10/06 12:29:47 neil Exp $
 30  *
 31  * Author  :  Legacy TINA
 32  *
 33  * Notes : Cholesky
 34  *
 35  *********
 36 */
 37 /** 
 38  *  @file
 39  *  @brief  Cholesky decomposition of a square matrix.  
 40  *
 41  *  Functions include;
 42  *
 43  *  - cholesky decomposition and backsubstitution for an array.
 44  *  - Higher level Matrix calls to these functions.
 45  *
 46  *  See Numerical Recipes (3rd ed) Ch2.9 pp 96-98 for more details regarding the method.
 47 */
 48 
 49 #include "mathNum_cholesky.h"
 50 
 51 #if HAVE_CONFIG_H
 52 #include <config.h>
 53 #endif
 54 
 55 #include <math.h>
 56 #include <tina/sys/sysDef.h>
 57 #include <tina/sys/sysPro.h>
 58 #include <tina/math/math_MatrPro.h>
 59 #include <tina/math/math_MatrDef.h>
 60 #include <tina/math/math_VecDef.h>
 61 #include <tina/math/math_VecPro.h>
 62 
 63 #define TINY 1.0e-20
 64 
 65 static int      cholesky1(double **A, int n)
 66 {
 67         int             i, j, k;
 68 
 69         for (i = 0; i < n; ++i)
 70         {
 71                 for (j = 0; j <= i; ++j)
 72                 {
 73                         double          sum = A[i][j];
 74 
 75                         A[j][i] = 0;
 76 
 77                         for (k = 0; k < j; ++k)
 78                                 sum -= A[i][k] * A[j][k];       /* computed previously */
 79 
 80                         if (i == j)
 81                         {
 82                                 if (sum < 0)
 83                                         return (0);
 84                                 sum = sqrt(sum);
 85                                 if (fabs(sum) < TINY)
 86                                         return (0);
 87                                 A[i][j] = sum;
 88                         } else
 89                                 A[i][j] = sum / A[j][j];
 90                 }
 91         }
 92         return (1);
 93 }
 94 
 95 #undef TINY
 96 
 97 static void     cholbksb(double **A, int n, double *x, double *b)
 98 {
 99         int             i, j;
100         double          sum;
101 
102         for (i = 0; i < n; i++)
103         {
104                 sum = b[i];
105                 for (j = 0; j < i; ++j)
106                         sum -= x[j] * A[i][j];
107                 x[i] = sum / A[i][i];
108         }
109 
110         for (i = n - 1; i >= 0; i--)
111         {
112                 sum = x[i];
113                 for (j = i + 1; j < n; j++)
114                         sum -= x[j] * A[j][i];
115                 x[i] = sum / A[i][i];
116         }
117 }
118 
119 Matrix         *matrix_cholesky_decomp(Matrix * A)
120 {
121         if (A == NULL || A->n != A->m)
122                 return (NULL);
123 
124         A = matrix_cast_fill(A, double_v);
125         if (!cholesky1(A->el.double_v, A->n))
126         {
127                 matrix_free((Matrix *) A);
128                 return (NULL);
129         }
130         return (A);
131 }
132 
133 Vector         *matrix_cholesky_back_sub(Matrix * A, Vector * b)
134 {
135         int             n;
136         Vector         *x = NULL;
137         Vector         *c = NULL;
138 
139         if (A == NULL || A->n != A->m)
140                 return (NULL);
141 
142         n = A->n;
143         if (b == NULL || b->n != n)
144                 return (NULL);
145 
146         A = matrix_cast_fill(A, double_v);
147         c = vector_cast(b, double_v);
148         x = vector_alloc(n, double_v);
149         cholbksb(A->el.double_v, n, (double *) x->data, (double *) c->data);
150         vector_free(c);
151         matrix_free((Matrix *) A);
152 
153         return (x);
154 }
155 
156 Vector         *matrix_cholesky_sol(Matrix * A, Vector * b)
157 {
158         int             n;
159         int             sucess;
160         Vector         *x = NULL;
161         Vector         *c = NULL;
162 
163         if (A == NULL || A->n != A->m)
164                 return (NULL);
165 
166         n = A->n;
167         if (b == NULL || b->n != n)
168                 return (NULL);
169 
170         A = matrix_cast_fill(A, double_v);
171         sucess = cholesky1(A->el.double_v, n);
172 
173         if (sucess)
174         {
175                 c = vector_cast(b, double_v);
176                 x = vector_alloc(n, double_v);
177                 cholbksb(A->el.double_v, n, (double *) x->data, (double *) c->data);
178                 vector_free(c);
179         }
180         matrix_free((Matrix *) A);
181         return (x);
182 }
183 
184 Vector         *matrix_cholesky_weighted_least_square(Matrix * A, Matrix * W, Vector * b)
185 {
186         Matrix         *At;
187         Matrix         *M;
188         Matrix         *temp;
189         Vector         *x;
190 
191         if (A == NULL || A->m != b->n || W == NULL || W->n != W->m || W->m != A->m)
192                 return (NULL);
193 
194         A = matrix_cast_fill(A, double_v);      /* to maintain precision */
195         At = matrix_transp(A);
196         x = matrix_vprod(W, b);
197         b = matrix_vprod(At, x);
198         vector_free(x);
199         temp = matrix_prod(W, A);
200         M = matrix_prod(At, temp);
201         matrix_free((Matrix *) temp);
202         x = matrix_cholesky_sol(M, b);
203         vector_free(b);
204         matrix_free((Matrix *) A);
205         matrix_free((Matrix *) At);
206         matrix_free((Matrix *) M);
207         return (x);
208 }
209 
210 Vector         *matrix_cholesky_least_square(Matrix * A, Vector * b)
211 {
212         Matrix         *At;
213         Matrix         *M;
214         Vector         *x;
215 
216         if (A == NULL || b == NULL || A->m != b->n)
217                 return (NULL);
218 
219         A = matrix_cast_fill(A, double_v);      /* to keep precision */
220         At = matrix_transp(A);
221         b = matrix_vprod(At, b);
222         M = matrix_prod(At, A);
223         x = matrix_cholesky_sol(M, b);
224         vector_free(b);
225         matrix_free((Matrix *) A);
226         matrix_free((Matrix *) At);
227         matrix_free((Matrix *) M);
228         return (x);
229 }
230 

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