[GRASS-SVN] r38888 - in grass/branches/develbranch_6: include/Make
lib/external lib/external/ccmath
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Aug 27 16:50:28 EDT 2009
Author: huhabla
Date: 2009-08-27 16:50:27 -0400 (Thu, 27 Aug 2009)
New Revision: 38888
Added:
grass/branches/develbranch_6/lib/external/ccmath/
grass/branches/develbranch_6/lib/external/ccmath/C01-matrix
grass/branches/develbranch_6/lib/external/ccmath/Makefile
grass/branches/develbranch_6/lib/external/ccmath/README
grass/branches/develbranch_6/lib/external/ccmath/atou1.c
grass/branches/develbranch_6/lib/external/ccmath/atovm.c
grass/branches/develbranch_6/lib/external/ccmath/ccmath.h
grass/branches/develbranch_6/lib/external/ccmath/ccmath_grass.h
grass/branches/develbranch_6/lib/external/ccmath/chouse.c
grass/branches/develbranch_6/lib/external/ccmath/chousv.c
grass/branches/develbranch_6/lib/external/ccmath/cmattr.c
grass/branches/develbranch_6/lib/external/ccmath/cmcpy.c
grass/branches/develbranch_6/lib/external/ccmath/cminv.c
grass/branches/develbranch_6/lib/external/ccmath/cmmul.c
grass/branches/develbranch_6/lib/external/ccmath/cmmult.c
grass/branches/develbranch_6/lib/external/ccmath/cmprt.c
grass/branches/develbranch_6/lib/external/ccmath/csolv.c
grass/branches/develbranch_6/lib/external/ccmath/cvmul.c
grass/branches/develbranch_6/lib/external/ccmath/eigen.c
grass/branches/develbranch_6/lib/external/ccmath/eigval.c
grass/branches/develbranch_6/lib/external/ccmath/evmax.c
grass/branches/develbranch_6/lib/external/ccmath/hconj.c
grass/branches/develbranch_6/lib/external/ccmath/heigval.c
grass/branches/develbranch_6/lib/external/ccmath/heigvec.c
grass/branches/develbranch_6/lib/external/ccmath/hevmax.c
grass/branches/develbranch_6/lib/external/ccmath/hmgen.c
grass/branches/develbranch_6/lib/external/ccmath/house.c
grass/branches/develbranch_6/lib/external/ccmath/housev.c
grass/branches/develbranch_6/lib/external/ccmath/ldumat.c
grass/branches/develbranch_6/lib/external/ccmath/ldvmat.c
grass/branches/develbranch_6/lib/external/ccmath/lgpl.license
grass/branches/develbranch_6/lib/external/ccmath/matprt.c
grass/branches/develbranch_6/lib/external/ccmath/mattr.c
grass/branches/develbranch_6/lib/external/ccmath/mcopy.c
grass/branches/develbranch_6/lib/external/ccmath/minv.c
grass/branches/develbranch_6/lib/external/ccmath/mmul.c
grass/branches/develbranch_6/lib/external/ccmath/ortho.c
grass/branches/develbranch_6/lib/external/ccmath/otrma.c
grass/branches/develbranch_6/lib/external/ccmath/otrsm.c
grass/branches/develbranch_6/lib/external/ccmath/psinv.c
grass/branches/develbranch_6/lib/external/ccmath/qrbdi.c
grass/branches/develbranch_6/lib/external/ccmath/qrbdu1.c
grass/branches/develbranch_6/lib/external/ccmath/qrbdv.c
grass/branches/develbranch_6/lib/external/ccmath/qrecvc.c
grass/branches/develbranch_6/lib/external/ccmath/qreval.c
grass/branches/develbranch_6/lib/external/ccmath/qrevec.c
grass/branches/develbranch_6/lib/external/ccmath/rmmult.c
grass/branches/develbranch_6/lib/external/ccmath/ruinv.c
grass/branches/develbranch_6/lib/external/ccmath/smgen.c
grass/branches/develbranch_6/lib/external/ccmath/solv.c
grass/branches/develbranch_6/lib/external/ccmath/solv.s
grass/branches/develbranch_6/lib/external/ccmath/solvps.c
grass/branches/develbranch_6/lib/external/ccmath/solvru.c
grass/branches/develbranch_6/lib/external/ccmath/solvtd.c
grass/branches/develbranch_6/lib/external/ccmath/sv2u1v.c
grass/branches/develbranch_6/lib/external/ccmath/sv2uv.c
grass/branches/develbranch_6/lib/external/ccmath/sv2val.c
grass/branches/develbranch_6/lib/external/ccmath/svdu1v.c
grass/branches/develbranch_6/lib/external/ccmath/svduv.c
grass/branches/develbranch_6/lib/external/ccmath/svdval.c
grass/branches/develbranch_6/lib/external/ccmath/trncm.c
grass/branches/develbranch_6/lib/external/ccmath/trnm.c
grass/branches/develbranch_6/lib/external/ccmath/unfl.c
grass/branches/develbranch_6/lib/external/ccmath/unitary.c
grass/branches/develbranch_6/lib/external/ccmath/utrncm.c
grass/branches/develbranch_6/lib/external/ccmath/utrnhm.c
grass/branches/develbranch_6/lib/external/ccmath/vmul.c
Modified:
grass/branches/develbranch_6/include/Make/Grass.make.in
grass/branches/develbranch_6/lib/external/Makefile
Log:
New mathematical library ccmath in lib/external to replace the
numerical recipes code in gmath library
Modified: grass/branches/develbranch_6/include/Make/Grass.make.in
===================================================================
--- grass/branches/develbranch_6/include/Make/Grass.make.in 2009-08-27 20:04:06 UTC (rev 38887)
+++ grass/branches/develbranch_6/include/Make/Grass.make.in 2009-08-27 20:50:27 UTC (rev 38888)
@@ -94,6 +94,7 @@
ARRAYSTATS_LIBNAME = grass_arraystats
BITMAP_LIBNAME = grass_bitmap
BTREE_LIBNAME = grass_btree
+CCMATH_LIBNAME = grass_ccmath
CLUSTER_LIBNAME = grass_cluster
COORCNV_LIBNAME = grass_coorcnv
DATETIME_LIBNAME = grass_datetime
@@ -192,6 +193,7 @@
ARRAYSTATSLIB = -l$(ARRAYSTATS_LIBNAME) $(GISLIB)
BITMAPLIB = -l$(BITMAP_LIBNAME) $(LINKMLIB)
BTREELIB = -l$(BTREE_LIBNAME)
+CCMATHLIB = -l$(CCMATH_LIBNAME)
CLUSTERLIB = -l$(CLUSTER_LIBNAME) $(IMAGERYLIB) $(GISLIB)
COORCNVLIB = -l$(COORCNV_LIBNAME)
DATETIMELIB = -l$(DATETIME_LIBNAME)
@@ -206,7 +208,7 @@
EDITLIB = -l$(EDIT_LIBNAME) $(GISLIB) $(VASKLIB)
G3DLIB = -l$(G3D_LIBNAME) $(GISLIB)
GISLIB = -l$(GIS_LIBNAME) $(DATETIMELIB) $(XDRLIB) $(SOCKLIB) $(INTLLIB)
-GMATHLIB = -l$(GMATH_LIBNAME) $(GISLIB)
+GMATHLIB = -l$(GMATH_LIBNAME) $(GISLIB) $(CCMATHLIB)
GPDELIB = -l$(GPDE_LIBNAME) $(GISLIB) $(G3DLIB)
GPROJLIB = -l$(GPROJ_LIBNAME) $(GISLIB) $(PROJLIB) $(GDALLIBS)
IBTREELIB = -l$(IBTREE_LIBNAME)
@@ -292,6 +294,7 @@
ARRAYSTATSDEP = $(ARCH_LIBDIR)/$(LIB_PREFIX)$(ARRAYSTATS_LIBNAME)$(LIB_SUFFIX)
BITMAPDEP = $(ARCH_LIBDIR)/$(LIB_PREFIX)$(BITMAP_LIBNAME)$(LIB_SUFFIX)
BTREEDEP = $(ARCH_LIBDIR)/$(LIB_PREFIX)$(BTREE_LIBNAME)$(LIB_SUFFIX)
+CCMATHDEP = $(ARCH_LIBDIR)/$(LIB_PREFIX)$(CCMATH_LIBNAME)$(LIB_SUFFIX)
COORCNVDEP = $(ARCH_LIBDIR)/$(LIB_PREFIX)$(COORCNV_LIBNAME)$(LIB_SUFFIX)
CLUSTERDEP = $(ARCH_LIBDIR)/$(LIB_PREFIX)$(CLUSTER_LIBNAME)$(LIB_SUFFIX)
DATETIMEDEP = $(ARCH_LIBDIR)/$(LIB_PREFIX)$(DATETIME_LIBNAME)$(LIB_SUFFIX)
Modified: grass/branches/develbranch_6/lib/external/Makefile
===================================================================
--- grass/branches/develbranch_6/lib/external/Makefile 2009-08-27 20:04:06 UTC (rev 38887)
+++ grass/branches/develbranch_6/lib/external/Makefile 2009-08-27 20:50:27 UTC (rev 38888)
@@ -2,6 +2,7 @@
MODULE_TOPDIR = ../..
SUBDIRS = \
+ ccmath \
bwidget \
shapelib
Added: grass/branches/develbranch_6/lib/external/ccmath/C01-matrix
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/C01-matrix (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/C01-matrix 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,1280 @@
+ Chapter 1
+
+ LINEAR ALGEBRA
+
+ Summary
+
+ The matrix algebra library contains functions that
+ perform the standard computations of linear algebra.
+ General areas covered are:
+
+ o Solution of Linear Systems
+ o Matrix Inversion
+ o Eigensystem Analysis
+ o Matrix Utility Operations
+ o Singular Value Decomposition
+
+ The operations covered here are fundamental to many
+ areas of mathematics and statistics. Thus, functions
+ in this library segment are called by other library
+ functions. Both real and complex valued matrices
+ are covered by functions in the first four of these
+ categories.
+
+
+ Notes on Contents
+
+ Functions in this library segment provide the basic operations of
+ numerical linear algebra and some useful utility functions for operations on
+ vectors and matrices. The following list describes the functions available for
+ operations with real-valued matrices.
+
+
+ o Solving and Inverting Linear Systems:
+
+ solv --------- solve a general system of real linear equations.
+ solvps ------- solve a real symmetric linear system.
+ solvru ------- solve a real right upper triangular linear system.
+ solvtd ------- solve a tridiagonal real linear system.
+
+ minv --------- invert a general real square matrix.
+ psinv -------- invert a real symmetric matrix.
+ ruinv -------- invert a right upper triangular matrix.
+
+
+ The solution of a general linear system and efficient algorithms for
+ solving special systems with symmetric and tridiagonal matrices are provided
+ by these functions. The general solution function employs a LU factorization
+ with partial pivoting and it is very robust. It will work efficiently on any
+ problem that is not ill-conditioned. The symmetric matrix solution is based
+ on a modified Cholesky factorization. It is best used on positive definite
+ matrices that do not require pivoting for numeric stability. Tridiagonal
+ solvers require order-N operations (N = dimension). Thus, they are highly
+ recommended for this important class of sparse systems. Two matrix inversion
+ routines are provided. The general inversion function is again LU based. It
+ is suitable for use on any stable (ie. well-conditioned) problem. The
+ Cholesky based symmetric matrix inversion is efficient and safe for use on
+ matrices known to be positive definite, such as the variance matrices
+ encountered in statistical computations. Both the solver and the inverse
+ functions are designed to enhance data locality. They are very effective
+ on modern microprocessors.
+
+
+ o Eigensystem Analysis:
+
+ eigen ------ extract all eigen values and vectors of a real
+ symmetric matrix.
+ eigval ----- extract the eigen values of a real symmetric matrix.
+ evmax ------ compute the eigen value of maximum absolute magnitude
+ and its corresponding vector for a symmetric matrix.
+
+
+ Eigensystem functions operate on real symmetric matrices. Two forms of
+ the general eigen routine are provided because the computation of eigen values
+ only is much faster when vectors are not required. The basic algorithms use
+ a Householder reduction to tridiagonal form followed by QR iterations with
+ shifts to enhance convergence. This has become the accepted standard for
+ symmetric eigensystem computation. The evmax function uses an efficient
+ iterative power method algorithm to extract the eigen value of maximum
+ absolute size and the corresponding eigenvector.
+
+
+ o Singular Value Decomposition:
+
+ svdval ----- compute the singular values of a m by n real matrix.
+ sv2val ----- compute the singular values of a real matrix
+ efficiently for m >> n.
+ svduv ------ compute the singular values and the transformation
+ matrices u and v for a real m by n matrix.
+ sv2uv ------ compute the singular values and transformation
+ matrices efficiently for m >> n.
+ svdu1v ----- compute the singular values and transformation
+ matrices u1 and v, where u1 overloads the input
+ with the first n column vectors of u.
+ sv2u1v ----- compute the singular values and the transformation
+ matrices u1 and v efficiently for m >> n.
+
+
+ Singular value decomposition is extremely useful when dealing with linear
+ systems that may be singular. Singular values with values near zero are flags
+ of a potential rank deficiency in the system matrix. They can be used to
+ identify the presence of an ill-conditioned problem and, in some cases, to
+ deal with the potential instability. They are applied to the linear least
+ squares problem in this library. Singular values also define some important
+ matrix norm parameters such as the 2-norm and the condition value. A complete
+ decomposition provides both singular values and an orthogonal decomposition of
+ vector spaces related to the matrix identifying the range and null-space.
+ Fortunately, a highly stable algorithm based on Householder reduction to
+ bidiagonal form and QR rotations can be used to implement the decomposition.
+ The library provides two forms with one more efficient when the dimensions
+ satisfy m > (3/2)n.
+
+
+ o Real Matrix Utilities:
+
+ rmmult ----- multiply two compatible real matrices.
+ mmul ------- multiply two real square matrices.
+ vmul ------- multiply a vector by a square matrix (transform).
+ mattr ------ compute the transpose of a general matrix.
+ trnm ------- transpose a real square matrix in place.
+ otrma ------ compute orthogonal conjugate of a square matrix.
+ otrsm ------ compute orthogonal conjugate of a symmetric matrix.
+ smgen ------ construct a symmetric matrix from its eigen values
+ and vectors.
+ ortho ------ generate a general orthogonal matrix (uses random
+ rotation generator).
+ mcopy ------ make a copy of an array.
+ matprt ----- print a matrix in row order with a specified
+ format for elements to stdout or a file (fmatprt).
+
+
+ The utility functions perform simple matrix operations such as matrix
+ multiplication, transposition, matrix conjugation, and linear transformation
+ of a vector. They are used to facilitate the rapid production of test and
+ application code requiring these operations. The function call overhead
+ associated with use of these utilities becomes negligible as the dimension of
+ the matrices increases. The implementation of these routines is designed to
+ enhance data locality and thus execution performance. A generator for
+ orthogonal matrices can be used to generate test matrices.
+
+
+ Linear Algebra with Complex Matrices
+
+ The complex section of the linear algebra library contains functions that
+ operate on general complex valued matrices and on Hermitian matrices.
+ Hermitian matrices are the complex analog of symmetric matrices. Complex
+ arithmetic is performed in-line in these functions to ensure efficient
+ execution.
+
+
+ o Solving and Inverting Complex Linear Systems:
+
+ csolv ------ solve a general complex system of linear equations.
+
+ cminv ------ invert a general complex matrix.
+
+
+ Both these functions are based on the robust LU factorization algorithm
+ with row pivots. They can be expected to solve or invert any system that is
+ well conditioned.
+
+
+ o Hermitian Eigensystem Analysis:
+
+ heigvec ---- extract the eigen values and vectors of a Hermitian
+ matrix.
+ heigval ---- compute the eigenvalues of a Hermitian matrix.
+ hevmax ----- compute the eigen value of largest absolute magnitude
+ and the corresponding vector of a Hermitian matrix.
+
+
+ The algorithms used for complex eigensystems are complex generalizations
+ of those employed in the real systems. The eigen values of Hermitian matrices
+ are real and their eigenvectors form a unitary matrix. As in the real case,
+ the function for eigen values only is provided as a time saver. These
+ routines have important application to many quantum mechanical problems.
+
+
+ o Complex Matrix Utilities:
+
+ cmmult ---- multiply two general, size compatible, complex matrices.
+ cmmul ----- multiply two square complex matrices.
+ cvmul ----- multiply a complex vector by a complex square matrix.
+ cmattr ---- transpose a general complex matrix.
+ trncm ----- transpose a complex square matrix in place.
+ hconj ----- transform a square complex matrix to its Hermitian
+ conjugate in place.
+ utrncm ---- compute the unitary transform of a complex matrix.
+ utrnhm ---- compute the unitary transform of a Hermitian matrix.
+ hmgen ----- generate a general Hermitian from its eigen values
+ and vectors.
+ unitary --- generate a general unitary matrix (uses a random
+ rotation generator).
+ cmcpy ----- copy a complex array.
+ cmprt ----- print a complex matrix in row order with a specified
+ format for matrix elements.
+
+
+ These utility operations replicate the utilities available for real
+ matrices. Matrix computations implemented in a manner that enhances data
+ locality. This ensures their efficiency on modern computers with a memory
+ hierarchy.
+
+-------------------------------------------------------------------------------
+
+ General Technical Comments
+
+ Efficient computation with matrices on modern processors must be
+ adapted to the storage scheme employed for matrix elements. The functions
+ of this library segment do not employ the multidimensional array intrinsic
+ of the C language. Access to elements employs the simple row-major scheme
+ described here.
+
+ Matrices are modeled by the library functions as arrays with elements
+ stored in row order. Thus, the element in the jth row and kth column of
+ the n by n matrix M, stored in the array mat[], is addressed by
+
+ M[j,k] = mat[n*j+k] , with 0 =< j,k <= n-1 .
+
+ (Remember that C employs zero as the starting index.) The storage order has
+ important implications for data locality.
+
+ The algorithms employed here all have excellent numerical stability, and
+ the default double precision arithmetic of C enhances this. Thus, any
+ problems encountered in using the matrix algebra functions will almost
+ certainly be due to an ill-conditioned matrix. (The Hilbert matrices,
+
+ H[i,j] = 1/(1+i+j) for i,j < n
+
+ form a good example of such ill-conditioned systems.) We remind the reader
+ that the appropriate response to such ill-conditioning is to seek an
+ alternative approach to the problem. The option of increasing precision has
+ already been exploited. Modification of the linear algebra algorithm code is
+ not normally effective in an ill-conditioned problem.
+
+------------------------------------------------------------------------------
+ FUNCTION SYNOPSES
+------------------------------------------------------------------------------
+
+ Linear System Solutions:
+-----------------------------------------------------------------------------
+
+solv
+
+ Solve a general linear system A*x = b.
+
+ int solv(double a[],double b[],int n)
+ a = array containing system matrix A in row order
+ (altered to L-U factored form by computation)
+ b = array containing system vector b at entry and
+ solution vector x at exit
+ n = dimension of system
+ return: 0 -> normal exit
+ -1 -> singular input
+
+ -----------------------------------------------------------
+
+solvps
+
+ Solve a symmetric positive definite linear system S*x = b.
+
+ int solvps(double a[],double b[],int n)
+ a = array containing system matrix S (altered to
+ Cholesky upper right factor by computation)
+ b = array containing system vector b as input and
+ solution vector x as output
+ n = dimension of system
+ return: 0 -> normal exit
+ 1 -> input matrix not positive definite
+
+ --------------------------------------------------------------
+
+solvtd
+
+ Solve a tridiagonal linear system M*x = y.
+
+ void solvtd(double a[],double b[],double c[],double x[],int m)
+ a = array containing m+1 diagonal elements of M
+ b = array of m elements below the main diagonal of M
+ c = array of m elements above the main diagonal
+ x = array containing the system vector y initially, and
+ the solution vector at exit (m+1 elements)
+ m = dimension parameter ( M is (m+1)x(m+1) )
+
+ --------------------------------------------------------------
+
+solvru
+
+ Solve an upper right triangular linear system T*x = b.
+
+ int solvru(double *a,double *b,int n)
+ a = pointer to array of upper right triangular matrix T
+ b = pointer to array of system vector
+ The computation overloads this with the
+ solution vector x.
+ n = dimension (dim(a)=n*n,dim(b)=n)
+ return value: f = status flag, with 0 -> normal exit
+ -1 -> system singular
+
+------------------------------------------------------------------------------
+
+ Matrix Inversion:
+------------------------------------------------------------------------------
+
+minv
+
+ Invert (in place) a general real matrix A -> Inv(A).
+
+ int minv(double a[],int n)
+ a = array containing the input matrix A
+ This is converted to the inverse matrix.
+ n = dimension of the system (i.e. A is n x n )
+ return: 0 -> normal exit
+ 1 -> singular input matrix
+
+ --------------------------------------------------------------
+
+psinv
+
+ Invert (in place) a symmetric real matrix, V -> Inv(V).
+
+ int psinv(double v[],int n)
+ v = array containing a symmetric input matrix
+ This is converted to the inverse matrix.
+ n = dimension of the system (dim(v)=n*n)
+ return: 0 -> normal exit
+ 1 -> input matrix not positive definite
+
+ The input matrix V is symmetric (V[i,j] = V[j,i]).
+
+ --------------------------------------------------------------
+
+ruinv
+
+ Invert an upper right triangular matrix T -> Inv(T).
+
+ int ruinv(double *a,int n)
+ a = pointer to array of upper right triangular matrix
+ This is replaced by the inverse matrix.
+ n = dimension (dim(a)=n*n)
+ return value: status flag, with 0 -> matrix inverted
+ -1 -> matrix singular
+
+
+-----------------------------------------------------------------------------
+
+ Symmetric Eigensystem Analysis:
+-----------------------------------------------------------------------------
+
+eigval
+
+ Compute the eigenvalues of a real symmetric matrix A.
+
+ void eigval(double *a,double *ev,int n)
+ a = pointer to array of symmetric n by n input
+ matrix A. The computation alters these values.
+ ev = pointer to array of the output eigenvalues
+ n = dimension parameter (dim(a)= n*n, dim(ev)= n)
+
+ --------------------------------------------------------------
+
+eigen
+
+ Compute the eigenvalues and eigenvectors of a real symmetric
+ matrix A.
+
+ void eigen(double *a,double *ev,int n)
+ double *a,*ev; int n;
+ a = pointer to store for symmetric n by n input
+ matrix A. The computation overloads this with an
+ orthogonal matrix of eigenvectors E.
+ ev = pointer to the array of the output eigenvalues
+ n = dimension parameter (dim(a)= n*n, dim(ev)= n)
+
+ The input and output matrices are related by
+
+ A = E*D*E~ where D is the diagonal matrix of eigenvalues
+ D[i,j] = ev[i] if i=j and 0 otherwise.
+
+ The columns of E are the eigenvectors.
+
+ ---------------------------------------------------------------
+
+evmax
+
+ Compute the maximum (absolute) eigenvalue and corresponding
+ eigenvector of a real symmetric matrix A.
+
+ double evmax(double a[],double u[],int n)
+ double a[],u[]; int n;
+ a = array containing symmetric input matrix A
+ u = array containing the n components of the eigenvector
+ at exit (vector normalized to 1)
+ n = dimension of system
+ return: ev = eigenvalue of A with maximum absolute value
+ HUGE -> convergence failure
+
+
+-----------------------------------------------------------------------------
+ Eigensystem Auxiliaries:
+-----------------------------------------------------------------------------
+
+ The following routines are used by the eigensystem functions.
+ They are not normally called by the user.
+
+house
+
+ Transform a real symmetric matrix to tridiagonal form.
+
+ void house(double *a,double *d,double *dp,int n)
+ a = pointer to array of the symmetric input matrix A
+ These values are altered by the computation.
+ d = pointer to array of output diagonal elements
+ dp = pointer to array of n-1 elements neighboring the
+ diagonal in the symmetric transformed matrix
+ n = dimension (dim(a)= n*n, dim(d)=dim(dp)=n)
+
+ The output arrays are related to the tridiagonal matrix T by
+
+ T[i,i+1] = T[i+1,i] = dp[i] for i=0 to n-2, and
+ T[i,i] = d[i] for i=0 to n-1.
+ --------------------------------------------------------------
+
+housev
+
+ Transform a real symmetric matrix to tridiagonal form and
+ compute the orthogonal matrix of this transformation.
+
+ void housev(double *a,double *d,double *dp,int n)
+ a = pointer to array of symmetric input matrix A
+ The computation overloads this array with the
+ orthogonal transformation matrix O.
+ d = pointer to array of diagonal output elements
+ dp = pointer to array of n-1 elements neighboring the
+ diagonal in the symmetric transformed matrix
+ n = dimension (dim(a)= n*n, dim(d)=dim(dp)=n)
+
+ The orthogonal transformation matrix O satisfies O~*T*O = A.
+ ----------------------------------------------------------------
+
+qreval
+
+ Perform a QR reduction of a real symmetric tridiagonal
+ matrix to diagonal form.
+
+ int qreval(double *ev,double *dp,int n)
+ ev = pointer to array of input diagonal elements
+ The computation overloads this array with
+ eigenvalues.
+ dp = pointer to array input elements neighboring the
+ diagonal. This array is altered by the computation.
+ n = dimension (dim(ev)=dim(dp)= n)
+
+ ---------------------------------------------------------------
+
+qrevec
+
+ Perform a QR reduction of a real symmetric tridiagonal matrix
+ to diagonal form and update an orthogonal transformation matrix.
+
+ int qrevec(double *ev,double *evec,double *dp,int n)
+ ev = pointer to array of input diagonal elements
+ that the computation overloads with eigenvalues
+ evec = pointer array of orthogonal input matrix
+ This is updated by the computation to a matrix
+ of eigenvectors.
+ dp = pointer to array input elements neighboring the
+ diagonal. This array is altered by the computation.
+ n = dimension (dim(ev)=dim(dp)= n)
+
+ This function operates on the output of 'housev'.
+
+------------------------------------------------------------------------------
+
+ Matrix Utilities:
+------------------------------------------------------------------------------
+
+mmul
+
+ Multiply two real square matrices C = A * B.
+
+ void mmul(double *c,double *a,double *b,int n)
+ double *a,*b,*c; int n;
+ a = pointer to store for left product matrix
+ b = pointer to store for right product matrix
+ c = pointer to store for output matrix
+ n = dimension (dim(a)=dim(b)=dim(c)=n*n)
+
+ -------------------------------------------------------
+
+rmmult
+
+ Multiply two matrices Mat = A*B.
+
+ void rmmult(double *mat,double *a,double *b,int m,int k,int n)
+ double mat[],a[],b[]; int m,k,n;
+ mat = array containing m by n product matrix at exit
+ a = input array containing m by k matrix
+ b = input array containing k by n matrix
+ (all matrices stored in row order)
+ m,k,n = dimension parameters of arrays
+
+ ----------------------------------------------------------
+
+vmul
+
+ Multiply a vector by a matrix Vp = Mat*V.
+
+ void vmul(double *vp,double *mat,double *v,int n)
+ vp = pointer to array containing output vector
+ mat = pointer to array containing input matrix in row order
+ v = pointer to array containing input vector
+ n = dimension of vectors (mat is n by n)
+
+ ----------------------------------------------------------------
+
+vnrm
+
+ Compute the inner product of two real vectors, p = u~*v.
+
+ double vnrm(double *u,double *v,int n)
+ u = pointer to array of input vector u
+ v = pointer to array of input vector v
+ n = dimension (dim(u)=dim(v)=n)
+ return value: p = u~*v (dot product of u and v)
+
+ -----------------------------------------------------
+
+trnm
+
+ Transpose a real square matrix in place A -> A~.
+
+ void trnm(double *a,int n)
+ a = pointer to array of n by n input matrix A
+ This is overloaded by the transpose of A.
+ n = dimension (dim(a)=n*n)
+
+ ---------------------------------------------------------
+
+mattr
+
+ Transpose an m by n matrix A = B~.
+
+ void mattr(double *a,double *b,int m,int n)
+ a = pointer to array containing output n by m matrix
+ b = pointer to array containing input m by n matrix
+ (matrices stored in row order)
+ m,n = dimension parameters (dim(a)=dim(b)=n*m)
+
+ ------------------------------------------------------------
+
+otrma
+
+ Perform an orthogonal similarity transform C = A*B*A~.
+
+ void otrma(double *c,double *a,double *b,int n)
+ c = pointer to array of output matrix C
+ a = pointer to array of transformation A
+ b = pointer to array of input matrix B
+ n = dimension (dim(a)=dim(b)=dim(c)=n*n)
+
+ -----------------------------------------------------------
+
+otrsm
+
+ Perform a similarity transform on a symmetric matrix S = A*B*A~.
+
+ void otrsm(double *sm,double *a,double *b,int n)
+ sm = pointer to array of output matrix S
+ a = pointer to array of transformation matrix A
+ b = pointer to array of symmetric input matrix B
+ n = dimension (dim(a)=dim(b)=dim(sm)=n*n)
+
+ ---------------------------------------------------------------
+
+smgen
+
+ Construct a symmetric matrix from specified eigenvalues and
+ eigenvectors.
+
+ void smgen(double *a,double *eval,double *evec,int n)
+ a = pointer to array containing output matrix
+ eval = pointer to array containing the n eigenvalues
+ evec = pointer to array containing eigenvectors
+ (n by n with kth column the vector corresponding
+ to the kth eigenvalue)
+ n = system dimension
+
+ If D is the diagonal matrix of eigenvalues
+ and E[i,j] = evec[j+n*i] , then A = E*D*E~.
+
+ ----------------------------------------------------------------
+
+ortho
+
+ Generate a general orthogonal transformation matrix, E~*E = I.
+
+ void ortho(double *e,int n)
+ e = pointer to array of orthogonal output matrix E
+ n = dimension of vector space (dim(e)=n*n)
+
+ This function calls on the uniform random generator 'unfl' to
+ produce random rotation angles. Therefore this random generator
+ should be initialized by a call of 'setunfl' before calling
+ ortho (see Chapter 7).
+
+ -----------------------------------------------------------------
+
+mcopy
+
+ Copy an array a = b.
+
+ void mcopy(double *a,double *b,int n)
+ a = array containing output values, identical to input
+ b at exit
+ b = input array
+ n = dimension of arrays
+
+ -----------------------------------------------------------
+
+matprt
+
+ Print an array in n rows of m columns to stdout.
+
+ void matprt(double *a,int n,int m,char *fmt)
+ a = pointer to input array stored in row order (size = n*m)
+ n = number of output rows
+ m = number of output columns
+ fmt= pointer to character array containing format string
+ (printf formats eg. " %f")
+
+ Long rows may overflow the print line.
+
+ ---------------------------------------------------------------
+
+fmatprt
+
+ Print formatted array output to a file.
+
+ void fmatprt(FILE *fp,double *a,int n,int m,char *fmt)
+ fp = pointer to file opened for writing
+ a = pointer to input array stored in row order (size = n*m)
+ n = number of output rows
+ m = number of output columns
+ fmt= pounter to character array containing format string
+ (printf formats eg. " %f")
+
+------------------------------------------------------------------------------
+
+ Singular Value Decomposition:
+------------------------------------------------------------------------------
+
+ A number of versions of the Singular Value Decomposition (SVD)
+ are implemented in the library. They support the efficient
+ computation of this important factorization for a real m by n
+ matrix A. The general form of the SVD is
+
+ A = U*S*V~ with S = | D |
+ | 0 |
+
+ where U is an m by m orthogonal matrix, V is an n by n orthogonal matrix,
+ D is the n by n diagonal matrix of singular value, and S is the singular
+ m by n matrix produced by the transformation.
+
+ The singular values computed by these functions provide important
+ information on the rank of the matrix A, and on several matrix
+ norms of A. The number of non-zero singular values d[i] in D
+ equal to the rank of A. The two norm of A is
+
+ ||A|| = max(d[i]) , and the condition number is
+
+ k(A) = max(d[i])/min(d[i]) .
+
+ The Frobenius norm of the matrix A is
+
+ Fn(A) = Sum(i=0 to n-1) d[i]^2 .
+
+ Singular values consistent with zero are easily recognized, since
+ the decomposition algorithms have excellent numerical stability.
+ The value of a 'zero' d[i] is no larger than a few times the
+ computational rounding error e.
+
+ The matrix U1 is formed from the first n orthonormal column vectors
+ of U. U1[i,j] = U[i,j] for i = 1 to m and j = 1 to n. A singular
+ value decomposition of A can also be expressed in terms of the m by\
+ n matrix U1, with
+
+ A = U1*D*V~ .
+
+ SVD functions with three forms of output are provided. The first
+ form computes only the singular values, while the second computes
+ the singular values and the U and V orthogonal transformation
+ matrices. The third form of output computes singular values, the
+ V matrix, and saves space by overloading the input array with
+ the U1 matrix.
+
+ Two forms of decomposition algorithm are available for each of the
+ three output types. One is computationally efficient when m ~ n.
+ The second, distinguished by the prefix 'sv2' in the function name,
+ employs a two stage Householder reduction to accelerate computation
+ when m substantially exceeds n. Use of functions of the second form
+ is recommended for m > 2n.
+
+ Singular value output from each of the six SVD functions satisfies
+
+ d[i] >= 0 for i = 0 to n-1.
+-------------------------------------------------------------------------------
+
+svdval
+
+ Compute the singular values of a real m by n matrix A.
+
+ int svdval(double *d,double *a,int m,int n)
+ d = pointer to double array of dimension n
+ (output = singular values of A)
+ a = pointer to store of the m by n input matrix A
+ (A is altered by the computation)
+ m = number of rows in A
+ n = number of columns in A (m>=n required)
+ return value: status flag with:
+ 0 -> success
+ -1 -> input error m < n
+
+ ------------------------------------------------------------
+
+sv2val
+
+ Compute singular values when m >> n.
+
+ int sv2val(double *d,double *a,int m,int n)
+ d = pointer to double array of dimension n
+ (output = singular values of A)
+ a = pointer to store of the m by n input matrix A
+ (A is altered by the computation)
+ m = number of rows in A
+ n = number of columns in A (m>=n required)
+ return value: status flag with:
+ 0 -> success
+ -1 -> input error m < n
+ --------------------------------------------------------------
+
+svduv
+
+ Compute the singular value transformation S = U~*A*V.
+
+ int svduv(double *d,double *a,double *u,int m,double *v,int n)
+ d = pointer to double array of dimension n
+ (output = singular values of A)
+ a = pointer to store of the m by n input matrix A
+ (A is altered by the computation)
+ u = pointer to store for m by m orthogonal matrix U
+ v = pointer to store for n by n orthogonal matrix V
+ m = number of rows in A
+ n = number of columns in A (m>=n required)
+ return value: status flag with:
+ 0 -> success
+ -1 -> input error m < n
+
+ --------------------------------------------------------------------
+
+sv2uv
+
+ Compute the singular value transformation when m >> n.
+
+ int sv2uv(double *d,double *a,double *u,int m,double *v,int n)
+ d = pointer to double array of dimension n
+ (output = singular values of A)
+ a = pointer to store of the m by n input matrix A
+ (A is altered by the computation)
+ u = pointer to store for m by m orthogonal matrix U
+ v = pointer to store for n by n orthogonal matrix V
+ m = number of rows in A
+ n = number of columns in A (m>=n required)
+ return value: status flag with:
+ 0 -> success
+ -1 -> input error m < n
+
+ ----------------------------------------------------------------
+
+svdu1v
+
+ Compute the singular value transformation with A overloaded by
+ the partial U-matrix.
+
+ int svdu1v(double *d,double *a,int m,double *v,int n)
+ d = pointer to double array of dimension n
+ (output = singular values of A)
+ a = pointer to store of the m by n input matrix A
+ (At output a is overloaded by the matrix U1
+ whose n columns are orthogonal vectors equal to
+ the first n columns of U.)
+ v = pointer to store for n by n orthogonal matrix V
+ m = number of rows in A
+ n = number of columns in A (m>=n required)
+ return value: status flag with:
+ 0 -> success
+ -1 -> input error m < n
+
+ ------------------------------------------------------------------
+
+sv2u1v
+
+ Compute the singular value transformation with partial U
+ matrix U1 efficiently for m >> n.
+
+ #include <math.h>
+ int sv2u1v(d,a,m,v,n)
+ double *d,*a,*v; int m,n;
+ d = pointer to double array of dimension n
+ (output = singular values of A)
+ a = pointer to store of the m by n input matrix A
+ (At output a is overloaded by the matrix U1
+ whose n columns are orthogonal vectors equal to
+ the first n columns of U.)
+ v = pointer to store for n by n orthogonal matrix V
+ m = number of rows in A
+ n = number of columns in A (m>=n required)
+ return value: status flag with:
+ 0 -> success
+ -1 -> input error m < n
+
+-------------------------------------------------------------------------------
+ Auxiliary Functions used in SVD Computation:
+-------------------------------------------------------------------------------
+
+ The following routines are used by the singular value decomposition
+ functions. They are not normally called by the user.
+
+qrbdi
+
+ Perform a QR reduction of a bidiagonal matrix.
+
+ int qrbdi(double *d,double *e,int m)
+ d = pointer to n-dimensional array of diagonal values
+ (overloaded by diagonal elements of reduced matrix)
+ e = pointer to store of superdiagonal values (loaded in
+ first m-1 elements of the array). Values are altered
+ by the computation.
+ m = dimension of the d and e arrays
+ return value: N = number of QR iterations required
+
+ -------------------------------------------------------------
+
+qrbdv
+
+ Perform a QR reduction of a bidiagonal matrix and update the
+ the orthogonal transformation matrices U and V.
+
+ int qrbdv(double *d,double *e,double *u,int m,double *v,int n)
+ d = pointer to n-dimensional array of diagonal values
+ (overloaded by diagonal elements of reduced matrix)
+ e = pointer to store of superdiagonal values (loaded in
+ first m-1 elements of the array). Values are altered
+ by the computation.
+ u = pointer to store of m by m orthogonal matrix U updated
+ by the computation
+ v = pointer to store of n by n orthogonal matrix V updated
+ by the computation
+ m = dimension parameter of the U matrix
+ n = size of the d and e arrays and the number of rows and
+ columns in the V matrix
+ return value: N = number of QR iterations required
+
+ ---------------------------------------------------------------
+
+qrbd1
+
+ Perform a QR reduction of a bidiagonal matrix and update the
+ transformation matrices U1 and V.
+
+ int qrbdu1(double *d,double *e,double *u1,int m,double *v,int n)
+ d = pointer to n-dimensional array of diagonal values
+ (overloaded by diagonal elements of reduced matrix)
+ e = pointer to store of superdiagonal values (loaded in
+ first m-1 elements of the array). Values are altered
+ by the computation.
+ u1 = pointer to store of m by n transformation matrix U1
+ updated by the computation
+ v = pointer to store of n by n orthogonal matrix V updated
+ by the computation
+ m = number of rows in the U1 matrix
+ n = size of the d and e arrays, number of columns in the U1
+ matrix, and the number of rows and columns in the V matrix.
+ return value: N = number of QR iterations required
+
+ ------------------------------------------------------------------
+
+ldumat
+
+ Compute a left Householder transform matrix U from the vectors
+ specifying the Householder reflections.
+
+ void ldumat(double *a,double *u,int m,int n)
+ a = pointer to store of m by n input matrix A. Elements
+ of A on and below the main diagonal specify the
+ vectors of n Householder reflections (see note below).
+ u = pointer to store for the m by m orthogonal output
+ matrix U.
+ m = number of rows in A and U, and number of columns in U.
+ n = number of columns in A
+
+ --------------------------------------------------------------
+
+ldvmat
+
+ Compute a right Householder transform matrix from the vectors
+ specifying the Householder reflections.
+
+ void ldvmat(double *a,double *v,int n)
+ a = pointer to store of n by n input matrix A. Elements
+ of A on and above the superdiagonal specify vectors
+ of a sequence of Householder reflections (see note below).
+ v = pointer to store for the n by n orthogonal output
+ matrix V
+ n = number of rows and columns in A and V
+
+ -----------------------------------------------------------------
+
+atou1
+
+ Overload a Householder left-factored matrix A with the first
+ n columns of the Householder orthogonal matrix.
+
+ void atou1(double *a,int m,int n)
+ a = pointer to store of m by n input matrix A. Elements
+ of A on and below the main diagonal specify the
+ vectors of n Householder reflections (see note below).
+ This array is overloaded by the first n columns of the
+ Householder transformation matrix.
+ m = number of rows in A
+ n = number of columns in A
+
+ ---------------------------------------------------------------
+
+atovm
+
+ Overload a Householder right-factored square matrix A with the
+ Householder transformation matrix V.
+
+ void atovm(double *v,int n)
+ v = pointer to store for the n by n orthogonal
+ output matrix V
+ n = number of rows and columns in V
+
+-------------------------------------------------------------------------------
+
+ Individual Householder reflections are specified by a vector h.
+ The corresponding orthogonal reflection matrix is given by
+
+ H = I - c* h~ .
+
+ Input matrices store the vector, normalized to have its leading
+ coefficient equal to one, and the normalization factor
+
+ c = 2/(h~*h) .
+
+ Storage for the vectors is by column starting at the diagonal for
+ a left transform, and by row starting at the superdiagonal for a
+ right transform. The first location holds c followed by
+ components 2 to k of the vector.
+
+
+-------------------------------------------------------------------------------
+ Complex Linear Algebra
+-------------------------------------------------------------------------------
+
+ Solution and Inverse:
+-------------------------------------------------------------------------------
+
+csolv
+
+ Solve a complex linear system A*x = b.
+
+ int csolv(Cpx *a,Cpx *b,int n)
+ a = pointer to array of n by n system matrix A
+ The computation alters this array to a LU factorization.
+ b = pointer to input array of system vector b
+ This is replaced by the solution vector b -> x.
+ n = dimension of system (dim(a)=n*n, dim(b)=n)
+ return value: status flag with: 0 -> valid solution
+ 1 -> system singular
+
+ ---------------------------------------------------------------
+
+cminv
+
+ Invert a general complex matrix in place A -> Inv(A).
+
+ int cminv(Cpx *a,int n)
+ a = pointer to input array of complex n by n matrix A
+ The computation replaces A by its inverse.
+ n = dimension of system (dim(a)=n*n)
+ return value: status flag with: 0 -> valid solution
+ 1 -> system singular
+
+------------------------------------------------------------------------------
+
+ Hermitian Eigensystems:
+------------------------------------------------------------------------------
+
+heigval
+
+ Compute the eigenvalues of a Hermitian matrix.
+
+ void heigval(Cpx *a,double *ev,int n)
+ a = pointer to array for the Hermitian matrix H
+ These values are altered by the computation.
+ ev = pointer to array that is loaded with the
+ eigenvalues of H by the computation
+ n = dimension (dim(a)=n*n, dim(ev)=n)
+
+ --------------------------------------------------------
+
+heigvec
+
+ Compute the eigenvalues and eigenvectors of a Hermitian
+ matrix.
+
+ void heigvec(Cpx *a,double *ev,int n)
+ a = pointer to array for the hermitian matrix H
+ This array is loaded with a unitary matrix of
+ eigenvectors E by the computation.
+ ev = pointer to array that is loaded with the
+ eigenvalues of H by the computation
+ n = dimension (dim(a)=n*n, dim(ev)=n)
+
+ The eigen vector matrix output E satisfies
+
+ E^*E = I and A = E*D*E^
+
+ where D[i,j] = ev[i] for i=j and 0 otherwise
+ and E^ is the Hermitian conjugate of E.
+ The columns of E are the eigenvectors of A.
+
+ ----------------------------------------------------------
+
+hevmax
+
+ Compute the eigenvalue of maximum absolute value and
+ the corresponding eigenvector of a Hermitian matrix.
+
+ double hevmax(Cpx *a,Cpx *u,int n)
+ Cpx *a,*u; int n;
+ a = pointer to array for the Hermitian matrix H
+ u = pointer to array for the eigenvector umax
+ n = dimension (dim(a)=n*n, dim(u)=n)
+ return value: emax = eigenvalue of H with largest
+ absolute value
+
+ The eigenvector u and eigenvalue emax are related by u^*A*u = emax.
+
+------------------------------------------------------------------------------
+ Hermitian Eigensystem Auxiliaries:
+------------------------------------------------------------------------------
+
+ The following routines are called by the Hermitian eigensystem
+ functions. They are not normally called by the user.
+
+chouse
+
+ Transform a Hermitian matrix H to real symmetric tridiagonal
+ form.
+
+ void chouse(Cpx *a,double *d,double *dp,int n)
+ a = pointer to input array of complex matrix elements of H
+ This array is altered by the computation
+ d = pointer to output array of real diagonal elements
+ dp = pointer to output array of real superdiagonal elements
+ n = system dimension, with:
+ dim(a) = n * n, dim(d) = dim(dn) = n;
+
+ -----------------------------------------------------------------
+
+chousv
+
+ Transform a Hermitian matrix H to real symmetric tridiagonal
+ form, and compute the unitary matrix of this transformation.
+
+ void chousv(Cpx *a,double *d,double *dp,int n)
+ a = pointer to input array of complex matrix elements of H
+ The computation replaces this with the unitary matrix
+ U of the transformation.
+ d = pointer to output array of real diagonal elements
+ dp = pointer to output array of real superdiagonal elements
+ n = system dimension, with:
+ dim(a) = n*n, dim(d) = dim(dn) = n;
+
+ The matrix U satisfies
+
+ A = U^*T*U where T is real and tridiagonal,
+ with T[i,i+1] = T[i+1,i] = dp[i] and T[i,i] = d[i].
+
+ ------------------------------------------------------------
+
+qrecvc
+
+ Use QR transformations to reduce a real symmetric tridiagonal
+ matrix to diagonal form, and update a unitary transformation
+ matrix.
+
+ void qrecvc(double *ev,Cpx *evec,double *dp,int n)
+ ev = pointer to input array of diagonal elements
+ The computation transforms these to eigenvalues
+ of the input matrix.
+ evec = pointer to input array of a unitary transformation
+ matrix U. The computation applies the QR rotations
+ to this matrix.
+ dp = pointer to input array of elements neighboring the
+ diagonal. These values are altered by the computation.
+ n = dimension parameter (dim(ev)=dim(dp)=n, dim(evec)=n*n)
+
+ This function operates on the output of 'chousv'.
+
+-------------------------------------------------------------------------------
+
+ Complex Matrix Utilities:
+-------------------------------------------------------------------------------
+
+cvmul
+
+ Transform a complex vector u = A*v.
+
+ void cvmul(Cpx *u,Cpx *a,Cpx *v,int n)
+ u = pointer to array of output vector u.
+ a = pointer to array of transform matrix A.
+ v = pointer to array of input vector v.
+ n = dimension (dim(u)=dim(v)=n, dim(a)=n*n)
+
+ -----------------------------------------------------------
+
+cvnrm
+
+ Compute a Hermitian inner product s = u^*v.
+
+ Cpx cvnrm(Cpx *u,Cpx *v,int n)
+ u = pointer to array of first vector u
+ v = pointer to array of second vector v
+ n = dimension (dim(u)=dim(v)=n)
+ return value: s = complex value of inner product
+
+ -----------------------------------------------------------
+
+cmmul
+
+ Multiply two square complex matrices C = A * B.
+
+ void cmmul(Cpx *c,Cpx *a,Cpx *b,int n)
+ a = pointer to input array of left matrix factor A
+ b = pointer to input array of right matrix factor B
+ c = pointer to array of output product matrix C
+ n = dimension parameter (dim(c)=dim(a)=dim(b)=n*n)
+
+ -------------------------------------------------------------
+
+cmmult
+
+ Multiply two complex matrices C = A * B.
+
+ void cmmult(Cpx *c,Cpx *a,Cpx *b,int n,int m,int l)
+ a = pointer to input array of right n by m factor matrix A
+ b = pointer to input array of left m by l factor matrix B
+ c = pointer to store for n by l output matrix C
+ n,m,l = system dimension parameters, with
+ (dim(c)=n*l, dim(a)=n*m, dim(b)=m*l)
+
+ ----------------------------------------------------------------
+
+hconj
+
+ Compute the Hermitian conjugate in place, A -> A^.
+
+ void hconj(Cpx *a,int n)
+ a = pointer to input array for the complex matrix A
+ This is converted to the Hermitian conjugate A^.
+ n = dimension (dim(a)=n*n)
+
+ ----------------------------------------------------------
+
+utrncm
+
+ Perform a unitary similarity transformation C = T*B*T^.
+
+ void utrncm(Cpx *cm,Cpx *a,Cpx *b,int n)
+ a = pointer to the array of the transform matrix T
+ b = pointer to the array of the input matrix B
+ cm = pointer to output array of the transformed matrix C
+ n = dimension (dim(cm)=dim(a)=dim(b)=n*n)
+
+ ---------------------------------------------------------------
+
+utrnhm
+
+ Perform a unitary similarity transformation on a Hermitian
+ matrix H' = T*H*T^.
+
+ void utrnhm(Cpx *hm,Cpx *a,Cpx *b,int n)
+ a = pointer to the array of the transform matrix T
+ b = pointer to the array of the Hermitian input matrix H
+ hm = pointer to array containing Hermitian output matrix H'
+ n = dimension (dim(cm)=dim(a)=dim(b)=n*n)
+
+ -----------------------------------------------------------------
+
+trncm
+
+ Transpose a complex square matrix in place A -> A~.
+
+ void trncm(Cpx *a,int n)
+ a = pointer to array of n by n complex matrix A
+ The computation replaces A by its transpose
+ n = dimension (dim(a)=n*n)
+
+ ---------------------------------------------------------
+
+cmattr
+
+ Compute the transpose A = B~ of a complex m by n matrix.
+
+ void cmattr(Cpx *a,Cpx *b,int m,int n)
+ a = pointer to output array of matrix A
+ b = pointer to input array of matrix B
+ m, n = matrix dimensions, with B m by n and A n by m
+ (dim(a)=dim(b)= m*n)
+
+ -----------------------------------------------------------------
+
+hmgen
+
+ Generate a Hermitian matrix with specified eigen values and
+ eigenvectors.
+
+ void hmgen(Cpx *h,double *ev,Cpx *u,int n)
+ h = pointer to complex array of output matrix H
+ ev = pointer to real array of input eigen values
+ u = pointer to complex array of unitary matrix U
+ n = dimension (dim(h)=dim(u)=n*n, dim(ev)=n)
+
+ If D is a diagonal matrix with D[i,j] = ev[i] for i=j and 0
+ otherwise. H = U*D*U^. The columns of U are eigenvectors.
+
+ -----------------------------------------------------------------
+
+unitary
+
+ Generate a random unitary transformation U.
+
+ void unitary(Cpx *u,int n)
+ u = pointer to complex output array for U
+ n = dimension (dim(u)=n*n)
+
+
+ This function calls on the uniform random generator 'unfl' to
+ produce random rotation angles. Therefore this random generator
+ should be initialized by a call of 'setunfl' before calling
+ 'unitary' (see Chapter 7).
+
+ ---------------------------------------------------------------------
+
+cmcpy
+
+ Copy a complex array A = B.
+
+ void cmcpy(Cpx *a,Cpx *b,int n)
+ a = pointer to store for output array
+ b = pointer to store for input array
+ n = dimension of complex arrays A and B
+
+ -----------------------------------------------------------
+
+cmprt
+
+ Print rows of a complex matrix in a specified format.
+
+ void cmprt(Cpx *a,int m,int n,char *f)
+ a = pointer to array of complex m by n matrix
+ m = number of columns
+ n = number of rows
+ f = character array holding format for complex number
+ output (ie., "%f, %f ")
+
+ Long rows may overflow the print line.
Added: grass/branches/develbranch_6/lib/external/ccmath/Makefile
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/Makefile (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/Makefile 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,13 @@
+MODULE_TOPDIR = ../../..
+
+LIB_NAME = $(CCMATH_LIBNAME)
+include $(MODULE_TOPDIR)/include/Make/Lib.make
+
+default: $(ARCH_INCDIR)/ccmath_grass.h
+ $(MAKE) lib
+
+$(ARCH_INCDIR)/ccmath_grass.h: ccmath_grass.h
+ $(INSTALL_DATA) ccmath_grass.h $(ARCH_INCDIR)/ccmath_grass.h
+
+#doxygen:
+DOXNAME=ccmath
Added: grass/branches/develbranch_6/lib/external/ccmath/README
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/README (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/README 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,5 @@
+The code in this directory is a part of the
+ccmath library version 2.2.1.
+
+This code is licensed under the terms of the LGPL.
+See the lgpl.license file for details.
Added: grass/branches/develbranch_6/lib/external/ccmath/atou1.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/atou1.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/atou1.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,51 @@
+/* atou1.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include <stdlib.h>
+void atou1(double *a, int m, int n)
+{
+ double *p0, *p, *q, *w;
+
+ int i, j, k, mm;
+
+ double s, h;
+
+ w = (double *)calloc(m, sizeof(double));
+ p0 = a + n * n - 1;
+ i = n - 1;
+ mm = m - n;
+ if (mm == 0) {
+ *p0 = 1.;
+ p0 -= n + 1;
+ --i;
+ ++mm;
+ }
+ for (; i >= 0; --i, ++mm, p0 -= n + 1) {
+ if (*p0 != 0.) {
+ for (j = 0, p = p0 + n; j < mm; p += n)
+ w[j++] = *p;
+ h = *p0;
+ *p0 = 1. - h;
+ for (j = 0, p = p0 + n; j < mm; p += n)
+ *p = -h * w[j++];
+ for (k = i + 1, q = p0 + 1; k < n; ++k) {
+ for (j = 0, p = q + n, s = 0.; j < mm; p += n)
+ s += w[j++] * *p;
+ s *= h;
+ for (j = 0, p = q + n; j < mm; p += n)
+ *p -= s * w[j++];
+ *q++ = -s;
+ }
+ }
+ else {
+ *p0 = 1.;
+ for (j = 0, p = p0 + n, q = p0 + 1; j < mm; ++j, p += n)
+ *p = *q++ = 0.;
+ }
+ }
+ free(w);
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/atovm.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/atovm.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/atovm.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,43 @@
+/* atovm.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+void atovm(double *v, int n)
+{
+ double *p0, *q0, *p, *q, *qq;
+
+ double h, s;
+
+ int i, j, k, mm;
+
+ q0 = v + n * n - 1;
+ *q0 = 1.;
+ q0 -= n + 1;
+ p0 = v + n * n - n - n - 1;
+ for (i = n - 2, mm = 1; i >= 0; --i, p0 -= n + 1, q0 -= n + 1, ++mm) {
+ if (i && *(p0 - 1) != 0.) {
+ for (j = 0, p = p0, h = 1.; j < mm; ++j, ++p)
+ h += *p * *p;
+ h = *(p0 - 1);
+ *q0 = 1. - h;
+ for (j = 0, q = q0 + n, p = p0; j < mm; ++j, q += n)
+ *q = -h * *p++;
+ for (k = i + 1, q = q0 + 1; k < n; ++k) {
+ for (j = 0, qq = q + n, p = p0, s = 0.; j < mm; ++j, qq += n)
+ s += *qq * *p++;
+ s *= h;
+ for (j = 0, qq = q + n, p = p0; j < mm; ++j, qq += n)
+ *qq -= s * *p++;
+ *q++ = -s;
+ }
+ }
+ else {
+ *q0 = 1.;
+ for (j = 0, p = q0 + 1, q = q0 + n; j < mm; ++j, q += n)
+ *q = *p++ = 0.;
+ }
+ }
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/ccmath.h
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/ccmath.h (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/ccmath.h 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,181 @@
+/* ccmath.h CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ *
+ * Modified by Soeren gebbert 2009/01/08
+ * Removed al unued functions in grass. Only the linear algebra
+ * functions are used.
+ * ------------------------------------------------------------------------
+ */
+/*
+ CCM
+
+ Numerical Analysis Toolkit Header File
+ ELF Shared Library Version
+*/
+ /* Required for Shared Library */
+#ifndef _CCMATH_H_
+#define _CCMATH_H_
+#define XMATH 1
+
+ /* Define File Pointers and Standard Library */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+ /* Definitions of Types */
+
+#ifndef NULL
+#define NULL ((void *)0
+#endif
+
+ /* Complex Types */
+
+#ifndef CPX
+struct complex {double re,im;};
+typedef struct complex Cpx;
+#define CPX 1
+#endif
+
+/* Linear Algebra */
+
+
+ /* Real Linear Systems */
+
+
+ int minv(double *a,int n) ;
+
+ int psinv(double *v,int n) ;
+
+ int ruinv(double *a,int n) ;
+
+ int solv(double *a,double *b,int n) ;
+
+ int solvps(double *s,double *x,int n) ;
+
+ int solvru(double *a,double *b,int n) ;
+
+ void solvtd(double *a,double *b,double *c,double *x,int m) ;
+
+ void eigen(double *a,double *eval,int n) ;
+
+ void eigval(double *a,double *eval,int n) ;
+
+ double evmax(double *a,double *u,int n) ;
+
+ int svdval(double *d,double *a,int m,int n) ;
+
+ int sv2val(double *d,double *a,int m,int n) ;
+
+ int svduv(double *d,double *a,double *u,int m,double *v,int n) ;
+
+ int sv2uv(double *d,double *a,double *u,int m,double *v,int n) ;
+
+ int svdu1v(double *d,double *a,int m,double *v,int n) ;
+
+ int sv2u1v(double *d,double *a,int m,double *v,int n) ;
+
+ void mmul(double *mat,double *a,double *b,int n) ;
+
+ void rmmult(double *mat,double *a,double *b,int m,int k,int n) ;
+
+ void vmul(double *vp,double *mat,double *v,int n) ;
+
+ double vnrm(double *u,double *v,int n) ;
+
+ void matprt(double *a,int n,int m,char *fmt) ;
+
+ void fmatprt(FILE *fp,double *a,int n,int m,char *fmt) ;
+
+ void trnm(double *a,int n) ;
+
+ void mattr(double *a,double *b,int m,int n) ;
+
+ void otrma(double *at,double *u,double *a,int n) ;
+
+ void otrsm(double *st,double *u,double *s0,int n) ;
+
+ void mcopy(double *a,double *b,int m) ;
+
+ void ortho(double *evc,int n) ;
+
+ void smgen(double *a,double *eval,double *evec,int n) ;
+
+ /* utility routines for real symmertic eigensystems */
+
+ void house(double *a,double *d,double *ud,int n) ;
+
+ void housev(double *a,double *d,double *ud,int n) ;
+
+ int qreval(double *eval,double *ud,int n) ;
+
+ int qrevec(double *eval,double *evec,double *dp,int n) ;
+
+ /* utility routines for singular value decomposition */
+
+ int qrbdi(double *d, double *e,int n) ;
+
+ int qrbdv(double *d, double *e,double *u,int m,double *v,int n) ;
+
+ int qrbdu1(double *d, double *e,double *u,int m,double *v,int n) ;
+
+ void ldumat(double *a,double *u,int m,int n) ;
+
+ void ldvmat(double *a,double *v,int n) ;
+
+ void atou1(double *a,int m,int n) ;
+
+ void atovm(double *v,int n) ;
+
+
+ /* Complex Matrix Algebra */
+
+
+ int cminv(Cpx *a,int n) ;
+
+ int csolv(Cpx *a,Cpx *b,int n) ;
+
+ void heigvec(Cpx *a,double *eval,int n) ;
+
+ void heigval(Cpx *a,double *eval,int n) ;
+
+ double hevmax(Cpx *a,Cpx *u,int n) ;
+
+ void cmmul(Cpx *c,Cpx *a,Cpx *b,int n) ;
+
+ void cmmult(Cpx *c,Cpx *a,Cpx *b,int m,int k,int n) ;
+
+ void cvmul(Cpx *vp,Cpx *mat,Cpx *v,int n) ;
+
+ Cpx cvnrm(Cpx *u,Cpx *v,int n) ;
+
+ void cmprt(Cpx *a,int n,int m,char *fmt) ;
+
+ void trncm(Cpx *a,int n) ;
+
+ void hconj(Cpx *u,int n) ;
+
+ void cmattr(Cpx *a,Cpx *b,int m,int n) ;
+
+ void utrncm(Cpx *at,Cpx *u,Cpx *a,int n) ;
+
+ void utrnhm(Cpx *ht,Cpx *u,Cpx *h0,int n) ;
+
+ void cmcpy(Cpx *a,Cpx *b,int n) ;
+
+ void unitary(Cpx *u,int n) ;
+
+ void hmgen(Cpx *h,double *eval,Cpx *u,int n) ;
+
+
+ /* utility routines for hermitian eigen problems */
+
+ void chouse(Cpx *a,double *d,double *ud,int n) ;
+
+ void chousv(Cpx *a,double *d,double *ud,int n) ;
+
+ void qrecvc(double *eval,Cpx *evec,double *ud,int n) ;
+#endif
Added: grass/branches/develbranch_6/lib/external/ccmath/ccmath_grass.h
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/ccmath_grass.h (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/ccmath_grass.h 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,190 @@
+/* ccmath.h CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ *
+ * Modified by Soeren gebbert 2009/01/08.
+ *
+ * \attention Do not include this header file directly. Use the gmath.h file which
+ * includes ccmath_grass_wrapper.h header file, to wrap the function names.
+ * ------------------------------------------------------------------------
+ */
+/*
+ CCM
+
+ Numerical Analysis Toolkit Header File
+ ELF Shared Library Version
+*/
+ /* Required for Shared Library */
+#ifndef _CCMATH_H_
+#define _CCMATH_H_
+#define XMATH 1
+
+ /* Define File Pointers and Standard Library */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
+ /* Definitions of Types */
+
+#ifndef NULL
+#define NULL ((void *)0
+#endif
+
+ /* Complex Types */
+
+#ifndef CPX
+struct complex {double re,im;};
+typedef struct complex Cpx;
+#define CPX 1
+#endif
+
+/* Linear Algebra */
+
+
+ /* Real Linear Systems */
+
+
+ int minv(double *a,int n) ;
+
+ int psinv(double *v,int n) ;
+
+ int ruinv(double *a,int n) ;
+
+ int solv(double *a,double *b,int n) ;
+
+ int solvps(double *s,double *x,int n) ;
+
+ int solvru(double *a,double *b,int n) ;
+
+ void solvtd(double *a,double *b,double *c,double *x,int m) ;
+
+ void eigenvect(double *a,double *eval,int n) ; /*temporary commented out to avoid conflicts with the existing
+ eigen implementation of gmath lib. Soeren gebbert 10.08.2009*/
+
+ void eigval(double *a,double *eval,int n) ;
+
+ double evmax(double *a,double *u,int n) ;
+
+ int svdval(double *d,double *a,int m,int n) ;
+
+ int sv2val(double *d,double *a,int m,int n) ;
+
+ int svduv(double *d,double *a,double *u,int m,double *v,int n) ;
+
+ int sv2uv(double *d,double *a,double *u,int m,double *v,int n) ;
+
+ int svdu1v(double *d,double *a,int m,double *v,int n) ;
+
+ int sv2u1v(double *d,double *a,int m,double *v,int n) ;
+
+ void mmul(double *mat,double *a,double *b,int n) ;
+
+ void rmmult(double *mat,double *a,double *b,int m,int k,int n) ;
+
+ void vmul(double *vp,double *mat,double *v,int n) ;
+
+ double vnrm(double *u,double *v,int n) ;
+
+ void matprt(double *a,int n,int m,char *fmt) ;
+
+ void fmatprt(FILE *fp,double *a,int n,int m,char *fmt) ;
+
+ void trnm(double *a,int n) ;
+
+ void mattr(double *a,double *b,int m,int n) ;
+
+ void otrma(double *at,double *u,double *a,int n) ;
+
+ void otrsm(double *st,double *u,double *s0,int n) ;
+
+ void mcopy(double *a,double *b,int m) ;
+
+ void ortho(double *evc,int n) ;
+
+ void smgen(double *a,double *eval,double *evec,int n) ;
+
+ /* utility routines for real symmertic eigensystems */
+
+ void house(double *a,double *d,double *ud,int n) ;
+
+ void housev(double *a,double *d,double *ud,int n) ;
+
+ int qreval(double *eval,double *ud,int n) ;
+
+ int qrevec(double *eval,double *evec,double *dp,int n) ;
+
+ /* utility routines for singular value decomposition */
+
+ int qrbdi(double *d, double *e,int n) ;
+
+ int qrbdv(double *d, double *e,double *u,int m,double *v,int n) ;
+
+ int qrbdu1(double *d, double *e,double *u,int m,double *v,int n) ;
+
+ void ldumat(double *a,double *u,int m,int n) ;
+
+ void ldvmat(double *a,double *v,int n) ;
+
+ void atou1(double *a,int m,int n) ;
+
+ void atovm(double *v,int n) ;
+
+
+ /* Complex Matrix Algebra */
+
+
+ int cminv(Cpx *a,int n) ;
+
+ int csolv(Cpx *a,Cpx *b,int n) ;
+
+ void heigvec(Cpx *a,double *eval,int n) ;
+
+ void heigval(Cpx *a,double *eval,int n) ;
+
+ double hevmax(Cpx *a,Cpx *u,int n) ;
+
+ void cmmul(Cpx *c,Cpx *a,Cpx *b,int n) ;
+
+ void cmmult(Cpx *c,Cpx *a,Cpx *b,int m,int k,int n) ;
+
+ void cvmul(Cpx *vp,Cpx *mat,Cpx *v,int n) ;
+
+ Cpx cvnrm(Cpx *u,Cpx *v,int n) ;
+
+ void cmprt(Cpx *a,int n,int m,char *fmt) ;
+
+ void trncm(Cpx *a,int n) ;
+
+ void hconj(Cpx *u,int n) ;
+
+ void cmattr(Cpx *a,Cpx *b,int m,int n) ;
+
+ void utrncm(Cpx *at,Cpx *u,Cpx *a,int n) ;
+
+ void utrnhm(Cpx *ht,Cpx *u,Cpx *h0,int n) ;
+
+ void cmcpy(Cpx *a,Cpx *b,int n) ;
+
+ void unitary(Cpx *u,int n) ;
+
+ void hmgen(Cpx *h,double *eval,Cpx *u,int n) ;
+
+
+ /* utility routines for hermitian eigen problems */
+
+ void chouse(Cpx *a,double *d,double *ud,int n) ;
+
+ void chousv(Cpx *a,double *d,double *ud,int n) ;
+
+ void qrecvc(double *eval,Cpx *evec,double *ud,int n) ;
+
+ /* Simulation support */
+
+ double unfl() ;
+
+ void setunfl(unsigned int seed) ;
+
+#endif
Added: grass/branches/develbranch_6/lib/external/ccmath/chouse.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/chouse.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/chouse.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,96 @@
+/* chouse.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include <stdlib.h>
+#include "ccmath.h"
+void chouse(Cpx * a, double *d, double *dp, int n)
+{
+ double sc, x, y;
+
+ Cpx cc, u, *q0;
+
+ int i, j, k, m, e;
+
+ Cpx *qw, *pc, *p;
+
+ q0 = (Cpx *) calloc(2 * n, sizeof(Cpx));
+ for (i = 0, p = q0 + n, pc = a; i < n; ++i, pc += n + 1)
+ *p++ = *pc;
+ for (j = 0, pc = a; j < n - 2; ++j, pc += n + 1) {
+ m = n - j - 1;
+ for (i = 1, sc = 0.; i <= m; ++i)
+ sc += pc[i].re * pc[i].re + pc[i].im * pc[i].im;
+ if (sc > 0.) {
+ sc = sqrt(sc);
+ p = pc + 1;
+ y = sc + (x = sqrt(p->re * p->re + p->im * p->im));
+ if (x > 0.) {
+ cc.re = p->re / x;
+ cc.im = p->im / x;
+ }
+ else {
+ cc.re = 1.;
+ cc.im = 0.;
+ }
+ x = 1. / sqrt(2. * sc * y);
+ y *= x;
+ for (i = 0, qw = pc + 1; i < m; ++i) {
+ q0[i].re = q0[i].im = 0.;
+ if (i) {
+ qw[i].re *= x;
+ qw[i].im *= -x;
+ }
+ else {
+ qw[0].re = y * cc.re;
+ qw[0].im = -y * cc.im;
+ }
+ }
+ for (i = 0, e = j + 2, p = pc + n + 1, y = 0.; i < m;
+ ++i, p += e++) {
+ q0[i].re += (u.re = qw[i].re) * p->re - (u.im =
+ qw[i].im) * p->im;
+ q0[i].im += u.re * p->im + u.im * p->re;
+ ++p;
+ for (k = i + 1; k < m; ++k, ++p) {
+ q0[i].re += qw[k].re * p->re - qw[k].im * p->im;
+ q0[i].im += qw[k].im * p->re + qw[k].re * p->im;
+ q0[k].re += u.re * p->re + u.im * p->im;
+ q0[k].im += u.im * p->re - u.re * p->im;
+ }
+ y += u.re * q0[i].re + u.im * q0[i].im;
+ }
+ for (i = 0; i < m; ++i) {
+ q0[i].re -= y * qw[i].re;
+ q0[i].re += q0[i].re;
+ q0[i].im -= y * qw[i].im;
+ q0[i].im += q0[i].im;
+ }
+ for (i = 0, e = j + 2, p = pc + n + 1; i < m; ++i, p += e++) {
+ for (k = i; k < m; ++k, ++p) {
+ p->re -= qw[i].re * q0[k].re + qw[i].im * q0[k].im
+ + q0[i].re * qw[k].re + q0[i].im * qw[k].im;
+ p->im -= qw[i].im * q0[k].re - qw[i].re * q0[k].im
+ + q0[i].im * qw[k].re - q0[i].re * qw[k].im;
+ }
+ }
+ }
+ d[j] = pc->re;
+ dp[j] = sc;
+ }
+ d[j] = pc->re;
+ d[j + 1] = (pc + n + 1)->re;
+ u = *(pc + 1);
+ dp[j] = sqrt(u.re * u.re + u.im * u.im);
+ for (j = 0, pc = a, qw = q0 + n; j < n; ++j, pc += n + 1) {
+ *pc = qw[j];
+ for (i = 1, p = pc + n; i < n - j; ++i, p += n) {
+ pc[i].re = p->re;
+ pc[i].im = -p->im;
+ }
+ }
+ free(q0);
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/chousv.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/chousv.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/chousv.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,123 @@
+/* chousv.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include <stdlib.h>
+#include "ccmath.h"
+void chousv(Cpx * a, double *d, double *dp, int n)
+{
+ double sc, x, y;
+
+ Cpx cc, u, *qs;
+
+ int i, j, k, m, e;
+
+ Cpx *qw, *pc, *p, *q;
+
+ qs = (Cpx *) calloc(2 * n, sizeof(Cpx));
+ q = qs + n;
+ for (j = 0, pc = a; j < n - 2; ++j, pc += n + 1, ++q) {
+ m = n - j - 1;
+ for (i = 1, sc = 0.; i <= m; ++i)
+ sc += pc[i].re * pc[i].re + pc[i].im * pc[i].im;
+ if (sc > 0.) {
+ sc = sqrt(sc);
+ p = pc + 1;
+ y = sc + (x = sqrt(p->re * p->re + p->im * p->im));
+ if (x > 0.) {
+ cc.re = p->re / x;
+ cc.im = p->im / x;
+ }
+ else {
+ cc.re = 1.;
+ cc.im = 0.;
+ }
+ q->re = -cc.re;
+ q->im = -cc.im;
+ x = 1. / sqrt(2. * sc * y);
+ y *= x;
+ for (i = 0, qw = pc + 1; i < m; ++i) {
+ qs[i].re = qs[i].im = 0.;
+ if (i) {
+ qw[i].re *= x;
+ qw[i].im *= -x;
+ }
+ else {
+ qw[0].re = y * cc.re;
+ qw[0].im = -y * cc.im;
+ }
+ }
+ for (i = 0, e = j + 2, p = pc + n + 1, y = 0.; i < m;
+ ++i, p += e++) {
+ qs[i].re += (u.re = qw[i].re) * p->re - (u.im =
+ qw[i].im) * p->im;
+ qs[i].im += u.re * p->im + u.im * p->re;
+ ++p;
+ for (k = i + 1; k < m; ++k, ++p) {
+ qs[i].re += qw[k].re * p->re - qw[k].im * p->im;
+ qs[i].im += qw[k].im * p->re + qw[k].re * p->im;
+ qs[k].re += u.re * p->re + u.im * p->im;
+ qs[k].im += u.im * p->re - u.re * p->im;
+ }
+ y += u.re * qs[i].re + u.im * qs[i].im;
+ }
+ for (i = 0; i < m; ++i) {
+ qs[i].re -= y * qw[i].re;
+ qs[i].re += qs[i].re;
+ qs[i].im -= y * qw[i].im;
+ qs[i].im += qs[i].im;
+ }
+ for (i = 0, e = j + 2, p = pc + n + 1; i < m; ++i, p += e++) {
+ for (k = i; k < m; ++k, ++p) {
+ p->re -= qw[i].re * qs[k].re + qw[i].im * qs[k].im
+ + qs[i].re * qw[k].re + qs[i].im * qw[k].im;
+ p->im -= qw[i].im * qs[k].re - qw[i].re * qs[k].im
+ + qs[i].im * qw[k].re - qs[i].re * qw[k].im;
+ }
+ }
+ }
+ d[j] = pc->re;
+ dp[j] = sc;
+ }
+ d[j] = pc->re;
+ cc = *(pc + 1);
+ d[j + 1] = (pc += n + 1)->re;
+ dp[j] = sc = sqrt(cc.re * cc.re + cc.im * cc.im);
+ q->re = cc.re /= sc;
+ q->im = cc.im /= sc;
+ for (i = 0, m = n + n, p = pc; i < m; ++i, --p)
+ p->re = p->im = 0.;
+ pc->re = 1.;
+ (pc -= n + 1)->re = 1.;
+ qw = pc - n;
+ for (m = 2; m < n; ++m, qw -= n + 1) {
+ for (j = 0, p = pc, pc->re = 1.; j < m; ++j, p += n) {
+ for (i = 0, q = p, u.re = u.im = 0.; i < m; ++i, ++q) {
+ u.re += qw[i].re * q->re - qw[i].im * q->im;
+ u.im += qw[i].re * q->im + qw[i].im * q->re;
+ }
+ for (i = 0, q = p, u.re += u.re, u.im += u.im; i < m; ++i, ++q) {
+ q->re -= u.re * qw[i].re + u.im * qw[i].im;
+ q->im -= u.im * qw[i].re - u.re * qw[i].im;
+ }
+ }
+ for (i = 0, p = qw + m - 1; i < n; ++i, --p)
+ p->re = p->im = 0.;
+ (pc -= n + 1)->re = 1.;
+ }
+ for (j = 1, p = a + n + 1, q = qs + n, u.re = 1., u.im = 0.; j < n;
+ ++j, ++p, ++q) {
+ sc = u.re * q->re - u.im * q->im;
+ u.im = u.im * q->re + u.re * q->im;
+ u.re = sc;
+ for (i = 1; i < n; ++i, ++p) {
+ sc = u.re * p->re - u.im * p->im;
+ p->im = u.re * p->im + u.im * p->re;
+ p->re = sc;
+ }
+ }
+ free(qs);
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/cmattr.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/cmattr.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/cmattr.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,18 @@
+/* cmattr.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include "ccmath.h"
+void cmattr(Cpx * a, Cpx * b, int m, int n)
+{
+ Cpx *p;
+
+ int i, j;
+
+ for (i = 0; i < n; ++i, ++b)
+ for (j = 0, p = b; j < m; ++j, p += n)
+ *a++ = *p;
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/cmcpy.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/cmcpy.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/cmcpy.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,15 @@
+/* cmcpy.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include "ccmath.h"
+void cmcpy(Cpx * a, Cpx * b, int n)
+{
+ int i;
+
+ for (i = 0; i < n; ++i)
+ *a++ = *b++;
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/cminv.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/cminv.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/cminv.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,149 @@
+/* cminv.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include <stdlib.h>
+#include "ccmath.h"
+int cminv(Cpx * a, int n)
+{
+ int i, j, k, m, lc, *le;
+
+ Cpx *ps, *p, *q, *pa, *pd;
+
+ Cpx z, h, *q0;
+
+ double s, t, tq = 0., zr = 1.e-15;
+
+ le = (int *)calloc(n, sizeof(int));
+ q0 = (Cpx *) calloc(n, sizeof(Cpx));
+ pa = pd = a;
+ for (j = 0; j < n; ++j, ++pa, pd += n + 1) {
+ if (j > 0) {
+ for (i = 0, p = pa, q = q0; i < n; ++i, p += n)
+ *q++ = *p;
+ for (i = 1; i < n; ++i) {
+ lc = i < j ? i : j;
+ z.re = z.im = 0.;
+ for (k = 0, p = pa + i * n - j, q = q0; k < lc; ++k, ++q, ++p) {
+ z.re += p->re * q->re - p->im * q->im;
+ z.im += p->im * q->re + p->re * q->im;
+ }
+ q0[i].re -= z.re;
+ q0[i].im -= z.im;
+ }
+ for (i = 0, p = pa, q = q0; i < n; ++i, p += n)
+ *p = *q++;
+ }
+ s = fabs(pd->re) + fabs(pd->im);
+ lc = j;
+ for (k = j + 1, ps = pd; k < n; ++k) {
+ ps += n;
+ if ((t = fabs(ps->re) + fabs(ps->im)) > s) {
+ s = t;
+ lc = k;
+ }
+ }
+ tq = tq > s ? tq : s;
+ if (s < zr * tq) {
+ free(le - j);
+ free(q0);
+ return -1;
+ }
+ *le++ = lc;
+ if (lc != j) {
+ p = a + n * j;
+ q = a + n * lc;
+ for (k = 0; k < n; ++k, ++p, ++q) {
+ h = *p;
+ *p = *q;
+ *q = h;
+ }
+ }
+ t = pd->re * pd->re + pd->im * pd->im;
+ h.re = pd->re / t;
+ h.im = -(pd->im) / t;
+ for (k = j + 1, ps = pd + n; k < n; ++k, ps += n) {
+ z.re = ps->re * h.re - ps->im * h.im;
+ z.im = ps->im * h.re + ps->re * h.im;
+ *ps = z;
+ }
+ *pd = h;
+ }
+ for (j = 1, pd = ps = a; j < n; ++j) {
+ for (k = 0, pd += n + 1, q = ++ps; k < j; ++k, q += n) {
+ z.re = q->re * pd->re - q->im * pd->im;
+ z.im = q->im * pd->re + q->re * pd->im;
+ *q = z;
+ }
+ }
+ for (j = 1, pa = a; j < n; ++j) {
+ ++pa;
+ for (i = 0, q = q0, p = pa; i < j; ++i, p += n)
+ *q++ = *p;
+ for (k = 0; k < j; ++k) {
+ h.re = h.im = 0.;
+ for (i = k, p = pa + k * n + k - j, q = q0 + k; i < j; ++i) {
+ h.re -= p->re * q->re - p->im * q->im;
+ h.im -= p->im * q->re + p->re * q->im;
+ ++p;
+ ++q;
+ }
+ q0[k] = h;
+ }
+ for (i = 0, q = q0, p = pa; i < j; ++i, p += n)
+ *p = *q++;
+ }
+ for (j = n - 2, pd = pa = a + n * n - 1; j >= 0; --j) {
+ --pa;
+ pd -= n + 1;
+ for (i = 0, m = n - j - 1, q = q0, p = pd + n; i < m; ++i, p += n)
+ *q++ = *p;
+ for (k = n - 1, ps = pa; k > j; --k, ps -= n) {
+ z.re = -ps->re;
+ z.im = -ps->im;
+ for (i = j + 1, p = ps + 1, q = q0; i < k; ++i, ++p, ++q) {
+ z.re -= p->re * q->re - p->im * q->im;
+ z.im -= p->im * q->re + p->re * q->im;
+ }
+ q0[--m] = z;
+ }
+ for (i = 0, m = n - j - 1, q = q0, p = pd + n; i < m; ++i, p += n)
+ *p = *q++;
+ }
+ for (k = 0, pa = a; k < n - 1; ++k, ++pa) {
+ for (i = 0, q = q0, p = pa; i < n; ++i, p += n)
+ *q++ = *p;
+ for (j = 0, ps = a; j < n; ++j, ps += n) {
+ if (j > k) {
+ h.re = h.im = 0.;
+ p = ps + j;
+ i = j;
+ }
+ else {
+ h = q0[j];
+ p = ps + k + 1;
+ i = k + 1;
+ }
+ for (; i < n; ++i, ++p) {
+ h.re += p->re * q0[i].re - p->im * q0[i].im;
+ h.im += p->im * q0[i].re + p->re * q0[i].im;
+ }
+ q0[j] = h;
+ }
+ for (i = 0, q = q0, p = pa; i < n; ++i, p += n)
+ *p = *q++;
+ }
+ for (j = n - 2, le--; j >= 0; --j) {
+ for (k = 0, p = a + j, q = a + *(--le); k < n; ++k, p += n, q += n) {
+ h = *p;
+ *p = *q;
+ *q = h;
+ }
+ }
+ free(le);
+ free(q0);
+ return 0;
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/cmmul.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/cmmul.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/cmmul.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,28 @@
+/* cmmul.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include "ccmath.h"
+void cmmul(Cpx * c, Cpx * a, Cpx * b, int n)
+{
+ Cpx s, *p, *q;
+
+ int i, j, k;
+
+ trncm(b, n);
+ for (i = 0; i < n; ++i, a += n) {
+ for (j = 0, q = b; j < n; ++j) {
+ for (k = 0, p = a, s.re = s.im = 0.; k < n; ++k) {
+ s.re += p->re * q->re - p->im * q->im;
+ s.im += p->im * q->re + p->re * q->im;
+ ++p;
+ ++q;
+ }
+ *c++ = s;
+ }
+ }
+ trncm(b, n);
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/cmmult.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/cmmult.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/cmmult.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,29 @@
+/* cmmult.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include <stdlib.h>
+#include "ccmath.h"
+void cmmult(Cpx * cm, Cpx * a, Cpx * b, int n, int m, int l)
+{
+ Cpx z, *q0, *p, *q;
+
+ int i, j, k;
+
+ q0 = (Cpx *) calloc(m, sizeof(Cpx));
+ for (i = 0; i < l; ++i, ++cm) {
+ for (k = 0, p = b + i; k < m; p += l)
+ q0[k++] = *p;
+ for (j = 0, p = a, q = cm; j < n; ++j, q += l) {
+ for (k = 0, z.re = z.im = 0.; k < m; ++k, ++p) {
+ z.re += p->re * q0[k].re - p->im * q0[k].im;
+ z.im += p->im * q0[k].re + p->re * q0[k].im;
+ }
+ *q = z;
+ }
+ }
+ free(q0);
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/cmprt.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/cmprt.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/cmprt.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,20 @@
+/* cmprt.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include "ccmath.h"
+void cmprt(Cpx * a, int m, int n, char *f)
+{
+ int i, j;
+
+ Cpx *p;
+
+ for (i = 0, p = a; i < m; ++i) {
+ for (j = 0; j < n; ++j, ++p)
+ printf(f, p->re, p->im);
+ printf("\n");
+ }
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/csolv.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/csolv.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/csolv.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,102 @@
+/* csolv.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include <stdlib.h>
+#include "ccmath.h"
+int csolv(Cpx * a, Cpx * b, int n)
+{
+ int i, j, k, lc;
+
+ Cpx *ps, *p, *q, *pa, *pd;
+
+ Cpx z, h, *q0;
+
+ double s, t, tq = 0., zr = 1.e-15;
+
+ q0 = (Cpx *) calloc(n, sizeof(Cpx));
+ pa = a;
+ pd = a;
+ for (j = 0; j < n; ++j, ++pa, pd += n + 1) {
+ if (j > 0) {
+ for (i = 0, p = pa, q = q0; i < n; ++i, p += n)
+ *q++ = *p;
+ for (i = 1; i < n; ++i) {
+ lc = i < j ? i : j;
+ z.re = z.im = 0.;
+ for (k = 0, p = pa + i * n - j, q = q0; k < lc; ++k, ++q, ++p) {
+ z.re += p->re * q->re - p->im * q->im;
+ z.im += p->im * q->re + p->re * q->im;
+ }
+ q0[i].re -= z.re;
+ q0[i].im -= z.im;
+ }
+ for (i = 0, p = pa, q = q0; i < n; ++i, p += n)
+ *p = *q++;
+ }
+ s = fabs(pd->re) + fabs(pd->im);
+ lc = j;
+ for (k = j + 1, ps = pd; k < n; ++k) {
+ ps += n;
+ if ((t = fabs(ps->re) + fabs(ps->im)) > s) {
+ s = t;
+ lc = k;
+ }
+ }
+ tq = tq > s ? tq : s;
+ if (s < zr * tq) {
+ free(q0);
+ return -1;
+ }
+ if (lc != j) {
+ h = b[j];
+ b[j] = b[lc];
+ b[lc] = h;
+ p = a + n * j;
+ q = a + n * lc;
+ for (k = 0; k < n; ++k) {
+ h = *p;
+ *p++ = *q;
+ *q++ = h;
+ }
+ }
+ t = pd->re * pd->re + pd->im * pd->im;
+ h.re = pd->re / t;
+ h.im = -(pd->im) / t;
+ for (k = j + 1, ps = pd + n; k < n; ++k, ps += n) {
+ z.re = ps->re * h.re - ps->im * h.im;
+ z.im = ps->im * h.re + ps->re * h.im;
+ *ps = z;
+ }
+ }
+ for (j = 1, ps = b + 1; j < n; ++j, ++ps) {
+ for (k = 0, p = a + n * j, q = b, z.re = z.im = 0.; k < j; ++k) {
+ z.re += p->re * q->re - p->im * q->im;
+ z.im += p->im * q->re + p->re * q->im;
+ ++p;
+ ++q;
+ }
+ ps->re -= z.re;
+ ps->im -= z.im;
+ }
+ for (j = n - 1, --ps, pd = a + n * n - 1; j >= 0; --j, pd -= n + 1) {
+ for (k = j + 1, p = pd + 1, q = b + j + 1, z.re = z.im = 0.; k < n;
+ ++k) {
+ z.re += p->re * q->re - p->im * q->im;
+ z.im += p->im * q->re + p->re * q->im;
+ ++p;
+ ++q;
+ }
+ h.re = ps->re - z.re;
+ h.im = ps->im - z.im;
+ t = pd->re * pd->re + pd->im * pd->im;
+ ps->re = (h.re * pd->re + h.im * pd->im) / t;
+ ps->im = (h.im * pd->re - h.re * pd->im) / t;
+ --ps;
+ }
+ free(q0);
+ return 0;
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/cvmul.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/cvmul.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/cvmul.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,36 @@
+/* cvmul.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include "ccmath.h"
+void cvmul(Cpx * u, Cpx * a, Cpx * v, int n)
+{
+ Cpx *q;
+
+ int i, j;
+
+ for (i = 0; i < n; ++i, ++u) {
+ u->re = u->im = 0.;
+ for (j = 0, q = v; j < n; ++j, ++a, ++q) {
+ u->re += a->re * q->re - a->im * q->im;
+ u->im += a->im * q->re + a->re * q->im;
+ }
+ }
+}
+
+Cpx cvnrm(Cpx * u, Cpx * v, int n)
+{
+ int k;
+
+ Cpx z;
+
+ z.re = z.im = 0.;
+ for (k = 0; k < n; ++k, ++u, ++v) {
+ z.re += u->re * v->re + u->im * v->im;
+ z.im += u->re * v->im - u->im * v->re;
+ }
+ return z;
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/eigen.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/eigen.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/eigen.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,19 @@
+/* eigen.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include <stdlib.h>
+#include "ccmath.h"
+void eigen(double *a, double *ev, int n)
+{
+ double *dp;
+
+ dp = (double *)calloc(n, sizeof(double));
+ housev(a, ev, dp, n);
+ qrevec(ev, a, dp, n);
+ trnm(a, n);
+ free(dp);
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/eigval.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/eigval.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/eigval.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,18 @@
+/* eigval.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include <stdlib.h>
+#include "ccmath.h"
+void eigval(double *a, double *ev, int n)
+{
+ double *dp;
+
+ dp = (double *)calloc(n, sizeof(double));
+ house(a, ev, dp, n);
+ qreval(ev, dp, n);
+ free(dp);
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/evmax.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/evmax.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/evmax.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,47 @@
+/* evmax.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include <stdlib.h>
+#include "ccmath.h"
+double evmax(double *a, double *u, int n)
+{
+ double *p, *q, *qm, *r, *s, *t;
+
+ double ev, evm, c, h;
+
+ int kc;
+
+ q = (double *)calloc(n, sizeof(double));
+ qm = q + n;
+ *(qm - 1) = 1.;
+ ev = 0.;
+ for (kc = 0; kc < 200; ++kc) {
+ h = c = 0.;
+ evm = ev;
+ for (p = u, r = a, s = q; s < qm;) {
+ *p = 0.;
+ for (t = q; t < qm;)
+ *p += *r++ * *t++;
+ c += *p * *p;
+ h += *p++ * *s++;
+ }
+ ev = c / h;
+ c = sqrt(c);
+ for (p = u, s = q; s < qm;) {
+ *p /= c;
+ *s++ = *p++;
+ }
+ if (((c = ev - evm) < 0. ? -c : c) < 1.e-16 * (ev < 0. ? -ev : ev)) {
+ free(q);
+ return ev;
+ }
+ }
+ free(q);
+ for (kc = 0; kc < n;)
+ u[kc++] = 0.;
+ return 0.;
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/hconj.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/hconj.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/hconj.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,26 @@
+/* hconj.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include "ccmath.h"
+void hconj(Cpx * a, int n)
+{
+ Cpx s, *p, *q;
+
+ int i, j, e;
+
+ for (i = 0, e = n - 1; i < n; ++i, --e, a += n + 1) {
+ for (j = 0, p = a + 1, q = a + n; j < e; ++j) {
+ s = *p;
+ s.im = -s.im;
+ p->re = q->re;
+ (p++)->im = -q->im;
+ *q = s;
+ q += n;
+ }
+ a->im = -a->im;
+ }
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/heigval.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/heigval.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/heigval.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,18 @@
+/* heigval.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include <stdlib.h>
+#include "ccmath.h"
+void heigval(Cpx * a, double *ev, int n)
+{
+ double *dp;
+
+ dp = (double *)calloc(n, sizeof(double));
+ chouse(a, ev, dp, n);
+ qreval(ev, dp, n);
+ free(dp);
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/heigvec.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/heigvec.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/heigvec.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,19 @@
+/* heigvec.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include <stdlib.h>
+#include "ccmath.h"
+void heigvec(Cpx * a, double *ev, int n)
+{
+ double *dp;
+
+ dp = (double *)calloc(n, sizeof(double));
+ chousv(a, ev, dp, n);
+ qrecvc(ev, a, dp, n);
+ hconj(a, n);
+ free(dp);
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/hevmax.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/hevmax.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/hevmax.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,42 @@
+/* hevmax.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include <stdlib.h>
+#include "ccmath.h"
+double hevmax(Cpx * a, Cpx * u, int n)
+{
+ Cpx *x, *p, h;
+
+ double e, ep, s, t, te = 1.e-12;
+
+ int k, j;
+
+ x = (Cpx *) calloc(n, sizeof(Cpx));
+ x[0].re = 1.;
+ e = 0.;
+ do {
+ for (k = 0, p = a, s = t = 0.; k < n; ++k) {
+ for (j = 0, h.re = h.im = 0.; j < n; ++j, ++p) {
+ h.re += p->re * x[j].re - p->im * x[j].im;
+ h.im += p->im * x[j].re + p->re * x[j].im;
+ }
+ s += h.re * h.re + h.im * h.im;
+ t += h.re * x[k].re + h.im * x[k].im;
+ u[k] = h;
+ }
+ ep = e;
+ e = s / t;
+ s = 1. / sqrt(s);
+ for (k = 0; k < n; ++k) {
+ u[k].re *= s;
+ u[k].im *= s;
+ x[k] = u[k];
+ }
+ } while (fabs(e - ep) > fabs(te * e));
+ free(x);
+ return e;
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/hmgen.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/hmgen.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/hmgen.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,29 @@
+/* hmgen.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include <stdlib.h>
+#include "ccmath.h"
+void hmgen(Cpx * h, double *ev, Cpx * u, int n)
+{
+ Cpx *v, *p;
+
+ int i, j;
+
+ double e;
+
+ v = (Cpx *) calloc(n * n, sizeof(Cpx));
+ cmcpy(v, u, n * n);
+ hconj(v, n);
+ for (i = 0, p = v; i < n; ++i) {
+ for (j = 0, e = ev[i]; j < n; ++j, ++p) {
+ p->re *= e;
+ p->im *= e;
+ }
+ }
+ cmmul(h, u, v, n);
+ free(v);
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/house.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/house.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/house.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,73 @@
+/* house.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include <stdlib.h>
+#include "ccmath.h"
+void house(double *a, double *d, double *dp, int n)
+{
+ double sc, x, y, h;
+
+ int i, j, k, m, e;
+
+ double *qw, *qs, *pc, *p;
+
+ qs = (double *)calloc(2 * n, sizeof(double));
+ for (j = 0, qw = qs + n, pc = a; j < n; pc += n + 1)
+ qw[j++] = *pc;
+ for (j = 0, pc = a; j < n - 2; ++j, pc += n + 1) {
+ m = n - j - 1;
+ for (i = 1, sc = 0.; i <= m; ++i)
+ sc += pc[i] * pc[i];
+ if (sc > 0.) {
+ sc = sqrt(sc);
+ if ((x = *(pc + 1)) < 0.) {
+ y = x - sc;
+ h = 1. / sqrt(-2. * sc * y);
+ }
+ else {
+ y = x + sc;
+ h = 1. / sqrt(2. * sc * y);
+ sc = -sc;
+ }
+ for (i = 0, qw = pc + 1; i < m; ++i) {
+ qs[i] = 0.;
+ if (i)
+ qw[i] *= h;
+ else
+ qw[i] = y * h;
+ }
+ for (i = 0, e = j + 2, p = pc + n + 1, h = 0.; i < m;
+ ++i, p += e++) {
+ qs[i] += (y = qw[i]) * *p++;
+ for (k = i + 1; k < m; ++k) {
+ qs[i] += qw[k] * *p;
+ qs[k] += y * *p++;
+ }
+ h += y * qs[i];
+ }
+ for (i = 0; i < m; ++i) {
+ qs[i] -= h * qw[i];
+ qs[i] += qs[i];
+ }
+ for (i = 0, e = j + 2, p = pc + n + 1; i < m; ++i, p += e++) {
+ for (k = i; k < m; ++k)
+ *p++ -= qw[i] * qs[k] + qs[i] * qw[k];
+ }
+ }
+ d[j] = *pc;
+ dp[j] = sc;
+ }
+ d[j] = *pc;
+ dp[j] = *(pc + 1);
+ d[j + 1] = *(pc + n + 1);
+ for (j = 0, pc = a, qw = qs + n; j < n; ++j, pc += n + 1) {
+ *pc = qw[j];
+ for (i = 1, p = pc + n; i < n - j; p += n)
+ pc[i++] = *p;
+ }
+ free(qs);
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/housev.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/housev.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/housev.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,82 @@
+/* housev.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include <stdlib.h>
+#include "ccmath.h"
+void housev(double *a, double *d, double *dp, int n)
+{
+ double sc, x, y, h;
+
+ int i, j, k, m, e;
+
+ double *qw, *qs, *pc, *p;
+
+ qs = (double *)calloc(n, sizeof(double));
+ for (j = 0, pc = a; j < n - 2; ++j, pc += n + 1) {
+ m = n - j - 1;
+ for (i = 1, sc = 0.; i <= m; ++i)
+ sc += pc[i] * pc[i];
+ if (sc > 0.) {
+ sc = sqrt(sc);
+ if ((x = *(pc + 1)) < 0.) {
+ y = x - sc;
+ h = 1. / sqrt(-2. * sc * y);
+ }
+ else {
+ y = x + sc;
+ h = 1. / sqrt(2. * sc * y);
+ sc = -sc;
+ }
+ for (i = 0, qw = pc + 1; i < m; ++i) {
+ qs[i] = 0.;
+ if (i)
+ qw[i] *= h;
+ else
+ qw[i] = y * h;
+ }
+ for (i = 0, e = j + 2, p = pc + n + 1, h = 0.; i < m;
+ ++i, p += e++) {
+ qs[i] += (y = qw[i]) * *p++;
+ for (k = i + 1; k < m; ++k) {
+ qs[i] += qw[k] * *p;
+ qs[k] += y * *p++;
+ }
+ h += y * qs[i];
+ }
+ for (i = 0; i < m; ++i) {
+ qs[i] -= h * qw[i];
+ qs[i] += qs[i];
+ }
+ for (i = 0, e = j + 2, p = pc + n + 1; i < m; ++i, p += e++) {
+ for (k = i; k < m; ++k)
+ *p++ -= qw[i] * qs[k] + qs[i] * qw[k];
+ }
+ }
+ d[j] = *pc;
+ dp[j] = sc;
+ }
+ d[j] = *pc;
+ dp[j] = *(pc + 1);
+ d[j + 1] = *(pc += n + 1);
+ free(qs);
+ for (i = 0, m = n + n, p = pc; i < m; ++i)
+ *p-- = 0.;
+ *pc = 1.;
+ *(pc -= n + 1) = 1.;
+ qw = pc - n;
+ for (m = 2; m < n; ++m, qw -= n + 1) {
+ for (j = 0, p = pc, *pc = 1.; j < m; ++j, p += n) {
+ for (i = 0, qs = p, h = 0.; i < m;)
+ h += qw[i++] * *qs++;
+ for (i = 0, qs = p, h += h; i < m;)
+ *qs++ -= h * qw[i++];
+ }
+ for (i = 0, p = qw + m; i < n; ++i)
+ *(--p) = 0.;
+ *(pc -= n + 1) = 1.;
+ }
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/ldumat.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/ldumat.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/ldumat.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,57 @@
+/* ldumat.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include <stdlib.h>
+void ldumat(double *a, double *u, int m, int n)
+{
+ double *p0, *q0, *p, *q, *w;
+
+ int i, j, k, mm;
+
+ double s, h;
+
+ w = (double *)calloc(m, sizeof(double));
+ for (i = 0, mm = m * m, q = u; i < mm; ++i)
+ *q++ = 0.;
+ p0 = a + n * n - 1;
+ q0 = u + m * m - 1;
+ mm = m - n;
+ i = n - 1;
+ for (j = 0; j < mm; ++j, q0 -= m + 1)
+ *q0 = 1.;
+ if (mm == 0) {
+ p0 -= n + 1;
+ *q0 = 1.;
+ q0 -= m + 1;
+ --i;
+ ++mm;
+ }
+ for (; i >= 0; --i, ++mm, p0 -= n + 1, q0 -= m + 1) {
+ if (*p0 != 0.) {
+ for (j = 0, p = p0 + n, h = 1.; j < mm; p += n)
+ w[j++] = *p;
+ h = *p0;
+ *q0 = 1. - h;
+ for (j = 0, q = q0 + m; j < mm; q += m)
+ *q = -h * w[j++];
+ for (k = i + 1, q = q0 + 1; k < m; ++k) {
+ for (j = 0, p = q + m, s = 0.; j < mm; p += m)
+ s += w[j++] * *p;
+ s *= h;
+ for (j = 0, p = q + m; j < mm; p += m)
+ *p -= s * w[j++];
+ *q++ = -s;
+ }
+ }
+ else {
+ *q0 = 1.;
+ for (j = 0, p = q0 + 1, q = q0 + m; j < mm; ++j, q += m)
+ *q = *p++ = 0.;
+ }
+ }
+ free(w);
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/ldvmat.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/ldvmat.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/ldvmat.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,46 @@
+/* ldvmat.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+void ldvmat(double *a, double *v, int n)
+{
+ double *p0, *q0, *p, *q, *qq;
+
+ double h, s;
+
+ int i, j, k, mm;
+
+ for (i = 0, mm = n * n, q = v; i < mm; ++i)
+ *q++ = 0.;
+ *v = 1.;
+ q0 = v + n * n - 1;
+ *q0 = 1.;
+ q0 -= n + 1;
+ p0 = a + n * n - n - n - 1;
+ for (i = n - 2, mm = 1; i > 0; --i, p0 -= n + 1, q0 -= n + 1, ++mm) {
+ if (*(p0 - 1) != 0.) {
+ for (j = 0, p = p0, h = 1.; j < mm; ++j, ++p)
+ h += *p * *p;
+ h = *(p0 - 1);
+ *q0 = 1. - h;
+ for (j = 0, q = q0 + n, p = p0; j < mm; ++j, q += n)
+ *q = -h * *p++;
+ for (k = i + 1, q = q0 + 1; k < n; ++k) {
+ for (j = 0, qq = q + n, p = p0, s = 0.; j < mm; ++j, qq += n)
+ s += *qq * *p++;
+ s *= h;
+ for (j = 0, qq = q + n, p = p0; j < mm; ++j, qq += n)
+ *qq -= s * *p++;
+ *q++ = -s;
+ }
+ }
+ else {
+ *q0 = 1.;
+ for (j = 0, p = q0 + 1, q = q0 + n; j < mm; ++j, q += n)
+ *q = *p++ = 0.;
+ }
+ }
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/lgpl.license
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/lgpl.license (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/lgpl.license 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,513 @@
+
+ GNU LESSER GENERAL PUBLIC LICENSE
+ Version 2.1, February 1999
+
+ Copyright (C) 1991, 1999 Free Software Foundation, Inc.
+ 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+[This is the first released version of the Lesser GPL. It also counts
+ as the successor of the GNU Library Public License, version 2, hence
+ the version number 2.1.]
+
+ Preamble
+
+ The licenses for most software are designed to take away your
+freedom to share and change it. By contrast, the GNU General Public
+Licenses are intended to guarantee your freedom to share and change
+free software--to make sure the software is free for all its users.
+
+ This license, the Lesser General Public License, applies to some
+specially designated software packages--typically libraries--of the
+Free Software Foundation and other authors who decide to use it. You
+can use it too, but we suggest you first think carefully about whether
+this license or the ordinary General Public License is the better
+strategy to use in any particular case, based on the explanations
+below.
+
+ When we speak of free software, we are referring to freedom of use,
+not price. Our General Public Licenses are designed to make sure that
+you have the freedom to distribute copies of free software (and charge
+for this service if you wish); that you receive source code or can get
+it if you want it; that you can change the software and use pieces of
+it in new free programs; and that you are informed that you can do
+these things.
+
+ To protect your rights, we need to make restrictions that forbid
+distributors to deny you these rights or to ask you to surrender these
+rights. These restrictions translate to certain responsibilities for
+you if you distribute copies of the library or if you modify it.
+
+ For example, if you distribute copies of the library, whether gratis
+or for a fee, you must give the recipients all the rights that we gave
+you. You must make sure that they, too, receive or can get the source
+code. If you link other code with the library, you must provide
+complete object files to the recipients, so that they can relink them
+with the library after making changes to the library and recompiling
+it. And you must show them these terms so they know their rights.
+
+ We protect your rights with a two-step method: (1) we copyright the
+library, and (2) we offer you this license, which gives you legal
+permission to copy, distribute and/or modify the library.
+
+ To protect each distributor, we want to make it very clear that
+there is no warranty for the free library. Also, if the library is
+modified by someone else and passed on, the recipients should know
+that what they have is not the original version, so that the original
+author's reputation will not be affected by problems that might be
+introduced by others.
+^L
+ Finally, software patents pose a constant threat to the existence of
+any free program. We wish to make sure that a company cannot
+effectively restrict the users of a free program by obtaining a
+restrictive license from a patent holder. Therefore, we insist that
+any patent license obtained for a version of the library must be
+consistent with the full freedom of use specified in this license.
+
+ Most GNU software, including some libraries, is covered by the
+ordinary GNU General Public License. This license, the GNU Lesser
+General Public License, applies to certain designated libraries, and
+is quite different from the ordinary General Public License. We use
+this license for certain libraries in order to permit linking those
+libraries into non-free programs.
+
+ When a program is linked with a library, whether statically or using
+a shared library, the combination of the two is legally speaking a
+combined work, a derivative of the original library. The ordinary
+General Public License therefore permits such linking only if the
+entire combination fits its criteria of freedom. The Lesser General
+Public License permits more lax criteria for linking other code with
+the library.
+
+ We call this license the "Lesser" General Public License because it
+does Less to protect the user's freedom than the ordinary General
+Public License. It also provides other free software developers Less
+of an advantage over competing non-free programs. These disadvantages
+are the reason we use the ordinary General Public License for many
+libraries. However, the Lesser license provides advantages in certain
+special circumstances.
+
+ For example, on rare occasions, there may be a special need to
+encourage the widest possible use of a certain library, so that it
+becomes a de-facto standard. To achieve this, non-free programs must be
+allowed to use the library. A more frequent case is that a free
+library does the same job as widely used non-free libraries. In this
+case, there is little to gain by limiting the free library to free
+software only, so we use the Lesser General Public License.
+
+ In other cases, permission to use a particular library in non-free
+programs enables a greater number of people to use a large body of
+free software. For example, permission to use the GNU C Library in
+non-free programs enables many more people to use the whole GNU
+operating system, as well as its variant, the GNU/Linux operating
+system.
+
+ Although the Lesser General Public License is Less protective of the
+users' freedom, it does ensure that the user of a program that is
+linked with the Library has the freedom and the wherewithal to run
+that program using a modified version of the Library.
+
+ The precise terms and conditions for copying, distribution and
+modification follow. Pay close attention to the difference between a
+"work based on the library" and a "work that uses the library". The
+former contains code derived from the library, whereas the latter must
+be combined with the library in order to run.
+^L
+ GNU LESSER GENERAL PUBLIC LICENSE
+ TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
+
+ 0. This License Agreement applies to any software library or other
+program which contains a notice placed by the copyright holder or
+other authorized party saying it may be distributed under the terms of
+this Lesser General Public License (also called "this License").
+Each licensee is addressed as "you".
+
+ A "library" means a collection of software functions and/or data
+prepared so as to be conveniently linked with application programs
+(which use some of those functions and data) to form executables.
+
+ The "Library", below, refers to any such software library or work
+which has been distributed under these terms. A "work based on the
+Library" means either the Library or any derivative work under
+copyright law: that is to say, a work containing the Library or a
+portion of it, either verbatim or with modifications and/or translated
+straightforwardly into another language. (Hereinafter, translation is
+included without limitation in the term "modification".)
+
+ "Source code" for a work means the preferred form of the work for
+making modifications to it. For a library, complete source code means
+all the source code for all modules it contains, plus any associated
+interface definition files, plus the scripts used to control
+compilation and installation of the library.
+
+ Activities other than copying, distribution and modification are not
+covered by this License; they are outside its scope. The act of
+running a program using the Library is not restricted, and output from
+such a program is covered only if its contents constitute a work based
+on the Library (independent of the use of the Library in a tool for
+writing it). Whether that is true depends on what the Library does
+and what the program that uses the Library does.
+
+ 1. You may copy and distribute verbatim copies of the Library's
+complete source code as you receive it, in any medium, provided that
+you conspicuously and appropriately publish on each copy an
+appropriate copyright notice and disclaimer of warranty; keep intact
+all the notices that refer to this License and to the absence of any
+warranty; and distribute a copy of this License along with the
+Library.
+
+ You may charge a fee for the physical act of transferring a copy,
+and you may at your option offer warranty protection in exchange for a
+fee.
+
+ 2. You may modify your copy or copies of the Library or any portion
+of it, thus forming a work based on the Library, and copy and
+distribute such modifications or work under the terms of Section 1
+above, provided that you also meet all of these conditions:
+
+ a) The modified work must itself be a software library.
+
+ b) You must cause the files modified to carry prominent notices
+ stating that you changed the files and the date of any change.
+
+ c) You must cause the whole of the work to be licensed at no
+ charge to all third parties under the terms of this License.
+
+ d) If a facility in the modified Library refers to a function or a
+ table of data to be supplied by an application program that uses
+ the facility, other than as an argument passed when the facility
+ is invoked, then you must make a good faith effort to ensure that,
+ in the event an application does not supply such function or
+ table, the facility still operates, and performs whatever part of
+ its purpose remains meaningful.
+
+ (For example, a function in a library to compute square roots has
+ a purpose that is entirely well-defined independent of the
+ application. Therefore, Subsection 2d requires that any
+ application-supplied function or table used by this function must
+ be optional: if the application does not supply it, the square
+ root function must still compute square roots.)
+
+These requirements apply to the modified work as a whole. If
+identifiable sections of that work are not derived from the Library,
+and can be reasonably considered independent and separate works in
+themselves, then this License, and its terms, do not apply to those
+sections when you distribute them as separate works. But when you
+distribute the same sections as part of a whole which is a work based
+on the Library, the distribution of the whole must be on the terms of
+this License, whose permissions for other licensees extend to the
+entire whole, and thus to each and every part regardless of who wrote
+it.
+
+Thus, it is not the intent of this section to claim rights or contest
+your rights to work written entirely by you; rather, the intent is to
+exercise the right to control the distribution of derivative or
+collective works based on the Library.
+
+In addition, mere aggregation of another work not based on the Library
+with the Library (or with a work based on the Library) on a volume of
+a storage or distribution medium does not bring the other work under
+the scope of this License.
+
+ 3. You may opt to apply the terms of the ordinary GNU General Public
+License instead of this License to a given copy of the Library. To do
+this, you must alter all the notices that refer to this License, so
+that they refer to the ordinary GNU General Public License, version 2,
+instead of to this License. (If a newer version than version 2 of the
+ordinary GNU General Public License has appeared, then you can specify
+that version instead if you wish.) Do not make any other change in
+these notices.
+^L
+ Once this change is made in a given copy, it is irreversible for
+that copy, so the ordinary GNU General Public License applies to all
+subsequent copies and derivative works made from that copy.
+
+ This option is useful when you wish to copy part of the code of
+the Library into a program that is not a library.
+
+ 4. You may copy and distribute the Library (or a portion or
+derivative of it, under Section 2) in object code or executable form
+under the terms of Sections 1 and 2 above provided that you accompany
+it with the complete corresponding machine-readable source code, which
+must be distributed under the terms of Sections 1 and 2 above on a
+medium customarily used for software interchange.
+
+ If distribution of object code is made by offering access to copy
+from a designated place, then offering equivalent access to copy the
+source code from the same place satisfies the requirement to
+distribute the source code, even though third parties are not
+compelled to copy the source along with the object code.
+
+ 5. A program that contains no derivative of any portion of the
+Library, but is designed to work with the Library by being compiled or
+linked with it, is called a "work that uses the Library". Such a
+work, in isolation, is not a derivative work of the Library, and
+therefore falls outside the scope of this License.
+
+ However, linking a "work that uses the Library" with the Library
+creates an executable that is a derivative of the Library (because it
+contains portions of the Library), rather than a "work that uses the
+library". The executable is therefore covered by this License.
+Section 6 states terms for distribution of such executables.
+
+ When a "work that uses the Library" uses material from a header file
+that is part of the Library, the object code for the work may be a
+derivative work of the Library even though the source code is not.
+Whether this is true is especially significant if the work can be
+linked without the Library, or if the work is itself a library. The
+threshold for this to be true is not precisely defined by law.
+
+ If such an object file uses only numerical parameters, data
+structure layouts and accessors, and small macros and small inline
+functions (ten lines or less in length), then the use of the object
+file is unrestricted, regardless of whether it is legally a derivative
+work. (Executables containing this object code plus portions of the
+Library will still fall under Section 6.)
+
+ Otherwise, if the work is a derivative of the Library, you may
+distribute the object code for the work under the terms of Section 6.
+Any executables containing that work also fall under Section 6,
+whether or not they are linked directly with the Library itself.
+^L
+ 6. As an exception to the Sections above, you may also combine or
+link a "work that uses the Library" with the Library to produce a
+work containing portions of the Library, and distribute that work
+under terms of your choice, provided that the terms permit
+modification of the work for the customer's own use and reverse
+engineering for debugging such modifications.
+
+ You must give prominent notice with each copy of the work that the
+Library is used in it and that the Library and its use are covered by
+this License. You must supply a copy of this License. If the work
+during execution displays copyright notices, you must include the
+copyright notice for the Library among them, as well as a reference
+directing the user to the copy of this License. Also, you must do one
+of these things:
+
+ a) Accompany the work with the complete corresponding
+ machine-readable source code for the Library including whatever
+ changes were used in the work (which must be distributed under
+ Sections 1 and 2 above); and, if the work is an executable linked
+ with the Library, with the complete machine-readable "work that
+ uses the Library", as object code and/or source code, so that the
+ user can modify the Library and then relink to produce a modified
+ executable containing the modified Library. (It is understood
+ that the user who changes the contents of definitions files in the
+ Library will not necessarily be able to recompile the application
+ to use the modified definitions.)
+
+ b) Use a suitable shared library mechanism for linking with the
+ Library. A suitable mechanism is one that (1) uses at run time a
+ copy of the library already present on the user's computer system,
+ rather than copying library functions into the executable, and (2)
+ will operate properly with a modified version of the library, if
+ the user installs one, as long as the modified version is
+ interface-compatible with the version that the work was made with.
+
+ c) Accompany the work with a written offer, valid for at
+ least three years, to give the same user the materials
+ specified in Subsection 6a, above, for a charge no more
+ than the cost of performing this distribution.
+
+ d) If distribution of the work is made by offering access to copy
+ from a designated place, offer equivalent access to copy the above
+ specified materials from the same place.
+
+ e) Verify that the user has already received a copy of these
+ materials or that you have already sent this user a copy.
+
+ For an executable, the required form of the "work that uses the
+Library" must include any data and utility programs needed for
+reproducing the executable from it. However, as a special exception,
+the materials to be distributed need not include anything that is
+normally distributed (in either source or binary form) with the major
+components (compiler, kernel, and so on) of the operating system on
+which the executable runs, unless that component itself accompanies
+the executable.
+
+ It may happen that this requirement contradicts the license
+restrictions of other proprietary libraries that do not normally
+accompany the operating system. Such a contradiction means you cannot
+use both them and the Library together in an executable that you
+distribute.
+^L
+ 7. You may place library facilities that are a work based on the
+Library side-by-side in a single library together with other library
+facilities not covered by this License, and distribute such a combined
+library, provided that the separate distribution of the work based on
+the Library and of the other library facilities is otherwise
+permitted, and provided that you do these two things:
+
+ a) Accompany the combined library with a copy of the same work
+ based on the Library, uncombined with any other library
+ facilities. This must be distributed under the terms of the
+ Sections above.
+
+ b) Give prominent notice with the combined library of the fact
+ that part of it is a work based on the Library, and explaining
+ where to find the accompanying uncombined form of the same work.
+
+ 8. You may not copy, modify, sublicense, link with, or distribute
+the Library except as expressly provided under this License. Any
+attempt otherwise to copy, modify, sublicense, link with, or
+distribute the Library is void, and will automatically terminate your
+rights under this License. However, parties who have received copies,
+or rights, from you under this License will not have their licenses
+terminated so long as such parties remain in full compliance.
+
+ 9. You are not required to accept this License, since you have not
+signed it. However, nothing else grants you permission to modify or
+distribute the Library or its derivative works. These actions are
+prohibited by law if you do not accept this License. Therefore, by
+modifying or distributing the Library (or any work based on the
+Library), you indicate your acceptance of this License to do so, and
+all its terms and conditions for copying, distributing or modifying
+the Library or works based on it.
+
+ 10. Each time you redistribute the Library (or any work based on the
+Library), the recipient automatically receives a license from the
+original licensor to copy, distribute, link with or modify the Library
+subject to these terms and conditions. You may not impose any further
+restrictions on the recipients' exercise of the rights granted herein.
+You are not responsible for enforcing compliance by third parties with
+this License.
+^L
+ 11. If, as a consequence of a court judgment or allegation of patent
+infringement or for any other reason (not limited to patent issues),
+conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License. If you cannot
+distribute so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you
+may not distribute the Library at all. For example, if a patent
+license would not permit royalty-free redistribution of the Library by
+all those who receive copies directly or indirectly through you, then
+the only way you could satisfy both it and this License would be to
+refrain entirely from distribution of the Library.
+
+If any portion of this section is held invalid or unenforceable under
+any particular circumstance, the balance of the section is intended to
+apply, and the section as a whole is intended to apply in other
+circumstances.
+
+It is not the purpose of this section to induce you to infringe any
+patents or other property right claims or to contest validity of any
+such claims; this section has the sole purpose of protecting the
+integrity of the free software distribution system which is
+implemented by public license practices. Many people have made
+generous contributions to the wide range of software distributed
+through that system in reliance on consistent application of that
+system; it is up to the author/donor to decide if he or she is willing
+to distribute software through any other system and a licensee cannot
+impose that choice.
+
+This section is intended to make thoroughly clear what is believed to
+be a consequence of the rest of this License.
+
+ 12. If the distribution and/or use of the Library is restricted in
+certain countries either by patents or by copyrighted interfaces, the
+original copyright holder who places the Library under this License
+may add an explicit geographical distribution limitation excluding those
+countries, so that distribution is permitted only in or among
+countries not thus excluded. In such case, this License incorporates
+the limitation as if written in the body of this License.
+
+ 13. The Free Software Foundation may publish revised and/or new
+versions of the Lesser General Public License from time to time.
+Such new versions will be similar in spirit to the present version,
+but may differ in detail to address new problems or concerns.
+
+Each version is given a distinguishing version number. If the Library
+specifies a version number of this License which applies to it and
+"any later version", you have the option of following the terms and
+conditions either of that version or of any later version published by
+the Free Software Foundation. If the Library does not specify a
+license version number, you may choose any version ever published by
+the Free Software Foundation.
+^L
+ 14. If you wish to incorporate parts of the Library into other free
+programs whose distribution conditions are incompatible with these,
+write to the author to ask for permission. For software which is
+copyrighted by the Free Software Foundation, write to the Free
+Software Foundation; we sometimes make exceptions for this. Our
+decision will be guided by the two goals of preserving the free status
+of all derivatives of our free software and of promoting the sharing
+and reuse of software generally.
+
+ NO WARRANTY
+
+ 15. BECAUSE THE LIBRARY IS LICENSED FREE OF CHARGE, THERE IS NO
+WARRANTY FOR THE LIBRARY, TO THE EXTENT PERMITTED BY APPLICABLE LAW.
+EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR
+OTHER PARTIES PROVIDE THE LIBRARY "AS IS" WITHOUT WARRANTY OF ANY
+KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE
+LIBRARY IS WITH YOU. SHOULD THE LIBRARY PROVE DEFECTIVE, YOU ASSUME
+THE COST OF ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
+
+ 16. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN
+WRITING WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY
+AND/OR REDISTRIBUTE THE LIBRARY AS PERMITTED ABOVE, BE LIABLE TO YOU
+FOR DAMAGES, INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR
+CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OR INABILITY TO USE THE
+LIBRARY (INCLUDING BUT NOT LIMITED TO LOSS OF DATA OR DATA BEING
+RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD PARTIES OR A
+FAILURE OF THE LIBRARY TO OPERATE WITH ANY OTHER SOFTWARE), EVEN IF
+SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH
+DAMAGES.
+
+ END OF TERMS AND CONDITIONS
+^L
+ How to Apply These Terms to Your New Libraries
+
+ If you develop a new library, and you want it to be of the greatest
+possible use to the public, we recommend making it free software that
+everyone can redistribute and change. You can do so by permitting
+redistribution under these terms (or, alternatively, under the terms
+of the ordinary General Public License).
+
+ To apply these terms, attach the following notices to the library.
+It is safest to attach them to the start of each source file to most
+effectively convey the exclusion of warranty; and each file should
+have at least the "copyright" line and a pointer to where the full
+notice is found.
+
+
+ <one line to give the library's name and a brief idea of what it
+does.>
+ Copyright (C) <year> <name of author>
+
+ This library is free software; you can redistribute it and/or
+ modify it under the terms of the GNU Lesser 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
+ Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser 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
+
+Also add information on how to contact you by electronic and paper
+mail.
+
+You should also get your employer (if you work as a programmer) or
+your
+school, if any, to sign a "copyright disclaimer" for the library, if
+necessary. Here is a sample; alter the names:
+
+ Yoyodyne, Inc., hereby disclaims all copyright interest in the
+ library `Frob' (a library for tweaking knobs) written by James
+Random Hacker.
+
+ <signature of Ty Coon>, 1 April 1990
+ Ty Coon, President of Vice
+
+That's all there is to it!
+
+
Added: grass/branches/develbranch_6/lib/external/ccmath/matprt.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/matprt.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/matprt.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,33 @@
+/* matprt.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include <stdio.h>
+void matprt(double *a, int n, int m, char *fmt)
+{
+ int i, j;
+
+ double *p;
+
+ for (i = 0, p = a; i < n; ++i) {
+ for (j = 0; j < m; ++j)
+ printf(fmt, *p++);
+ printf("\n");
+ }
+}
+
+void fmatprt(FILE * fp, double *a, int n, int m, char *fmt)
+{
+ int i, j;
+
+ double *p;
+
+ for (i = 0, p = a; i < n; ++i) {
+ for (j = 0; j < m; ++j)
+ fprintf(fp, fmt, *p++);
+ fprintf(fp, "\n");
+ }
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/mattr.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/mattr.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/mattr.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,17 @@
+/* mattr.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+void mattr(double *a, double *b, int m, int n)
+{
+ double *p;
+
+ int i, j;
+
+ for (i = 0; i < n; ++i, ++b)
+ for (j = 0, p = b; j < m; ++j, p += n)
+ *a++ = *p;
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/mcopy.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/mcopy.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/mcopy.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,16 @@
+/* mcopy.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+void mcopy(double *a, double *b, int m)
+{
+ double *p, *q;
+
+ int k;
+
+ for (p = a, q = b, k = 0; k < m; ++k)
+ *p++ = *q++;
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/minv.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/minv.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/minv.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,123 @@
+/* minv.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include <stdlib.h>
+#include "ccmath.h"
+int minv(double *a, int n)
+{
+ int lc, *le;
+
+ double s, t, tq = 0., zr = 1.e-15;
+
+ double *pa, *pd, *ps, *p, *q, *q0;
+
+ int i, j, k, m;
+
+ le = (int *)malloc(n * sizeof(int));
+ q0 = (double *)malloc(n * sizeof(double));
+ for (j = 0, pa = pd = a; j < n; ++j, ++pa, pd += n + 1) {
+ if (j > 0) {
+ for (i = 0, q = q0, p = pa; i < n; ++i, p += n)
+ *q++ = *p;
+ for (i = 1; i < n; ++i) {
+ lc = i < j ? i : j;
+ for (k = 0, p = pa + i * n - j, q = q0, t = 0.; k < lc; ++k)
+ t += *p++ * *q++;
+ q0[i] -= t;
+ }
+ for (i = 0, q = q0, p = pa; i < n; ++i, p += n)
+ *p = *q++;
+ }
+ s = fabs(*pd);
+ lc = j;
+ for (k = j + 1, ps = pd; k < n; ++k) {
+ if ((t = fabs(*(ps += n))) > s) {
+ s = t;
+ lc = k;
+ }
+ }
+ tq = tq > s ? tq : s;
+ if (s < zr * tq) {
+ free(le - j);
+ free(q0);
+ return -1;
+ }
+ *le++ = lc;
+ if (lc != j) {
+ for (k = 0, p = a + n * j, q = a + n * lc; k < n; ++k) {
+ t = *p;
+ *p++ = *q;
+ *q++ = t;
+ }
+ }
+ for (k = j + 1, ps = pd, t = 1. / *pd; k < n; ++k)
+ *(ps += n) *= t;
+ *pd = t;
+ }
+ for (j = 1, pd = ps = a; j < n; ++j) {
+ for (k = 0, pd += n + 1, q = ++ps; k < j; ++k, q += n)
+ *q *= *pd;
+ }
+ for (j = 1, pa = a; j < n; ++j) {
+ ++pa;
+ for (i = 0, q = q0, p = pa; i < j; ++i, p += n)
+ *q++ = *p;
+ for (k = 0; k < j; ++k) {
+ t = 0.;
+ for (i = k, p = pa + k * n + k - j, q = q0 + k; i < j; ++i)
+ t -= *p++ * *q++;
+ q0[k] = t;
+ }
+ for (i = 0, q = q0, p = pa; i < j; ++i, p += n)
+ *p = *q++;
+ }
+ for (j = n - 2, pd = pa = a + n * n - 1; j >= 0; --j) {
+ --pa;
+ pd -= n + 1;
+ for (i = 0, m = n - j - 1, q = q0, p = pd + n; i < m; ++i, p += n)
+ *q++ = *p;
+ for (k = n - 1, ps = pa; k > j; --k, ps -= n) {
+ t = -(*ps);
+ for (i = j + 1, p = ps, q = q0; i < k; ++i)
+ t -= *++p * *q++;
+ q0[--m] = t;
+ }
+ for (i = 0, m = n - j - 1, q = q0, p = pd + n; i < m; ++i, p += n)
+ *p = *q++;
+ }
+ for (k = 0, pa = a; k < n - 1; ++k, ++pa) {
+ for (i = 0, q = q0, p = pa; i < n; ++i, p += n)
+ *q++ = *p;
+ for (j = 0, ps = a; j < n; ++j, ps += n) {
+ if (j > k) {
+ t = 0.;
+ p = ps + j;
+ i = j;
+ }
+ else {
+ t = q0[j];
+ p = ps + k + 1;
+ i = k + 1;
+ }
+ for (; i < n;)
+ t += *p++ * q0[i++];
+ q0[j] = t;
+ }
+ for (i = 0, q = q0, p = pa; i < n; ++i, p += n)
+ *p = *q++;
+ }
+ for (j = n - 2, le--; j >= 0; --j) {
+ for (k = 0, p = a + j, q = a + *(--le); k < n; ++k, p += n, q += n) {
+ t = *p;
+ *p = *q;
+ *q = t;
+ }
+ }
+ free(le);
+ free(q0);
+ return 0;
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/mmul.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/mmul.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/mmul.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,24 @@
+/* mmul.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include "ccmath.h"
+void mmul(double *c, double *a, double *b, int n)
+{
+ double *p, *q, s;
+
+ int i, j, k;
+
+ trnm(b, n);
+ for (i = 0; i < n; ++i, a += n) {
+ for (j = 0, q = b; j < n; ++j) {
+ for (k = 0, p = a, s = 0.; k < n; ++k)
+ s += *p++ * *q++;
+ *c++ = s;
+ }
+ }
+ trnm(b, n);
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/ortho.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/ortho.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/ortho.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,40 @@
+/* ortho.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include "ccmath.h"
+static double tpi = 6.28318530717958647;
+
+void ortho(double *e, int n)
+{
+ int i, j, k, m;
+
+ double *p, *q, c, s, a, unfl();
+
+ for (i = 0, p = e; i < n; ++i) {
+ for (j = 0; j < n; ++j) {
+ if (i == j)
+ *p++ = 1.;
+ else
+ *p++ = 0.;
+ }
+ }
+ for (i = 0, m = n - 1; i < m; ++i) {
+ for (j = i + 1; j < n; ++j) {
+ a = tpi * unfl();
+ c = cos(a);
+ s = sin(a);
+ p = e + n * i;
+ q = e + n * j;
+ for (k = 0; k < n; ++k) {
+ a = *p * c + *q * s;
+ *q = *q * c - *p * s;
+ *p++ = a;
+ ++q;
+ }
+ }
+ }
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/otrma.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/otrma.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/otrma.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,29 @@
+/* otrma.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include <stdlib.h>
+void otrma(double *c, double *a, double *b, int n)
+{
+ double z, *q0, *p, *s, *t;
+
+ int i, j, k;
+
+ q0 = (double *)calloc(n, sizeof(double));
+ for (i = 0; i < n; ++i, ++c) {
+ for (j = 0, t = b; j < n; ++j) {
+ for (k = 0, s = a + i * n, z = 0.; k < n; ++k)
+ z += *t++ * *s++;
+ q0[j] = z;
+ }
+ for (j = 0, p = c, t = a; j < n; ++j, p += n) {
+ for (k = 0, s = q0, z = 0.; k < n; ++k)
+ z += *t++ * *s++;
+ *p = z;
+ }
+ }
+ free(q0);
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/otrsm.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/otrsm.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/otrsm.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,31 @@
+/* otrsm.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include <stdlib.h>
+void otrsm(double *sm, double *a, double *b, int n)
+{
+ double z, *q0, *p, *s, *t;
+
+ int i, j, k;
+
+ q0 = (double *)calloc(n, sizeof(double));
+ for (i = 0; i < n; ++i) {
+ for (j = 0, t = b; j < n; ++j) {
+ for (k = 0, s = a + i * n, z = 0.; k < n; ++k)
+ z += *t++ * *s++;
+ q0[j] = z;
+ }
+ for (j = 0, p = sm + i, t = a; j <= i; ++j, p += n) {
+ for (k = 0, s = q0, z = 0.; k < n; ++k)
+ z += *t++ * *s++;
+ *p = z;
+ if (j < i)
+ sm[i * n + j] = z;
+ }
+ }
+ free(q0);
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/psinv.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/psinv.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/psinv.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,45 @@
+/* psinv.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include "ccmath.h"
+int psinv(double *v, int n)
+{
+ double z, *p, *q, *r, *s, *t;
+
+ int j, k;
+
+ for (j = 0, p = v; j < n; ++j, p += n + 1) {
+ for (q = v + j * n; q < p; ++q)
+ *p -= *q * *q;
+ if (*p <= 0.)
+ return -1;
+ *p = sqrt(*p);
+ for (k = j + 1, q = p + n; k < n; ++k, q += n) {
+ for (r = v + j * n, s = v + k * n, z = 0.; r < p;)
+ z += *r++ * *s++;
+ *q -= z;
+ *q /= *p;
+ }
+ }
+ trnm(v, n);
+ for (j = 0, p = v; j < n; ++j, p += n + 1) {
+ *p = 1. / *p;
+ for (q = v + j, t = v; q < p; t += n + 1, q += n) {
+ for (s = q, r = t, z = 0.; s < p; s += n)
+ z -= *s * *r++;
+ *q = z * *p;
+ }
+ }
+ for (j = 0, p = v; j < n; ++j, p += n + 1) {
+ for (q = v + j, t = p - j; q <= p; q += n) {
+ for (k = j, r = p, s = q, z = 0.; k < n; ++k)
+ z += *r++ * *s++;
+ *t++ = (*q = z);
+ }
+ }
+ return 0;
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/qrbdi.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/qrbdi.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/qrbdi.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,77 @@
+/* qrbdi.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include "ccmath.h"
+int qrbdi(double *dm, double *em, int m)
+{
+ int i, j, k, n;
+
+ double u, x, y, a, b, c, s, t;
+
+ for (j = 1, t = fabs(dm[0]); j < m; ++j)
+ if ((s = fabs(dm[j]) + fabs(em[j - 1])) > t)
+ t = s;
+ t *= 1.e-15;
+ n = 100 * m;
+ for (j = 0; m > 1 && j < n; ++j) {
+ for (k = m - 1; k > 0; --k) {
+ if (fabs(em[k - 1]) < t)
+ break;
+ if (fabs(dm[k - 1]) < t) {
+ for (i = k, s = 1., c = 0.; i < m; ++i) {
+ a = s * em[i - 1];
+ b = dm[i];
+ em[i - 1] *= c;
+ dm[i] = u = sqrt(a * a + b * b);
+ s = -a / u;
+ c = b / u;
+ }
+ break;
+ }
+ }
+ y = dm[k];
+ x = dm[m - 1];
+ u = em[m - 2];
+ a = (y + x) * (y - x) - u * u;
+ s = y * em[k];
+ b = s + s;
+ u = sqrt(a * a + b * b);
+ if (u > 0.) {
+ c = sqrt((u + a) / (u + u));
+ if (c != 0.)
+ s /= (c * u);
+ else
+ s = 1.;
+ for (i = k; i < m - 1; ++i) {
+ b = em[i];
+ if (i > k) {
+ a = s * em[i];
+ b *= c;
+ em[i - 1] = u = sqrt(x * x + a * a);
+ c = x / u;
+ s = a / u;
+ }
+ a = c * y + s * b;
+ b = c * b - s * y;
+ s *= dm[i + 1];
+ dm[i] = u = sqrt(a * a + s * s);
+ y = c * dm[i + 1];
+ c = a / u;
+ s /= u;
+ x = c * b + s * y;
+ y = c * y - s * b;
+ }
+ }
+ em[m - 2] = x;
+ dm[m - 1] = y;
+ if (fabs(x) < t)
+ --m;
+ if (m == k + 1)
+ --m;
+ }
+ return j;
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/qrbdu1.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/qrbdu1.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/qrbdu1.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,94 @@
+/* qrbdu1.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include "ccmath.h"
+int qrbdu1(double *dm, double *em, double *um, int mm, double *vm, int m)
+{
+ int i, j, k, n, jj, nm;
+
+ double u, x, y, a, b, c, s, t, w, *p, *q;
+
+ for (j = 1, t = fabs(dm[0]); j < m; ++j)
+ if ((s = fabs(dm[j]) + fabs(em[j - 1])) > t)
+ t = s;
+ t *= 1.e-15;
+ n = 100 * m;
+ nm = m;
+ for (j = 0; m > 1 && j < n; ++j) {
+ for (k = m - 1; k > 0; --k) {
+ if (fabs(em[k - 1]) < t)
+ break;
+ if (fabs(dm[k - 1]) < t) {
+ for (i = k, s = 1., c = 0.; i < m; ++i) {
+ a = s * em[i - 1];
+ b = dm[i];
+ em[i - 1] *= c;
+ dm[i] = u = sqrt(a * a + b * b);
+ s = -a / u;
+ c = b / u;
+ for (jj = 0, p = um + k - 1; jj < mm; ++jj, p += nm) {
+ q = p + i - k + 1;
+ w = c * *p + s * *q;
+ *q = c * *q - s * *p;
+ *p = w;
+ }
+ }
+ break;
+ }
+ }
+ y = dm[k];
+ x = dm[m - 1];
+ u = em[m - 2];
+ a = (y + x) * (y - x) - u * u;
+ s = y * em[k];
+ b = s + s;
+ u = sqrt(a * a + b * b);
+ if (u > 0.) {
+ c = sqrt((u + a) / (u + u));
+ if (c != 0.)
+ s /= (c * u);
+ else
+ s = 1.;
+ for (i = k; i < m - 1; ++i) {
+ b = em[i];
+ if (i > k) {
+ a = s * em[i];
+ b *= c;
+ em[i - 1] = u = sqrt(x * x + a * a);
+ c = x / u;
+ s = a / u;
+ }
+ a = c * y + s * b;
+ b = c * b - s * y;
+ for (jj = 0, p = vm + i; jj < nm; ++jj, p += nm) {
+ w = c * *p + s * *(p + 1);
+ *(p + 1) = c * *(p + 1) - s * *p;
+ *p = w;
+ }
+ s *= dm[i + 1];
+ dm[i] = u = sqrt(a * a + s * s);
+ y = c * dm[i + 1];
+ c = a / u;
+ s /= u;
+ x = c * b + s * y;
+ y = c * y - s * b;
+ for (jj = 0, p = um + i; jj < mm; ++jj, p += nm) {
+ w = c * *p + s * *(p + 1);
+ *(p + 1) = c * *(p + 1) - s * *p;
+ *p = w;
+ }
+ }
+ }
+ em[m - 2] = x;
+ dm[m - 1] = y;
+ if (fabs(x) < t)
+ --m;
+ if (m == k + 1)
+ --m;
+ }
+ return j;
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/qrbdv.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/qrbdv.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/qrbdv.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,94 @@
+/* qrbdv.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include "ccmath.h"
+int qrbdv(double *dm, double *em, double *um, int mm, double *vm, int m)
+{
+ int i, j, k, n, jj, nm;
+
+ double u, x, y, a, b, c, s, t, w, *p, *q;
+
+ for (j = 1, t = fabs(dm[0]); j < m; ++j)
+ if ((s = fabs(dm[j]) + fabs(em[j - 1])) > t)
+ t = s;
+ t *= 1.e-15;
+ n = 100 * m;
+ nm = m;
+ for (j = 0; m > 1 && j < n; ++j) {
+ for (k = m - 1; k > 0; --k) {
+ if (fabs(em[k - 1]) < t)
+ break;
+ if (fabs(dm[k - 1]) < t) {
+ for (i = k, s = 1., c = 0.; i < m; ++i) {
+ a = s * em[i - 1];
+ b = dm[i];
+ em[i - 1] *= c;
+ dm[i] = u = sqrt(a * a + b * b);
+ s = -a / u;
+ c = b / u;
+ for (jj = 0, p = um + k - 1; jj < mm; ++jj, p += mm) {
+ q = p + i - k + 1;
+ w = c * *p + s * *q;
+ *q = c * *q - s * *p;
+ *p = w;
+ }
+ }
+ break;
+ }
+ }
+ y = dm[k];
+ x = dm[m - 1];
+ u = em[m - 2];
+ a = (y + x) * (y - x) - u * u;
+ s = y * em[k];
+ b = s + s;
+ u = sqrt(a * a + b * b);
+ if (u != 0.) {
+ c = sqrt((u + a) / (u + u));
+ if (c != 0.)
+ s /= (c * u);
+ else
+ s = 1.;
+ for (i = k; i < m - 1; ++i) {
+ b = em[i];
+ if (i > k) {
+ a = s * em[i];
+ b *= c;
+ em[i - 1] = u = sqrt(x * x + a * a);
+ c = x / u;
+ s = a / u;
+ }
+ a = c * y + s * b;
+ b = c * b - s * y;
+ for (jj = 0, p = vm + i; jj < nm; ++jj, p += nm) {
+ w = c * *p + s * *(p + 1);
+ *(p + 1) = c * *(p + 1) - s * *p;
+ *p = w;
+ }
+ s *= dm[i + 1];
+ dm[i] = u = sqrt(a * a + s * s);
+ y = c * dm[i + 1];
+ c = a / u;
+ s /= u;
+ x = c * b + s * y;
+ y = c * y - s * b;
+ for (jj = 0, p = um + i; jj < mm; ++jj, p += mm) {
+ w = c * *p + s * *(p + 1);
+ *(p + 1) = c * *(p + 1) - s * *p;
+ *p = w;
+ }
+ }
+ }
+ em[m - 2] = x;
+ dm[m - 1] = y;
+ if (fabs(x) < t)
+ --m;
+ if (m == k + 1)
+ --m;
+ }
+ return j;
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/qrecvc.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/qrecvc.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/qrecvc.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,78 @@
+/* qrecvc.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include "ccmath.h"
+void qrecvc(double *ev, Cpx * evec, double *dp, int n)
+{
+ double cc, sc = 0.0, d, x = 0.0, y, h = 0.0, tzr = 1.e-15;
+
+ int i, j, k, m, nqr = 50 * n;
+
+ Cpx *p;
+
+ for (j = 0, m = n - 1; j < nqr; ++j) {
+ while (1) {
+ if (m < 1)
+ break;
+ k = m - 1;
+ if (fabs(dp[k]) <= fabs(ev[m]) * tzr)
+ --m;
+ else {
+ x = (ev[k] - ev[m]) / 2.;
+ h = sqrt(x * x + dp[k] * dp[k]);
+ if (m > 1 && fabs(dp[m - 2]) > fabs(ev[k]) * tzr)
+ break;
+ if ((cc = sqrt((1. + x / h) / 2.)) != 0.)
+ sc = dp[k] / (2. * cc * h);
+ else
+ sc = 1.;
+ x += ev[m];
+ ev[m--] = x - h;
+ ev[m--] = x + h;
+ for (i = 0, p = evec + n * (m + 1); i < n; ++i, ++p) {
+ h = p[0].re;
+ p[0].re = cc * h + sc * p[n].re;
+ p[n].re = cc * p[n].re - sc * h;
+ h = p[0].im;
+ p[0].im = cc * h + sc * p[n].im;
+ p[n].im = cc * p[n].im - sc * h;
+ }
+ }
+ }
+ if (x > 0.)
+ d = ev[m] + x - h;
+ else
+ d = ev[m] + x + h;
+ cc = 1.;
+ y = 0.;
+ ev[0] -= d;
+ for (k = 0; k < m; ++k) {
+ x = ev[k] * cc - y;
+ y = dp[k] * cc;
+ h = sqrt(x * x + dp[k] * dp[k]);
+ if (k > 0)
+ dp[k - 1] = sc * h;
+ ev[k] = cc * h;
+ cc = x / h;
+ sc = dp[k] / h;
+ ev[k + 1] -= d;
+ y *= sc;
+ ev[k] = cc * (ev[k] + y) + ev[k + 1] * sc * sc + d;
+ for (i = 0, p = evec + n * k; i < n; ++i, ++p) {
+ h = p[0].re;
+ p[0].re = cc * h + sc * p[n].re;
+ p[n].re = cc * p[n].re - sc * h;
+ h = p[0].im;
+ p[0].im = cc * h + sc * p[n].im;
+ p[n].im = cc * p[n].im - sc * h;
+ }
+ }
+ ev[k] = ev[k] * cc - y;
+ dp[k - 1] = ev[k] * sc;
+ ev[k] = ev[k] * cc + d;
+ }
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/qreval.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/qreval.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/qreval.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,59 @@
+/* qreval.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include "ccmath.h"
+int qreval(double *ev, double *dp, int n)
+{
+ double cc, sc = 0.0, d, x, y, h, tzr = 1.e-15;
+
+ int j, k, m, mqr = 8 * n;
+
+ for (j = 0, m = n - 1;; ++j) {
+ while (1) {
+ if (m < 1)
+ return 0;
+ k = m - 1;
+ if (fabs(dp[k]) <= fabs(ev[m]) * tzr)
+ --m;
+ else {
+ x = (ev[k] - ev[m]) / 2.;
+ h = sqrt(x * x + dp[k] * dp[k]);
+ if (m > 1 && fabs(dp[m - 2]) > fabs(ev[k]) * tzr)
+ break;
+ x += ev[m];
+ ev[m--] = x - h;
+ ev[m--] = x + h;
+ }
+ }
+ if (j > mqr)
+ return -1;
+ if (x > 0.)
+ d = ev[m] + x - h;
+ else
+ d = ev[m] + x + h;
+ cc = 1.;
+ y = 0.;
+ ev[0] -= d;
+ for (k = 0; k < m; ++k) {
+ x = ev[k] * cc - y;
+ y = dp[k] * cc;
+ h = sqrt(x * x + dp[k] * dp[k]);
+ if (k > 0)
+ dp[k - 1] = sc * h;
+ ev[k] = cc * h;
+ cc = x / h;
+ sc = dp[k] / h;
+ ev[k + 1] -= d;
+ y *= sc;
+ ev[k] = cc * (ev[k] + y) + ev[k + 1] * sc * sc + d;
+ }
+ ev[k] = ev[k] * cc - y;
+ dp[k - 1] = ev[k] * sc;
+ ev[k] = ev[k] * cc + d;
+ }
+ return 0;
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/qrevec.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/qrevec.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/qrevec.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,75 @@
+/* qrevec.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include <math.h>
+int qrevec(double *ev, double *evec, double *dp, int n)
+{
+ double cc, sc = 0.0, d, x, y, h, tzr = 1.e-15;
+
+ int i, j, k, m, mqr = 8 * n;
+
+ double *p;
+
+ for (j = 0, m = n - 1;; ++j) {
+ while (1) {
+ if (m < 1)
+ return 0;
+ k = m - 1;
+ if (fabs(dp[k]) <= fabs(ev[m]) * tzr)
+ --m;
+ else {
+ x = (ev[k] - ev[m]) / 2.;
+ h = sqrt(x * x + dp[k] * dp[k]);
+ if (m > 1 && fabs(dp[m - 2]) > fabs(ev[k]) * tzr)
+ break;
+ if ((cc = sqrt((1. + x / h) / 2.)) != 0.)
+ sc = dp[k] / (2. * cc * h);
+ else
+ sc = 1.;
+ x += ev[m];
+ ev[m--] = x - h;
+ ev[m--] = x + h;
+ for (i = 0, p = evec + n * (m + 1); i < n; ++i, ++p) {
+ h = p[0];
+ p[0] = cc * h + sc * p[n];
+ p[n] = cc * p[n] - sc * h;
+ }
+ }
+ }
+ if (j > mqr)
+ return -1;
+ if (x > 0.)
+ d = ev[m] + x - h;
+ else
+ d = ev[m] + x + h;
+ cc = 1.;
+ y = 0.;
+ ev[0] -= d;
+ for (k = 0; k < m; ++k) {
+ x = ev[k] * cc - y;
+ y = dp[k] * cc;
+ h = sqrt(x * x + dp[k] * dp[k]);
+ if (k > 0)
+ dp[k - 1] = sc * h;
+ ev[k] = cc * h;
+ cc = x / h;
+ sc = dp[k] / h;
+ ev[k + 1] -= d;
+ y *= sc;
+ ev[k] = cc * (ev[k] + y) + ev[k + 1] * sc * sc + d;
+ for (i = 0, p = evec + n * k; i < n; ++i, ++p) {
+ h = p[0];
+ p[0] = cc * h + sc * p[n];
+ p[n] = cc * p[n] - sc * h;
+ }
+ }
+ ev[k] = ev[k] * cc - y;
+ dp[k - 1] = ev[k] * sc;
+ ev[k] = ev[k] * cc + d;
+ }
+ return 0;
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/rmmult.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/rmmult.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/rmmult.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,26 @@
+/* rmmult.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include <stdlib.h>
+void rmmult(double *rm, double *a, double *b, int n, int m, int l)
+{
+ double z, *q0, *p, *q;
+
+ int i, j, k;
+
+ q0 = (double *)calloc(m, sizeof(double));
+ for (i = 0; i < l; ++i, ++rm) {
+ for (k = 0, p = b + i; k < m; p += l)
+ q0[k++] = *p;
+ for (j = 0, p = a, q = rm; j < n; ++j, q += l) {
+ for (k = 0, z = 0.; k < m;)
+ z += *p++ * q0[k++];
+ *q = z;
+ }
+ }
+ free(q0);
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/ruinv.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/ruinv.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/ruinv.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,31 @@
+/* ruinv.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+int ruinv(double *a, int n)
+{
+ int j;
+
+ double fabs();
+
+ double tt, z, *p, *q, *r, *s, *t;
+
+ for (j = 0, tt = 0., p = a; j < n; ++j, p += n + 1)
+ if ((z = fabs(*p)) > tt)
+ tt = z;
+ tt *= 1.e-16;
+ for (j = 0, p = a; j < n; ++j, p += n + 1) {
+ if (fabs(*p) < tt)
+ return -1;
+ *p = 1. / *p;
+ for (q = a + j, t = a; q < p; t += n + 1, q += n) {
+ for (s = q, r = t, z = 0.; s < p; s += n)
+ z -= *s * *r++;
+ *q = z * *p;
+ }
+ }
+ return 0;
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/smgen.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/smgen.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/smgen.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,19 @@
+/* smgen.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+void smgen(double *a, double *eval, double *evec, int n)
+{
+ double *p, *q, *ps, *r, *s, *t, *v = evec + n * n;
+
+ for (ps = a, p = evec; p < v; p += n) {
+ for (q = evec; q < v; q += n, ++ps) {
+ *ps = 0.;
+ for (r = eval, s = p, t = q; r < eval + n;)
+ *ps += *r++ * *s++ * *t++;
+ }
+ }
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/solv.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/solv.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/solv.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,71 @@
+/* solv.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU general
+ * public license. ( See the gpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include <stdlib.h>
+#include "ccmath.h"
+int solv(double *a, double *b, int n)
+{
+ int i, j, k, lc;
+
+ double *ps, *p, *q, *pa, *pd;
+
+ double *q0, s, t, tq = 0., zr = 1.e-15;
+
+ q0 = (double *)calloc(n, sizeof(double));
+ for (j = 0, pa = a, pd = a; j < n; ++j, ++pa, pd += n + 1) {
+ if (j) {
+ for (i = 0, q = q0, p = pa; i < n; ++i, p += n)
+ *q++ = *p;
+ for (i = 1; i < n; ++i) {
+ lc = i < j ? i : j;
+ for (k = 0, p = pa + i * n - j, q = q0, t = 0.; k < lc; ++k)
+ t += *p++ * *q++;
+ q0[i] -= t;
+ }
+ for (i = 0, q = q0, p = pa; i < n; ++i, p += n)
+ *p = *q++;
+ }
+ s = fabs(*pd);
+ lc = j;
+ for (k = j + 1, ps = pd; k < n; ++k) {
+ if ((t = fabs(*(ps += n))) > s) {
+ s = t;
+ lc = k;
+ }
+ }
+ tq = tq > s ? tq : s;
+ if (s < zr * tq) {
+ free(q0);
+ return -1;
+ }
+ if (lc != j) {
+ t = b[j];
+ b[j] = b[lc];
+ b[lc] = t;
+ for (k = 0, p = a + n * j, q = a + n * lc; k < n; ++k) {
+ t = *p;
+ *p++ = *q;
+ *q++ = t;
+ }
+ }
+ for (k = j + 1, ps = pd, t = 1. / *pd; k < n; ++k)
+ *(ps += n) *= t;
+ }
+ for (j = 1, ps = b + 1; j < n; ++j) {
+ for (k = 0, p = a + n * j, q = b, t = 0.; k < j; ++k)
+ t += *p++ * *q++;
+ *ps++ -= t;
+ }
+ for (j = n - 1, --ps, pd = a + n * n - 1; j >= 0; --j, pd -= n + 1) {
+ for (k = j + 1, p = pd, q = b + j, t = 0.; k < n; ++k)
+ t += *++p * *++q;
+ *ps -= t;
+ *ps-- /= *pd;
+ }
+ free(q0);
+ return 0;
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/solv.s
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/solv.s (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/solv.s 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,387 @@
+ .file "solv2.c"
+ .version "01.01"
+gcc2_compiled.:
+.section .rodata
+ .align 4
+.LC0:
+ .long 0x9ee75616,0x3cd203af
+.text
+ .align 4
+.globl solv
+ .type solv, at function
+solv:
+ pushl %ebp
+ movl %esp,%ebp
+ subl $72,%esp
+ pushl %edi
+ pushl %esi
+ pushl %ebx
+ fldz
+ pushl $8
+ movl 16(%ebp),%edx
+ pushl %edx
+ fstpt -60(%ebp)
+ call calloc
+ movl %eax,-20(%ebp)
+ movl $0,-4(%ebp)
+ movl 8(%ebp),%ecx
+ movl %ecx,-12(%ebp)
+ movl %ecx,-16(%ebp)
+ addl $8,%esp
+ fldt -60(%ebp)
+ movl 16(%ebp),%edi
+ cmpl %edi,-4(%ebp)
+ jge .L72
+ leal 0(,%edi,8),%edx
+ movl %edx,-24(%ebp)
+ addl $8,%edx
+ movl %edx,-32(%ebp)
+ movl $0,-40(%ebp)
+ movl 12(%ebp),%ecx
+ movl %ecx,-44(%ebp)
+ movl $0,-48(%ebp)
+ .align 4
+.L7:
+ cmpl $0,-4(%ebp)
+ je .L8
+ movl $0,-64(%ebp)
+ movl -20(%ebp),%edi
+ movl %edi,-72(%ebp)
+ movl -12(%ebp),%ebx
+ movl 16(%ebp),%edx
+ cmpl %edx,-64(%ebp)
+ jge .L10
+ .align 4
+.L12:
+ movl -72(%ebp),%ecx
+ movl (%ebx),%eax
+ movl %eax,(%ecx)
+ movl 4(%ebx),%eax
+ movl %eax,4(%ecx)
+ addl $8,%ecx
+ movl %ecx,-72(%ebp)
+ incl -64(%ebp)
+ addl -24(%ebp),%ebx
+ movl 16(%ebp),%edi
+ cmpl %edi,-64(%ebp)
+ jl .L12
+.L10:
+ movl $1,-64(%ebp)
+ movl 16(%ebp),%edx
+ cmpl %edx,-64(%ebp)
+ jge .L15
+ movl -48(%ebp),%ecx
+ movl %ecx,-28(%ebp)
+ movl -20(%ebp),%edi
+ addl $8,%edi
+ movl %edi,-68(%ebp)
+ movl %edx,-36(%ebp)
+ .align 4
+.L17:
+ movl -64(%ebp),%edx
+ movl %edx,-8(%ebp)
+ movl -4(%ebp),%ecx
+ cmpl %ecx,%edx
+ jle .L18
+ movl %ecx,-8(%ebp)
+.L18:
+ xorl %esi,%esi
+ movl -36(%ebp),%edi
+ movl -12(%ebp),%edx
+ leal (%edx,%edi,8),%eax
+ movl %eax,%ebx
+ subl -28(%ebp),%ebx
+ movl -20(%ebp),%ecx
+ movl %ecx,%edi
+ movl %ecx,-72(%ebp)
+ movl -8(%ebp),%ecx
+ fldz
+ cmpl %esi,%ecx
+ jle .L20
+ .align 4
+.L22:
+ fldl (%ebx)
+ fmull (%edi)
+ faddp %st,%st(1)
+ addl $8,%edi
+ addl $8,%ebx
+ incl %esi
+ cmpl %esi,%ecx
+ jg .L22
+.L20:
+ movl -72(%ebp),%ecx
+ movl -68(%ebp),%edx
+ fldl (%edx)
+ fsubp %st,%st(1)
+ fstpl (%edx)
+ addl $8,%edx
+ movl %edx,-68(%ebp)
+ movl 16(%ebp),%ecx
+ addl %ecx,-36(%ebp)
+ incl -64(%ebp)
+ cmpl %ecx,-64(%ebp)
+ jl .L17
+.L15:
+ movl $0,-64(%ebp)
+ movl -20(%ebp),%edi
+ movl %edi,-72(%ebp)
+ movl -12(%ebp),%ebx
+ movl 16(%ebp),%edx
+ cmpl %edx,-64(%ebp)
+ jge .L8
+ .align 4
+.L28:
+ movl -72(%ebp),%ecx
+ movl (%ecx),%eax
+ movl %eax,(%ebx)
+ movl 4(%ecx),%eax
+ movl %eax,4(%ebx)
+ addl $8,%ecx
+ movl %ecx,-72(%ebp)
+ incl -64(%ebp)
+ addl -24(%ebp),%ebx
+ movl 16(%ebp),%edi
+ cmpl %edi,-64(%ebp)
+ jl .L28
+.L8:
+ movl -16(%ebp),%edx
+ fldl (%edx)
+ fabs
+ movl -4(%ebp),%ecx
+ movl %ecx,-8(%ebp)
+ movl %ecx,%esi
+ incl %esi
+ movl %edx,-68(%ebp)
+ cmpl %esi,16(%ebp)
+ jle .L31
+ .align 4
+.L33:
+ movl -24(%ebp),%edi
+ addl %edi,-68(%ebp)
+ movl -68(%ebp),%edx
+ fldl (%edx)
+ fabs
+ fcom %st(1)
+ fnstsw %ax
+ andb $69,%ah
+ jne .L73
+ fstp %st(1)
+ movl %esi,-8(%ebp)
+ jmp .L32
+ .align 4
+.L73:
+ fstp %st(0)
+.L32:
+ incl %esi
+ cmpl %esi,16(%ebp)
+ jg .L33
+.L31:
+ fld %st(0)
+ fxch %st(2)
+ fcom %st(1)
+ fnstsw %ax
+ andb $69,%ah
+ jne .L74
+ fstp %st(2)
+ jmp .L36
+ .align 4
+.L74:
+ fstp %st(0)
+.L36:
+ fldl .LC0
+ fmul %st(2),%st
+ fcompp
+ fnstsw %ax
+ andb $69,%ah
+ jne .L38
+ fstp %st(0)
+ movl -20(%ebp),%ecx
+ pushl %ecx
+ call free
+ movl $-1,%eax
+ jmp .L71
+ .align 4
+.L38:
+ movl -4(%ebp),%edi
+ cmpl %edi,-8(%ebp)
+ je .L39
+ movl -44(%ebp),%edx
+ fldl (%edx)
+ movl -8(%ebp),%ecx
+ movl 12(%ebp),%edi
+ movl (%edi,%ecx,8),%eax
+ movl %eax,(%edx)
+ movl 4(%edi,%ecx,8),%eax
+ movl %eax,4(%edx)
+ fstpl (%edi,%ecx,8)
+ xorl %esi,%esi
+ movl -40(%ebp),%edx
+ movl 8(%ebp),%ecx
+ leal (%ecx,%edx,8),%ebx
+ movl 16(%ebp),%eax
+ imull -8(%ebp),%eax
+ leal (%ecx,%eax,8),%eax
+ movl %eax,-72(%ebp)
+ cmpl %esi,16(%ebp)
+ jle .L39
+ .align 4
+.L43:
+ fldl (%ebx)
+ movl -72(%ebp),%edi
+ movl (%edi),%eax
+ movl %eax,(%ebx)
+ movl 4(%edi),%eax
+ movl %eax,4(%ebx)
+ addl $8,%ebx
+ fstpl (%edi)
+ addl $8,%edi
+ movl %edi,-72(%ebp)
+ incl %esi
+ cmpl %esi,16(%ebp)
+ jg .L43
+.L39:
+ movl -4(%ebp),%esi
+ incl %esi
+ movl -16(%ebp),%edx
+ movl %edx,-68(%ebp)
+ fld1
+ fdivl (%edx)
+ cmpl %esi,16(%ebp)
+ jle .L75
+ .align 4
+.L48:
+ movl -24(%ebp),%ecx
+ addl %ecx,-68(%ebp)
+ movl -68(%ebp),%edi
+ fldl (%edi)
+ fmul %st(1),%st
+ fstpl (%edi)
+ incl %esi
+ cmpl %esi,16(%ebp)
+ jg .L48
+.L75:
+ fstp %st(0)
+ movl 16(%ebp),%edx
+ addl %edx,-40(%ebp)
+ addl $8,-44(%ebp)
+ addl $8,-48(%ebp)
+ incl -4(%ebp)
+ addl $8,-12(%ebp)
+ movl -32(%ebp),%ecx
+ addl %ecx,-16(%ebp)
+ cmpl %edx,-4(%ebp)
+ jl .L7
+.L72:
+ fstp %st(0)
+ movl $1,-4(%ebp)
+ movl 12(%ebp),%edi
+ addl $8,%edi
+ movl %edi,-68(%ebp)
+ movl 16(%ebp),%edx
+ cmpl %edx,-4(%ebp)
+ jge .L52
+ movl 16(%ebp),%eax
+ .align 4
+.L54:
+ xorl %esi,%esi
+ movl 8(%ebp),%ecx
+ leal (%ecx,%eax,8),%ebx
+ movl 12(%ebp),%edi
+ movl %edi,-72(%ebp)
+ fldz
+ cmpl %esi,-4(%ebp)
+ jle .L56
+ .align 4
+.L58:
+ fldl (%ebx)
+ movl -72(%ebp),%edx
+ fmull (%edx)
+ faddp %st,%st(1)
+ addl $8,%edx
+ movl %edx,-72(%ebp)
+ addl $8,%ebx
+ incl %esi
+ cmpl %esi,-4(%ebp)
+ jg .L58
+.L56:
+ movl -68(%ebp),%ecx
+ fldl (%ecx)
+ fsubp %st,%st(1)
+ fstpl (%ecx)
+ addl $8,%ecx
+ movl %ecx,-68(%ebp)
+ addl 16(%ebp),%eax
+ incl -4(%ebp)
+ movl 16(%ebp),%edi
+ cmpl %edi,-4(%ebp)
+ jl .L54
+.L52:
+ movl 16(%ebp),%edx
+ decl %edx
+ movl %edx,-4(%ebp)
+ addl $-8,-68(%ebp)
+ movl 16(%ebp),%eax
+ imull %eax,%eax
+ movl 8(%ebp),%ecx
+ leal -8(%ecx,%eax,8),%eax
+ movl %eax,-16(%ebp)
+ testl %edx,%edx
+ jl .L62
+ movl 16(%ebp),%edi
+ leal 8(,%edi,8),%edi
+ movl %edi,-64(%ebp)
+ leal 0(,%edx,8),%eax
+ .align 4
+.L64:
+ movl -4(%ebp),%esi
+ incl %esi
+ movl -16(%ebp),%ebx
+ movl 12(%ebp),%edx
+ addl %eax,%edx
+ movl %edx,-72(%ebp)
+ fldz
+ cmpl %esi,16(%ebp)
+ jle .L66
+ .align 4
+.L68:
+ addl $8,%ebx
+ addl $8,-72(%ebp)
+ fldl (%ebx)
+ movl -72(%ebp),%ecx
+ fmull (%ecx)
+ faddp %st,%st(1)
+ incl %esi
+ cmpl %esi,16(%ebp)
+ jg .L68
+.L66:
+ movl -68(%ebp),%edi
+ fldl (%edi)
+ fsubp %st,%st(1)
+ fstl (%edi)
+ movl -16(%ebp),%edx
+ fdivl (%edx)
+ fstpl (%edi)
+ addl $-8,%edi
+ movl %edi,-68(%ebp)
+ addl $-8,%eax
+ movl -64(%ebp),%ecx
+ subl %ecx,%edx
+ movl %edx,-16(%ebp)
+ decl -4(%ebp)
+ jns .L64
+.L62:
+ movl -20(%ebp),%edi
+ pushl %edi
+ call free
+ xorl %eax,%eax
+.L71:
+ leal -84(%ebp),%esp
+ popl %ebx
+ popl %esi
+ popl %edi
+ movl %ebp,%esp
+ popl %ebp
+ ret
+.Lfe1:
+ .size solv,.Lfe1-solv
+ .ident "GCC: (GNU) 2.7.2"
Added: grass/branches/develbranch_6/lib/external/ccmath/solvps.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/solvps.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/solvps.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,39 @@
+/* solvps.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include "ccmath.h"
+int solvps(double *a, double *b, int n)
+{
+ double *p, *q, *r, *s, t;
+
+ int j, k;
+
+ for (j = 0, p = a; j < n; ++j, p += n + 1) {
+ for (q = a + j * n; q < p; ++q)
+ *p -= *q * *q;
+ if (*p <= 0.)
+ return -1;
+ *p = sqrt(*p);
+ for (k = j + 1, q = p + n; k < n; ++k, q += n) {
+ for (r = a + j * n, s = a + k * n, t = 0.; r < p;)
+ t += *r++ * *s++;
+ *q -= t;
+ *q /= *p;
+ }
+ }
+ for (j = 0, p = a; j < n; ++j, p += n + 1) {
+ for (k = 0, q = a + j * n; k < j;)
+ b[j] -= b[k++] * *q++;
+ b[j] /= *p;
+ }
+ for (j = n - 1, p = a + n * n - 1; j >= 0; --j, p -= n + 1) {
+ for (k = j + 1, q = p + n; k < n; q += n)
+ b[j] -= b[k++] * *q;
+ b[j] /= *p;
+ }
+ return 0;
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/solvru.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/solvru.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/solvru.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,28 @@
+/* solvru.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+int solvru(double *a, double *b, int n)
+{
+ int j, k;
+
+ double fabs();
+
+ double s, t, *p, *q;
+
+ for (j = 0, s = 0., p = a; j < n; ++j, p += n + 1)
+ if ((t = fabs(*p)) > s)
+ s = t;
+ s *= 1.e-16;
+ for (j = n - 1, p = a + n * n - 1; j >= 0; --j, p -= n + 1) {
+ for (k = j + 1, q = p + 1; k < n;)
+ b[j] -= b[k++] * *q++;
+ if (fabs(*p) < s)
+ return -1;
+ b[j] /= *p;
+ }
+ return 0;
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/solvtd.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/solvtd.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/solvtd.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,23 @@
+/* solvtd.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+void solvtd(double *a, double *b, double *c, double *x, int m)
+{
+ double s;
+
+ int j;
+
+ for (j = 0; j < m; ++j) {
+ s = b[j] / a[j];
+ a[j + 1] -= s * c[j];
+ x[j + 1] -= s * x[j];
+ }
+ for (j = m, s = 0.; j >= 0; --j) {
+ x[j] -= s * c[j];
+ s = (x[j] /= a[j]);
+ }
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/sv2u1v.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/sv2u1v.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/sv2u1v.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,136 @@
+/* sv2u1v.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include <stdlib.h>
+#include "ccmath.h"
+int sv2u1v(double *d, double *a, int m, double *v, int n)
+{
+ double *p, *p1, *q, *pp, *w, *e;
+
+ double s, t, h, r, sv;
+
+ int i, j, k, mm, nm, ms;
+
+ if (m < n)
+ return -1;
+ w = (double *)calloc(m + n, sizeof(double));
+ e = w + m;
+ for (i = 0, mm = m, p = a; i < n; ++i, --mm, p += n + 1) {
+ if (mm > 1) {
+ sv = h = 0.;
+ for (j = 0, q = p, s = 0.; j < mm; ++j, q += n) {
+ w[j] = *q;
+ s += *q * *q;
+ }
+ if (s > 0.) {
+ h = sqrt(s);
+ if (*p < 0.)
+ h = -h;
+ s += *p * h;
+ s = 1. / s;
+ t = 1. / (w[0] += h);
+ sv = 1. + fabs(*p / h);
+ for (k = 1, ms = n - i; k < ms; ++k) {
+ for (j = 0, q = p + k, r = 0.; j < mm; q += n)
+ r += w[j++] * *q;
+ r = r * s;
+ for (j = 0, q = p + k; j < mm; q += n)
+ *q -= r * w[j++];
+ }
+ for (j = 1, q = p; j < mm;)
+ *(q += n) = w[j++] * t;
+ }
+ *p = sv;
+ d[i] = -h;
+ }
+ if (mm == 1)
+ d[i] = *p;
+ }
+ for (i = 0, q = v, p = a; i < n; ++i) {
+ for (j = 0; j < n; ++j, ++q, ++p) {
+ if (j < i)
+ *q = 0.;
+ else if (j == i)
+ *q = d[i];
+ else
+ *q = *p;
+ }
+ }
+ atou1(a, m, n);
+ for (i = 0, mm = n, nm = n - 1, p = v; i < n; ++i, --mm, --nm, p += n + 1) {
+ if (i && mm > 1) {
+ sv = h = 0.;
+ for (j = 0, q = p, s = 0.; j < mm; ++j, q += n) {
+ w[j] = *q;
+ s += *q * *q;
+ }
+ if (s > 0.) {
+ h = sqrt(s);
+ if (*p < 0.)
+ h = -h;
+ s += *p * h;
+ s = 1. / s;
+ t = 1. / (w[0] += h);
+ sv = 1. + fabs(*p / h);
+ for (k = 1, ms = n - i; k < ms; ++k) {
+ for (j = 0, q = p + k, r = 0.; j < mm; q += n)
+ r += w[j++] * *q;
+ for (j = 0, q = p + k, r *= s; j < mm; q += n)
+ *q -= r * w[j++];
+ }
+ for (k = 0, p1 = a + i; k < m; ++k, p1 += n) {
+ for (j = 0, q = p1, r = 0.; j < mm;)
+ r += w[j++] * *q++;
+ for (j = 0, q = p1, r *= s; j < mm;)
+ *q++ -= r * w[j++];
+ }
+ }
+ *p = sv;
+ d[i] = -h;
+ }
+ if (mm == 1)
+ d[i] = *p;
+ p1 = p + 1;
+ if (nm > 1) {
+ sv = h = 0.;
+ for (j = 0, q = p1, s = 0.; j < nm; ++j, ++q)
+ s += *q * *q;
+ if (s > 0.) {
+ h = sqrt(s);
+ if (*p1 < 0.)
+ h = -h;
+ sv = 1. + fabs(*p1 / h);
+ s += *p1 * h;
+ s = 1. / s;
+ t = 1. / (*p1 += h);
+ for (k = n, ms = n * (n - i); k < ms; k += n) {
+ for (j = 0, q = p1, pp = p1 + k, r = 0.; j < nm; ++j)
+ r += *q++ * *pp++;
+ for (j = 0, q = p1, pp = p1 + k, r *= s; j < nm; ++j)
+ *pp++ -= r * *q++;
+ }
+ for (j = 1, q = p1 + 1; j < nm; ++j)
+ *q++ *= t;
+ }
+ *p1 = sv;
+ e[i] = -h;
+ }
+ if (nm == 1)
+ e[i] = *p1;
+ }
+ atovm(v, n);
+ qrbdu1(d, e, a, m, v, n);
+ for (i = 0; i < n; ++i) {
+ if (d[i] < 0.) {
+ d[i] = -d[i];
+ for (j = 0, p = v + i; j < n; ++j, p += n)
+ *p = -*p;
+ }
+ }
+ free(w);
+ return 0;
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/sv2uv.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/sv2uv.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/sv2uv.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,134 @@
+/* sv2uv.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include <stdlib.h>
+#include "ccmath.h"
+int sv2uv(double *d, double *a, double *u, int m, double *v, int n)
+{
+ double *p, *p1, *q, *pp, *w, *e;
+
+ double s, t, h, r, sv;
+
+ int i, j, k, mm, nm, ms;
+
+ if (m < n)
+ return -1;
+ w = (double *)calloc(m + n, sizeof(double));
+ e = w + m;
+ for (i = 0, mm = m, p = a; i < n; ++i, --mm, p += n + 1) {
+ if (mm > 1) {
+ sv = h = 0.;
+ for (j = 0, q = p, s = 0.; j < mm; ++j, q += n) {
+ w[j] = *q;
+ s += *q * *q;
+ }
+ if (s > 0.) {
+ h = sqrt(s);
+ if (*p < 0.)
+ h = -h;
+ s += *p * h;
+ s = 1. / s;
+ t = 1. / (w[0] += h);
+ sv = 1. + fabs(*p / h);
+ for (k = 1, ms = n - i; k < ms; ++k) {
+ for (j = 0, q = p + k, r = 0.; j < mm; q += n)
+ r += w[j++] * *q;
+ r = r * s;
+ for (j = 0, q = p + k; j < mm; q += n)
+ *q -= r * w[j++];
+ }
+ for (j = 1, q = p; j < mm;)
+ *(q += n) = w[j++] * t;
+ }
+ *p = sv;
+ d[i] = -h;
+ }
+ if (mm == 1)
+ d[i] = *p;
+ }
+ ldumat(a, u, m, n);
+ for (i = 0, q = a; i < n; ++i) {
+ for (j = 0; j < n; ++j, ++q) {
+ if (j < i)
+ *q = 0.;
+ else if (j == i)
+ *q = d[i];
+ }
+ }
+ for (i = 0, mm = n, nm = n - 1, p = a; i < n; ++i, --mm, --nm, p += n + 1) {
+ if (i && mm > 1) {
+ sv = h = 0.;
+ for (j = 0, q = p, s = 0.; j < mm; ++j, q += n) {
+ w[j] = *q;
+ s += *q * *q;
+ }
+ if (s > 0.) {
+ h = sqrt(s);
+ if (*p < 0.)
+ h = -h;
+ s += *p * h;
+ s = 1. / s;
+ t = 1. / (w[0] += h);
+ sv = 1. + fabs(*p / h);
+ for (k = 1, ms = n - i; k < ms; ++k) {
+ for (j = 0, q = p + k, r = 0.; j < mm; q += n)
+ r += w[j++] * *q;
+ for (j = 0, q = p + k, r *= s; j < mm; q += n)
+ *q -= r * w[j++];
+ }
+ for (k = 0, p1 = u + i; k < m; ++k, p1 += m) {
+ for (j = 0, q = p1, r = 0.; j < mm;)
+ r += w[j++] * *q++;
+ for (j = 0, q = p1, r *= s; j < mm;)
+ *q++ -= r * w[j++];
+ }
+ }
+ *p = sv;
+ d[i] = -h;
+ }
+ if (mm == 1)
+ d[i] = *p;
+ p1 = p + 1;
+ if (nm > 1) {
+ sv = h = 0.;
+ for (j = 0, q = p1, s = 0.; j < nm; ++j, ++q)
+ s += *q * *q;
+ if (s > 0.) {
+ h = sqrt(s);
+ if (*p1 < 0.)
+ h = -h;
+ sv = 1. + fabs(*p1 / h);
+ s += *p1 * h;
+ s = 1. / s;
+ t = 1. / (*p1 += h);
+ for (k = n, ms = n * (n - i); k < ms; k += n) {
+ for (j = 0, q = p1, pp = p1 + k, r = 0.; j < nm; ++j)
+ r += *q++ * *pp++;
+ for (j = 0, q = p1, pp = p1 + k, r *= s; j < nm; ++j)
+ *pp++ -= r * *q++;
+ }
+ for (j = 1, q = p1 + 1; j < nm; ++j)
+ *q++ *= t;
+ }
+ *p1 = sv;
+ e[i] = -h;
+ }
+ if (nm == 1)
+ e[i] = *p1;
+ }
+ ldvmat(a, v, n);
+ qrbdv(d, e, u, m, v, n);
+ for (i = 0; i < n; ++i) {
+ if (d[i] < 0.) {
+ d[i] = -d[i];
+ for (j = 0, p = v + i; j < n; ++j, p += n)
+ *p = -*p;
+ }
+ }
+ free(w);
+ return 0;
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/sv2val.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/sv2val.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/sv2val.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,105 @@
+/* sv2val.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include <stdlib.h>
+#include "ccmath.h"
+int sv2val(double *d, double *a, int m, int n)
+{
+ double *p, *p1, *q, *w, *v;
+
+ double s, h, u;
+
+ int i, j, k, mm, nm, ms;
+
+ if (m < n)
+ return -1;
+ w = (double *)calloc(m, sizeof(double));
+ for (i = 0, mm = m, p = a; i < n && mm > 1; ++i, --mm, p += n + 1) {
+ for (j = 0, q = p, s = 0.; j < mm; ++j, q += n) {
+ w[j] = *q;
+ s += *q * *q;
+ }
+ if (s > 0.) {
+ h = sqrt(s);
+ if (*p < 0.)
+ h = -h;
+ s += *p * h;
+ s = 1. / s;
+ w[0] += h;
+ for (k = 1, ms = n - i; k < ms; ++k) {
+ for (j = 0, q = p + k, u = 0.; j < mm; q += n)
+ u += w[j++] * *q;
+ u = u * s;
+ for (j = 0, q = p + k; j < mm; q += n)
+ *q -= u * w[j++];
+ }
+ *p = -h;
+ }
+ }
+ for (i = 0, p = a; i < n; ++i, p += n) {
+ for (j = 0, q = p; j < i; ++j)
+ *q++ = 0.;
+ }
+ for (i = 0, mm = n, nm = n - 1, p = a; i < n; ++i, --mm, --nm, p += n + 1) {
+ if (i && mm > 1) {
+ for (j = 0, q = p, s = 0.; j < mm; ++j, q += n) {
+ w[j] = *q;
+ s += *q * *q;
+ }
+ if (s > 0.) {
+ h = sqrt(s);
+ if (*p < 0.)
+ h = -h;
+ s += *p * h;
+ s = 1. / s;
+ w[0] += h;
+ for (k = 1, ms = n - i; k < ms; ++k) {
+ for (j = 0, q = p + k, u = 0.; j < mm; q += n)
+ u += w[j++] * *q;
+ u *= s;
+ for (j = 0, q = p + k; j < mm; q += n)
+ *q -= u * w[j++];
+ }
+ *p = -h;
+ }
+ }
+ p1 = p + 1;
+ if (nm > 1) {
+ for (j = 0, q = p1, s = 0.; j < nm; ++j, ++q)
+ s += *q * *q;
+ if (s > 0.) {
+ h = sqrt(s);
+ if (*p1 < 0.)
+ h = -h;
+ s += *p1 * h;
+ s = 1. / s;
+ *p1 += h;
+ for (k = n, ms = n * (m - i); k < ms; k += n) {
+ for (j = 0, q = p1, v = p1 + k, u = 0.; j < nm; ++j)
+ u += *q++ * *v++;
+ u *= s;
+ for (j = 0, q = p1, v = p1 + k; j < nm; ++j)
+ *v++ -= u * *q++;
+ }
+ *p1 = -h;
+ }
+ }
+ }
+ for (j = 0, p = a; j < n; ++j, p += n + 1) {
+ d[j] = *p;
+ if (j < n - 1)
+ w[j] = *(p + 1);
+ else
+ w[j] = 0.;
+ }
+ qrbdi(d, w, n);
+ for (i = 0; i < n; ++i)
+ if (d[i] < 0.)
+ d[i] = -d[i];
+ free(w);
+ return 0;
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/svdu1v.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/svdu1v.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/svdu1v.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,93 @@
+/* svdu1v.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include <stdlib.h>
+#include "ccmath.h"
+int svdu1v(double *d, double *a, int m, double *v, int n)
+{
+ double *p, *p1, *q, *pp, *w, *e;
+
+ double s, h, r, t, sv;
+
+ int i, j, k, mm, nm, ms;
+
+ if (m < n)
+ return -1;
+ w = (double *)calloc(m + n, sizeof(double));
+ e = w + m;
+ for (i = 0, mm = m, nm = n - 1, p = a; i < n; ++i, --mm, --nm, p += n + 1) {
+ if (mm > 1) {
+ sv = h = 0.;
+ for (j = 0, q = p, s = 0.; j < mm; ++j, q += n) {
+ w[j] = *q;
+ s += *q * *q;
+ }
+ if (s > 0.) {
+ h = sqrt(s);
+ if (*p < 0.)
+ h = -h;
+ s += *p * h;
+ s = 1. / s;
+ t = 1. / (w[0] += h);
+ sv = 1. + fabs(*p / h);
+ for (k = 1, ms = n - i; k < ms; ++k) {
+ for (j = 0, q = p + k, r = 0.; j < mm; q += n)
+ r += w[j++] * *q;
+ r *= s;
+ for (j = 0, q = p + k; j < mm; q += n)
+ *q -= r * w[j++];
+ }
+ for (j = 1, q = p; j < mm;)
+ *(q += n) = t * w[j++];
+ }
+ *p = sv;
+ d[i] = -h;
+ }
+ if (mm == 1)
+ d[i] = *p;
+ p1 = p + 1;
+ sv = h = 0.;
+ if (nm > 1) {
+ for (j = 0, q = p1, s = 0.; j < nm; ++j, ++q)
+ s += *q * *q;
+ if (s > 0.) {
+ h = sqrt(s);
+ if (*p1 < 0.)
+ h = -h;
+ sv = 1. + fabs(*p1 / h);
+ s += *p1 * h;
+ s = 1. / s;
+ t = 1. / (*p1 += h);
+ for (k = n, ms = n * (m - i); k < ms; k += n) {
+ for (j = 0, q = p1, pp = p1 + k, r = 0.; j < nm; ++j)
+ r += *q++ * *pp++;
+ r *= s;
+ for (j = 0, q = p1, pp = p1 + k; j < nm; ++j)
+ *pp++ -= r * *q++;
+ }
+ for (j = 1, q = p1 + 1; j < nm; ++j)
+ *q++ *= t;
+ }
+ *p1 = sv;
+ e[i] = -h;
+ }
+ if (nm == 1)
+ e[i] = *p1;
+ }
+ ldvmat(a, v, n);
+ atou1(a, m, n);
+ qrbdu1(d, e, a, m, v, n);
+ for (i = 0; i < n; ++i) {
+ if (d[i] < 0.) {
+ d[i] = -d[i];
+ for (j = 0, p = v + i; j < n; ++j, p += n)
+ *p = -*p;
+ }
+ }
+ free(w);
+ return 0;
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/svduv.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/svduv.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/svduv.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,93 @@
+/* svduv.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include <stdlib.h>
+#include "ccmath.h"
+int svduv(double *d, double *a, double *u, int m, double *v, int n)
+{
+ double *p, *p1, *q, *pp, *w, *e;
+
+ double s, h, r, t, sv;
+
+ int i, j, k, mm, nm, ms;
+
+ if (m < n)
+ return -1;
+ w = (double *)calloc(m + n, sizeof(double));
+ e = w + m;
+ for (i = 0, mm = m, nm = n - 1, p = a; i < n; ++i, --mm, --nm, p += n + 1) {
+ if (mm > 1) {
+ sv = h = 0.;
+ for (j = 0, q = p, s = 0.; j < mm; ++j, q += n) {
+ w[j] = *q;
+ s += *q * *q;
+ }
+ if (s > 0.) {
+ h = sqrt(s);
+ if (*p < 0.)
+ h = -h;
+ s += *p * h;
+ s = 1. / s;
+ t = 1. / (w[0] += h);
+ sv = 1. + fabs(*p / h);
+ for (k = 1, ms = n - i; k < ms; ++k) {
+ for (j = 0, q = p + k, r = 0.; j < mm; q += n)
+ r += w[j++] * *q;
+ r *= s;
+ for (j = 0, q = p + k; j < mm; q += n)
+ *q -= r * w[j++];
+ }
+ for (j = 1, q = p; j < mm;)
+ *(q += n) = t * w[j++];
+ }
+ *p = sv;
+ d[i] = -h;
+ }
+ if (mm == 1)
+ d[i] = *p;
+ p1 = p + 1;
+ sv = h = 0.;
+ if (nm > 1) {
+ for (j = 0, q = p1, s = 0.; j < nm; ++j, ++q)
+ s += *q * *q;
+ if (s > 0.) {
+ h = sqrt(s);
+ if (*p1 < 0.)
+ h = -h;
+ sv = 1. + fabs(*p1 / h);
+ s += *p1 * h;
+ s = 1. / s;
+ t = 1. / (*p1 += h);
+ for (k = n, ms = n * (m - i); k < ms; k += n) {
+ for (j = 0, q = p1, pp = p1 + k, r = 0.; j < nm; ++j)
+ r += *q++ * *pp++;
+ r *= s;
+ for (j = 0, q = p1, pp = p1 + k; j < nm; ++j)
+ *pp++ -= r * *q++;
+ }
+ for (j = 1, q = p1 + 1; j < nm; ++j)
+ *q++ *= t;
+ }
+ *p1 = sv;
+ e[i] = -h;
+ }
+ if (nm == 1)
+ e[i] = *p1;
+ }
+ ldvmat(a, v, n);
+ ldumat(a, u, m, n);
+ qrbdv(d, e, u, m, v, n);
+ for (i = 0; i < n; ++i) {
+ if (d[i] < 0.) {
+ d[i] = -d[i];
+ for (j = 0, p = v + i; j < n; ++j, p += n)
+ *p = -*p;
+ }
+ }
+ free(w);
+ return 0;
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/svdval.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/svdval.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/svdval.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,80 @@
+/* svdval.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include <stdlib.h>
+#include "ccmath.h"
+int svdval(double *d, double *a, int m, int n)
+{
+ double *p, *p1, *q, *w, *v;
+
+ double s, h, u;
+
+ int i, j, k, mm, nm, ms;
+
+ if (m < n)
+ return -1;
+ w = (double *)calloc(m, sizeof(double));
+ for (i = 0, mm = m, nm = n - 1, p = a; i < n; ++i, --mm, --nm, p += n + 1) {
+ if (mm > 1) {
+ for (j = 0, q = p, s = 0.; j < mm; ++j, q += n) {
+ w[j] = *q;
+ s += *q * *q;
+ }
+ if (s > 0.) {
+ h = sqrt(s);
+ if (*p < 0.)
+ h = -h;
+ s += *p * h;
+ s = 1. / s;
+ w[0] += h;
+ for (k = 1, ms = n - i; k < ms; ++k) {
+ for (j = 0, q = p + k, u = 0.; j < mm; q += n)
+ u += w[j++] * *q;
+ u *= s;
+ for (j = 0, q = p + k; j < mm; q += n)
+ *q -= u * w[j++];
+ }
+ *p = -h;
+ }
+ }
+ p1 = p + 1;
+ if (nm > 1) {
+ for (j = 0, q = p1, s = 0.; j < nm; ++j, ++q)
+ s += *q * *q;
+ if (s > 0.) {
+ h = sqrt(s);
+ if (*p1 < 0.)
+ h = -h;
+ s += *p1 * h;
+ s = 1. / s;
+ *p1 += h;
+ for (k = n, ms = n * (m - i); k < ms; k += n) {
+ for (j = 0, q = p1, v = p1 + k, u = 0.; j < nm; ++j)
+ u += *q++ * *v++;
+ u *= s;
+ for (j = 0, q = p1, v = p1 + k; j < nm; ++j)
+ *v++ -= u * *q++;
+ }
+ *p1 = -h;
+ }
+ }
+ }
+
+ for (j = 0, p = a; j < n; ++j, p += n + 1) {
+ d[j] = *p;
+ if (j != n - 1)
+ w[j] = *(p + 1);
+ else
+ w[j] = 0.;
+ }
+ qrbdi(d, w, n);
+ for (i = 0; i < n; ++i)
+ if (d[i] < 0.)
+ d[i] = -d[i];
+ free(w);
+ return 0;
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/trncm.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/trncm.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/trncm.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,23 @@
+/* trncm.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include "ccmath.h"
+void trncm(Cpx * a, int n)
+{
+ Cpx s, *p, *q;
+
+ int i, j, e;
+
+ for (i = 0, e = n - 1; i < n - 1; ++i, --e, a += n + 1) {
+ for (j = 0, p = a + 1, q = a + n; j < e; ++j) {
+ s = *p;
+ *p++ = *q;
+ *q = s;
+ q += n;
+ }
+ }
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/trnm.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/trnm.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/trnm.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,22 @@
+/* trnm.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+void trnm(double *a, int n)
+{
+ double s, *p, *q;
+
+ int i, j, e;
+
+ for (i = 0, e = n - 1; i < n - 1; ++i, --e, a += n + 1) {
+ for (p = a + 1, q = a + n, j = 0; j < e; ++j) {
+ s = *p;
+ *p++ = *q;
+ *q = s;
+ q += n;
+ }
+ }
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/unfl.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/unfl.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/unfl.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,22 @@
+/* unfl.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+static unsigned int a=69069U,c=244045795U;
+static unsigned int s,h,sbuf[256];
+double unfl()
+{ int i;
+ i=(int)(s>>24); s=sbuf[i];
+ h=a*h+c; sbuf[i]=h;
+ return s*2.328306436538696e-10;
+}
+void setunfl(unsigned int k)
+{ int j;
+ for(h=k,j=0; j<=256 ;++j){
+ h=a*h+c;
+ if(j<256) sbuf[j]=h; else s=h;
+ }
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/unitary.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/unitary.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/unitary.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,99 @@
+/* unitary.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include <stdlib.h>
+#include "ccmath.h"
+static double tpi = 6.283185307179586;
+
+static void uortho(double *g, int n);
+
+double unfl();
+
+void unitary(Cpx * u, int n)
+{
+ int i, j, k, m;
+
+ Cpx h, *v, *e, *p, *r;
+
+ double *g, *q, a;
+
+ m = n * n;
+ g = (double *)calloc(n * n, sizeof(double));
+ v = (Cpx *) calloc(m + n, sizeof(Cpx));
+ e = v + m;
+ h.re = 1.;
+ h.im = 0.;
+ for (i = 0; i < n; ++i) {
+ a = tpi * unfl();
+ e[i].re = cos(a);
+ e[i].im = sin(a);
+ a = h.re * e[i].re - h.im * e[i].im;
+ h.im = h.im * e[i].re + h.re * e[i].im;
+ h.re = a;
+ }
+ h.im = -h.im;
+ for (i = 0; i < n; ++i) {
+ a = e[i].re * h.re - e[i].im * h.im;
+ e[i].im = e[i].re * h.im + e[i].im * h.re;
+ e[i].re = a;
+ }
+ uortho(g, n);
+ for (i = 0, p = v, q = g; i < n; ++i) {
+ for (j = 0; j < n; ++j)
+ (p++)->re = *q++;
+ }
+ for (i = 0, p = v; i < n; ++i) {
+ for (j = 0, h = e[i]; j < n; ++j, ++p) {
+ a = h.re * p->re - h.im * p->im;
+ p->im = h.im * p->re + h.re * p->im;
+ p->re = a;
+ }
+ }
+ uortho(g, n);
+ for (i = m = 0, p = u; i < n; ++i, m += n) {
+ for (j = 0; j < n; ++j, ++p) {
+ p->re = p->im = 0.;
+ for (k = 0, q = g + m, r = v + j; k < n; ++k, r += n) {
+ p->re += *q * r->re;
+ p->im += *q++ * r->im;
+ }
+ }
+ }
+ free(g);
+ free(v);
+}
+
+static void uortho(double *g, int n)
+{
+ int i, j, k, m;
+
+ double *p, *q, c, s, a;
+
+ for (i = 0, p = g; i < n; ++i) {
+ for (j = 0; j < n; ++j) {
+ if (i == j)
+ *p++ = 1.;
+ else
+ *p++ = 0.;
+ }
+ }
+ for (i = 0, m = n - 1; i < m; ++i) {
+ for (j = i + 1; j < n; ++j) {
+ a = tpi * unfl();
+ c = cos(a);
+ s = sin(a);
+ p = g + n * i;
+ q = g + n * j;
+ for (k = 0; k < n; ++k) {
+ a = *p * c + *q * s;
+ *q = *q * c - *p * s;
+ *p++ = a;
+ ++q;
+ }
+ }
+ }
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/utrncm.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/utrncm.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/utrncm.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,36 @@
+/* utrncm.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include <stdlib.h>
+#include "ccmath.h"
+void utrncm(Cpx * cm, Cpx * a, Cpx * b, int n)
+{
+ Cpx z, *q0, *p, *s, *t;
+
+ int i, j, k;
+
+ q0 = (Cpx *) calloc(n, sizeof(Cpx));
+ for (i = 0; i < n; ++i, ++cm) {
+ for (j = 0, t = b; j < n; ++j) {
+ z.re = z.im = 0.;
+ for (k = 0, s = a + i * n; k < n; ++k, ++s, ++t) {
+ z.re += t->re * s->re + t->im * s->im;
+ z.im += t->im * s->re - t->re * s->im;
+ }
+ q0[j] = z;
+ }
+ for (j = 0, p = cm, t = a; j < n; ++j, p += n) {
+ z.re = z.im = 0.;
+ for (k = 0, s = q0; k < n; ++k, ++t, ++s) {
+ z.re += t->re * s->re - t->im * s->im;
+ z.im += t->im * s->re + t->re * s->im;
+ }
+ *p = z;
+ }
+ }
+ free(q0);
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/utrnhm.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/utrnhm.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/utrnhm.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,40 @@
+/* utrnhm.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+#include <stdlib.h>
+#include "ccmath.h"
+void utrnhm(Cpx * hm, Cpx * a, Cpx * b, int n)
+{
+ Cpx z, *q0, *p, *s, *t;
+
+ int i, j, k;
+
+ q0 = (Cpx *) calloc(n, sizeof(Cpx));
+ for (i = 0; i < n; ++i) {
+ for (j = 0, t = b; j < n; ++j) {
+ z.re = z.im = 0.;
+ for (k = 0, s = a + i * n; k < n; ++k, ++s, ++t) {
+ z.re += t->re * s->re + t->im * s->im;
+ z.im += t->im * s->re - t->re * s->im;
+ }
+ q0[j] = z;
+ }
+ for (j = 0, p = hm + i, t = a; j <= i; ++j, p += n) {
+ z.re = z.im = 0.;
+ for (k = 0, s = q0; k < n; ++k, ++t, ++s) {
+ z.re += t->re * s->re - t->im * s->im;
+ z.im += t->im * s->re + t->re * s->im;
+ }
+ *p = z;
+ if (j < i) {
+ z.im = -z.im;
+ hm[i * n + j] = z;
+ }
+ }
+ }
+ free(q0);
+}
Added: grass/branches/develbranch_6/lib/external/ccmath/vmul.c
===================================================================
--- grass/branches/develbranch_6/lib/external/ccmath/vmul.c (rev 0)
+++ grass/branches/develbranch_6/lib/external/ccmath/vmul.c 2009-08-27 20:50:27 UTC (rev 38888)
@@ -0,0 +1,30 @@
+/* vmul.c CCMATH mathematics library source code.
+ *
+ * Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
+ * This code may be redistributed under the terms of the GNU library
+ * public license (LGPL). ( See the lgpl.license file for details.)
+ * ------------------------------------------------------------------------
+ */
+void vmul(double *vp, double *mat, double *v, int n)
+{
+ double s, *q;
+
+ int k, i;
+
+ for (k = 0; k < n; ++k) {
+ for (i = 0, q = v, s = 0.; i < n; ++i)
+ s += *mat++ * *q++;
+ *vp++ = s;
+ }
+}
+
+double vnrm(double *u, double *v, int n)
+{
+ double s;
+
+ int i;
+
+ for (i = 0, s = 0.; i < n; ++i)
+ s += *u++ * *v++;
+ return s;
+}
More information about the grass-commit
mailing list