Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members  

guma_lineq.h

Go to the documentation of this file.
00001 /* LIBGUL - Geometry Utility Library
00002  * Copyright (C) 1998-1999 Norbert Irmer
00003  *
00004  * This library is free software; you can redistribute it and/or
00005  * modify it under the terms of the GNU Library General Public
00006  * License as published by the Free Software Foundation; either
00007  * version 2 of the License, or (at your option) any later version.
00008  *
00009  * This library is distributed in the hope that it will be useful,
00010  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00011  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00012  * Library General Public License for more details.
00013  *
00014  * You should have received a copy of the GNU Library General Public
00015  * License along with this library; if not, write to the
00016  * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
00017  * Boston, MA 02111-1307, USA.
00018  */
00019  
00020 #ifndef GUMA_LINEQ_H
00021 #define GUMA_LINEQ_H
00022 
00023 namespace guma {
00024 
00025 using gul::Ptr;
00026 using gul::rtr;
00027 using gul::Square;
00028 using gul::point;
00029 
00030 /*-----------------------------------------------------------------------*//**
00031   Decompose a matrix into lower and upper triangular matrix. the function
00032   throws an exception when an singular matrix is encountered (for
00033   convenience, so that not only the calling program can react to this,
00034   but also functions higher in the call stack)                              */
00035 /*--------------------------------------------------------------------------*/
00036 template< class T >
00037 GULAPI bool LUDecomposition( int nRowsCols, Ptr< Ptr<T> >& A,
00038                              int *perm_sign, Ptr<int>& index );
00039 
00040 /*---------------------------------------------------------------*//**
00041   Solve system of linear equations with help of LU-Decomposition
00042   (given in A).                                                     */
00043 /*------------------------------------------------------------------*/
00044 template< class T >
00045 GULAPI void LUBackSubstitution( int nRowsCols, Ptr<int>& index, 
00046                                 Ptr< Ptr<T> >& A, Ptr<T>& b );
00047 
00048 /*-------------------------------------------------------------------------*//**
00049   Multiply a band matrix and a vector: b = A*x, m1 = number of subdiagonal
00050   elements, m2 = number of superdiagonal elements, n = number of rows.
00051   It demonstrates the indexing scheme for band matrices                       */
00052 /*----------------------------------------------------------------------------*/
00053 template< class T >
00054 GULAPI void BandMVProduct( int n, int m1, int m2, 
00055                            Ptr< Ptr<T> > a, Ptr<T> x, Ptr<T> b);
00056 
00057 /*-----------------------------------------------------------------------*//**
00058   Decompose a band matrix into lower and upper triangular matrix: A = L*U.
00059   m1 = number of subdiagonal elements, m2 = number of superdiagonal elements,
00060   n  = number of rows and columns. L is returned in L[0..n-1][0..m1-1],
00061   U in the input matrix A; the rowwise permutation from pivoting is returned
00062   in 'index', its sign in 'perm_sign'.                                      */
00063 /*--------------------------------------------------------------------------*/
00064 template< class T >
00065 GULAPI void BandDecomposition( 
00066                         int n, int m1, int m2, Ptr< Ptr<T> >& A,
00067                         Ptr< Ptr<T> >& L, int *perm_sign, Ptr<int>& index );
00068 
00069 /*-----------------------------------------------------------------------*//**
00070   Solve system of linear equations with help of LU-Decomposition of
00071   a band matrix (given by L and A).                                                     */
00072 /*--------------------------------------------------------------------------*/
00073 template< class T >
00074 GULAPI void BandBackSubstitution( 
00075                    int n, int m1, int m2, const Ptr< Ptr<T> >& A,
00076                    const Ptr< Ptr<T> >& L, const Ptr<int>& index, Ptr<T>& b );
00077 
00078 
00079 /*----------------------------------------------------------------*//**
00080   Solve systems of linear equations with Gauss Jordan Elimination.
00081 
00082     (A)nn * (b')nm = (b)nm
00083     (A)nn * (A')nn = (1)nn 
00084 
00085   After that, matrix 'A' contains the inverse matrix of 'A',
00086   and the columns of 'b' contain the solution vectors of the
00087   different right sides in 'b'. (see "Numerical Recipes in C")       */
00088 /*-------------------------------------------------------------------*/
00089 template< class T >
00090 GULAPI bool GaussJordan( int nRowsCols, int nRightSides, 
00091                          Ptr< Ptr<T> >& A, Ptr< Ptr<T> > &b );
00092 
00093 
00094 /*---------------------------------------------------------------------*//**
00095   Calculate (a^2 + b^2)^{1/2}, with precautions to avoid destructive
00096   over- or underflow.                                                     */
00097 /*------------------------------------------------------------------------*/
00098 template< class T >
00099 inline T Pythagoras( const T& a, const T& b )
00100 {
00101   T absa,absb;
00102   
00103   absa = rtr<T>::fabs(a);
00104   absb = rtr<T>::fabs(b);
00105   if( absa > absb ) 
00106     return absa * rtr<T>::sqrt((T)1+Square(absb/absa));
00107   
00108   return absb == (T)0 ? (T)0 : absb * rtr<T>::sqrt((T)1+Square(absa/absb));
00109 }
00110 
00111 /*------------------------------------------------------------------------*//**
00112   Singular Value Decomposition.
00113   =============================
00114 
00115   A = U * W * V^{T} with
00116 
00117     A = m x n Matrix     (m >= n)
00118     U = m x n Matrix
00119     W = n x n Matrix
00120     V = n x n Matrix
00121         
00122     U,V orthogonal ( i.e. U^{T} * U = (1), V^{T} * V = (1) )
00123     W diagonal matrix, containing the singular values.                       */
00124 /*---------------------------------------------------------------------------*/
00125 template< class T >
00126 GULAPI void SVDecomposition( int m, int n, Ptr< Ptr<T> >& a,
00127                           Ptr<T>& w, Ptr< Ptr<T> >& v );
00128 
00129 /*-------------------------------------------------------------------*//**
00130   Solve system of linear equations, when a Singular Value Decomposition
00131   is given                                                              */
00132 /*----------------------------------------------------------------------*/  
00133 template< class T >            
00134 GULAPI void SVBackSubstitution( int m, int n,
00135                          Ptr< Ptr<T> >& u, Ptr<T>& w, 
00136                          Ptr< Ptr<T> >& v, Ptr<T>& b,
00137                          Ptr<T>& x );
00138 template< class T >            
00139 GULAPI void SVBackSubstitution( int m, int n,
00140                          Ptr< Ptr<T> >& u, Ptr<T>& w, 
00141                          Ptr< Ptr<T> >& v, Ptr<T>& b );
00142 
00143 /*---------------------------------------------------------------------*//**
00144   Sets very small singular values to zero. This avoids disastrous
00145   rounding errors, when the coefficient matrix of a system of
00146   linear equations is ill-conditioned.
00147   If some singular values were zero, or had to be set to zero, the function
00148   returns false, else it returns true                                     */
00149 /*------------------------------------------------------------------------*/
00150 template< class T >
00151 GULAPI bool SVZero( int n, Ptr<T>& w );
00152 
00153 
00154 /*-----------------------------------------------------------------------*//**
00155   Intersect two lines in 3-dimensions, using singular decomposition.
00156   (remark: this is a overkill, if you alredy know that the two lines 
00157   intersect (but on the other hand this is the first oppurtinity
00158   to do something senseful with the above SVD algorithm:)))))
00159 
00160   retcode = 0  => intersection point found
00161   retcode = 1  => lines are parallel
00162   retcode = 2  => lines are not parallel, but don't intersect either
00163                   (what's the english word for "windschief" ?)
00164   In the last two cases the solution minimizes |x|^2 or |Ax + b|,
00165   and the midst of the two points on the two lines is returned.              */
00166 /*---------------------------------------------------------------------------*/
00167 template< class T >
00168 int SVDIntersectLines( point<T> P1, point<T> T1,
00169                        point<T> P2, point<T> T2,
00170                        T *alpha1, T *alpha2, point<T> *P );
00171 
00172 }
00173 
00174 #endif
00175 
00176 
00177 
00178 
00179 
00180 
00181 
00182 
00183 
00184 
00185 
00186 
00187 
00188 
00189 
00190 
00191 
00192 
00193 
00194 
00195 
00196 
00197 
00198 

Generated on Mon Jan 21 04:17:34 2002 for GUL 0.6 - Geometry Utility Library by doxygen1.2.13.1 written by Dimitri van Heesch, © 1997-2001