Algebra

The Algebra project provides structures, tools and methods of the usual algebra techniques. More...

Namespaces
namespace	STK::lapack
	namespace enclosing the wrappers of the lapack routines.

Classes
class	STK::CG< MultFunctor, ColVector, InitFunctor >
	The conjugate gradient method is an algorithm for the numerical solution of particular systems of linear equations, namely those whose matrix is symmetric and positive-definite. More...
class	STK::PCG< MultFunctor, CondFunctor, ColVector, InitFunctor >
	In most cases, preconditioning is necessary to ensure fast convergence of the conjugate gradient method. More...
class	STK::ILeastSquare< Derived >
	The class ILeastSquare is an interface class for the methods solving the least-square problem. More...
class	STK::InvertMatrix< Matrix, Size_ >
	The InvertMatrix class is a functor class allowing to compute the inverse of a symmetric matrix. More...
class	STK::IQr< Derived >
	The class IQr is an interface class for the method computing the QR Decomposition of a ArrayXX. More...
class	STK::ISvd< Derived >
	Compute the Singular Value Decomposition of an array. More...
class	STK::ISymEigen< Derived >
	The class ISymEigen is an interface class for the method computing the eigenvalue Decomposition of a symmetric ArrayXX. More...
class	STK::lapack::MultiLeastSquare< ArrayB, ArrayA >
	The class MultiLeastSQquare solve the least square problem when the response b is multidimensional. More...
class	STK::MultiLeastSquare< ArrayB, ArrayA >
	The class MultiLeastSQquare solve the least square problem when the response b is multidimensional. More...
class	STK::Qr
	The class Qr perform the QR decomposition of an ArrayXX. More...
class	STK::Svd< Array >
	The class Svd compute the Singular Value Decomposition of a Array with the Golub-Reinsch Algorithm. More...
class	STK::SymEigen< SquareArray >
	The class SymEigen compute the eigenvalue Decomposition of a symmetric ArrayXX. More...
class	STK::lapack::Qr
	{ More...
class	STK::lapack::Svd
	{ More...
class	STK::lapack::SymEigen< SquareArray >
	{ More...

Functions
template<class Lhs >
bool	STK::cholesky (ExprBase< Lhs > const &A, Array2DDiagonal< typename Lhs::Type > &D, Array2DLowerTriangular< typename Lhs::Type > &L)
	Compute the Cholesky decomposition of a symmetric matrix.
template<class Type >
Real	STK::compGivens (Type const &y, Type const &z, Type &cosinus, Type &sinus)
	Compute Givens rotation.
template<class TContainer2D >
void	STK::rightGivens (ArrayBase< TContainer2D > &M, int j1, int j2, typename TContainer2D::Type const &cosinus, typename TContainer2D::Type const &sinus)
	Apply Givens rotation.
template<class TContainer2D >
void	STK::leftGivens (ArrayBase< TContainer2D > &M, int i1, int i2, typename TContainer2D::Type const &cosinus, typename TContainer2D::Type const &sinus)
	left multiplication by a Givens ArrayXX.
template<class TContainer2D >
void	STK::gramSchmidt (ArrayBase< TContainer2D > &A)
	The gramSchmidt method perform the Gram Schmidt orthonormalization of an Array of Real.
template<class Vector >
Real	STK::house (ArrayBase< Vector > &x)
	Compute the Householder vector v of a vector x.
template<class Lhs , class Rhs >
Real	STK::dotHouse (ExprBase< Lhs > const &x, ExprBase< Rhs > const &v)
	dot product with a Householder vector.
template<class Lhs , class Rhs >
void	STK::applyLeftHouseholderVector (ArrayBase< Lhs > const &M, ExprBase< Rhs > const &v)
	left multiplication by a Householder vector.
template<class Lhs , class Rhs >
void	STK::applyLeftHouseholderArray (ArrayBase< Lhs > const &M, ArrayBase< Rhs > const &H)
	left multiplication by a Householder array.
template<class Lhs , class Rhs >
void	STK::applyRightHouseholderVector (ArrayBase< Lhs > const &M, ExprBase< Rhs > const &v)
	right multiplication by a Householder vector.
template<class TContainer2D , class Rhs >
void	STK::applyRightHouseholderArray (ArrayBase< TContainer2D > const &M, ArrayBase< Rhs > const &H)
	left multiplication by a Householder ArrayXX.
template<class Matrix >
hidden::AlgebraTraits< InvertMatrix< Matrix, hidden::Traits< Matrix >::sizeRows_ > >::Result	STK::invert (Matrix const &mat)
	Utility function allowing to compute the inverse of a matrix.
int	STK::lapack::gelsd (int m, int n, int nrhs, Real *a, int lda, Real *b, int ldb, Real *s, Real *rcond, int *rank, Real *work, int lWork, int *iwork)
	wrapper of the LAPACK GELSD routine: computes the minimum-norm solution to a real linear least squares problem.
int	STK::lapack::geqrf (int m, int n, Real *a, int lda, Real *tau, Real *work, int lWork)
	wrapper of the LAPACK DGEQRF routine: computes the Qr decomposition of a matrix.
int	STK::lapack::gesdd (char jobz, int m, int n, Real *a, int lda, Real *s, Real *u, int ldu, Real *vt, int ldvt, Real *work, int lWork, int *iWork)
	wrapper of the LAPACK DGESDD routine.
int	STK::lapack::syevr (char jobz, char range, char uplo, int n, Real *a, int lda, Real vl, Real vu, int il, int iu, Real abstol, int *m, Real *w, Real *z, int ldz, int *isuppz, Real *work, int lWork, int *iwork, int liwork)
	wrapper of the LAPACK SYEVR routine.
template<class TContainer2D >
TContainer2D &	STK::transpose (ArrayBase< TContainer2D > &A)
	The transpose method allow to transpose an array in place.
void	STK::blas::dgemm (int m, int n, int nrhs, Real *a, int lda, Real *b, int ldb, Real *s, Real *rcond, int *rank, Real *work, int lWork, int *iwork)
	wrapper of the BLAS gemm routine: performs one of the matrix-matrix operations
void	STK::blas::dgesv_ (int *n, int *nrhs, double *A, int *lda, int *ipiv, double *B, int *ldb, int *info)
	wrapper of the BLAS gesv routine:

Detailed Description

The Algebra project provides structures, tools and methods of the usual algebra techniques.

The Algebra project proposes some set of template function for computing dot product, weighted dot product, vector norm, weighted vector norm and so on... It implement also linear algebra methods for the ArrayXX and the ArraySquareX classes:

• The Qr decomposition of an arbitrary matrix of Real,

  See also

  Qr

• The svd decomposition of an arbitrary matrix of real,

  See also

  Svd

• An Eigenvalue decomposition for symmetric (square) matrices and a generalized inverse method for such matrices,

  See also

  SymEigen.

• methods for solving the linear least square problem.

It proposes also some set of method for performing

• Givens rotations on a matrix
• GramSchmidt orthogonalization of the columns of a matrix/array
• Householder rotations of an array/matrix.
• Cholesky decomposition.

Function Documentation

applyLeftHouseholderArray()

template<class Lhs , class Rhs >

void STK::applyLeftHouseholderArray	(	ArrayBase< Lhs > const &	M,
		ArrayBase< Rhs > const &	H
	)

left multiplication by a Householder array.

Perform a left multiplication of the Array M by a Householder Matrix H. M <- HM with H = I + WZ'. The matrix H is represented as a product of elementary reflectors H = H(1) H(2) . . . H(k)

Parameters

M	the matrix to multiply
H	the matrix with the Householder vectors

Definition at line 142 of file STK_Householder.h.

{
  typedef typename hidden::Traits<Rhs>::Col ColVector;
  // compute the number of iterations
  int first = H.beginCols(), last = std::min( H.lastIdxCols(), H.lastIdxRows());
  // get range of the first Householder vector
  Range range_ve(last, H.lastIdxRows(), 0);
  // iterations
  for (int j=last; j>=first; j--)
  {
    // apply left Householder vector to M
    ColVector v(H.asDerived(), range_ve, j);
    applyLeftHouseholderVector(M, v);
    // decrease range of the Householder vector
    range_ve.decFirst(1);
  }
}

References STK::applyLeftHouseholderVector().

Referenced by main().

applyLeftHouseholderVector()

template<class Lhs , class Rhs >

void STK::applyLeftHouseholderVector	(	ArrayBase< Lhs > const &	M,
		ExprBase< Rhs > const &	v
	)

left multiplication by a Householder vector.

Perform a left multiplication of the matrix M with a Householder matrix . Overwrite M with HM.

Parameters

M	the matrix to multiply (input/output)
v	the Householder vector (input)

Definition at line 109 of file STK_Householder.h.

{
  typedef typename hidden::Traits<Lhs>::Col ColVector;
  // get beta
  Real beta = v.front();
  if (beta)
  {
    // Multiplication of the cols by P=I+beta vv'
    for (int j=M.beginCols(); j<M.endCols(); j++)
    {
      // a ref on the jth column of M
      ColVector Mj(M.asDerived(), v.range(), j);
      // Computation of aux=beta* <v,M^j>
      Real aux =  dotHouse( Mj, v) * beta;
      // updating columns
      Mj.front() += aux;
      for(int i = v.begin()+1; i < v.end(); ++i) Mj[i] += v[i] * aux;
    }
  }
}

References STK::dotHouse().

Referenced by STK::applyLeftHouseholderArray(), STK::bidiag(), STK::IQr< Derived >::compQ(), and STK::Qr::qr().

applyRightHouseholderArray()

template<class TContainer2D , class Rhs >

void STK::applyRightHouseholderArray	(	ArrayBase< TContainer2D > const &	M,
		ArrayBase< Rhs > const &	H
	)

left multiplication by a Householder ArrayXX.

Perform a right multiplication of the matrix M with a Householder Marix H. M <- MP with H = I + WZ'. The Householder vectors are stored in the rows of H.

Parameters

M	the Array to multiply
H	the Householder ArrayXX

Definition at line 202 of file STK_Householder.h.

{
  // compute the number of iterations
  int first = H.beginCols(), last = std::min( H.lastIdxCols(), H.lastIdxRows());
  // get range of the first Householder vector
  Range range_ve(last, H.lastIdxRows());
  // iterations
  for (int j=last; j>=first; j--)
  {
    // apply left Householder vector to M
    typename Rhs::Col v(H.asDerived(), range_ve, j);
    applyRightHouseholderVector(M, v);
    // decrease range of the Householder vector
    range_ve.decFirst(1);
  }
}

References STK::applyRightHouseholderVector().

applyRightHouseholderVector()

template<class Lhs , class Rhs >

void STK::applyRightHouseholderVector	(	ArrayBase< Lhs > const &	M,
		ExprBase< Rhs > const &	v
	)

right multiplication by a Householder vector.

Perform a right multiplication of the matrix M with a Householder matrix . Overwrite M with MH.

Parameters

M	the matrix to multiply (input/output)
v	the Householder vector (input)

Definition at line 170 of file STK_Householder.h.

{
  // get beta
  Real beta = v.front();
  if (beta)
  {
    // Multiplication of the cols by P=I+beta vv'
    for (int i=M.beginRows(); i<M.endRows(); i++)
    {
      // a ref on the ith row of M
      typename hidden::Traits<Lhs>::Row Mi(M.asDerived(), v.range(), i);
      // Computation of aux=beta* <v,M_i>
      Real aux =  dotHouse( Mi, v) * beta;
      // updating column
      Mi.front() += aux;
      // Computation of M_i + beta <v,M_i>  v = M_i + aux v'
      for (int i=v.begin()+1; i<v.end(); i++) Mi[i] +=  v[i] * aux;
    }
  }
}

References STK::dotHouse().

Referenced by STK::applyRightHouseholderArray(), and STK::bidiag().

cholesky()

template<class Lhs >

bool STK::cholesky	(	ExprBase< Lhs > const &	A,
		Array2DDiagonal< typename Lhs::Type > &	D,
		Array2DLowerTriangular< typename Lhs::Type > &	L
	)

Compute the Cholesky decomposition of a symmetric matrix.

Compute a closely related variant of the classical Cholesky decomposition, the LDL decomposition where L is a lower unit triangular matrix and D is a diagonal matrix.

Parameters

A	the square and symmetric matrix to decompose, only the lower part is used
D	diagonal matrix of the decomposition
L	lower triangular matrix of the decomposition

Returns

true if no error, false otherwise (results are misleading)

Definition at line 63 of file STK_Cholesky.h.

{
  typedef typename Lhs::Type Type;
  if(A.rows() != A.cols())
  { STKRUNTIME_ERROR_NO_ARG(cholesky,A must be square symmetric matrix);}
  // start
  D.resize(A.rows(), A.cols());
  L.resize(A.rows(), A.cols());
  bool nozero = true;
  for (int j=A.beginCols(); j<A.endCols(); ++j )
  {
    Type sum1 = A(j,j);
    for (int k=A.beginCols(); k<j; ++k) { sum1 -= L(j,k) * L(j,k) * D[k];}
    D[j] = sum1;
    nozero &= (sum1!=0);
    L(j,j) = Type(1);
    for (int i=j+1; i<A.endRows(); ++i)
    {
      Type sum2 = A(i,j);
      for (int k=A.beginCols(); k<j; ++k) { sum2 -= L(i,k)*L(j,k)* D[k];}
      L(i,j) = sum1 ? sum2/sum1 : sum2;
    }
  }
  return nozero;
}

References STK::cholesky(), and STKRUNTIME_ERROR_NO_ARG.

Referenced by STK::cholesky(), and main().

compGivens()

template<class Type >

Real STK::compGivens	(	Type const &	y,
		Type const &	z,
		Type &	cosinus,
		Type &	sinus
	)

Compute Givens rotation.

Compute the Givens rotation

in order to eliminate the coefficient z and return the value r of the rotated element.

See also

http://en.wikipedia.org/wiki/Givens_rotation

Parameters

y	The coefficient to rotate (input)
z	the coefficient to eliminate (input)
cosinus	the cosinus of the Givens rotation (output)
sinus	the sinus of the Givens rotation rotation (output)

Definition at line 75 of file STK_Givens.h.

{
  // trivial case
  if (z == 0)
  {
    sinus   = 0.0;
    cosinus = 1.0;
    return y;
  }
  // compute Givens rotation avoiding underflow and overflow
  if (std::abs(z) > std::abs(y))
  {
    Type aux = y/z;
    Type t   = sign(z, sqrt(1.0+aux*aux));
    sinus    = 1.0/t;
    cosinus  = sinus * aux;
    return t*z;
  }
  else
  {
    Type aux = z/y;
    Type t   = sign(y, sqrt(1.0+aux*aux));
    cosinus  = 1.0/t;
    sinus    = cosinus * aux;
    return t*y;
  }
}

References STK::sign().

Referenced by STK::IQr< Derived >::eraseCol(), STK::IQr< Derived >::insertCol(), and STK::IQr< Derived >::pushBackCol().

dgemm()

void STK::blas::dgemm ( int  m, int  n, int  nrhs, Real *  a, int  lda, Real *  b, int  ldb, Real *  s, Real *  rcond, int *  rank, Real *  work, int  lWork, int *  iwork  )

inline

wrapper of the BLAS gemm routine: performs one of the matrix-matrix operations

    C := alpha*op( A )*op( B ) + beta*C,

 where  op( X ) is one of

    op( X ) = X   or   op( X ) = X**T,
@verbatim
alpha and beta are scalars, and A, B and C are matrices, with op( A )
an m by k matrix,  op( B )  a  k by n matrix and  C an m by n matrix.

Parameters

[in]	TRANSA	CHARACTER*1 On entry, TRANSA specifies the form of op( A ) to be used in the matrix multiplication as follows: TRANSA = 'N' or 'n', op( A ) = A. TRANSA = 'T' or 't', op( A ) = A**T. TRANSA = 'C' or 'c', op( A ) = A**T.
[in]	TRANSB	Char[0] On entry, TRANSB specifies the form of op( B ) to be used in the matrix multiplication as follows: TRANSB = 'N' or 'n', op( B ) = B. TRANSB = 'T' or 't', op( B ) = B**T. TRANSB = 'C' or 'c', op( B ) = B**T.
[in]	M	int On entry, M specifies the number of rows of the matrix op( A ) and of the matrix C. M must be at least zero.
[in]	N	int On entry, N specifies the number of columns of the matrix op( B ) and the number of columns of the matrix C. N must be at least zero.
[in]	K	int On entry, K specifies the number of columns of the matrix op( A ) and the number of rows of the matrix op( B ). K must be at least zero.
[in]	ALPHA	Real. On entry, ALPHA specifies the scalar alpha.
[in]	A	a Real array of dimension (LDA, ka), where ka is k when TRANSA = 'N' or 'n', and is m otherwise. Before entry with TRANSA = 'N' or 'n', the leading m by k part of the array A must contain the matrix A, otherwise the leading k by m part of the array A must contain the matrix A.
[in]	LDA	int On entry, LDA specifies the first dimension of A as declared in the calling (sub) program. When TRANSA = 'N' or 'n' then LDA must be at least max( 1, m ), otherwise LDA must be at least max( 1, k ).
[in]	B	Real array of dimension (LDB, kb), where kb is n when TRANSB = 'N' or 'n', and is k otherwise. Before entry with TRANSB = 'N' or 'n', the leading k by n part of the array B must contain the matrix B, otherwise the leading n by k part of the array B must contain the matrix B.
[in]	LDB	INTEGER On entry, LDB specifies the first dimension of B as declared in the calling (sub) program. When TRANSB = 'N' or 'n' then LDB must be at least max( 1, k ), otherwise LDB must be at least max( 1, n ).
[in]	BETA	DOUBLE PRECISION. On entry, BETA specifies the scalar beta. When BETA is supplied as zero then C need not be set on input.
[in,out]	C	Real array of dimension (LDC, N) Before entry, the leading m by n part of the array C must contain the matrix C, except when beta is zero, in which case C need not be set on entry. On exit, the array C is overwritten by the m by n matrix ( alpha*op( A )*op( B ) + beta*C ).
[in]	LDC	int On entry, LDC specifies the first dimension of C as declared in the calling (sub) program. LDC must be at least max( 1, m ).

Definition at line 148 of file STK_blas_util.h.

{
  int info = 1;
#ifdef STKUSEBLAS
#ifdef STKREALAREFLOAT
  sgemm_( &m, &n, &nrhs, a, &lda, b, &ldb, s, rcond, rank, work, &lWork, iwork, &info);
#else
  dgemm_( &m, &n, &nrhs
          , a, &lda
          , b, &ldb
          , s
          , rcond, rank
          , work, &lWork, iwork, &info);
#endif
#endif
  return info;
}

dgesv_()

void STK::blas::dgesv_	(	int *	n,
		int *	nrhs,
		double *	A,
		int *	lda,
		int *	ipiv,
		double *	B,
		int *	ldb,
		int *	info
	)

wrapper of the BLAS gesv routine:

dotHouse()

template<class Lhs , class Rhs >

Real STK::dotHouse	(	ExprBase< Lhs > const &	x,
		ExprBase< Rhs > const &	v
	)

dot product with a Householder vector.

Scalar product of an 1D container with a Householder vector d = < x,v>. The first componant of a Householder vector is 1.0

Parameters

x	first vector
v	the Householder vector

Definition at line 90 of file STK_Householder.h.

{
  // compute the product
  Real sum = x[ v.begin()] /* *1.0 */;
  for (int i= v.begin()+1; i<v.end(); i++) sum += x[i] * v[i];
  // return <x,v>
  return(sum);
}

References STK::sum().

Referenced by STK::applyLeftHouseholderVector(), and STK::applyRightHouseholderVector().

gelsd()

int STK::lapack::gelsd ( int  m, int  n, int  nrhs, Real *  a, int  lda, Real *  b, int  ldb, Real *  s, Real *  rcond, int *  rank, Real *  work, int  lWork, int *  iwork  )

inline

wrapper of the LAPACK GELSD routine: computes the minimum-norm solution to a real linear least squares problem.

GELSD computes the minimum-norm solution to a real linear least squares problem: minimize 2-norm(| b - A*x |) using the singular value decomposition (SVD) of A. A is an M-by-N matrix which may be rank-deficient.

Parameters

[in]	m	The number of rows of A. m >= 0 .
[in]	n	The number of columns of A. n>= 0 .
[in]	nrhs	The number of right hand sides, i.e., the number of columns of the matrices B and X. nrhs >= 0.
[in]	a	On entry, the M-by-N matrix A. On exit, A has been destroyed.
[in]	lda	The leading dimension of the array A. lda >= max(1,m).
[in,out]	b	On entry, the M-by-NRHS right hand side matrix B. * On exit, B is overwritten by the N-by-NRHS solution matrix X. * If m >= n and RANK = n, the residual sum-of-squares for the solution in * the i-th column is given by the sum of squares of elements n+1:m in that * column. *
[in]	ldb	The leading dimension of the array B. ldb >= max(1,max(m,n)) .
[out]	s	The singular values of A in decreasing order. * The condition number of A in the 2-norm = s(1)/s(min(m,n)). *
[out]	rcond	used to determine the effective rank of A. * Singular values s(i) <= RCOND*s(1) are treated as zero. * If RCOND < 0, machine precision is used instead. *
[out]	rank	The effective rank of A, i.e., the number of singular values which are greater than RCOND*S(1).
[out]	work	On exit, if info = 0, work(1) returns the optimal lWork.
[in]	lWork	The dimension of the array work. lWork must be at least 1. * The exact minimum amount of workspace needed depends on M, N and NRHS. * As long as lWork is at least 12*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2, * if M is greater than or equal to N or 12*M + 2*M*SMLSIZ + 8*M*NLVL + M*NRHS + (SMLSIZ+1)**2, * if M is less than N, the code will execute correctly. * SMLSIZ is returned by ILAENV and is equal to the maximum size of * the subproblems at the bottom of the computation tree (usually about 25), * and NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 ) * For good performance, lWork should generally be larger. * If lWork = -1, then a workspace query is assumed; the routine only * calculates the optimal size of the work array, returns this value as the * first entry of the work array, and no error message related to lWork is * issued by XERBLA. *
	iwork	array of integer of dimension (MAX(1,LiWork)) * LiWork >= 3 * MINMN * NLVL + 11 * MINMN, where MINMN = MIN( M,N ). * On exit, if info = 0, iWork(1) returns the minimum LiWork *

Returns

info

*    = 0:  successful exit
*    < 0:  if info = -i, the i-th argument had an illegal value.
*    > 0:  the algorithm for computing the SVD failed to converge;
*    if info = i, i off-diagonal elements of an intermediate bidiagonal
*    form did not converge to zero.
*  

*    Further Details
*    ===============
*      Based on contributions by
*      Ming Gu and Ren-Cang Li, Computer Science Division, University of
*      California at Berkeley, USA
*      Osni Marques, LBNL/NERSC, USA
* 

Definition at line 171 of file STK_lapack_Util.h.

{
  int info = 1;
#ifdef STKUSELAPACK
#ifdef STKREALAREFLOAT
  sgelsd_( &m, &n, &nrhs, a, &lda, b, &ldb, s, rcond, rank, work, &lWork, iwork, &info);
#else
  dgelsd_( &m, &n, &nrhs, a, &lda, b, &ldb, s, rcond, rank, work, &lWork, iwork, &info);
#endif
#endif
  return info;
}

Referenced by STK::lapack::MultiLeastSquare< ArrayB, ArrayA >::computeLS().

geqrf()

int STK::lapack::geqrf ( int  m, int  n, Real *  a, int  lda, Real *  tau, Real *  work, int  lWork  )

inline

wrapper of the LAPACK DGEQRF routine: computes the Qr decomposition of a matrix.

Parameters

[in]	m	The number of rows of the matrix A. M >= 0.
[in]	n	The number of columns of the matrix A. N >= 0.
[in,out]	a	Real array, dimension (lda, N) * On entry, the M-by-N matrix A. * On exit, the elements on and above the diagonal of the array * contain the min(M,N)-by-N upper trapezoidal matrix R (R is * upper triangular if m >= n); the elements below the diagonal, * with the array TAU, represent the orthogonal matrix Q as a * product of min(m,n) elementary reflectors (see Further Details). *
[in]	lda	The leading dimension of the array A. lda >= max(1,M).
[out]	tau	Real array, dimension min(M,N) The scalar factors of the elementary reflectors (see Further Details).
[in,out]	work	Real array, dimension (MAX(1,lWork)) * On exit, if info = 0, work(1) returns the optimal lWork. *
[in]	lWork	The dimension of the array work * lWork >= max(1,N). * For optimum performance lWork >= N*NB, where NB is the optimal blocksize. * * If lWork = -1, then a workspace query is assumed; the routine * only calculates the optimal size of the work array, returns * this value as the first entry of the work array, and no error * message related to lWork is issued by XERBLA. *

Returns

info

*  = 0:  successful exit
*  < 0:  if info = -i, the i-th argument had an illegal value
* 

*  Further Details
*  ===============
*
*  The matrix Q is represented as a product of elementary reflectors
*
*     Q = H(1) H(2) . . . H(k), where k = min(m,n).
*
*  Each H(i) has the form
*
*     H(i) = I - tau * v * v'
*
*  where tau is a real scalar, and v is a real vector with
*  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i),
*  and tau in TAU(i).
* 

Definition at line 250 of file STK_lapack_Util.h.

{
  int info = 1;
 
#ifdef STKUSELAPACK
#ifdef STKREALAREFLOAT
  sgeqrf_(&m, &n, a, &lda, tau, work, &lWork, &info);
#else
  dgeqrf_(&m, &n, a, &lda, tau, work, &lWork, &info);
#endif
#endif
 
  return info;
}

Referenced by STK::lapack::Qr::computeQr().

gesdd()

int STK::lapack::gesdd ( char  jobz, int  m, int  n, Real *  a, int  lda, Real *  s, Real *  u, int  ldu, Real *  vt, int  ldvt, Real *  work, int  lWork, int *  iWork  )

inline

wrapper of the LAPACK DGESDD routine.

DGESDD computes the Svd decomposition of a matrix.

*   DGESDD computes the singular value decomposition (SVD) of a real
*   m-by-n matrix a, optionally computing the left and right singular
*   vectors.  If singular vectors are desired, it uses a
*   divide-and-conquer algorithm.
*
*   The SVD is written
*
*      a = u * SIGMA * transpose(v)
*
*   where SIGMA is an m-by-n matrix which is zero except for its
*   min(m,n) diagonal elements, u is an m-by-m orthogonal matrix, and
*   v is an n-by-n orthogonal matrix. The diagonal elements of SIGMA
*   are the singular values of a; they are real and non-negative, and
*   are returned in descending order. The first min(m,n) columns of
*   u and v are the left and right singular vectors of a.
*  

Note

routine returns vt = v**T, not v.

Arguments:

Parameters

[in]	jobz	* jobz is Char*1 * Specifies options for computing all or part of the matrix u: * = 'A': all m columns of u and all n rows of v**T are * returned in the arrays u and vt; * = 'S': the first min(m,n) columns of u and the first * min(m,n) rows of v**T are returned in the arrays U * and vt; * = 'O': If m >= n, the first n columns of u are overwritten * on the array a and all rows of v**T are returned in * the array vt; * otherwise, all columns of u are returned in the * array u and the first m rows of v**T are overwritten * in the array a; * = 'N': no columns of u or rows of v**T are computed. *
[in]	m	* m is Integer * The number of rows of the input matrix a. m >= 0. *
[in]	n	* n is Integer * The number of columns of the input matrix a. n >= 0. *
[in,out]	a	* a is STK::Real array, dimension (lda,n) * On entry, the m-by-n matrix a. * On exit, * if jobz = 'O', a is overwritten with the first n columns * of u (the left singular vectors, stored * columnwise) if m >= n; * a is overwritten with the first m rows * of v**T (the right singular vectors, stored * rowwise) otherwise. * if jobz .ne. 'O', the contents of a are destroyed. *
[in]	lda	* lda is Integer * The leading dimension of the array a. lda >= max(1,m). *
[out]	s	* s is STK::Real array, dimension (min(m,n)) * The singular values of a, sorted so that s[i] >= s[i+1]. *
[out]	u	* u is STK::Real array, dimension (ldu,ucol) * ucol = m if jobz = 'A' or jobz = 'O' and m < n; * ucol = min(m,n) if jobz = 'S'. * If jobz = 'A' or jobz = 'O' and m < n, u contains the m-by-m * orthogonal matrix u; * if jobz = 'S', u contains the first min(m,n) columns of u * (the left singular vectors, stored columnwise); * if jobz = 'O' and m >= n, or jobz = 'N', u is not referenced. *
[in]	ldu	* ldu is Integer * The leading dimension of the array U. ldu >= 1; if * jobz = 'S' or 'A' or jobz = 'O' and m < n, ldu >= m. *
[out]	vt	* vt is STK::Real array, dimension (ldvt,n) * If jobz = 'A' or jobz = 'O' and m >= n, vt contains the * N-by-N orthogonal matrix v**T; * if jobz = 'S', vt contains the first min(m,n) rows of * v**T (the right singular vectors, stored rowwise); * if jobz = 'O' and m < n, or jobz = 'N', vt is not referenced. *
[in]	ldvt	* ldvt is Integer * The leading dimension of the array vt. ldvt >= 1; if * jobz = 'A' or jobz = 'O' and m >= n, ldvt >= n; * if jobz = 'S', ldvt >= min(m,n). *
[out]	work	* work is STK::Real array, dimension (MAX(1,lWork)) * On exit, if info = 0, work(1) returns the optimal lWork; *
[in]	lWork	* lWork is Integer * The dimension of the array work. lWork >= 1. * If jobz = 'N', * lWork >= 3*min(m,n) + max(max(m,n),7*min(m,n)). * If jobz = 'O', * lWork >= 3*min(m,n) + * max(max(m,n),5*min(m,n)*min(m,n)+4*min(m,n)). * If jobz = 'S' or 'A' * lWork >= min(m,n)*(6+4*min(m,n))+max(m,n) * For good performance, lWork should generally be larger. * If lWork = -1 but other input arguments are legal, work(1) * returns the optimal lWork. *
[out]	iWork	* iWork is Integer array, dimension (8*min(m,n)) *

Returns

info

*          info is Integer
*          = 0:  successful exit.
*          < 0:  if info = -i, the i-th argument had an illegal value.
*          > 0:  DBDSDC did not converge, updating process failed.
* 

Authors:

Author

Univ. of Tennessee

Univ. of California Berkeley

Univ. of Colorado Denver

NAG Ltd.

Contributors:

Ming Gu and Huan Ren, Computer Science Division, University of
California at Berkeley, USA

Definition at line 436 of file STK_lapack_Util.h.

{
  int info = 1;
#ifdef STKUSELAPACK
#ifdef STKREALAREFLOAT
  sgesdd_( &jobz, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work, &lWork, iWork, &info);
#else
  dgesdd_( &jobz, &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work, &lWork, iWork, &info);
#endif
#endif
  return info;
}

Referenced by STK::lapack::Svd::computeSvd().

gramSchmidt()

template<class TContainer2D >

void STK::gramSchmidt

(

ArrayBase< TContainer2D > &

A

)

The gramSchmidt method perform the Gram Schmidt orthonormalization of an Array of Real.

Parameters

A	the container to orthnormalize

Definition at line 52 of file STK_GramSchmidt.h.

{
  // orthonormalize
  for (int j= A.beginCols(); j< A.endCols(); j++)
  {
    for( int i= A.beginCols(); i < j; i++)
    {
      const Real dotij = dot(A.asDerived().col(i), A.asDerived().col(j));
      A.asDerived().col(j) -= A.asDerived().col(i) * dotij;
    }
    // normalize
    const Real norm = A.asDerived().col(j).norm();
    if (norm)
    {
      A.asDerived().col(j)/=norm;
    }
  }
}

References STK::dot(), and STK::norm().

Referenced by STK::LinearAAModel< Array >::simul().

house()

template<class Vector >

Real STK::house

(

ArrayBase< Vector > &

x

)

Compute the Householder vector v of a vector x.

Given a vector x, compute the vector v of the matrix of Householder such that . The vector v is of the form : and is stored in x. The value 1 is skipped and is stored in front of v. The method return the value v1.

Parameters

x	the vector to rotate, it is overwritten by v

Definition at line 56 of file STK_Householder.h.

{
  // compute L^{\infty} norm of X, declare result and norm2 of X
  Real scale  = x.normInf(), v1, norm2 = 0.0;
  if (scale)  // if not 0.0
  {
    norm2 = (x.asDerived().sub(_R(x.begin()+1,x.lastIdx()))/=scale).norm2();
  }
  // check if the lower part is significative
  if (norm2 < Arithmetic<Real>::epsilon())
  { v1 = x.front(); x.front() = 0.0; }
  else
  {
    Real s, aux1 = x.front()/scale;
    // compute v1 = P_v X and beta of the Householder vector
    v1 =  (norm2 = sign(aux1, sqrt(aux1*aux1+norm2))) * scale;
    // compute and save beta
    x.front() = (s = aux1-norm2)/norm2;
    // comp v and save it
    x.asDerived().sub(_R(x.begin()+1,x.lastIdx())) /= s;
  }
  return v1;
}

References _R, and STK::sign().

Referenced by STK::bidiag(), STK::Qr::qr(), and STK::SymEigen< SquareArray >::tridiagonalize().

invert()

template<class Matrix >

hidden::AlgebraTraits< InvertMatrix< Matrix, hidden::Traits< Matrix >::sizeRows_ > >::Result STK::invert

(

Matrix const &

mat

)

Utility function allowing to compute the inverse of a matrix.

Definition at line 197 of file STK_InvertMatrix.h.

{
  enum
  {
    sizeRows_ = hidden::Traits<Matrix>::sizeRows_,
    sizeCols_ = hidden::Traits<Matrix>::sizeCols_,
    size_ = (sizeRows_<sizeCols_) ? sizeRows_ : sizeCols_
  };
  return InvertMatrix<Matrix, sizeRows_>(mat)();
}

Referenced by main(), STK::BSplineRegression< YArray, XVector, Weights >::regressionStep(), STK::MultidimRegression< Array, Weight >::regressionStep(), STK::MultidimRegression< Array, Weight >::regressionStep(), and STK::BSplineRegression< YArray, XVector, Weights >::regressionStep().

leftGivens()

template<class TContainer2D >

void STK::leftGivens	(	ArrayBase< TContainer2D > &	M,
		int	i1,
		int	i2,
		typename TContainer2D::Type const &	cosinus,
		typename TContainer2D::Type const &	sinus
	)

left multiplication by a Givens ArrayXX.

Perform a left multiplication of the matrix M with a Givens ArrayXX on the row1 and row2. row1 should be less than row2. The Array M is passed as const as we are using reference on the two rows we want to rotate.

See also

http://en.wikipedia.org/wiki/Givens_rotation

Parameters

M	the matix to multiply
i1,i2	the first and second rows
cosinus,sinus	the cosinus and sinus of the givens rotation

Definition at line 148 of file STK_Givens.h.

{
  typedef typename TContainer2D::Type Type;
  // apply left Givens rotation
  for (int j = M.beginCols(); j< M.endCols(); j++)
  {
    const Type aux1 = M.elt(i1, j), aux2 = M.elt(i2, j);
    M.elt(i1, j) = cosinus * aux1 + sinus * aux2;
    M.elt(i2, j) = cosinus * aux2 - sinus * aux1;
  }
}

Referenced by STK::IQr< Derived >::eraseCol(), and STK::IQr< Derived >::insertCol().

rightGivens()

template<class TContainer2D >

void STK::rightGivens	(	ArrayBase< TContainer2D > &	M,
		int	j1,
		int	j2,
		typename TContainer2D::Type const &	cosinus,
		typename TContainer2D::Type const &	sinus
	)

Apply Givens rotation.

Perform a right multiplication of the Container M with a Givens Array on the col1 and col2. col1 should be less than col2. The Array M is passed as const as we are using reference on the two cols we want to rotate.

See also

http://en.wikipedia.org/wiki/Givens_rotation

Parameters

M	the Container to multiply
j1,j2	the first and second columns
cosinus,sinus	the cosinus and sinus of the givens rotation

Definition at line 119 of file STK_Givens.h.

{
  typedef typename TContainer2D::Type Type;
  // Apply givens rotation
  for (int i = M.beginRows(); i < M.endRows(); i++)
  {
    const Type aux1 = M.elt(i, j1), aux2 = M.elt(i, j2);
    M.elt(i, j1) = cosinus * aux1 + sinus * aux2;
    M.elt(i, j2) = cosinus * aux2 - sinus * aux1;
  }
}

Referenced by STK::Svd< Array >::diag(), STK::SymEigen< SquareArray >::diagonalize(), STK::IQr< Derived >::eraseCol(), STK::IQr< Derived >::insertCol(), STK::leftEliminate(), STK::IQr< Derived >::pushBackCol(), and STK::Svd< Array >::rightEliminate().

syevr()

int STK::lapack::syevr ( char  jobz, char  range, char  uplo, int  n, Real *  a, int  lda, Real  vl, Real  vu, int  il, int  iu, Real  abstol, int *  m, Real *  w, Real *  z, int  ldz, int *  isuppz, Real *  work, int  lWork, int *  iwork, int  liwork  )

inline

wrapper of the LAPACK SYEVR routine.

SYEVR computes the eigenvalues of a symmetric square matrix.

Parameters

[in]	jobz	* CHARACTER*1 * = 'N': Compute eigenvalues only; * = 'V': Compute eigenvalues and eigenvectors. *
[in]	range	* CHARACTER*1 * = 'A': all eigenvalues will be found. * = 'V': all eigenvalues in the half-open interval (VL_,VU_] will be found. * = 'I': the IL_-th through IU_-th eigenvalues will be found. *
[in]	uplo	* CHARACTER*1 * = 'U': Upper triangle of A is stored; * = 'L': Lower triangle of A is stored. *
[in]	n	The order of the matrix A. N >= 0.
[in,out]	a	Real array, dimension (lda, N) * On entry, the symmetric matrix A. If UPLO_ = 'U', the leading * N-by-N upper triangular part of A contains the upper triangular part * of the matrix A. If UPLO_ = 'L', the leading N-by-N lower triangular * part of A contains the lower triangular part of the matrix A. On * exit, the lower triangle (if UPLO_='L') or the upper triangle * (if UPLO_='U') of A, including the diagonal, is destroyed. *
[in]	lda	The leading dimension of the array A. lda >= max(1,N).
[in]	vl,vu	* Real If RANGE_='V', the lower and upper bounds * of the interval to be searched for eigenvalues. VL_ < VU_. Not * referenced if RANGE_ = 'A' or 'I'. *
[in]	il,iu	* Integer If RANGE_='I', the indices (in ascending order) * of the smallest and largest eigenvalues to be returned. * 1 <= IL_ <= IU_ <= NL, if NL > 0; IL_ = 1 and IU_ = 0 if NL = 0. Not * referenced if RANGE_ = 'A' or 'V'. *
[in]	abstol	* The absolute error tolerance for the eigenvalues. An approximate * eigenvalue is accepted as converged when it is determined * to lie in an interval [a,b] of width less than or equal to * ABSTOL + EPS * max( |a|,|b| ) , * where EPS is the machine precision. If ABSTOL is less than or * equal to zero, then EPS*|T| will be used in its place, where * |T| is the 1-norm of the tridiagonal matrix obtained by reducing A * to tridiagonal form. * If high relative accuracy is important, set ABSTOL to SLAMCH( * 'Safe minimum' ). Doing so will guarantee that eigenvalues are * computed to high relative accuracy when possible in future * releases. The current code does not make any guarantees about * high relative accuracy, but future releases will. See J. Barlow * and J. Demmel, "Computing Accurate Eigensystems of Scaled Diagonally * Dominant Matrices", LAPACK Working Note #7, for a discussion of * which matrices define their eigenvalues to high relative accuracy. *

See also

"Computing Small Singular Values of Bidiagonal Matrices with Guaranteed High Relative Accuracy," by Demmel and Kahan, LAPACK Working Note #3.

Parameters

[out]	m	* The total number of eigenvalues found. 0 <= M <= NL. If RANGE_ * = 'A', M = NL, and if RANGE_ = 'I', M = IU_-IL_+1. *
[out]	w	* array, dimension (NL) * The first M elements contain the selected eigenvalues in * ascending order. *
[out]	z	* array, dimension (LDZ, max(1,M)) * If jobz_ = 'V', then if info = 0, the first M columns of Z contain * the orthonormal eigenvectors of the matrix A corresponding * to the selected eigenvalues, with the i-th column of Z holding * the eigenvector associated with W(i). If jobz_ = 'N', then Z is * not referenced. Note: the user must ensure that at least * max(1,M) columns are supplied in the array Z; if RANGE_ = 'V', * the exact value of M is not known in advance and an upper bound * must be used. Supplying N columns is always safe. *
[in]	ldz	* The leading dimension of the array Z. LDZ >= 1, and if jobz_ = * 'V', LDZ >= max(1,N). *
[out]	isuppz	array, dimension ( 2*max(1,M) ) * The support of the eigenvectors in Z, i.e., the indices indicating * the nonzero elements in Z. The i-th eigenvector is * nonzero only in elements ISUPPZ( 2*i-1 ) through ISUPPZ( 2*i ). *
[in,out]	work	Real array, dimension (MAX(1,lWork)) * On exit, if info = 0, work(1) returns the optimal lWork. *
[in]	lWork	The dimension of the array work * lWork >= max(1,26*N). For optimal efficiency, lWork >= (NB+6)*N, where * NB is the max of the blocksize for SSYTRD and SORMTR returned by ILAENV. * If lWork = -1, then a workspace query is assumed; the routine only * calculates the optimal sizes of the work and iWork arrays, * returns these values as the first entries of the work and iWork * arrays, and no error message related to lWork or LiWork is * issued by XERBLA. *
[in,out]	iwork	array, dimension (MAX(1,LiWork)) * On exit, if info = 0, iWork(1) returns the optimal lWork. *
[in]	liwork	The dimension of the array iWork. * LiWork >= max(1,10*N). If LiWork = -1, then a workspace query is * assumed; the routine only calculates the optimal sizes of the work and * iWork arrays, returns these values as the first entries of the work and * iWork arrays, and no error message related to lWork or LiWork is * issued by XERBLA. *

Returns

info

*  = 0:  successful exit
*  < 0:  if info = -i, the i-th argument had an illegal value
*  > 0:  Internal error
* 

Definition at line 606 of file STK_lapack_Util.h.

{
  int info = 0;
#ifdef STKUSELAPACK
#ifdef STKREALAREFLOAT
  ssyevr_(&jobz, &range, &uplo, &n, a, &lda, &vl, &vu, &il,
          &iu, &abstol, m, w, z, &ldz, isuppz, work,
          &lWork, iwork, &liwork, &info);
#else
  dsyevr_(&jobz, &range, &uplo, &n, a, &lda, &vl, &vu, &il,
          &iu, &abstol, m, w, z, &ldz, isuppz, work,
          &lWork, iwork, &liwork, &info);
#endif
#endif
  return info;
}

Referenced by STK::lapack::SymEigen< SquareArray >::runImpl().

transpose()

template<class TContainer2D >

TContainer2D & STK::transpose

(

ArrayBase< TContainer2D > &

A

)

The transpose method allow to transpose an array in place.

Parameters

A	the container to transpose

Definition at line 51 of file STK_transpose.h.

{
  if(A.rows()!=A.cols())
  {
    TContainer2D Aux = A.transpose();
    A = Aux;
  }
  else // transpose in place
  {
    for (int i= A.beginRows(); i< A.endRows(); i++)
      for (int j= i+1; j< A.endCols(); j++) { std::swap(A(i,j),A(j,i));}
  }
  return A.asDerived();
}

Referenced by STK::lapack::Svd::runImpl().