---

guma_lineq.h

Go to the documentation of this file.

/* LIBGUL - Geometry Utility Library
 * Copyright (C) 1998-1999 Norbert Irmer
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Library General Public License for more details.
 *
 * You should have received a copy of the GNU Library General Public
 * License along with this library; if not, write to the
 * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 * Boston, MA 02111-1307, USA.
 */
 
#ifndef GUMA_LINEQ_H
#define GUMA_LINEQ_H

namespace guma {

using gul::Ptr;
using gul::rtr;
using gul::Square;
using gul::point;

/*-----------------------------------------------------------------------*//**
  Decompose a matrix into lower and upper triangular matrix. the function
  throws an exception when an singular matrix is encountered (for
  convenience, so that not only the calling program can react to this,
  but also functions higher in the call stack)                              */
/*--------------------------------------------------------------------------*/
template< class T >
GULAPI bool LUDecomposition( int nRowsCols, Ptr< Ptr<T> >& A,
                             int *perm_sign, Ptr<int>& index );

/*---------------------------------------------------------------*//**
  Solve system of linear equations with help of LU-Decomposition
  (given in A).                                                     */
/*------------------------------------------------------------------*/
template< class T >
GULAPI void LUBackSubstitution( int nRowsCols, Ptr<int>& index, 
                                Ptr< Ptr<T> >& A, Ptr<T>& b );

/*-------------------------------------------------------------------------*//**
  Multiply a band matrix and a vector: b = A*x, m1 = number of subdiagonal
  elements, m2 = number of superdiagonal elements, n = number of rows.
  It demonstrates the indexing scheme for band matrices                       */
/*----------------------------------------------------------------------------*/
template< class T >
GULAPI void BandMVProduct( int n, int m1, int m2, 
                           Ptr< Ptr<T> > a, Ptr<T> x, Ptr<T> b);

/*-----------------------------------------------------------------------*//**
  Decompose a band matrix into lower and upper triangular matrix: A = L*U.
  m1 = number of subdiagonal elements, m2 = number of superdiagonal elements,
  n  = number of rows and columns. L is returned in L[0..n-1][0..m1-1],
  U in the input matrix A; the rowwise permutation from pivoting is returned
  in 'index', its sign in 'perm_sign'.                                      */
/*--------------------------------------------------------------------------*/
template< class T >
GULAPI void BandDecomposition( 
                        int n, int m1, int m2, Ptr< Ptr<T> >& A,
                        Ptr< Ptr<T> >& L, int *perm_sign, Ptr<int>& index );

/*-----------------------------------------------------------------------*//**
  Solve system of linear equations with help of LU-Decomposition of
  a band matrix (given by L and A).                                                     */
/*--------------------------------------------------------------------------*/
template< class T >
GULAPI void BandBackSubstitution( 
                   int n, int m1, int m2, const Ptr< Ptr<T> >& A,
                   const Ptr< Ptr<T> >& L, const Ptr<int>& index, Ptr<T>& b );


/*----------------------------------------------------------------*//**
  Solve systems of linear equations with Gauss Jordan Elimination.

    (A)nn * (b')nm = (b)nm
    (A)nn * (A')nn = (1)nn 

  After that, matrix 'A' contains the inverse matrix of 'A',
  and the columns of 'b' contain the solution vectors of the
  different right sides in 'b'. (see "Numerical Recipes in C")       */
/*-------------------------------------------------------------------*/
template< class T >
GULAPI bool GaussJordan( int nRowsCols, int nRightSides, 
                         Ptr< Ptr<T> >& A, Ptr< Ptr<T> > &b );


/*---------------------------------------------------------------------*//**
  Calculate (a^2 + b^2)^{1/2}, with precautions to avoid destructive
  over- or underflow.                                                     */
/*------------------------------------------------------------------------*/
template< class T >
inline T Pythagoras( const T& a, const T& b )
{
  T absa,absb;
  
  absa = rtr<T>::fabs(a);
  absb = rtr<T>::fabs(b);
  if( absa > absb ) 
    return absa * rtr<T>::sqrt((T)1+Square(absb/absa));
  
  return absb == (T)0 ? (T)0 : absb * rtr<T>::sqrt((T)1+Square(absa/absb));
}

/*------------------------------------------------------------------------*//**
  Singular Value Decomposition.
  =============================

  A = U * W * V^{T} with

    A = m x n Matrix     (m >= n)
    U = m x n Matrix
    W = n x n Matrix
    V = n x n Matrix
        
    U,V orthogonal ( i.e. U^{T} * U = (1), V^{T} * V = (1) )
    W diagonal matrix, containing the singular values.                       */
/*---------------------------------------------------------------------------*/
template< class T >
GULAPI void SVDecomposition( int m, int n, Ptr< Ptr<T> >& a,
                          Ptr<T>& w, Ptr< Ptr<T> >& v );

/*-------------------------------------------------------------------*//**
  Solve system of linear equations, when a Singular Value Decomposition
  is given                                                              */
/*----------------------------------------------------------------------*/  
template< class T >            
GULAPI void SVBackSubstitution( int m, int n,
                         Ptr< Ptr<T> >& u, Ptr<T>& w, 
                         Ptr< Ptr<T> >& v, Ptr<T>& b,
                         Ptr<T>& x );
template< class T >            
GULAPI void SVBackSubstitution( int m, int n,
                         Ptr< Ptr<T> >& u, Ptr<T>& w, 
                         Ptr< Ptr<T> >& v, Ptr<T>& b );

/*---------------------------------------------------------------------*//**
  Sets very small singular values to zero. This avoids disastrous
  rounding errors, when the coefficient matrix of a system of
  linear equations is ill-conditioned.
  If some singular values were zero, or had to be set to zero, the function
  returns false, else it returns true                                     */
/*------------------------------------------------------------------------*/
template< class T >
GULAPI bool SVZero( int n, Ptr<T>& w );


/*-----------------------------------------------------------------------*//**
  Intersect two lines in 3-dimensions, using singular decomposition.
  (remark: this is a overkill, if you alredy know that the two lines 
  intersect (but on the other hand this is the first oppurtinity
  to do something senseful with the above SVD algorithm:)))))

  retcode = 0  => intersection point found
  retcode = 1  => lines are parallel
  retcode = 2  => lines are not parallel, but don't intersect either
                  (what's the english word for "windschief" ?)
  In the last two cases the solution minimizes |x|^2 or |Ax + b|,
  and the midst of the two points on the two lines is returned.              */
/*---------------------------------------------------------------------------*/
template< class T >
int SVDIntersectLines( point<T> P1, point<T> T1,
                       point<T> P2, point<T> T2,
                       T *alpha1, T *alpha2, point<T> *P );

}

#endif

---