lu.h

Go to the documentation of this file.
00001 // -*- c-basic-offset: 4 -*-
00002 
00029 /* The LU decomposition function was taken from the rocketworkbench project:
00030  * 
00031  * lu.c  -  PA = LU factorisation with pivoting
00032  * $Id$
00033  * Copyright (C) 2000
00034  *    Antoine Lefebvre <antoine.lefebvre@polymtl.ca>
00035  *
00036  *
00037  * Licensed under the GPLv2
00038  */
00039 
00040 
00041 #ifndef _HUGIN_MATH_LU_H
00042 #define _HUGIN_MATH_LU_H
00043 
00044 
00045 /* NOTE on matrix representation
00046  * -----------------------------
00047  * All matrix should be allocate the following way
00048  *
00049  * matrix = (double *) malloc (sizeof(double) * line * column)
00050  *
00051  * to access the value at line 2, column 4, you do it like that
00052  *
00053  * matrix[2 + 4*line]
00054  */
00055 
00056 
00057 /* Find the solution of linear system of equation using the
00058  * LU factorisation method as explain in
00059  * 'Advanced engineering mathematics' bye Erwin Kreyszig.
00060  *
00061  * This algorithm will also do column permutation to find
00062  * the larger pivot.
00063  *
00064  * ARGUMENTS
00065  * ---------
00066  * matrix: the augmented matrix of coefficient in the system
00067  *         with right hand side value.
00068  *
00069  * solution: the solution vector
00070  *
00071  * neq: number of equation in the system
00072  *
00073  * Antoine Lefebvre
00074  *    february 6, 2000 Initial version
00075  *    october 20, 2000 revision of the permutation method
00076  */
00077 
00078 extern "C"
00079 {
00080 int math_lu_solve(double *matrix, double *solution, int neq);
00081 }
00082 
00083 
00084 namespace hugin_utils
00085 {
00086 
00088     class LMS_Solver
00089     {
00090     public:
00095         LMS_Solver(unsigned nEq)
00096         {
00097             m_nEq = nEq;
00098             m_AtA = new double[nEq*(nEq+1)];
00099             for (unsigned i=0; i < nEq*(nEq+1); i++) m_AtA[i] = 0;
00100         }
00101 
00102         ~LMS_Solver()
00103         {
00104             delete[] m_AtA;
00105         }
00106 
00108         template <class Iter>
00109         void addRow(Iter Arow, double b)
00110         {
00111             for( unsigned i=0; i<m_nEq; ++i)
00112             {
00113                 // calculate  Atb
00114                 m_AtA[i + m_nEq*m_nEq]+=Arow[i]*b;
00115                 for( unsigned j=0; j<m_nEq; ++j)
00116                 {
00117                     m_AtA[i + j*m_nEq] += Arow[i]*Arow[j];
00118                 }
00119             }
00120         }
00121 
00123         template <class Vector>
00124         bool solve(Vector & x)
00125         {
00126             double * solution = new double[m_nEq];
00127             bool ret = math_lu_solve(m_AtA, solution, m_nEq) != 0;
00128             for (unsigned i=0; i < m_nEq; i++) {
00129                 x[i] = solution[i];
00130             }
00131             delete[] solution;
00132             return ret;
00133         }
00134 
00135     protected:
00136         unsigned m_nEq;
00137         // matrix that contains [AtA  Atb]
00138         // uses fortran array access conventions
00139         double * m_AtA;
00140     };
00141 
00142 } // namespace
00143 
00144 #endif // _H
00145 
00146 
00147 

Generated on Mon Sep 1 01:25:39 2014 for Hugintrunk by  doxygen 1.3.9.1