---

guma Namespace Reference

Compounds
struct	guma::bbt_node
class	guma::kdnnrec
class	guma::kdrec
class	guma::kdtnode
class	guma::kdtsearch
class	guma::kdtsnode
class	guma::newton_info
class	guma::random_des
class	guma::rkdtree
class	guma::rkdtree_base
Typedefs
typedef guma::bbt_node	BBTNode
Functions
template<class T> GULAPI bool	GaussJordan (int nRowsCols, int nRightSides, Ptr< Ptr< T > > &A, Ptr< Ptr< T > > &b)
	Solve systems of linear equations with Gauss Jordan Elimination. More...
template GULAPI bool	GaussJordan (int nRowsCols, int m, Ptr< Ptr< float > > &A, Ptr< Ptr< float > > &b)
template GULAPI bool	GaussJordan (int nRowsCols, int m, Ptr< Ptr< double > > &A, Ptr< Ptr< double > > &b)
template<class T> GULAPI bool	LUDecomposition (int nRowsCols, Ptr< Ptr< T > > &A, int *perm_sign, Ptr< int > &index)
	Decompose a matrix into lower and upper triangular matrix. More...
template GULAPI bool	LUDecomposition (int nRowsCols, Ptr< Ptr< float > > &A, int *perm_sign, Ptr< int > &index)
template GULAPI bool	LUDecomposition (int nRowsCols, Ptr< Ptr< double > > &A, int *perm_sign, Ptr< int > &index)
template<class T> GULAPI void	LUBackSubstitution (int nRowsCols, Ptr< int > &index, Ptr< Ptr< T > > &A, Ptr< T > &b)
	Solve system of linear equations with help of LU-Decomposition (given in A). More...
template GULAPI void	LUBackSubstitution (int nRowsCols, Ptr< int > &index, Ptr< Ptr< float > > &A, Ptr< float > &b)
template GULAPI void	LUBackSubstitution (int nRowsCols, Ptr< int > &index, Ptr< Ptr< double > > &A, Ptr< double > &b)
template<class T> GULAPI void	BandMVProduct (int n, int m1, int m2, Ptr< Ptr< T > > a, Ptr< T > x, 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. More...
template GULAPI void	BandMVProduct (int n, int m1, int m2, Ptr< Ptr< float > > a, Ptr< float > x, Ptr< float > b)
template GULAPI void	BandMVProduct (int n, int m1, int m2, Ptr< Ptr< double > > a, Ptr< double > x, Ptr< double > b)
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)
	Decompose a band matrix into lower and upper triangular matrix: A = L*U. More...
template GULAPI void	BandDecomposition (int n, int m1, int m2, Ptr< Ptr< float > > &A, Ptr< Ptr< float > > &L, int *perm_sign, Ptr< int > &index)
template GULAPI void	BandDecomposition (int n, int m1, int m2, Ptr< Ptr< double > > &A, Ptr< Ptr< double > > &L, int *perm_sign, Ptr< int > &index)
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 system of linear equations with help of LU-Decomposition of a band matrix (given by L and A). More...
template GULAPI void	BandBackSubstitution (int n, int m1, int m2, const Ptr< Ptr< float > > &A, const Ptr< Ptr< float > > &L, const Ptr< int > &index, Ptr< float > &b)
template GULAPI void	BandBackSubstitution (int n, int m1, int m2, const Ptr< Ptr< double > > &A, const Ptr< Ptr< double > > &L, const Ptr< int > &index, Ptr< double > &b)
template<class T> GULAPI void	SVDecomposition (int m, int n, Ptr< Ptr< T > > &a, Ptr< T > &w, Ptr< Ptr< T > > &v)
	Singular Value Decomposition. More...
template GULAPI void	SVDecomposition (int m, int n, Ptr< Ptr< float > > &a, Ptr< float > &w, Ptr< Ptr< float > > &v)
template GULAPI void	SVDecomposition (int m, int n, Ptr< Ptr< double > > &a, Ptr< double > &w, Ptr< Ptr< double > > &v)
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)
	Solve system of linear equations, when a Singular Value Decomposition is given. More...
template GULAPI void	SVBackSubstitution (int m, int n, Ptr< Ptr< float > > &u, Ptr< float > &w, Ptr< Ptr< float > > &v, Ptr< float > &b, Ptr< float > &x)
template GULAPI void	SVBackSubstitution (int m, int n, Ptr< Ptr< double > > &u, Ptr< double > &w, Ptr< Ptr< double > > &v, Ptr< double > &b, Ptr< double > &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)
	Solve system of linear equations, when a Singular Value Decomposition is given. More...
template GULAPI void	SVBackSubstitution (int m, int n, Ptr< Ptr< float > > &u, Ptr< float > &w, Ptr< Ptr< float > > &v, Ptr< float > &b)
template GULAPI void	SVBackSubstitution (int m, int n, Ptr< Ptr< double > > &u, Ptr< double > &w, Ptr< Ptr< double > > &v, Ptr< double > &b)
template<class T> GULAPI bool	SVZero (int n, Ptr< T > &w)
	Sets very small singular values to zero. More...
template GULAPI bool	SVZero (int n, Ptr< float > &w)
template GULAPI bool	SVZero (int n, Ptr< double > &w)
template<class T> int	SVDIntersectLines (point< T > P1, point< T > T1, point< T > P2, point< T > T2, T *alpha1, T *alpha2, point< T > *P)
	Intersect two lines in 3-dimensions, using singular decomposition. More...
template int	SVDIntersectLines (point< float > P1, point< float > T1, point< float > P2, point< float > T2, float *alpha1, float *alpha2, point< float > *P)
template int	SVDIntersectLines (point< double > P1, point< double > T1, point< double > P2, point< double > T2, double *alpha1, double *alpha2, point< double > *P)
template<class T> T	Pythagoras (const T &a, const T &b)
	Calculate (a^2 + b^2)^{1/2}, with precautions to avoid destructive over- or underflow. More...
template<class T> bool	GoldenSectionSearch (T ax, T bx, T cx, T(*func)(T, void *), void *fdat, T tol, int max_iterations, T *xmin, T *fmin)
	Searches the minimum of a function, defined over an intervall. More...
template bool	GoldenSectionSearch (float ax, float bx, float cx, float(*func)(float, void *), void *, float tol, int max_iterations, float *xmin, float *fmin)
template bool	GoldenSectionSearch (double ax, double bx, double cx, double(*func)(double, void *), void *, double tol, int max_iterations, double *xmin, double *fmin)
template<class T> bool	LinearSearch (const Ptr< T > &p, const Ptr< T > &xold, T fold, const Ptr< T > &g, const Ptr< T > &x1, const Ptr< T > &x2, newton_info< T > &I, bool &found, bool &check, Ptr< T > &x, T *fx)
template bool	LinearSearch (const Ptr< float > &p, const Ptr< float > &xold, float fold, const Ptr< float > &g, const Ptr< float > &x1, const Ptr< float > &x2, newton_info< float > &I, bool &found, bool &check, Ptr< float > &x, float *fx)
template<class T> bool	Newton (Ptr< T > &x, const Ptr< T > &x1, const Ptr< T > &x2, newton_info< T > &I, int maxits)
template bool	Newton (Ptr< float > &x, const Ptr< float > &x1, const Ptr< float > &x2, newton_info< float > &I, int maxits)
template bool	Newton (Ptr< double > &x, const Ptr< double > &x1, const Ptr< double > &x2, newton_info< double > &I, int maxits)
GULAPI uint32	GetSeed ()
GULAPI void	PseudoDes (uint32 *lword, uint32 *irword)
GULAPI void	HeapSort (size_t nelem, size_t elsize, void *elptr, int(*usrfunc)(void *, void *, void *), void *usrdata)
GULAPI void	PtrHeapSort (size_t nelem, void **elptr, int(*usrfunc)(void *, void *, void *), void *usrdata)
GULAPI void	heapsort (int nA, Ptr< void *> &A)
GULAPI void	heapsort (int nA, Ptr< int > &A)
GULAPI int	BinarySearch (int x, int nA, const gul::Ptr< int > &A)
GULAPI int	BinarySearch (void *x, int nA, const gul::Ptr< void *> &A)
GULAPI void	InitBBT (BBTNode *head)
GULAPI int	SearchElementInBBT (BBTNode *head, void *key, int compare(void *, void *), BBTNode **element)
GULAPI BBTNode *	SearchSuccessorInBBT (BBTNode *element, BBTNode *head)
GULAPI BBTNode *	SearchPredecessorInBBT (BBTNode *element, BBTNode *head)
GULAPI void	InsertElementIntoBBT (BBTNode *element, BBTNode *where, int side, BBTNode *head)
GULAPI void	DeleteElementFromBBT (BBTNode *element, BBTNode *head)
GULAPI void	DumpBBTTree (BBTNode *node, int level, void dumpkey_func(void *))
GULAPI void	DumpBBTNode (BBTNode *node, int level)
template<class T> void	_heapsort (int nA, gul::Ptr< T > &A)
template<class T> int	_BinarySearch (T x, int nA, const gul::Ptr< T > &A)
template<class T> int	_BinarySearch (T x, int nA, const gul::Ptr< T > &A, int *rleft, int *rright)
template<class U> void	CommonCoordinateSystem (const point< U > &B, const point< U > &P, const point< U > &Ab, const point< U > &Pb, Ptr< Ptr< U > > &T)
template void	CommonCoordinateSystem (const point< float > &B, const point< float > &P, const point< float > &Ab, const point< float > &Pb, Ptr< Ptr< float > > &T)
template void	CommonCoordinateSystem (const point< double > &B, const point< double > &P, const point< double > &Ab, const point< double > &Pb, Ptr< Ptr< double > > &T)
template<class U> void	CommonCoordinateSystem (const gul::point< U > &B, const gul::point< U > &P, const gul::point< U > &Ab, const gul::point< U > &Pb, gul::Ptr< gul::Ptr< U > > &T)

---

Typedef Documentation

typedef struct guma::bbt_node guma::BBTNode

Referenced by gul::RefMap::Node::bbtnode(), gul::Map::Node::bbtnode(), gugr::edge_interval_tree::dump_edge_interval(), gul::RefMap::First(), gul::Map< kdnnrec< U, V > >::First(), gugr::edge_interval_tree::FreeNode(), gugr::edge_interval_tree::Head(), gugr::edge_interval_tree::InsertNode(), gul::RefMap::Last(), gul::Map< kdnnrec< U, V > >::Last(), gul::Map< kdnnrec< U, V > >::Map(), gul::RefMap::NewNode(), gul::Map< kdnnrec< U, V > >::NewNode(), gugr::edge_interval_tree::NewNode(), gul::Map< kdnnrec< U, V > >::NewNodeV(), gul::RefMap::Node::Node(), gul::Map::Node::Node(), gul::RefMap::operator=(), gul::Map< kdnnrec< U, V > >::operator=(), gul::RefMap::RefMap(), gul::RefMap::RemoveNode(), gul::Map< kdnnrec< U, V > >::RemoveNode(), gugr::edge_interval_tree::RemoveNode(), gul::RefMap::SearchNode(), gul::Map< kdnnrec< U, V > >::SearchNode(), gugr::edge_interval_tree::SearchNode(), gul::Map< kdnnrec< U, V > >::~Map(), and gul::RefMap::~RefMap().

---

Function Documentation

template<class T> int _BinarySearch (  T    x, int    nA, const gul::Ptr< T > &    A, int *    rleft, int *    rright )  [inline]

Definition at line 219 of file guma_sorting.h. { int mid,right,left; if( nA == 0 ) { *rleft = -1; *rright = -1; return -1; } if( A[0] >= x ) { *rright = 0; if( A[0] > x ) { *rleft = -1; return -1; } else { *rleft = 0; return 0; } } if( A[nA-1] <= x ) { *rleft = nA-1; if( A[nA-1] < x ) { *rright = -1; return -1; } else { *rright = nA-1; return nA-1; } } left = 0; right = nA; mid = (left+right)/2; while( (mid!=left) && (mid!=right) ) { if( A[mid] == x ) break; else if( A[mid] > x ) { right = mid; mid = (left+right)>>1; } else { left = mid; mid = (left+right)>>1; } } if( A[mid] != x ) { *rleft = left; *rright = right; return -1; } else { *rleft = mid; *rright = mid; } return mid; }

template<class T> int _BinarySearch (  T    x, int    nA, const gul::Ptr< T > &    A )  [inline]

Definition at line 168 of file guma_sorting.h. { int mid,right,left; if( nA == 0 ) return -1; if( A[0] >= x ) { if( A[0] > x ) return -1; return 0; } if( A[nA-1] <= x ) { if( A[nA-1] < x ) return -1; return nA-1; } left = 0; right = nA; mid = (left+right)/2; while( (mid!=left) && (mid!=right) ) { if( A[mid] == x ) break; else if( A[mid] > x ) { right = mid; mid = (left+right)>>1; } else { left = mid; mid = (left+right)>>1; } } if( A[mid] != x ) return -1; return mid; }

template<class T> void _heapsort (  int    nA, gul::Ptr< T > &    A )  [inline]

Definition at line 109 of file guma_sorting.h. { gul::Ptr<T>& ra = A; T rra; int n = nA, l, ir, i, j; if( n < 2 ) return; l = (n >> 1) + 1; ir = n; for(;;) { if( l > 1 ) { rra = ra[l-2]; l--; } else { rra = ra[ir-1]; ra[ir-1] = ra[0]; if( --ir == 1 ) { ra[0] = rra; break; } } i = l; j = l+l; while( j <= ir ) { if( (j < ir) && (ra[j-1] < ra[j]) ) j++; if( rra < ra[j-1] ) { ra[i-1] = ra[j-1]; i = j; j <<= 1; } else j = ir + 1; } ra[i-1] = rra; } }

template GULAPI void BandBackSubstitution (  int    n, int    m1, int    m2, const Ptr< Ptr< double > > &    A, const Ptr< Ptr< double > > &    L, const Ptr< int > &    index, Ptr< double > &    b )

template GULAPI void BandBackSubstitution (  int    n, int    m1, int    m2, const Ptr< Ptr< float > > &    A, const Ptr< Ptr< float > > &    L, const Ptr< int > &    index, Ptr< float > &    b )

template<class T> GULAPI void guma::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 system of linear equations with help of LU-Decomposition of a band matrix (given by L and A). Definition at line 393 of file guma_lineq.cpp. { int i,k,l,mm; T dum; mm = m1 + m2 + 1; l = m1; for( k = 0; k < n; k++ ) /* Forward substitution */ { i = index[k]; if( i != k ) Swap( b[k], b[i] ); if( l < n ) l++; for( i = k+1; i < l; i++ ) b[i] -= L[k][i-k-1]*b[k]; } l = 1; for( i = n-1; i >= 0; i-- ) /* Forward substitution */ { dum = b[i]; for( k = 1; k < l; k++ ) dum -= A[i][k]*b[k+i]; b[i] = dum/A[i][0]; if( l < mm ) l++; } }

template GULAPI void BandDecomposition (  int    n, int    m1, int    m2, Ptr< Ptr< double > > &    A, Ptr< Ptr< double > > &    L, int *    perm_sign, Ptr< int > &    index )

template GULAPI void BandDecomposition (  int    n, int    m1, int    m2, Ptr< Ptr< float > > &    A, Ptr< Ptr< float > > &    L, int *    perm_sign, Ptr< int > &    index )

template<class T> GULAPI void guma::BandDecomposition (  int    n, int    m1, int    m2, Ptr< Ptr< T > > &    A, Ptr< Ptr< T > > &    L, int *    perm_sign, Ptr< int > &    index )

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'. Definition at line 320 of file guma_lineq.cpp. { int i,j,k,l,mm,d; T dum; mm = m1 + m2 + 1; l = m1; for( i = 0; i < m1; i++ ) { for( j = m1-i; j < mm; j++ ) A[i][j-l] = A[i][j]; l--; for( j = mm-l-1; j < mm; j++ ) A[i][j] = 0.0; } d = 1; l = m1; for( k = 0; k < n; k++ ) { dum = A[k][0]; i = k; if( l < n ) l++; for( j = k+1; j < l; j++ ) { if( rtr<T>::fabs(A[j][0]) > rtr<T>::fabs(dum) ) { dum = A[j][0]; i = j; } } index[k] = i; if( dum == 0.0 ) /* if algorithmically singular */ A[k][0] = rtr<T>::tiny(); if( i != k ) /* interchange rows */ { d = -d; for( j = 0; j < mm; j++ ) Swap( A[k][j], A[i][j] ); } for( i = k+1; i < l; i++ ) // i = i'+1 { dum = A[i][0]/A[k][0]; L[k][i-k-1] = dum; for( j = 1; j < mm; j++ ) A[i][j-1] = A[i][j] - dum*A[k][j]; A[i][mm-1] = 0.0; } } *perm_sign = d; }

template GULAPI void BandMVProduct (  int    n, int    m1, int    m2, Ptr< Ptr< double > >    a, Ptr< double >    x, Ptr< double >    b )

template GULAPI void BandMVProduct (  int    n, int    m1, int    m2, Ptr< Ptr< float > >    a, Ptr< float >    x, Ptr< float >    b )

template<class T> GULAPI void guma::BandMVProduct (  int    n, int    m1, int    m2, Ptr< Ptr< T > >    a, Ptr< T >    x, 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 Definition at line 286 of file guma_lineq.cpp. { int i,k,j,m; for( i = 0; i < n; i++ ) { k = i - m1; m = Min(m1+m2+1,n-k); b[i] = 0; for( j = Max(0,-k); j < m; j++ ) b[i] += a[i][j] * x[j+k]; } }

GULAPI int guma::BinarySearch (  void *    x, int    nA, const gul::Ptr< void *> &    A )

Definition at line 153 of file guma_sorting.cpp. { return _BinarySearch( x, nA, A); }

GULAPI int guma::BinarySearch (  int    x, int    nA, const gul::Ptr< int > &    A )

Definition at line 151 of file guma_sorting.cpp. { return _BinarySearch( x, nA, A); }

template<class U> void CommonCoordinateSystem (  const gul::point< U > &    B, const gul::point< U > &    P, const gul::point< U > &    Ab, const gul::point< U > &    Pb, gul::Ptr< gul::Ptr< U > > &    T )

template void CommonCoordinateSystem (  const point< double > &    B, const point< double > &    P, const point< double > &    Ab, const point< double > &    Pb, Ptr< Ptr< double > > &    T )

template void CommonCoordinateSystem (  const point< float > &    B, const point< float > &    P, const point< float > &    Ab, const point< float > &    Pb, Ptr< Ptr< float > > &    T )

template<class U> void CommonCoordinateSystem (  const point< U > &    B, const point< U > &    P, const point< U > &    Ab, const point< U > &    Pb, Ptr< Ptr< U > > &    T )

Definition at line 49 of file guma_transform.cpp. Referenced by GUL_CommonCoordinateSystem(). { point<U> a1,a2,a3,b1,b2,b3; Ptr< Ptr<U> > Ta,Tb; int i; Ta.reserve_pool(3); Tb.reserve_pool(3); for( i = 0; i < 3; i++ ) { Ta[i].reserve_pool(3); Tb[i].reserve_pool(3); } normalize( &a1, B ); normalize( &a2, P-((P*a1)*a1) ); normalize( &a3, cross_product(a1,a2) ); normalize( &b1, -Ab ); normalize( &b2, P-((P*b1)*b1) ); normalize( &b3, cross_product(b1,b2) ); Tb[0][0] = b1[0]; Tb[0][1] = b1[1]; Tb[0][2] = b1[2]; Tb[1][0] = b2[0]; Tb[1][1] = b2[1]; Tb[1][2] = b2[2]; Tb[2][0] = b3[0]; Tb[2][1] = b3[1]; Tb[2][2] = b3[2]; GaussJordan( 3, 0, Tb, Tb ); Ta[0][0] = a1[0]; Ta[0][1] = a1[1]; Ta[0][2] = a1[2]; Ta[1][0] = a2[0]; Ta[1][1] = a2[1]; Ta[1][2] = a2[2]; Ta[2][0] = a3[0]; Ta[2][1] = a3[1]; Ta[2][2] = a3[2]; MMProduct( 3, 3, 3, Tb, Ta, T ); T[3][0] = B[0]; T[3][1] = B[1]; T[3][2] = B[2]; T[3][3] = (U)1; T[0][3] = T[1][3] = T[2][3] = (U)0; }

GULAPI void guma::DeleteElementFromBBT (  BBTNode *    element, BBTNode *    head )

Definition at line 601 of file guma_sorting.cpp. Referenced by gul::RefMap::RemoveNode(), gul::Map< kdnnrec< U, V > >::RemoveNode(), and gugr::edge_interval_tree::RemoveNode(). { BBTNode *P,*Q,*Pk,*Pk_1,*A,*B,*X; int ak, ak_1, ak_ini; P = head->right; if( P == NULL ) return; P = Q = element; if( P->right == NULL ) /* Knoten suchen, bei dem rechter oder */ { /* linker Ast = NULL: */ Pk = P; ak = 1; Pk_1 = Pk->parent; ak_1 = ak_ini = Pk->F.l; } else if( P->left == NULL ) { Pk = P; ak = -1; Pk_1 = Pk->parent; ak_1 = ak_ini = Pk->F.l; } else { P = P->right; /* nach rechts gehen */ if( P->left == NULL ) { Pk = Q; if( Q->F.l == 1 ) { Q->parent->right = P; P->F.l = 1; } else { Q->parent->left = P; P->F.l = -1; } P->parent = Q->parent; P->F.b = Q->F.b; ak = 1; Pk_1 = P; ak_1 = -1; ak_ini = 1; } else { do /* suchen bis Nachfolger in symmetrischer */ { /* Ordnung gefunden */ Pk = P; P = P->left; } while( P != NULL ); ak = -1; Pk_1 = Pk->parent; ak_1 = ak_ini = Pk->F.l; } } /* --- Knoten loeschen: -------------------------------------------- */ if( ak_1 == -1 ) /* wenn Pk linker Sohn */ { Pk_1->left = (ak == -1 ? Pk->right : Pk->left); if( Pk_1->left != NULL ) { Pk_1->left->parent = Pk_1; Pk_1->left->F.l = -1; } } else { Pk_1->right = (ak == -1 ? Pk->right : Pk->left); if( Pk_1->right != NULL ) { Pk_1->right->parent = Pk_1; Pk_1->right->F.l = 1; } } if( Pk != Q ) { Pk->left = Q->left; /* Pk uebernimmt nun die Rolle von Q */ if( Pk->left != NULL ) Pk->left->parent = Pk; Pk->right = Q->right; if( Pk->right != NULL ) Pk->right->parent = Pk; Pk->F.b = Q->F.b; Pk->parent = Q->parent; Pk->F.l = Q->F.l; if( Pk->F.l == -1 ) Pk->parent->left = Pk; else Pk->parent->right = Pk; } /* Q kann jetzt freigeben werden */ /* --- Balance-Faktoren updaten ------------------------------------- */ ak = ak_ini; Pk = Pk_1; while( 1 ) { if( Pk == head ) /* Listenkopf */ { (*((int *)&head->left))--; /* Hoehe des Baums um 1 erniedrigen */ break; /*** fertig ***/ } else if( Pk->F.b == ak ) { Pk->F.b = 0; ak = Pk->F.l; /* k-- */ Pk = Pk->parent; } else if( Pk->F.b == 0 ) { Pk->F.b = -ak; break; /*** fertig ***/ } else /* rebalancieren erforderlich, */ { /* da (Pk->B == -ak) */ if( ak == -1 ) /*** Ausfuehrung fuer (ak == -1) ***/ { A = Pk; ak_1 = Pk->F.l; Pk_1 = Pk->parent; B = A->right; /* B = Zweig auf der anderen Seite */ if( B->F.b == 1 ) /* B->B == -ak => Single Rotation */ { A->right = B->left; /* LINK(A,-ak) = LINK(B,ak) */ if( A->right != NULL ) { A->right->parent = A; A->right->F.l = 1; } B->left = A; /* LINK(B,ak) = A */ if( B->left != NULL ) { B->left->parent = B; B->left->F.l = -1; } A->F.b = B->F.b = 0; ak = ak_1; /* k-- */ Pk = Pk_1; if( ak == -1 ) { Pk->left = B; if( Pk->left != NULL ) { Pk->left->parent = Pk; Pk->left->F.l = -1; } } else { Pk->right = B; if( Pk->right != NULL ) { Pk->right->parent = Pk; Pk->right->F.l = 1; } } } else if( B->F.b == -1 ) /* (B->B == ak) => Double Rotation */ { X = B->left; /* X = LINK(B,ak) */ A->right = X->left; /* LINK(A,-ak) = LINK(X,ak) */ if( A->right != NULL ) { A->right->parent = A; A->right->F.l = 1; } B->left = X->right; /* LINK(B,ak) = LINK(X,-ak) */ if( B->left != NULL ) { B->left->parent = B; B->left->F.l = -1; } X->left = A; /* LINK(X,ak) = A */ if( X->left != NULL ) { X->left->parent = X; X->left->F.l = -1; } X->right = B; /* LINK(X,-ak) = B */ if( X->right != NULL ) { X->right->parent = X; X->right->F.l = 1; } if( X->F.b == -1 ) /* (X->B == ak) */ { A->F.b = 0; /* A->B = 0 */ B->F.b = 1; /* B->B = -ak */ } else if( X->F.b == 1 ) /* (X->B == -ak) */ { A->F.b = -1; /* A->B = ak */ B->F.b = 0; /* B->B = 0 */ } else /* B->B = 0 */ { A->F.b = 0; B->F.b = 0; } X->F.b = 0; ak = ak_1; /* k-- */ Pk = Pk_1; if( ak == -1 ) { Pk->left = X; if( Pk->left != NULL ) { Pk->left->parent = Pk; Pk->left->F.l = -1; } } else { Pk->right = X; if( Pk->right != NULL ) { Pk->right->parent = Pk; Pk->right->F.l = 1; } } } else /* (B->B == 0) => Single Rotation */ { A->right = B->left; /* LINK(A,-ak) = LINK(B,ak) */ if( A->right != NULL ) { A->right->parent = A; A->right->F.l = 1; } B->left = A; /* LINK(B,ak) = A */ if( B->left != NULL ) { B->left->parent = B; B->left->F.l = -1; } B->F.b = -1; /* B->B = ak */ ak = ak_1; /* k-- */ Pk = Pk_1; if( ak == -1 ) { Pk->left = B; if( Pk->left != NULL ) { Pk->left->parent = Pk; Pk->left->F.l = -1; } } else { Pk->right = B; if( Pk->right != NULL ) { Pk->right->parent = Pk; Pk->right->F.l = 1; } } break; /*** fertig ***/ } } else /*** dasselbe nochmal fuer (ak == +1) ***/ { A = Pk; ak_1 = Pk->F.l; Pk_1 = Pk->parent; B = A->left; /* B = Zweig auf der anderen Seite */ if( B->F.b == -1 ) /* (B->B == -ak) => Single Rotation */ { A->left = B->right; /* LINK(A,-ak) = LINK(B,ak) */ if( A->left != NULL ) { A->left->parent = A; A->left->F.l = -1; } B->right = A; if( B->right != NULL ) { B->right->parent = B; B->right->F.l = 1; } A->F.b = B->F.b = 0; ak = ak_1; /* k-- */ Pk = Pk_1; if( ak == -1 ) { Pk->left = B; if( Pk->left != NULL ) { Pk->left->parent = Pk; Pk->left->F.l = -1; } } else { Pk->right = B; if( Pk->right != NULL ) { Pk->right->parent = Pk; Pk->right->F.l = 1; } } } else if( B->F.b == 1 ) /* (B->B == ak) => Double Rotation */ { X = B->right; /* X = LINK(B,ak) */ A->left = X->right; /* LINK(A,-ak) = LINK(X,ak) */ if( A->left != NULL ) { A->left->parent = A; A->left->F.l = -1; } B->right = X->left; /* LINK(B,ak) = LINK(X,-ak) */ if( B->right != NULL ) { B->right->parent = B; B->right->F.l = 1; } X->right = A; /* LINK(X,ak) = A */ if( X->right != NULL ) { X->right->parent = X; X->right->F.l = 1; } X->left = B; /* LINK(X,-ak) = B */ if( X->left != NULL ) { X->left->parent = X; X->left->F.l = -1; } if( X->F.b == 1 ) /* (X->B == ak) */ { A->F.b = 0; /* A->B = 0 */ B->F.b = -1; /* B->B = -ak */ } else if( X->F.b == -1 ) /* (X->B == -ak) */ { A->F.b = 1; /* A->B = ak */ B->F.b = 0; /* B->B = 0 */ } else /* B->B = 0 */ { A->F.b = 0; B->F.b = 0; } X->F.b = 0; ak = ak_1; /* k-- */ Pk = Pk_1; if( ak == -1 ) { Pk->left = X; if( Pk->left != NULL ) { Pk->left->parent = Pk; Pk->left->F.l = -1; } } else { Pk->right = X; if( Pk->right != NULL ) { Pk->right->parent = Pk; Pk->right->F.l = 1; } } } else /* (B->B == 0) => Single Rotation */ { A->left = B->right; /* LINK(A,-ak) = LINK(B,ak) */ if( A->left != NULL ) { A->left->parent = A; A->left->F.l = -1; } B->right = A; /* LINK(B,ak) = A */ if( B->right != NULL ) { B->right->parent = B; B->right->F.l = 1; } B->F.b = 1; /* B->B = ak */ ak = ak_1; /* k-- */ Pk = Pk_1; if( ak == -1 ) { Pk->left = B; if( Pk->left != NULL ) { Pk->left->parent = Pk; Pk->left->F.l = -1; } } else { Pk->right = B; if( Pk->right != NULL ) { Pk->right->parent = Pk; Pk->right->F.l = 1; } } break; /*** fertig ***/ } } } /*** Ende "while"-Schleife ***/ } return; }

GULAPI void guma::DumpBBTNode (  BBTNode *    node, int    level )

Definition at line 1097 of file guma_sorting.cpp. { int i; if( node == NULL ) return; if( node->left != NULL ) DumpBBTNode( node->left, level+4 ); for( i = 0; i < level; i++ ) { cout << " "; } cout << (void *)node << ":[b:" << node->F.b << ", l:" << node->F.l << "] (l=" << (void *)node->left << ", r=" << (void *)node->right << ", p=" << (void *)node->parent << ")\n"; if( node->right != NULL ) DumpBBTNode( node->right, level+4 ); }

GULAPI void guma::DumpBBTTree (  BBTNode *    node, int    level, void dumpkey_func(void *)    )

Definition at line 1069 of file guma_sorting.cpp. { int i; if( node == NULL ) return; if( node->left != NULL ) DumpBBTTree( node->left, level+4, dumpkey_func ); for( i = 0; i < level; i++ ) { cout << " "; } if( node->key != NULL ) dumpkey_func( node->key ); cout << "\n"; if( node->right != NULL ) DumpBBTTree( node->right, level+4, dumpkey_func ); }

template GULAPI bool GaussJordan (  int    nRowsCols, int    m, Ptr< Ptr< double > > &    A, Ptr< Ptr< double > > &    b )

template GULAPI bool GaussJordan (  int    nRowsCols, int    m, Ptr< Ptr< float > > &    A, Ptr< Ptr< float > > &    b )

template<class T> GULAPI bool guma::GaussJordan (  int    nRowsCols, int    nRightSides, Ptr< Ptr< T > > &    A, Ptr< 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") Definition at line 59 of file guma_lineq.cpp. { Ptr<int> indexc,indexr; Ptr<T> ipivot; int icol,irow, n,i,j,k,l,ll,m; T big,pivotinv,dummy; n = nRowsCols - 1; m = nRightSides - 1; indexc.reserve_pool(n+1); indexr.reserve_pool(n+1); ipivot.reserve_pool(n+1); for( j = 0; j <= n; j++ ) ipivot[j] = 0; for( i = 0; i <= n; i++ ) { big = 0.0; for( j = 0; j <= n; j++ ) { if( ipivot[j] != 1 ) { for( k = 0; k <= n; k++ ) { if( ipivot[k] == 0 ) { if( rtr<T>::fabs( A[j][k] ) >= big ) { big = rtr<T>::fabs( A[j][k] ); irow = j; icol = k; } } else if( ipivot[k] > 1 ) return false; } } } ipivot[icol]++; if( irow != icol ) { for( l = 0; l <= n; l++ ) Swap( A[irow][l], A[icol][l] ); for( l = 0; l <= m; l++ ) Swap( b[irow][l], b[icol][l] ); } indexr[i] = irow; indexc[i] = icol; if( A[icol][icol] == 0.0 ) return false; pivotinv = (T)1 / A[icol][icol]; A[icol][icol] = 1.0; for( l = 0; l <= n; l++ ) A[icol][l] *= pivotinv; for( l = 0; l <= m; l++ ) b[icol][l] *= pivotinv; for( ll = 0; ll <= n; ll++ ) { if( ll != icol ) { dummy = A[ll][icol]; A[ll][icol] = 0.0; for( l = 0; l <= n; l++ ) A[ll][l] -= A[icol][l] * dummy; for( l = 0; l <= m; l++ ) b[ll][l] -= b[icol][l] * dummy; } } } for( l = n; l >= 0; l-- ) { if( indexr[l] != indexc[l] ) { for( k = 0; k <= n; k++ ) Swap( A[k][indexr[l]], A[k][indexc[l]] ); } } return true; }

GULAPI gul::uint32 guma::GetSeed (    )

Definition at line 91 of file guma_random.cpp. Referenced by guma::random_des::Random(). { uint32 seed; int fd,n; fd = open( "/dev/random", O_RDONLY ); if( fd < 0 ) { cout << "GetSeed: Couldn't open \"/dev/random\"\n"; return (uint32)time(0); } n = read( fd, &seed, sizeof(seed) ); if( n != sizeof(seed) ) { close(fd); cout << "GetSeed: Couldn't read " << sizeof(seed) << "bytes\n"; return (uint32)time(0); } close(fd); return seed; }

template bool GoldenSectionSearch (  double    ax, double    bx, double    cx, double(*    func)(double, void *), void *   , double    tol, int    max_iterations, double *    xmin, double *    fmin )

template bool GoldenSectionSearch (  float    ax, float    bx, float    cx, float(*    func)(float, void *), void *   , float    tol, int    max_iterations, float *    xmin, float *    fmin )

template<class T> bool guma::GoldenSectionSearch (  T    ax, T    bx, T    cx, T(*    func)(T, void *), void *    fdat, T    tol, int    max_iterations, T *    xmin, T *    fmin )

Searches the minimum of a function, defined over an intervall. uses Golden Section search (a three point bracketing method), see "Numerical recipes in C", intervall: ax < bx < cx, function values: f(ax) > f(bx) < f(cx) Definition at line 35 of file guma_minimize.cpp. { T x0,x1,x2,x3,f1,f2; int iteration; x0 = ax; x3 = cx; if( rtr<T>::fabs(cx-bx) > rtr<T>::fabs(bx-ax) ) { x1 = bx; x2 = bx + rtr<T>::golden_c()*(cx-bx); } else { x2 = bx; x1 = bx - rtr<T>::golden_c()*(bx-ax); } f1 = func( x1, fdat ); f2 = func( x2, fdat ); for( iteration = 0; iteration < max_iterations; iteration++ ) { if( f2 < f1 ) { x0 = x1; x1 = x2; x2 = rtr<T>::golden_r() * x2 + rtr<T>::golden_c() * x3; f1 = f2; f2 = func( x2, fdat ); } else { x3 = x2; x2 = x1; x1 = rtr<T>::golden_r() * x2 + rtr<T>::golden_c() * x0; f2 = f1; f1 = func( x1, fdat ); } if( rtr<T>::fabs(x3-x0) <= tol*(rtr<T>::fabs(x1)+rtr<T>::fabs(x2)) ) break; } if( f1 < f2 ) { *xmin = x1; *fmin = f1; } else { *xmin = x2; *fmin = f2; } if( iteration == max_iterations ) return false; return true; }

GULAPI void guma::heapsort (  int    nA, gul::Ptr< int > &    A )

Definition at line 148 of file guma_sorting.cpp. { _heapsort( nA, A ); }

GULAPI void guma::heapsort (  int    nA, gul::Ptr< void *> &    A )

Definition at line 147 of file guma_sorting.cpp. { _heapsort( nA, A ); }

GULAPI void guma::HeapSort (  size_t    nelem, size_t    elsize, void *    elptr, int(*    usrfunc)(void *, void *, void *), void *    usrdata )

Definition at line 38 of file guma_sorting.cpp. { char *ra, *rra; size_t n = nelem, l, ir, i, j; rra = (char *)alloca( elsize ); if( rra == NULL ) return; ra = (char *)elptr; if( n < 2 ) return; l = (n >> 1) + 1; ir = n; for(;;) { if( l > 1 ) { memcpy( rra, &ra[(l-2)*elsize], elsize ); l--; } else { memcpy( rra, &ra[(ir-1)*elsize], elsize ); memcpy( &ra[(ir-1)*elsize], ra, elsize ); if( --ir == 1 ) { memcpy( ra, rra, elsize ); break; } } i = l; j = l+l; while( j <= ir ) { if( (j < ir) && (usrfunc( &ra[(j-1)*elsize], &ra[j*elsize], usrdata ) < 0) ) j++; if( usrfunc(rra, &ra[(j-1)*elsize], usrdata) < 0 ) { memcpy( &ra[(i-1)*elsize], &ra[(j-1)*elsize], elsize ); i = j; j <<= 1; } else j = ir + 1; } memcpy( &ra[(i-1)*elsize], rra, elsize ); } return; }

GULAPI void guma::InitBBT (  BBTNode *    head )

Definition at line 166 of file guma_sorting.cpp. Referenced by gul::RefMap::DeleteElems(), gul::Map< kdnnrec< U, V > >::DeleteElems(), gugr::edge_interval_tree::edge_interval_tree(), gul::Map< kdnnrec< U, V > >::Map(), and gul::RefMap::RefMap(). { head->F.b = 1; head->F.l = 0; (*((int *)(&head->left))) = 0; /* height = 0 */ head->right = NULL; }

GULAPI void guma::InsertElementIntoBBT (  BBTNode *    element, BBTNode *    where, int    side, BBTNode *    head )

Definition at line 281 of file guma_sorting.cpp. Referenced by gul::RefMap::InsertNode(), gul::Map< kdnnrec< U, V > >::InsertNode(), and gugr::edge_interval_tree::InsertNode(). { BBTNode *Q, *Pk, *Pk_1, *A, *B, *X; int ak, ak_1; Q = element; /* neuen Knoten initialisieren */ Q->left = Q->right = NULL; Q->F.b = 0; Pk = where; ak = side; if( ak == -1 ) { Pk->left = Q; if( Pk->left != NULL ) { Pk->left->parent = Pk; Pk->left->F.l = -1; } } else { Pk->right = Q; if( Pk->right != NULL ) { Pk->right->parent = Pk; Pk->right->F.l = 1; } } while( (Pk->F.b == 0) && (Pk != head) ) /* Balance Faktoren aktualisieren */ { Pk->F.b = ak; ak = Pk->F.l; /* k-- */ Pk = Pk->parent; } if( Pk == head ) { (*((int *)(&head->left)))++; /* Baumhoehe im Headerknoten erhoehen */ return; /*** fertig ***/ } if( ak == -1 ) /*** eventuell Rebalancieren, Ausfuehrung fuer (ak == -1) ***/ { A = Pk; Pk_1 = Pk->parent; ak_1 = Pk->F.l; if( A->F.b == 1 ) /* (A->B == -ak) */ { A->F.b = 0; return; /*** fertig ***/ } else /* A->B == ak */ { B = Pk->left; /* B = LINK(A,ak) */ if( B->F.b == -1 ) /* (B->B == ak) => Single Rotation */ { A->left = B->right; /* LINK(A,ak) = LINK(B,-ak) */ if( A->left != NULL ) { A->left->parent = A; A->left->F.l = -1; } B->right = A; /* LINK(B,-ak) = A */ if( B->right != NULL ) { B->right->parent = B; B->right->F.l = 1; } A->F.b = B->F.b = 0; ak = ak_1; /* k-- */ Pk = Pk_1; if( ak == -1 ) /* Wurzel Teilbaum = B */ { Pk->left = B; if( Pk->left != NULL ) { Pk->left->parent = Pk; Pk->left->F.l = -1; } } else { Pk->right = B; if( Pk->right != NULL ) { Pk->right->parent = Pk; Pk->right->F.l = 1; } } return; /*** fertig ***/ } else /* (B->B == -ak) => Double Rotation */ { X = B->right; /* X = LINK(B,-ak) */ A->left = X->right; /* LINK(A,ak) = LINK(X,-ak) */ if( A->left != NULL ) { A->left->parent = A; A->left->F.l = -1; } B->right = X->left; /* LINK(B,-ak) = LINK(X,ak) */ if( B->right != NULL ) { B->right->parent = B; B->right->F.l = 1; } X->left = B; /* LINK(X,ak) = B */ if( X->left != NULL ) { X->left->parent = X; X->left->F.l = -1; } X->right = A; /* LINK(X,-ak) = A */ if( X->right != NULL ) { X->right->parent = X; X->right->F.l = 1; } if( X->F.b == 1 ) /* (X.B == -ak) */ { B->F.b = -1; /* B->B = ak */ A->F.b = 0; /* A->B = 0 */ } else if( X->F.b == -1 ) /* (X.B == ak ) */ { B->F.b = 0; /* B->B = 0 */ A->F.b = 1; /* A->B = -ak */ } else /* (X.B == 0), (Anm: ist der Fall, wenn X == Q)*/ { B->F.b = 0; /* B->B = 0 */ A->F.b = 0; /* A->B = 0 */ } X->F.b = 0; /* X->B = 0 */ ak = ak_1; /* k-- */ Pk = Pk_1; if( ak == -1 ) /* Wurzel Teilbaum = X */ { Pk->left = X; if( Pk->left != NULL ) { Pk->left->parent = Pk; Pk->left->F.l = -1; } } else { Pk->right = X; if( Pk->right != NULL ) { Pk->right->parent = Pk; Pk->right->F.l = 1; } } return; /*** fertig ***/ } } } else /*** Ausfuehrung fuer (ak == 1) ***/ { A = Pk; Pk_1 = Pk->parent; ak_1 = Pk->F.l; if( A->F.b == -1 ) /* (A->B == -ak) */ { A->F.b = 0; return; /*** fertig ***/ } else /* A->B == ak */ { B = Pk->right; /* B = LINK(A,ak) */ if( B->F.b == 1 ) /* (B->B == ak) => Single Rotation */ { A->right = B->left; /* LINK(A,ak) = LINK(B,-ak) */ if( A->right != NULL ) { A->right->parent = A; A->right->F.l = 1; } B->left = A; /* LINK(B,-ak) = A */ if( B->left != NULL ) { B->left->parent = B; B->left->F.l = -1; } A->F.b = B->F.b = 0; ak = ak_1; /* k-- */ Pk = Pk_1; if( ak == -1 ) /* Wurzel Teilbaum = B */ { Pk->left = B; if( Pk->left != NULL ) { Pk->left->parent = Pk; Pk->left->F.l = -1; } } else { Pk->right = B; if( Pk->right != NULL ) { Pk->right->parent = Pk; Pk->right->F.l = 1; } } return; /*** fertig ***/ } else /* (B->B == -ak) => Double Rotation */ { X = B->left; /* X = LINK(B,-ak) */ A->right = X->left; /* LINK(A,ak) = LINK(X,-ak) */ if( A->right != NULL ) { A->right->parent = A; A->right->F.l = 1; } B->left = X->right; /* LINK(B,-ak) = LINK(X,ak) */ if( B->left != NULL ) { B->left->parent = B; B->left->F.l = -1; } X->right = B; /* LINK(X,ak) = B */ if( X->right != NULL ) { X->right->parent = X; X->right->F.l = 1; } X->left = A; /* LINK(X,-ak) = A */ if( X->left != NULL ) { X->left->parent = X; X->left->F.l = -1; } if( X->F.b == -1 ) /* (X.B == -ak) */ { B->F.b = 1; /* B->B = ak */ A->F.b = 0; /* A->B = 0 */ } else if( X->F.b == 1 ) /* (X.B == ak ) */ { B->F.b = 0; /* B->B = 0 */ A->F.b = -1; /* A->B = -ak */ } else /* (X.B == 0), (Anm: ist der Fall, wenn X == Q)*/ { B->F.b = 0; /* B->B = 0 */ A->F.b = 0; /* A->B = 0 */ } X->F.b = 0; /* X->B = 0 */ ak = ak_1; /* k-- */ Pk = Pk_1; if( ak == -1 ) /* Wurzel Teilbaum = X */ { Pk->left = X; if( Pk->left != NULL ) { Pk->left->parent = Pk; Pk->left->F.l = -1; } } else { Pk->right = X; if( Pk->right != NULL ) { Pk->right->parent = Pk; Pk->right->F.l = 1; } } return; /*** fertig ***/ } } } }

template bool LinearSearch (  const Ptr< float > &    p, const Ptr< float > &    xold, float    fold, const Ptr< float > &    g, const Ptr< float > &    x1, const Ptr< float > &    x2, newton_info< float > &    I, bool &    found, bool &    check, Ptr< float > &    x, float *    fx )

template<class T> bool LinearSearch (  const Ptr< T > &    p, const Ptr< T > &    xold, T    fold, const Ptr< T > &    g, const Ptr< T > &    x1, const Ptr< T > &    x2, newton_info< T > &    I, bool &    found, bool &    check, Ptr< T > &    x, T *    fx )

Definition at line 50 of file guma_newton.cpp. {i+1}-x_{i}| < TOLX) */ Ptr<T>& x, /* return: point in domain */ T *fx /* return: function value there */ ) { T sum, slope, test, temp, alamin, alam, alphamin, alpha, stpmax; T a, b, alam2, disc, f2, fold2, rhs1, rhs2, tmplam; int i, n; T f; found = false; n = I.dim; // maximale Schrittweite bestimmen alphamin = 1.0; for( i = 0; i < n; i++ ) { if( xold[i] + p[i] < x1[i] ) alpha = (x1[i] - xold[i]) / p[i]; else if( xold[i] + p[i] > x2[i] ) alpha = (x2[i] - xold[i]) / p[i]; else alpha = 1.0; if( rtr<T>::fabs(alpha) == 0.0 ) { p[i] = 0.0; alpha = 1.0; } if( alpha < alphamin ) alphamin = alpha; } for( sum = 0.0, i = 0; i < n; i++ ) sum += p[i] * p[i]; stpmax = rtr<T>::fabs(alphamin) * rtr<T>::sqrt(sum); /*-------------------------------------------------------------*/ for( sum = 0.0, i = 0; i < n; i++ ) sum += p[i] * p[i]; sum = rtr<T>::sqrt( sum ); if( sum > stpmax ) for( i = 0; i < n; i++ ) p[i] *= stpmax / sum; /* Skalieren, falls Schritt zu gross */ for( slope = 0.0, i = 0; i < n; i++ ) slope += g[i] * p[i]; test = 0.0; for( i = 0; i < n; i++ ) { temp = rtr<T>::fabs(p[i]) / Max(rtr<T>::fabs(xold[i]),(T)1); if( temp > test ) test = temp; } if( test == (T)0 ) test = rtr<T>::tiny(); alamin = I.tolx / test; alam = (T)1; /* Im ersten Durchgang volle Schrittweite */ for( ;; ) { for( i = 0; i < n; i++ ) { x[i] = xold[i] + alam * p[i]; if( x[i] < x1[i] ) x[i] = x1[i]; else if( x[i] > x2[i] ) x[i] = x2[i]; } found = I.calcf( x, &f ); if( found ) { check = false; break; } if( alam < alamin ) /* Konvergenz von dx, */ { /* aufrufendes Programm sollte Konvergenz */ for( i = 0; i < n; i++ ) /* verifizieren */ x[i] = xold[i]; check = true; break; } if( f <= fold + I.alf * alam * slope ) /* Funktionswert nimmt stark */ { check = false; break; } /* genug ab */ if( alam == (T)1 ) /* falls 1.Durchgang */ tmplam = -slope / ((T)2 * (f-fold-slope)); else { rhs1 = f - fold - alam * slope; rhs2 = f2 - fold2 - alam2 * slope; a = (rhs1/(alam*alam) - rhs2/(alam2*alam2)) / (alam-alam2); b = (-alam2*rhs1/(alam*alam) + alam*rhs2/(alam2*alam2)) / (alam-alam2); if( a == (T)0 ) tmplam = -slope/((T)2*b); else { disc = b*b - (T)3*a*slope; if( disc < (T)0 ) /* roundoff problem */ return false; tmplam = (-b + rtr<T>::sqrt(disc)) / ((T)3*a); } if( tmplam > (T)0.5 * alam ) tmplam = (T)0.5 * alam; /* lambda <= 0.5 * lambda1 */ } alam2 = alam; f2 = f; fold2 = fold; alam = Max( tmplam, (T)0.1*alam ); /* lambda >= 0.1 * lambda1 */ } *fx = f; return true; }

template GULAPI void LUBackSubstitution (  int    nRowsCols, Ptr< int > &    index, Ptr< Ptr< double > > &    A, Ptr< double > &    b )

template GULAPI void LUBackSubstitution (  int    nRowsCols, Ptr< int > &    index, Ptr< Ptr< float > > &    A, Ptr< float > &    b )

template<class T> GULAPI void guma::LUBackSubstitution (  int    nRowsCols, Ptr< int > &    index, Ptr< Ptr< T > > &    A, Ptr< T > &    b )

Solve system of linear equations with help of LU-Decomposition (given in A). Definition at line 237 of file guma_lineq.cpp. { int ii,i,j,ip,n; T sum; n = nRowsCols-1; ii = -1; for( i = 0; i <= n; i++ ) { ip = index[i]; sum = b[ip]; b[ip] = b[i]; if( ii >= 0 ) { for( j = ii; j <= i-1; j++ ) sum -= A[i][j] * b[j]; } else if( sum != 0.0 ) ii = i; b[i] = sum; } for( i = n; i >= 0; i-- ) { sum = b[i]; for( j = i+1; j <= n; j++ ) sum -= A[i][j] * b[j]; b[i] = sum / A[i][i]; } }

template GULAPI bool LUDecomposition (  int    nRowsCols, Ptr< Ptr< double > > &    A, int *    perm_sign, Ptr< int > &    index )

template GULAPI bool LUDecomposition (  int    nRowsCols, Ptr< Ptr< float > > &    A, int *    perm_sign, Ptr< int > &    index )

template<class T> GULAPI bool guma::LUDecomposition (  int    nRowsCols, Ptr< Ptr< T > > &    A, int *    perm_sign, Ptr< int > &    index )

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) Definition at line 155 of file guma_lineq.cpp. { Ptr<T> vv; T big,sum,dummy,temp; int i,j,k,imax,n,d; n = nRowsCols - 1; vv.reserve_pool(n+1); /* vv will contain pivot scale factor for each row */ d = 1; for( i = 0; i <= n; i++ ) { big = 0.0; for( j = 0; j <= n; j++ ) { if( (temp = rtr<T>::fabs(A[i][j])) > big ) big = temp; } if( big == 0.0 ) return false; vv[i] = (T)1 / big; } for( j = 0; j <= n; j++ ) /* loop over columns */ { for( i = 0; i < j; i++ ) { sum = A[i][j]; for( k = 0; k < i; k++ ) sum -= A[i][k] * A[k][j]; A[i][j] = sum; } big = 0.0; for( i = j; i <= n; i++ ) { sum = A[i][j]; for( k = 0; k < j; k++ ) sum -= A[i][k] * A[k][j]; A[i][j] = sum; if( (dummy = vv[i] * rtr<T>::fabs(sum)) >= big ) // search pivot element { big = dummy; imax = i; } } if( j != imax ) /* pivot element not on diagonal */ { for( k = 0; k <= n; k++ ) Swap( A[imax][k], A[j][k] ); /* swap rows */ d = -d; vv[imax] = vv[j]; /* swap scale factors */ } index[j] = imax; if( A[j][j] == 0.0 ) A[j][j] = rtr<T>::tiny(); if( j != n ) /* divide by pivot element */ { dummy = (T)1 / A[j][j]; for( i = j+1; i <= n; i++ ) A[i][j] *= dummy; } } /* next column */ *perm_sign = d; return true; }

template bool Newton (  Ptr< double > &    x, const Ptr< double > &    x1, const Ptr< double > &    x2, newton_info< double > &    I, int    maxits )

template bool Newton (  Ptr< float > &    x, const Ptr< float > &    x1, const Ptr< float > &    x2, newton_info< float > &    I, int    maxits )

template<class T> bool guma::Newton (  Ptr< T > &    x, const Ptr< T > &    x1, const Ptr< T > &    x2, newton_info< T > &    I, int    maxits )

Definition at line 215 of file guma_newton.cpp. { int its, i, j, n; T f, fold, sum, temp, test; bool result, found, check; Ptr<int> indx; Ptr<T> g,p,xold; Ptr< Ptr<T> > A; int sign; n = I.dim; A.reserve_pool(n); for( i = 0; i < n; i++ ) A[i].reserve_pool(n); indx.reserve_pool(n); g.reserve_pool(n); p.reserve_pool(n); xold.reserve_pool(n); /* ------- Hauptschleife: ---------------------------------------- */ found = I.calcf( x, &f ); if( found ) { return true; } for( its = 1; its < maxits; its++ ) { for( i = 0; i < n; i++ ) /* Gradient von f=1/2*F*F */ { /* berechnen (=> df = F*J */ for( sum = 0.0, j = 0; j < n; j++ ) sum += I.fvec[j] * I.fjac[j][i]; g[i] = sum; } for( i = 0; i < n; i++ ) xold[i] = x[i]; fold = f; for( i = 0; i < n; i++ ) p[i] = -I.fvec[i]; /* rechte Seite der Gleichung */ MCopy(n,n,I.fjac,A); result = LUDecomposition( n, A, &sign, indx ); if( !result ) return(false); LUBackSubstitution( n, indx, A, p ); /* ------------------------------------------------------------------- */ result = LinearSearch( p, xold, fold, g, x1, x2, I, found, check, x, &f ); // if( (!result) || check ) return false; if( !result ) return false; // convergence of dx may happen when there is no better solution then a // point on the border of the surface, so for the time beeing i am // simply ignoring the 'check' flag (i am not sure, if this is allright // in all cases) if( found ) { return true; } test = 0.0; for( i = 0; i < n; i++ ) { temp = (rtr<T>::fabs(x[i]-xold[i])) / Max(rtr<T>::fabs(x[i]),(T)1); if( temp > test ) test = temp; } if( test < I.tolx ) return true; } /* maximum iterations exceeded */ return false; }

GULAPI void guma::PseudoDes (  gul::uint32 *    lword, gul::uint32 *    irword )

Definition at line 119 of file guma_random.cpp. Referenced by guma::random_des::des_random(). { static uint32 c1[4] = { 0xbaa96887L, 0x1e17d32cL, 0x03bcdc3cL, 0x0f33d1b2L }; static uint32 c2[4] = { 0x4b0f3b58L, 0xe874f0c3L, 0x6955c5a6L, 0x55a7ca46L }; uint32 i,ia,ib,iswap,itmph=0,itmpl=0; for( i = 0; i < 4; i++ ) { ia = (iswap=(*irword)) ^ c1[i]; itmpl = ia & 0xffff; itmph = ia >> 16; ib = itmpl*itmpl + ~(itmph*itmph); ia = (ib>>16) | ((ib&0xffff)<<16); *irword = (*lword) ^ ((ia ^ c2[i]) + itmpl*itmph); *lword = iswap; } }

GULAPI void guma::PtrHeapSort (  size_t    nelem, void **    elptr, int(*    usrfunc)(void *, void *, void *), void *    usrdata )

Definition at line 93 of file guma_sorting.cpp. { void **ra, *rra; size_t n = nelem, l, ir, i, j; ra = elptr; if( n < 2 ) return; l = (n >> 1) + 1; ir = n; for(;;) { if( l > 1 ) { rra = ra[l-2]; l--; } else { rra = ra[ir-1]; ra[ir-1] = ra[0]; if( --ir == 1 ) { ra[0] = rra; break; } } i = l; j = l+l; while( j <= ir ) { if( (j < ir) && (usrfunc( ra[j-1], ra[j], usrdata ) < 0) ) j++; if( usrfunc( rra, ra[j-1], usrdata) < 0 ) { ra[i-1] = ra[j-1]; i = j; j <<= 1; } else j = ir + 1; } ra[i-1] = rra; } }

template<class T> T Pythagoras (  const T &    a, const T &    b )  [inline]

Calculate (a^2 + b^2)^{1/2}, with precautions to avoid destructive over- or underflow. Definition at line 99 of file guma_lineq.h. { 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)); }

GULAPI int guma::SearchElementInBBT (  BBTNode *    head, void *    key, int compare(void *, void *)   , BBTNode **    element )

Definition at line 183 of file guma_sorting.cpp. Referenced by gul::RefMap::SearchNode(), gul::Map< kdnnrec< U, V > >::SearchNode(), and gugr::edge_interval_tree::SearchNode(). { int result; BBTNode *P,*el; el = head; result = 1; P = head->right; while( P != NULL ) /* Knoten suchen */ { el = P; result = compare( key, P->key ); if( result < 0 ) { P = P->left; } else if( result > 0 ) { P = P->right; } else break; } *element = el; return(result); }

GULAPI BBTNode * guma::SearchPredecessorInBBT (  BBTNode *    element, BBTNode *    head )

Definition at line 249 of file guma_sorting.cpp. Referenced by gul::RefMap::Predecessor(), and gul::Map< kdnnrec< U, V > >::Predecessor(). { BBTNode *N = element, *R, *P; if( N->left ) { N = N->left; for(;;) { R = N->right; if( R == NULL ) return N; N = R; } } else { for(;;) { P = N->parent; if( P == head ) return NULL; if( N->F.l == +1 ) return P; N = P; } } }

GULAPI BBTNode * guma::SearchSuccessorInBBT (  BBTNode *    element, BBTNode *    head )

Definition at line 219 of file guma_sorting.cpp. Referenced by gul::RefMap::Successor(), and gul::Map< kdnnrec< U, V > >::Successor(). { BBTNode *N = element, *L, *P; if( N->right ) { N = N->right; for(;;) { L = N->left; if( L == NULL ) return N; N = L; } } else { for(;;) { if( N == head ) return NULL; P = N->parent; if( N->F.l == -1 ) return P; N = P; } } }

template GULAPI void SVBackSubstitution (  int    m, int    n, Ptr< Ptr< double > > &    u, Ptr< double > &    w, Ptr< Ptr< double > > &    v, Ptr< double > &    b )

template GULAPI void SVBackSubstitution (  int    m, int    n, Ptr< Ptr< float > > &    u, Ptr< float > &    w, Ptr< Ptr< float > > &    v, Ptr< float > &    b )

template<class T> GULAPI void guma::SVBackSubstitution (  int    m, int    n, Ptr< Ptr< T > > &    u, Ptr< T > &    w, Ptr< Ptr< T > > &    v, Ptr< T > &    b )

Solve system of linear equations, when a Singular Value Decomposition is given. Definition at line 764 of file guma_lineq.cpp. { int jj,j,i; T s; Ptr<T> tmp; tmp.reserve_pool(n); for( j = 1; j <= n; j++ ) { s = 0.0; if( w[j-1] != 0.0 ) { for( i = 1; i <= m; i++ ) s += u[i-1][j-1] * b[i-1]; s /= w[j-1]; } tmp[j-1] = s; } for( j = 1; j <= n; j++ ) { s = 0.0; for ( jj = 1; jj <= n; jj++ ) s += v[j-1][jj-1] * tmp[jj-1]; b[j-1] = s; } return; }

template GULAPI void SVBackSubstitution (  int    m, int    n, Ptr< Ptr< double > > &    u, Ptr< double > &    w, Ptr< Ptr< double > > &    v, Ptr< double > &    b, Ptr< double > &    x )

template GULAPI void SVBackSubstitution (  int    m, int    n, Ptr< Ptr< float > > &    u, Ptr< float > &    w, Ptr< Ptr< float > > &    v, Ptr< float > &    b, Ptr< float > &    x )

template<class T> GULAPI void guma::SVBackSubstitution (  int    m, int    n, Ptr< Ptr< T > > &    u, Ptr< T > &    w, Ptr< Ptr< T > > &    v, Ptr< T > &    b, Ptr< T > &    x )

Solve system of linear equations, when a Singular Value Decomposition is given. Definition at line 718 of file guma_lineq.cpp. { int jj,j,i; T s; Ptr<T> tmp; tmp.reserve_pool(n); for( j = 1; j <= n; j++ ) { s = 0.0; if( w[j-1] != 0.0 ) { for( i = 1; i <= m; i++ ) s += u[i-1][j-1] * b[i-1]; s /= w[j-1]; } tmp[j-1] = s; } for( j = 1; j <= n; j++ ) { s = 0.0; for ( jj = 1; jj <= n; jj++ ) s += v[j-1][jj-1] * tmp[jj-1]; x[j-1] = s; } return; }

template GULAPI void SVDecomposition (  int    m, int    n, Ptr< Ptr< double > > &    a, Ptr< double > &    w, Ptr< Ptr< double > > &    v )

template GULAPI void SVDecomposition (  int    m, int    n, Ptr< Ptr< float > > &    a, Ptr< float > &    w, Ptr< Ptr< float > > &    v )

template<class T> GULAPI void guma::SVDecomposition (  int    m, int    n, Ptr< Ptr< T > > &    a, Ptr< T > &    w, Ptr< Ptr< T > > &    v )

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. Definition at line 445 of file guma_lineq.cpp. { int flag,i,its,j,jj,k,l,nm; T anorm,c,f,g,h,s,scale,x,y,z; Ptr<T> rv1; rv1.reserve_pool(n+1); /* ----- Householder reduction to bidiagonal form ------------------- */ g = scale = anorm = 0.0; for( i = 1; i <= n; i++ ) { l = i+1; rv1[i-1] = scale * g; g = s = scale = 0.0; if( i <= m ) { for( k = i; k <=m; k++ ) scale += rtr<T>::fabs( a[k-1][i-1] ); if( scale != 0.0 ) { for( k = i; k <= m; k++ ) { a[k-1][i-1] /= scale; s += a[k-1][i-1] * a[k-1][i-1]; } f = a[i-1][i-1]; g = -Sign( rtr<T>::sqrt(s), f ); h = f*g - s; a[i-1][i-1] = f-g; for( j=l; j<=n; j++ ) { for( s=0.0,k=i; k<=m; k++ ) s += a[k-1][i-1] * a[k-1][j-1]; f = s / h; for( k = i; k <= m; k++ ) a[k-1][j-1] += f * a[k-1][i-1]; } for( k = i; k <=m; k++ ) a[k-1][i-1] *= scale; } } w[i-1] = scale * g; g = s = scale = 0.0; if( (i <= m) && (i != n) ) { for( k = l; k <= n; k++ ) scale += rtr<T>::fabs( a[i-1][k-1] ); if( scale != 0.0 ) { for( k = l; k <= n; k++ ) { a[i-1][k-1] /= scale; s += a[i-1][k-1] * a[i-1][k-1]; } f = a[i-1][l-1]; g = -Sign( rtr<T>::sqrt(s), f ); h = f*g - s; a[i-1][l-1] = f-g; for( k = l; k <= n; k++ ) rv1[k-1] = a[i-1][k-1] / h; for( j = l; j <= m; j++ ) { for( s = 0.0,k=l; k <= n; k++ ) s += a[j-1][k-1] * a[i-1][k-1]; for( k = l; k <= n; k++ ) a[j-1][k-1] += s * rv1[k-1]; } for( k = l; k <= n; k++ ) a[i-1][k-1] *= scale; } } anorm = Max( anorm, (rtr<T>::fabs(w[i-1]) + rtr<T>::fabs(rv1[i-1])) ); } /* ------- akkumulation of right side transformations ------------- */ for( i = n; i >= 1; i-- ) { if( i < n ) { if( g != 0.0 ) { for( j = l; j <= n; j++ ) v[j-1][i-1] = (a[i-1][j-1] / a[i-1][l-1]) / g; for( j = l; j <= n; j++ ) { for( s=0.0,k=l; k<=n; k++ ) s += a[i-1][k-1] * v[k-1][j-1]; for( k = l; k <= n; k++ ) v[k-1][j-1] += s * v[k-1][i-1]; } } for( j = l; j <= n; j++ ) v[i-1][j-1] = v[j-1][i-1] = 0.0; } v[i-1][i-1] = 1.0; g = rv1[i-1]; l = i; } /* -------- accumulation of left side transformations --------------- */ for( i = Min(m,n); i >=1; i-- ) { l = i+1; g = w[i-1]; for( j = l; j <= n; j++ ) a[i-1][j-1] = 0.0; if( g != 0.0 ) { g = (T)1 / g; for( j = l; j <= n; j++ ) { for( s=0.0,k=l; k <= m; k++ ) s += a[k-1][i-1] * a[k-1][j-1]; f = (s/a[i-1][i-1]) * g; for( k = i; k <= m; k++ ) a[k-1][j-1] += f * a[k-1][i-1]; } for( j = i; j <= m; j++ ) a[j-1][i-1] *= g; } else for( j = i; j <= m; j++ ) a[j-1][i-1] = 0.0; a[i-1][i-1] += 1.0; } /* -------- Diagonalize the bidiagonal form --------------------- */ for( k = n; k >= 1; k-- ) { for( its = 1; its <= 30; its++ ) { flag = 1; for( l = k; l >= 1; l-- ) { nm = l-1; if( rtr<T>::fabs(rv1[l-1]) + anorm == anorm ) { flag = 0; break; } if( rtr<T>::fabs(w[nm-1]) + anorm == anorm ) break; } if( flag != 0 ) { c = 0.0; s = 1.0; for( i = l; i <= k; i++ ) { f = s * rv1[i-1]; rv1[i-1] = c * rv1[i-1]; if( rtr<T>::fabs(f) + anorm == anorm ) break; g = w[i-1]; h = Pythagoras( f, g ); w[i-1] = h; h = (T)1 / h; c = g * h; s = -f * h; for( j = 1; j <= m; j++ ) { y = a[j-1][nm-1]; z = a[j-1][i-1]; a[j-1][nm-1] = y * c + z * s; a[j-1][i-1] = z * c - y * s; } } } z = w[k-1]; if( l == k ) { if( z < 0.0 ) { w[k-1] = -z; for( j = 1; j <= n; j++ ) v[j-1][k-1] = -v[j-1][k-1]; } break; } if( its == 30 ) throw InternalError(); x = w[l-1]; nm = k-1; y = w[nm-1]; g = rv1[nm-1]; h = rv1[k-1]; f = ((y-z)*(y+z)+(g-h)*(g+h))/((T)2*h*y); g = Pythagoras( f, (T)1 ); f = ((x-z)*(x+z) + h*((y/(f+Sign(g,f)))-h)) / x; c = s = 1.0; for( j = l; j <= nm; j++ ) { i = j+1; g = rv1[i-1]; y = w[i-1]; h = s * g; g = c * g; z = Pythagoras( f, h ); rv1[j-1] = z; c = f/z; s = h/z; f = x*c + g*s; g = g*c - x*s; h = y*s; y *= c; for( jj = 1; jj <= n; jj++ ) { x = v[jj-1][j-1]; z = v[jj-1][i-1]; v[jj-1][j-1] = x*c + z*s; v[jj-1][i-1] = z*c - x*s; } z = Pythagoras( f, h ); w[j-1] = z; if( z != 0.0 ) { z = (T)1 / z; c = f*z; s = h*z; } f = c*g + s*y; x = c*y - s*g; for( jj = 1; jj <= m; jj++ ) { y = a[jj-1][j-1]; z = a[jj-1][i-1]; a[jj-1][j-1] = y*c + z*s; a[jj-1][i-1] = z*c - y*s; } } rv1[l-1] = 0.0; rv1[k-1] = f; w[k-1] = x; } } }

template int SVDIntersectLines (  point< double >    P1, point< double >    T1, point< double >    P2, point< double >    T2, double *    alpha1, double *    alpha2, point< double > *    P )

template int SVDIntersectLines (  point< float >    P1, point< float >    T1, point< float >    P2, point< float >    T2, float *    alpha1, float *    alpha2, point< float > *    P )

template<class T> int guma::SVDIntersectLines (  point< T >    P1, point< T >    T1, point< T >    P2, point< T >    T2, T *    alpha1, T *    alpha2, point< T > *    P )

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. Definition at line 858 of file guma_lineq.cpp. { Ptr< Ptr<T> > u,v; Ptr<T> w, x, b; point<T> L1,L2; int code,i; u.reserve_pool(3); for( i = 0; i < 3; i++ ) u[i].reserve_pool(2); v.reserve_pool(2); for( i = 0; i < 2; i++ ) v[i].reserve_pool(2); w.reserve_pool(2); x.reserve_pool(2); b.reserve_pool(3); code = 0; u[0][0] = T1.x; u[0][1] = -T2.x; u[1][0] = T1.y; u[1][1] = -T2.y; u[2][0] = T1.z; u[2][1] = -T2.z; SVDecomposition( 3, 2, u, w, v ); if( SVZero( 2, w ) ) // are the lines parallel ? code = 1; b[0] = P2.x - P1.x; b[1] = P2.y - P1.y; b[2] = P2.z - P1.z; SVBackSubstitution( 3, 2, u, w, v, b, x ); L1.x = P1.x + x[0] * T1.x; L1.y = P1.y + x[0] * T1.y; L1.z = P1.z + x[0] * T1.z; L2.x = P2.x + x[1] * T2.x; L2.y = P2.y + x[1] * T2.y; L2.z = P2.z + x[1] * T2.z; // not parallel, but points too far away from each other ? if( !code && !rel_equal(L1, L2, rtr<T>::zero_tol() ) ) code = 2; P->x = (L1.x + L2.x) / (T)2; P->y = (L1.y + L2.y) / (T)2; P->z = (L1.z + L2.z) / (T)2; *alpha1 = x[0]; *alpha2 = x[1]; return code; }

template GULAPI bool SVZero (  int    n, Ptr< double > &    w )

template GULAPI bool SVZero (  int    n, Ptr< float > &    w )

template<class T> GULAPI bool guma::SVZero (  int    n, Ptr< T > &    w )

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 Definition at line 809 of file guma_lineq.cpp. { T wmax,wmin; int j; bool well_conditioned; T limit = (T)n * rtr<T>::epsilon(); wmax = wmin = rtr<T>::fabs(w[0]); for( j = 2; j <= n; j++ ) { if( rtr<T>::fabs(w[j-1]) > wmax ) wmax = rtr<T>::fabs(w[j-1]); else if( rtr<T>::fabs(w[j-1]) < wmin ) wmin = rtr<T>::fabs(w[j-1]); } // ill conditioned ? if( (wmax == 0.0) || (wmin / wmax < limit) ) well_conditioned = false; else well_conditioned = true; /* cout << "wmin/wmax = " << wmin / wmax << "\n"; */ wmin = wmax * limit; for( j = 1; j <= n; j++ ) if( wmin > rtr<T>::fabs(w[j-1]) ) w[j-1] = 0.0; return well_conditioned; }

---