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

Linux Cross Reference
Tina6/tina-libs/tina/math/mathMatv_chol.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_chol.c,v $
 27  * Date    :  $Date: 2003/10/06 12:29:47 $
 28  * Version :  $Revision: 1.5 $
 29  * CVS Id  :  $Id: mathMatv_chol.c,v 1.5 2003/10/06 12:29:47 neil Exp $
 30  *
 31  * Author  :  Legacy TINA
 32  *
 33  * Notes : cholesky.c Utilities for square root decomposition of a symmetric matrix.
 34  *
 35  *********
 36 */
 37 /** 
 38  *  @file
 39  *  @brief Utilities for square root decomposition of a symmetric matrix .      
 40  *
 41  *  Includes Cholesky and square-root free Cholesky decompositions. See Numerical Recipes
 42  *  (3rd ed) pp96 - 98 for more details of the method.
 43  */
 44 
 45 #include "mathMatv_chol.h"
 46 
 47 #if HAVE_CONFIG_H
 48 #include <config.h>
 49 #endif
 50 
 51 #include <stdio.h>
 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 
 58 /*
 59 Cholesky decomposition of n by n matrix s (s = L Lt).
 60 Overwrites L into lower triangle of s.
 61 */
 62 Bool            cholesky(int n, double **s)
 63 {
 64         int             i, j, k;
 65 
 66         for (k = 0; k < n; k++)
 67         {
 68                 if (s[k][k] < 0.0)
 69                         return (false);
 70                 if (s[k][k] < FLT_MIN)
 71                         s[k][k] = FLT_MIN;
 72                 s[k][k] = sqrt(s[k][k]);
 73                 for (i = k + 1; i < n; i++)
 74                 {
 75                         s[i][k] /= s[k][k];
 76                         for (j = k + 1; j <= i; j++)
 77                                 s[i][j] -= s[i][k] * s[j][k];
 78                 }
 79         }
 80         for (i = 0; i < n; i++)
 81                 for (j = i + 1; j < n; j++)
 82                         s[i][j] = 0.0;
 83         return (true);
 84 }
 85 
 86 void            mat_cholesky(Mat * s)
 87 {
 88         cholesky(s->n, s->el);
 89 }
 90 
 91 /*
 92 Solves ly=x in place. l n by n lower triangular matrix. Solution
 93 y overwrites x.
 94 */
 95 void            lower_solve(int n, double **l, double *x)
 96 {
 97         int             i, j;
 98 
 99         for (i = 0; i < n; i++)
100         {
101                 for (j = 0; j < i; j++)
102                         x[i] -= l[i][j] * x[j];
103                 x[i] /= l[i][i];
104         }
105 }
106 
107 /*
108 Square root free Cholesky decomposition s = l d lt.
109 l lower triangular, d stored on diagonal of l.
110 */
111 Bool            cholesky_ldl(int n, double **s, double **l, double *d)
112 {
113         int             i, j, k;
114 
115         for (i = 0; i < n; i++)
116         {
117                 for (j = 0; j <= i; j++)
118                         l[i][j] = s[i][j];
119                 for (j = i + 1; j < n; j++)
120                         l[i][j] = 0;
121         }
122         for (k = 0; k < n; k++)
123         {
124                 d[k] = l[k][k];
125                 l[k][k] = 1;
126                 for (i = k + 1; i < n; i++)
127                 {
128                         l[i][k] /= d[k];
129                         for (j = k + 1; j <= i; j++)
130                                 l[i][j] -= d[k] * l[i][k] * l[j][k];
131                 }
132         }
133         return (true);
134 }
135 
136 void            mat_cholesky_ldl(Mat * s, Mat * l, Vec * d)
137 {
138         cholesky_ldl(s->n, s->el, l->el, d->el);
139 }
140 

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