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

Classes
class	STK::EigenvaluesSymmetric
	The class EigenvaluesSymmetric compute the eigenvalue Decomposition of a symmetric Matrix. More...
class	STK::Array2D< Real >
	Specialization of the Array2D class for Real values. More...
class	STK::MatrixDiagonal
class	STK::MatrixLowerTriangular
	Declaration of the lower triangular matrix class. More...
class	STK::MatrixSquare
	Derivation of the MatrixSquare class for square arrays of Real. More...
class	STK::MatrixUpperTriangular
	Declaration of the upper triangular matrix class. More...
class	STK::ArrayHo< Real >
	Specialization of the templated class `ArrayHo` for Real. More...
class	STK::Qr
	The class Qr perform the QR decomposition of a Matrix. More...
class	STK::Svd
	The class Svd compute the Singular Value Decomposition of a Matrix with the Golub-Reinsch Algorithm. More...
class	STK::BinOpBase< Op, ExpLeft, ExpRight >
	Binary Operator bbase class. More...
class	STK::BinOpBase< Op, Real, ExpRight >
	Specialized Binary Operator Base class for the case Real Op Exp. More...
class	STK::BinOpBase< Op, ExpLeft, Real >
	Specialized Binary Operator Base class for the case Exp Op Real. More...
class	STK::BinOpBase< Mult, Matrix, ExpRight >
	Specialized Binary Operator Base class for the case Matrix * ExpRight. More...
class	STK::BinOpBase< Mult, ExpLeft, Matrix >
	Specialized Binary Operator Base class for the case ExpLeft * Matrix. More...
class	STK::BinOpBase< Mult, MatrixSquare, ExpRight >
	Specialized Binary Operator Base class for the case Matrix * ExpRight. More...
class	STK::BinOpBase< Mult, ExpLeft, MatrixSquare >
	Specialized Binary Operator Base class for the case ExpLeft * MatrixSquare. More...
class	STK::BinOpBase< Mult, MatrixUpperTriangular, ExpRight >
	Specialized Binary Operator Base class for the case MatrixUpperTriangular * ExpRight. More...
class	STK::BinOpBase< Mult, ExpLeft, MatrixUpperTriangular >
	Specialized Binary Operator Base class for the case ExpLeft * MatrixUpperTriangular. More...
class	STK::BinOpBase< Mult, MatrixLowerTriangular, ExpRight >
	Specialized Binary Operator Base class for the case MatrixLowerTriangular * ExpRight. More...
class	STK::BinOpBase< Mult, ExpLeft, MatrixLowerTriangular >
	Specialized Binary Operator Base class for the case ExpLeft * MatrixLowerTriangular. More...
class	STK::BinOp< Op, ExpLeft, ExpRight >
	Binary Operator class derived from BinOpBase class which is specialized. More...
class	STK::UnOp< Op, Exp >
	UnOp class for unary operators. More...
struct	STK::Plus
	Binary Operator Plus. More...
struct	STK::Minus
	Binary Operator Minus. More...
struct	STK::Mult
	Binary Operator Mult. More...
struct	STK::Div
	Binary Operator Div. More...
struct	STK::Uplus
	Unary Operator Plus. More...
struct	STK::Uminus
	Unary Operator Minus. More...
class	STK::Array1D< Real >
	Specialization of the templated class `Array1D` for Real. More...
Typedefs
typedef Array2D< Real >	STK::Matrix
	Specialization of the Array2D class for Real values.
typedef ArrayHo< Real >	STK::Point
	final class for a Real horizontal container.
typedef Array1D< Real >	STK::Vector
	Real one dimensional Vertical Array.
Functions
Real	STK::compGivens (Real const &y, Real const &z, Real &cosinus, Real &sinus)
	Compute Givens rotation.
template<class TContainerHo , class TContainerVe , class TContainer2D >
void	STK::rightGivens (const ITContainer2D< Real, TContainerHo, TContainerVe, TContainer2D > &M, Integer const &col1, Integer const &col2, Real const &cosinus, Real const &sinus)
	Apply Givens rotation.
template<class TContainerHo , class TContainerVe , class TContainer2D >
void	STK::leftGivens (const ITContainer2D< Real, TContainerHo, TContainerVe, TContainer2D > &M, Integer const &row1, Integer const &row2, Real const &cosinus, Real const &sinus)
	left multiplication by a Givens Matrix.
template<class TContainerHo , class TContainerVe , class TContainer2D >
void	STK::gramSchmidt (ITContainer2D< Real, TContainerHo, TContainerVe, TContainer2D > &A)
	The gramSchmidt method perform the Gram Schmidt orthonormalization of a ITContainer2D of Real.
template<class TContainer1D >
Real	STK::house (ITContainer1D< Real, TContainer1D > &x)
	Compute the Householder vector v of a vector x.
template<class TContainer1D_1 , class TContainer1D_2 >
Real	STK::dotHouse (const ITContainer1D< Real, TContainer1D_1 > &x, const ITContainer1D< Real, TContainer1D_2 > &v)
	dot product with a Householder vector.
template<class TContainerHo , class TContainerVe , class TContainer2D , class TContainer1D >
void	STK::leftHouseholder (const ITContainer2D< Real, TContainerHo, TContainerVe, TContainer2D > &M, const ITContainer1D< Real, TContainer1D > &v)
	left multiplication by a Householder vector.
template<class TContainerHo , class TContainerVe , class TContainer2D , class TContainer1D >
void	STK::rightHouseholder (const ITContainer2D< Real, TContainerHo, TContainerVe, TContainer2D > &M, const ITContainer1D< Real, TContainer1D > &v)
	right multiplication by a Householder vector.
template<class TContainerHo , class TContainerVe , class TContainer2D >
void	STK::leftHouseholder (const ITContainer2D< Real, TContainerHo, TContainerVe, TContainer2D > &M, Matrix const &H)
	left multiplication by a Householder Matrix.
template<class TContainerHo , class TContainerVe , class TContainer2D >
void	STK::rightHouseholder (const ITContainer2D< Real, TContainerHo, TContainerVe, TContainer2D > &M, Matrix const &H)
	left multiplication by a Householder Matrix.
template<class Container1D >
Real	STK::sum (ITContainer1D< Real, Container1D > const &x)
	Sum the element of the container.
template<class Container1D1 , class Container1D2 >
Real	STK::weightedSum (ITContainer1D< Real, Container1D1 > const &x, ITContainer1D< Real, Container1D2 > const &w)
	Weighted sum of the element of the container.
template<class Container1D >
Real	STK::normInf (ITContainer1D< Real, Container1D > const &x)
	Compute the infinite norm.
template<class Container1D1 , class Container1D2 >
Real	STK::weightedNormInf (ITContainer1D< Real, Container1D1 > const &x, ITContainer1D< Real, Container1D2 > const &w)
	Compute the weighted infinite norm.
template<class Container1D >
Real	STK::normTwo (ITContainer1D< Real, Container1D > const &x)
	compute the norm two
template<class Container1D1 , class Container1D2 >
Real	STK::weightedNormTwo (ITContainer1D< Real, Container1D1 > const &x, ITContainer1D< Real, Container1D2 > const &w)
	compute the weighted norm two
template<class Container1D >
Real	STK::normTwo2 (ITContainer1D< Real, Container1D > const &x)
	Compute the squared norm two.
template<class Container1D1 , class Container1D2 >
Real	STK::weightedNormTwo2 (ITContainer1D< Real, Container1D1 > const &x, ITContainer1D< Real, Container1D2 > const &w)
	Compute the squared norm two.
template<class Container1D1 , class Container1D2 >
Real	STK::dot (ITContainer1D< Real, Container1D1 > const &x, ITContainer1D< Real, Container1D2 > const &y)
	Compute the dot product.
template<class Container1D1 , class Container1D2 , class Container1D3 >
Real	STK::weightedDot (ITContainer1D< Real, Container1D1 > const &x, ITContainer1D< Real, Container1D2 > const &y, ITContainer1D< Real, Container1D3 > const &w)
	Compute the dot product.
template<class Container1D1 , class Container1D2 >
Real	STK::dist (ITContainer1D< Real, Container1D1 > const &x, ITContainer1D< Real, Container1D2 > const &y)
	Compute the distance between two vectors.
template<class Container1D1 , class Container1D2 , class Container1D3 >
Real	STK::weightedDist (ITContainer1D< Real, Container1D1 > const &x, ITContainer1D< Real, Container1D2 > const &y, ITContainer1D< Real, Container1D3 > const &w)
	Compute the weighted distance between two vectors.
Real	STK::trace (MatrixSquare const &A)
	trace of a square matrix
MatrixSquare *	STK::mult (MatrixSquare const &A, MatrixSquare const &B)
	MatrixSquare multiplication.
Matrix *	STK::mult (Matrix const &A, Matrix const &B)
	Matrix multiplication.
Matrix *	STK::mult (MatrixSquare const &A, Matrix const &B)
	Matrix multiplication.
template<class Container1D >
Vector *	STK::mult (Matrix const &A, ITContainer1D< Real, Container1D > const &X)
	Matrix multiplication by a Vector.
template<class TContainerHo , class TContainerVe , class TContainer2D , class TContainer1D >
void	STK::mult (ITContainer2D< Real, TContainerHo, TContainerVe, TContainer2D > const &A, ITContainer1D< Real, TContainer1D > const &X, ITContainer1D< Real, TContainer1D > &Y)
	Matrix multiplication by a Vector.
template<class Container1D >
void	STK::mult (MatrixSquare const &A, ITContainer1D< Real, Container1D > const &X, ITContainer1D< Real, Container1D > &Y)
	Square Matrix multiplication by a Vector.
template<class TContainerHo , class TContainerVe , class TContainer2D , class TContainer1D >
void	STK::multLefTranspose (ITContainer2D< Real, TContainerHo, TContainerVe, TContainer2D > const &A, ITContainer1D< Real, TContainer1D > const &X, ITContainer1D< Real, TContainer1D > &Y)
	Matrix multiplication by a Vector.
template<class Container1D >
void	STK::multLeftTranspose (MatrixSquare const &A, ITContainer1D< Real, Container1D > const &X, ITContainer1D< Real, Container1D > &Y)
	transposed square Matrix multiplication by a Vector.
Matrix *	STK::multLeftTranspose (Matrix const &A, Matrix const &B)
	Matrix multiplication.
template<class Container1D >
Vector *	STK::multLeftTranspose (Matrix const &A, ITContainer1D< Real, Container1D > const &X)
	Matrix by Vector multiplication.
Matrix *	STK::weightedMultLeftTranspose (Matrix const &A, Matrix const &B, Vector const &weights)
	Weighted matrix multiplication.
MatrixSquare *	STK::multLeftTranspose (Matrix const &A)
	Matrix multiplication.
MatrixSquare *	STK::multRightTranspose (Matrix const &A)
	Matrix multiplication.
MatrixSquare *	STK::weightedMultLeftTranspose (Matrix const &A, Vector const &weights)
	Weighted matrix multiplication.
template<class TContainerHo , class TContainerVe , class TContainer2D >
Real	STK::normInf (ITContainer2D< Real, TContainerHo, TContainerVe, TContainer2D > const &A)
	Compute the infinity norm of a 2D container.
template<class TContainerHo_1 , class TContainerVe_1 , class TContainer2D_1 , class TContainerHo_2 , class TContainerVe_2 , class TContainer2D_2 >
void	STK::transpose (ITContainer2D< Real, TContainerHo_1, TContainerVe_1, TContainer2D_1 > const &A, ITContainer2D< Real, TContainerHo_2, TContainerVe_2, TContainer2D_2 > &At)
	transpose a matrix
template<class TContainerHo , class TContainerVe , class TContainer2D >
ITContainer2D< Real, TContainerHo, TContainerVe, TContainer2D > &	STK::transpose (ITContainer2D< Real, TContainerHo, TContainerVe, TContainer2D > &Q)
	Transpose a Matrix.
template<class TContainerHo1 , class TContainerVe1 , class TContainer2D1 , class TContainerHo2 , class TContainerVe2 , class TContainer2D2 , class TContainerHo3 , class TContainerVe3 , class TContainer2D3 >
void	STK::mult (ITContainer2D< Real, TContainerHo1, TContainerVe1, TContainer2D1 > &C, const ITContainer2D< Real, TContainerHo2, TContainerVe2, TContainer2D2 > &A, const ITContainer2D< Real, TContainerHo3, TContainerVe3, TContainer2D3 > &B)
	Matrix multiplication.
Vector *	STK::multLeftTranspose (Matrix const &A, Vector const &X)
	*

Detailed Description

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

The main contribution of the Algebra project are

standard Real containers like Point, Vector, Matrix, ... using containers specialization
templated expressions evaluation for scientific computing
See also:
http://osl.iu.edu/~tveldhui/papers/techniques/
usual methods of linear algebra
QR and SVD decompositions for matrices
an Eigenvalue decomposition for symmetric matrices and a generalized inversion method for such matrices.

Typedef Documentation

typedef Array2D<Real> STK::Matrix

Specialization of the Array2D class for Real values.

A Matrix is a column oriented 2D container of Real.

Definition at line 52 of file STK_Matrix.h.

typedef ArrayHo<Real> STK::Point

final class for a Real horizontal container.

A Point is a row oriented 1D container of Real. It support templates expression for partial evaluation at compile-time.

Definition at line 51 of file STK_Point.h.

typedef Array1D<Real> STK::Vector

Real one dimensional Vertical Array.

A Vector is a column oriented 1D container of Real. It supports templates expression or partial evaluation at compile-time.

Definition at line 53 of file STK_Vector.h.

Function Documentation

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

Compute Givens rotation.

Compute the Givens rotation

$\begin{bmatrix} c & s \\ -s & c \end{bmatrix} \begin{bmatrix} a \\ b \end{bmatrix} = \begin{bmatrix} r \\ 0 \end{bmatrix}.$

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 44 of file STK_Givens.cpp.

References STK::abs(), and STK::sign().

Referenced by STK::Qr::eraseCol(), STK::Qr::insertCol(), and STK::Qr::pushBackCol().

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

template<class TContainerHo , class TContainerVe , class TContainer2D >

void STK::rightGivens	(	const ITContainer2D< Real, TContainerHo, TContainerVe, TContainer2D > &	M,
		Integer const &	col1,
		Integer const &	col2,
		Real const &	cosinus,
		Real const &	sinus
	)

Apply Givens rotation.

Perform a right multiplication of the Container M with a Givens Matrix on the col1 and col2. col1 should be less than col2. The Matrix 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
col1	the first col
col2	the second col
cosinus	the cosinus of the givens rotation
sinus	the sinus of the givens rotation

Definition at line 104 of file STK_Givens.h.

Referenced by STK::Svd::diag(), STK::EigenvaluesSymmetric::diagonalize(), STK::Qr::eraseCol(), STK::Qr::insertCol(), STK::Svd::leftEliminate(), STK::Qr::pushBackCol(), and STK::Svd::rightEliminate().

{
  // A ref on the col1 of M
  TContainerVe Mcol1(M.asLeaf(), M.rangeVe(), col1);
  // A ref on the col2 of M
  TContainerVe Mcol2(M.asLeaf(), M.rangeVe(), col2);
  // Apply givens rotation
  for (Integer i = M.firstRow(); i <=M.lastRow(); i++)
  {
    Real aux1 = Mcol1[i], aux2 = Mcol2[i];
    Mcol1[i] = cosinus * aux1 + sinus * aux2;
    Mcol2[i] = cosinus * aux2 - sinus * aux1;
  }
}

template<class TContainerHo , class TContainerVe , class TContainer2D >

void STK::leftGivens	(	const ITContainer2D< Real, TContainerHo, TContainerVe, TContainer2D > &	M,
		Integer const &	row1,
		Integer const &	row2,
		Real const &	cosinus,
		Real const &	sinus
	)

left multiplication by a Givens Matrix.

Perform a left multiplication of the matrix M with a Givens Matrix on the row1 and row2. row1 should be less than row2. The Matrix 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
row1	the first row
row2	the second row
cosinus	the cosinus of the givens rotation
sinus	the sinus of the givens rotation

Definition at line 148 of file STK_Givens.h.

Referenced by STK::Qr::eraseCol(), and STK::Qr::insertCol().

{
  // A ref on the row iter of R_
  TContainerHo Mrow1(M.asLeaf(), M.rangeHo(), row1);
  // A ref on the row iter1 of R_
  TContainerHo Mrow2(M.asLeaf(), M.rangeHo(), row2);
  // apply right Givens rotation
  for (Integer j = M.firstCol(); j<= M.lastCol(); j++)
  {
    Real aux1 = Mrow1[j], aux2 = Mrow2[j];
    Mrow1[j] = cosinus * aux1 + sinus * aux2;
    Mrow2[j] = cosinus * aux2 - sinus * aux1;
  }
}

template<class TContainerHo , class TContainerVe , class TContainer2D >

void STK::gramSchmidt ( ITContainer2D< Real, TContainerHo, TContainerVe, TContainer2D > & A )

The gramSchmidt method perform the Gram Schmidt orthonormalization of a ITContainer2D of Real.

Parameters:

A	the container to orthnormalize

Definition at line 56 of file STK_GramSchmidt.h.

References STK::dot(), STK::IContainer2D::firstCol(), STK::IContainer2D::lastCol(), STK::norm(), and STK::normTwo().

Referenced by STK::LinearAAModel::simul().

{
  // get dimensions
  const Integer first_col = A.firstCol(), last_col = A.lastCol();
  // orthonormalize
  for (Integer j= first_col; j<= last_col; j++)
  {
    for( Integer i= first_col; i < j; i++)
    {
      A[j] -= A[i] * dot(A[i], A[j]);
    }
    // normalize
    const Real norm = normTwo(A[j]);
    if (norm) A[j] /= norm;
  }
}

template<class TContainer1D >

Real STK::house ( ITContainer1D< Real, TContainer1D > & x )

Compute the Householder vector v of a vector x.

Given a vector x, compute the vector v of the matrix of Householder $ P=I-2vv'/(v'v) $ such that $ Px = v1 e_1 $ . The vector v is of the form : $ (1,x_2/s,...,x_n/s)' $ and is stored in x. The value 1 is skipped and $\beta = -2/(v'v)$ 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 64 of file STK_Householder.h.

References STK::ITContainer1D< TYPE, TContainer1D >::first(), STK::ITContainer1D< TYPE, TContainer1D >::front(), STK::ITContainer1D< TYPE, TContainer1D >::last(), STK::normInf(), and STK::sign().

Referenced by STK::Svd::bidiag(), STK::Qr::qr(), and STK::EigenvaluesSymmetric::tridiagonalize().

{
  // compute L^{\infty} norm of X
  Real scale  = normInf(x);
  // first and last index fo the essential Householder vector
  Integer first = x.first()+1, last = x.last();
  // result and norm2 of X
  Real v1, norm2 = 0.0;
  // normalize the vector 
  if (scale)  // if not 0.0
  {
    for (Integer i=first; i<=last; i++)
    { x[i] /= scale; norm2 += x[i]*x[i];}
  }
  // check if the lower part is significative
  if (norm2 < Arithmetic<Real>::epsilon())
  {
    v1 = x.front(); x.front() = 0.0; // beta = 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
    for (Integer i=first; i<=last; i++) x[i] /= s;
  }
  return v1;
}

template<class TContainer1D_1 , class TContainer1D_2 >

Real STK::dotHouse	(	const ITContainer1D< Real, TContainer1D_1 > &	x,
		const ITContainer1D< Real, TContainer1D_2 > &	v
	)

dot product with a Householder vector.

Scalar product of a TContainer1D with a Householder vector d = < x,v>. The first composant of a Householder vector is 1.0

Parameters:

x	first vector
v	the Householder vector

Definition at line 107 of file STK_Householder.h.

References STK::ITContainer1D< TYPE, TContainer1D >::first(), STK::ITContainer1D< TYPE, TContainer1D >::last(), and STK::sum().

Referenced by STK::leftHouseholder(), and STK::rightHouseholder().

{
  // first and last index fo the essential Householder vector
  const Integer first = v.first()+1, last = v.last();
  // compute the product
  Real sum = x[first-1] /* *1.0 */;
  for (Integer i=first; i<=last; i++) sum += x[i] * v[i];
  // return <x,v>
  return(sum);
}

template<class TContainerHo , class TContainerVe , class TContainer2D , class TContainer1D >

void STK::leftHouseholder	(	const ITContainer2D< Real, TContainerHo, TContainerVe, TContainer2D > &	M,
		const ITContainer1D< Real, TContainer1D > &	v
	)

left multiplication by a Householder vector.

Perform a left multiplication of the matrix M with a Householder matrix $ H=I+beta vv' $ . Overwrite M with HM.

Parameters:

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

Definition at line 133 of file STK_Householder.h.

References beta(), and STK::dotHouse().

Referenced by STK::Svd::bidiag(), and STK::leftHouseholder().

{
  // get beta
  Real beta = v.front();
  if (beta)
  {
    // get range of the Householder vector
    Range range_ve = v.range();
    // Multiplication of the cols by P=I+beta vv'
    for (Integer j=M.firstCol(); j<=M.lastCol(); j++)
    {
      // a ref on the jth column of M
      TContainerVe Mj(M.asLeaf(), range_ve, j);
      // Computation of aux=beta* <v,M^j>
      Real aux =  dotHouse( Mj, v) * beta;
      // updating row X.first()
      Mj.front() += aux;
      // essential range of v
      const Integer first =  v.first()+1, last = v.last();
      // Computation of M^j + beta <v,M^j>  v = M^j + aux v
      for (Integer i=first; i<=last; i++)
        Mj[i] +=  v[i] * aux;
    }
  }
}

template<class TContainerHo , class TContainerVe , class TContainer2D , class TContainer1D >

void STK::rightHouseholder	(	const ITContainer2D< Real, TContainerHo, TContainerVe, TContainer2D > &	M,
		const ITContainer1D< Real, TContainer1D > &	v
	)

right multiplication by a Householder vector.

Perform a right multiplication of the matrix M with a Householder matrix $ H=I+beta vv' $ . Overwrite M with MH.

Parameters:

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

Definition at line 181 of file STK_Householder.h.

References beta(), and STK::dotHouse().

Referenced by STK::Svd::bidiag(), and STK::rightHouseholder().

{
  // get beta
  Real beta = v.front();
  if (beta)
  {
    // Multiplication of the cols by P=I+beta vv'
    for (Integer i=M.firstRow(); i<=M.lastRow(); i++)
    {
      // a ref on the ith row of M
      TContainerHo Mi(M.asLeaf(), v.range(), i);
      // Computation of aux=beta* <v,M_i>
      Real aux =  dotHouse( Mi, v) * beta;
      // updating column X.first()
      Mi.front() += aux;
      // essential range of v
      const Integer first =  v.first()+1, last = v.last();
      // Computation of M_i + beta <v,M_i>  v = M_i + aux v'
      for (Integer i=first; i<=last; i++)
        Mi[i] +=  v[i] * aux;
    }
  }
}

template<class TContainerHo , class TContainerVe , class TContainer2D >

void STK::leftHouseholder	(	const ITContainer2D< Real, TContainerHo, TContainerVe, TContainer2D > &	M,
		Matrix const &	H
	)

left multiplication by a Householder Matrix.

Perform a left multiplication of the Matrix M with a Householder Marix H. M <- HM with H = I + WZ'. The Householder vectors are stored in the columns of H.

Parameters:

M	the matrix to multiply
H	the Householder Matrix

Definition at line 227 of file STK_Householder.h.

References STK::IContainer2D::firstCol(), STK::IContainer2D::lastCol(), STK::IContainer2D::lastRow(), STK::leftHouseholder(), and STK::min().

{
  // compute the number of iterations
  Integer first = H.firstCol()
        , last = min( H.lastCol(), H.lastRow());
  // get range of the first Householder vector
  Range range_ve(last, H.lastRow());
  // iterations
  for (Integer j=last; j>=first; j--)
  {
    // apply left Householder vector to M
    leftHouseholder(M, H(range_ve, j));
    // decrease range of the Householder vector
    range_ve.decFirst();
  }
}

template<class TContainerHo , class TContainerVe , class TContainer2D >

void STK::rightHouseholder	(	const ITContainer2D< Real, TContainerHo, TContainerVe, TContainer2D > &	M,
		Matrix const &	H
	)

left multiplication by a Householder Matrix.

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 Matrix to multiply
H	the Householder Matrix

Definition at line 264 of file STK_Householder.h.

References STK::IContainer2D::firstCol(), STK::IContainer2D::lastCol(), STK::IContainer2D::lastRow(), STK::min(), and STK::rightHouseholder().

{
  // compute the number of iterations
  Integer first = H.firstCol()
        , last = min( H.lastCol(), H.lastRow());
  // get range of the first Householder vector
  Range range_ve(last, H.lastRow());
  // iterations
  for (Integer j=last; j>=first; j--)
  {
    // apply left Householder vector to M
    rightHouseholder(M, H(range_ve, j));
    // decrease range of the Householder vector
    range_ve.decFirst();
  }
}

template<class Container1D >

Real STK::sum ( ITContainer1D< Real, Container1D > const & x )

Sum the element of the container.

Sum of the element of the Container1D x

$s= \sum_{i=1}^n x_i$

Parameters:

[in] x vector to treat

Returns:: the sum of the element of the container

Definition at line 57 of file STK_LinAlgebra1D.h.

References STK::ITContainer1D< TYPE, TContainer1D >::first(), and STK::ITContainer1D< TYPE, TContainer1D >::last().

Referenced by STK::Funct::betaRatio_cf(), STK::Funct::betaRatio_sr(), STK::Funct::coefs_even_se(), STK::Funct::coefs_odd_se(), STK::ITStatModel< Real, Point, Vector, Matrix >::compLogLikelihood(), STK::Stat::Univariate< Real, TContainer1D >::compStatistics(), STK::LocalVariance::computeCovarianceMatrices(), STK::LocalVariance::computeWeightedCovarianceMatrices(), STK::Stat::Univariate< Real, TContainer1D >::compWeightedStatistics(), STK::Funct::dev0(), STK::EigenvaluesSymmetric::diagonalize(), STK::dot(), STK::dotHouse(), STK::Funct::expm1(), STK::Funct::g0(), STK::Funct::gammaRatio_ae(), STK::Funct::gammaRatio_sr(), STK::EigenvaluesSymmetric::ginv(), STK::Funct::lanczosSerie(), STK::Law::MultivariateNormal< Container1D >::logLikelihood(), STK::Stat::mean(), STK::mult(), STK::multLeftTranspose(), STK::multRightTranspose(), STK::Funct::p1evl(), STK::Funct::polevl(), STK::Funct::serie_up(), STK::sumAlternateSerie(), STK::sumSerie(), STK::trace(), STK::Stat::variance(), STK::Stat::varianceWithFixedMean(), STK::weightedDot(), STK::weightedMultLeftTranspose(), and STK::weightedSum().

{
  // get dimensions
  const Integer first = x.first(), last = x.last();
  // compute the sum
  Real sum = 0.0;
  for (Integer i=first; i<=last; i++)
    sum += x[i];
  return (sum);
}

template<class Container1D1 , class Container1D2 >

Real STK::weightedSum	(	ITContainer1D< Real, Container1D1 > const &	x,
		ITContainer1D< Real, Container1D2 > const &	w
	)

Weighted sum of the element of the container.

Sum of the weighted elements of the Container1D x

$s= \sum_{i=1}^n w_i x_i$

Parameters:

[in]	x	vector to treat
	w	the weights to apply

Returns:: the weighted sum of the container

Definition at line 79 of file STK_LinAlgebra1D.h.

References STK::ITContainer1D< TYPE, TContainer1D >::first(), STK::Range::isIn(), STK::ITContainer1D< TYPE, TContainer1D >::last(), STK::ITContainer1D< TYPE, TContainer1D >::range(), and STK::sum().

{
#ifdef STK_DEBUG
    if (!x.range().isIn(w.range()))
      throw std::runtime_error("In weightedSum(x, w) "
                               "x.range() not include in w.range()");
#endif
  // dimensions
  const Integer first = x.first(), last = x.last();
  // compute the weighted sum
  Real sum = 0.0;
  for (Integer i=first; i<=last; i++)
    sum += w[i] * x[i];
  return (sum);
}

template<class Container1D >

Real STK::normInf ( ITContainer1D< Real, Container1D > const & x )

Compute the infinite norm.

Compute the maximal absolute value of the container x

$s= \max_{i=1}^n |x_i|$

Parameters:

[in] x vector to treat

Returns:: the infinite norm f the container

Definition at line 106 of file STK_LinAlgebra1D.h.

References STK::abs(), STK::ITContainer1D< TYPE, TContainer1D >::first(), STK::ITContainer1D< TYPE, TContainer1D >::last(), and STK::max().

Referenced by STK::house(), STK::normTwo(), and STK::normTwo2().

{
  // get dimensions
  const Integer first = x.first(), last = x.last();
  // compute the maxmal value
  Real scale = 0.0;
  for (Integer i=first; i<=last; i++)
    scale = max(scale, abs(x[i]));
  return (scale);
}

template<class Container1D1 , class Container1D2 >

Real STK::weightedNormInf	(	ITContainer1D< Real, Container1D1 > const &	x,
		ITContainer1D< Real, Container1D2 > const &	w
	)

Compute the weighted infinite norm.

Compute the maximal absolute weighted value of the container x

$s= \max_{i=1}^n |w_i x_i|$

Parameters:

[in]	x	vector to treat
	w	the weights to apply

Returns:: the weighted infinite norm

Definition at line 127 of file STK_LinAlgebra1D.h.

References STK::abs(), STK::ITContainer1D< TYPE, TContainer1D >::first(), STK::Range::isIn(), STK::ITContainer1D< TYPE, TContainer1D >::last(), STK::max(), and STK::ITContainer1D< TYPE, TContainer1D >::range().

Referenced by STK::weightedNormTwo(), and STK::weightedNormTwo2().

{
#ifdef STK_DEBUG
    if (!x.range().isIn(w.range()))
      throw std::runtime_error("In weightedNormInf(x, w) "
                               "x.range() not include in w.range()");
#endif
  // get dimensions
  const Integer first = x.first(), last = x.last();
  // compute weighted norm inf
  Real scale = 0.0;
  for (Integer i=first; i<=last; i++)
    scale = max(scale, abs(w[i]*x[i]));
  return (scale);
}

template<class Container1D >

Real STK::normTwo ( ITContainer1D< Real, Container1D > const & x )

compute the norm two

Compute the norm of the container x avoiding overflow

$\|x\| = \sqrt{\sum_{i=1}^n x^2_i }$

Parameters:

[in] x vector to treat

Returns:: the norm two of the container

Definition at line 154 of file STK_LinAlgebra1D.h.

References STK::ITContainer1D< TYPE, TContainer1D >::first(), STK::ITContainer1D< TYPE, TContainer1D >::last(), STK::norm(), and STK::normInf().

Referenced by STK::gramSchmidt().

{
  // compute the maximal value of x
  Real scale =normInf(x), norm =0.0;
  if (scale)
  {
    // get dimensions
    const Integer first = x.first(), last = x.last();
    // sum squared normalized values
    for (Integer i = first; i<=last; i++)
    {
      const Real aux = x[i]/scale;
      norm += aux * aux;
    }
  }
  // rescale sum
  return (Real(sqrt(double(norm)))*scale);
}

template<class Container1D1 , class Container1D2 >

Real STK::weightedNormTwo	(	ITContainer1D< Real, Container1D1 > const &	x,
		ITContainer1D< Real, Container1D2 > const &	w
	)

compute the weighted norm two

Compute the weighted norm of the container x avoiding overflow

$\|x\| = \sqrt{\sum_{i=1}^n w_i x^2_i }$

Parameters:

[in]	x	vector to treat
	w	the weights to apply

Returns:: the weighted two norm

Definition at line 183 of file STK_LinAlgebra1D.h.

References STK::ITContainer1D< TYPE, TContainer1D >::first(), STK::Range::isIn(), STK::ITContainer1D< TYPE, TContainer1D >::last(), STK::ITContainer1D< TYPE, TContainer1D >::range(), and STK::weightedNormInf().

{
#ifdef STK_DEBUG
    if (!x.range().isIn(w.range()))
      throw std::runtime_error("In weightedNormTwo(x, w) "
                               "x.range() not include in w.range()");
#endif
  // compute the maximal value of x
  Real scale = weightedNormInf(x, w), norm2 =0.0;
  if (scale)
  {
    // get dimensions
    const Integer first = x.first(), last = x.last();
    // compute norm2
    for (Integer i = first; i<=last; i++)
    {
      const Real aux = (w[i]*x[i])/scale;
      norm2 += aux * aux;
    }
  }
  // rescale sum
  return (Real(sqrt(double(norm2)))*scale);
}

template<class Container1D >

Real STK::normTwo2 ( ITContainer1D< Real, Container1D > const & x )

Compute the squared norm two.

Compute the square norm of the Container1D x avoiding overflow

$\|x\|^2 = \sum_{i=1}^n x^2_i$

Parameters:

[in] x vector to treat

Definition at line 218 of file STK_LinAlgebra1D.h.

References STK::ITContainer1D< TYPE, TContainer1D >::first(), STK::ITContainer1D< TYPE, TContainer1D >::last(), STK::norm(), and STK::normInf().

{
  Real scale =normInf(x), norm =0.0;
  if (scale)
  {
    // get dimensions
    const Integer first = x.first(), last = x.last();
    // sum squared normalized values
    for (Integer i = first; i<=last; i++)
    {
      Real aux = x[i]/scale;
      norm += aux * aux;
    }
  }
  // scale result
  return (norm*scale*scale);
}

template<class Container1D1 , class Container1D2 >

Real STK::weightedNormTwo2	(	ITContainer1D< Real, Container1D1 > const &	x,
		ITContainer1D< Real, Container1D2 > const &	w
	)

Compute the squared norm two.

Compute the weighted square norm two of the container x avoiding overflow

$\|x\|^2 = \sum_{i=1}^n w_i x^2_i$

Parameters:

[in]	x	vector to treat
	w	the weights to apply

Returns:: the weighted square norm two

Definition at line 247 of file STK_LinAlgebra1D.h.

References STK::ITContainer1D< TYPE, TContainer1D >::first(), STK::Range::isIn(), STK::ITContainer1D< TYPE, TContainer1D >::last(), STK::norm(), STK::ITContainer1D< TYPE, TContainer1D >::range(), and STK::weightedNormInf().

{
#ifdef STK_DEBUG
  if (!x.range().isIn(w.range()))
    throw std::runtime_error("In weightedNormTwo2(x, w) "
                             "x.range() not include in w.range()");
#endif
  Real scale =weightedNormInf(x, w), norm =0.0;
  if (scale)
  {
    // get dimensions
    const Integer first = x.first(), last = x.last();
    // compute norm2
    for (Integer i = first; i <= last; i++)
    {
      const Real aux = x[i]/scale;
      norm += w[i] * aux * aux;
    }
  }
  // scale result
  return (norm*scale*scale);
}

template<class Container1D1 , class Container1D2 >

Real STK::dot	(	ITContainer1D< Real, Container1D1 > const &	x,
		ITContainer1D< Real, Container1D2 > const &	y
	)

Compute the dot product.

Dot product of the vector x and the vector y: d = <x, y>.

$<x,y> = \sum_{i=1}^n x_i y_i$

The common range of the vectors is used. The value outside the range are thus interpreted as zero.

Parameters:

[in]	x	first vector
[in]	y	second vector

Returns:: the dot product of the two vectors

Definition at line 285 of file STK_LinAlgebra1D.h.

References STK::ITContainer1D< TYPE, TContainer1D >::first(), STK::isOdd(), STK::ITContainer1D< TYPE, TContainer1D >::last(), STK::max(), STK::min(), and STK::sum().

Referenced by STK::Qr::compQ(), STK::Svd::compU(), STK::Svd::compV(), STK::gramSchmidt(), STK::Qr::insertCol(), STK::mult(), STK::Qr::pushBackCol(), and STK::Qr::qr().

{
  // compute the valid range
  const Integer first = max(x.first(), y.first()) , last = min(x.last(), y.last());
  // compute the sum product
  Real sum=0.0;
  for (Integer i = first; i<last; i+=2)
    sum += x[i] * y[i] + x[i+1] * y[i+1];
  // check if last is odd
  if (isOdd(last)) sum+=x[last]*y[last];
  return (sum);
}

template<class Container1D1 , class Container1D2 , class Container1D3 >

Real STK::weightedDot	(	ITContainer1D< Real, Container1D1 > const &	x,
		ITContainer1D< Real, Container1D2 > const &	y,
		ITContainer1D< Real, Container1D3 > const &	w
	)

Compute the dot product.

Weighted dot product of the vector x and the vector y:

$<x,y> = \sum_{i=1}^n w_i x_i y_i.$

The common range of the vectors is used. The value outside the range are thus interpreted as zero.

Parameters:

[in]	x	first vector
[in]	y	second vector
[in]	w	the weights to apply

Returns:: the weighted dot product of the two vectors

Definition at line 314 of file STK_LinAlgebra1D.h.

References STK::ITContainer1D< TYPE, TContainer1D >::first(), STK::Range::isIn(), STK::isOdd(), STK::ITContainer1D< TYPE, TContainer1D >::last(), STK::max(), STK::min(), STK::ITContainer1D< TYPE, TContainer1D >::range(), and STK::sum().

{
  // compute the valid range
  const Integer first = max(x.first(), y.first()) , last = min(x.last(), y.last());
#ifdef STK_DEBUG
  if (!Range(first,last).isIn(w.range()))
    throw std::runtime_error("In weightedDot(x, w) "
                             "Range(first,last) not include in w.range()");
#endif
  // compute the sum product
  Real sum=0.0;
  for (Integer i = first; i<last; i+=2)
    sum += w[i]*x[i] * y[i] + w[i+1]*x[i+1] * y[i+1];
  // check if last is odd
  if (isOdd(last)) sum += w[last]*x[last]*y[last];
  return (sum);
}

template<class Container1D1 , class Container1D2 >

Real STK::dist	(	ITContainer1D< Real, Container1D1 > const &	x,
		ITContainer1D< Real, Container1D2 > const &	y
	)

Compute the distance between two vectors.

Compute the Euclidian distance between x and y without overflow.

$d(x,y) = || x - y|| = \sqrt{\sum_{i=1}^n |x_i - y_i|^2}.$

The common range of the vectors is used. The value outside the range are thus interpreted as equal.

Parameters:

[in]	x	first vector
[in]	y	second vector

Returns:: the Euclidian distance between x and y

Definition at line 349 of file STK_LinAlgebra1D.h.

References STK::abs(), STK::ITContainer1D< TYPE, TContainer1D >::first(), STK::ITContainer1D< TYPE, TContainer1D >::last(), STK::max(), and STK::min().

Referenced by STK::LocalVariance::minimalDistance(), and STK::LocalVariance::prim().

{
  // compute the valid range
  const Integer first = max(x.first(), y.first()) , last = min(x.last(), y.last());
  // compute the maximal difference
  Real scale = 0.;
  for (Integer i = first; i<=last; i++)
    scale = max(scale, abs(x[i] - y[i]));
  // Compute the norm
  Real norm2 = 0.;
  if (scale)
  { // comp the norm^2
    for (Integer i = first; i<=last; i++)
    {
      const Real aux = (x[i]-y[i])/scale;
      norm2 += aux * aux;
    }
  }
  // rescale sum
  return (Real(sqrt(double(norm2)))*scale);
}

template<class Container1D1 , class Container1D2 , class Container1D3 >

Real STK::weightedDist	(	ITContainer1D< Real, Container1D1 > const &	x,
		ITContainer1D< Real, Container1D2 > const &	y,
		ITContainer1D< Real, Container1D3 > const &	w
	)

Compute the weighted distance between two vectors.

Compute the weighted Euclidian distance between x and y without overflow.

$d(x,y) = || x - y|| = \sqrt{\sum_{i=1}^n w_i |x_i - y_i|^2}.$

The common range of the vectors is used. The value outside the range are thus interpreted as equal.

Parameters:

[in]	x	first vector
[in]	y	second vector
[in]	w	the weight of the data

Returns:: the weighted Euclidian distance between x and y

Definition at line 388 of file STK_LinAlgebra1D.h.

References STK::abs(), STK::ITContainer1D< TYPE, TContainer1D >::first(), STK::Range::isIn(), STK::ITContainer1D< TYPE, TContainer1D >::last(), STK::max(), STK::min(), and STK::ITContainer1D< TYPE, TContainer1D >::range().

{
  // compute the valid range
  const Integer first = max(x.first(), y.first()) , last = min(x.last(), y.last());
#ifdef STK_DEBUG
  if (!Range(first,last).isIn(w.range()))
    throw std::runtime_error("In weightedDist(x, w) "
                             "Range(first,last) not include in w.range()");
#endif
  // compute the maximal difference
  Real scale = 0., norm2= 0.;
  for (Integer i = first; i<=last; i++)
    scale = max(scale, abs(w[i]*(x[i] - y[i])));
  // Compute the norm
  if (scale)
  { // comp the norm^2
    for (Integer i = first; i<=last; i++)
    {
      const Real aux = (x[i]-y[i])/scale;
      norm2 += w[i]*aux * aux;
    }
  }
  // rescale sum
  return (Real(sqrt(double(norm2)))*scale);
}

Real STK::trace ( MatrixSquare const & A )

trace of a square matrix

Compute the trace of the matrix A.

Parameters:

[in] A the Matrix

Returns:: the sum of the diagonal elements

Definition at line 49 of file STK_LinAlgebra2D.cpp.

References STK::IContainer2D::firstCol(), STK::IContainer2D::lastCol(), and STK::sum().

{
  Real sum = 0.0;
  const Integer first = A.firstCol(), last = A.lastCol();
  for (Integer k = first; k<=last; k++)
    sum += A(k, k);
  return sum;
}

MatrixSquare * STK::mult	(	MatrixSquare const &	A,
		MatrixSquare const &	B
	)

MatrixSquare multiplication.

Perform the matrix product of the square matrix A with the square matrix B. The matrices A and B must have the the correct range.

The result is created on the stack.

Parameters:

[in]	A	Left operand
[in]	B	Right operand

Returns:: a pointer on the result as a MatrixSquare

Definition at line 70 of file STK_LinAlgebra2D.cpp.

References STK::IContainer2D::firstCol(), STK::IContainer2D::lastCol(), STK::MatrixSquare::range(), and STK::sum().

Referenced by STK::LocalVariance::computeAxis(), STK::MultidimRegression::prediction(), STK::BSplineRegression::prediction(), STK::ILinearReduct::projection(), STK::MultidimRegression::regression(), STK::BSplineRegression::regression(), STK::LinearAAModel::simul(), STK::MultidimRegression::wregression(), and STK::BSplineRegression::wregression().

{
#ifdef STK_DEBUG
  if (A.range() != B.range())
    throw std::runtime_error("In MatrixSquare* mult(A, B) "
                             "A.range() != B.range()");
#endif
  // create result
  MatrixSquare* res = new MatrixSquare(A.range());
  // indexes
  const Integer first = A.firstCol(), last = A.lastCol();
  for (Integer i=first; i<=last; i++)
  {
    for (Integer j=first; j<=last; j++)
    {
      Real sum = 0.0;
      for (Integer k = first; k<=last; k++)
        sum += A(i, k) * B(k, j);
      (*res)(i, j) = sum;
    }
  }
  return res;
}

Matrix * STK::mult	(	Matrix const &	A,
		Matrix const &	B
	)

Matrix multiplication.

Perform the matrix product of the matrix A with the matrix B. The matrix A and B must have the the correct range.

The result is created on the stack.

Parameters:

[in]	A	Left operand
[in]	B	Right operand

Returns:: a pointer on the result as a Matrix

Definition at line 186 of file STK_LinAlgebra2D.cpp.

References STK::IContainer2D::firstCol(), STK::IContainer2D::firstRow(), STK::IContainer2D::lastCol(), STK::IContainer2D::lastRow(), STK::IContainer2D::rangeHo(), STK::IContainer2D::rangeVe(), and STK::sum().

{
#ifdef STK_DEBUG
  if (A.rangeHo() != B.rangeVe())
    throw std::runtime_error("In Matrix* mult(A, B) "
                             "A.rangeHo() != B.rangeVe()");
#endif
  // create result
  Matrix* res = new Matrix(A.rangeVe(), B.rangeHo());
  // indexes
  const Integer first_row = A.firstRow(), last_row = A.lastRow();
  const Integer first_col = B.firstCol(), last_col = B.lastCol();
  const Integer first = A.firstCol(), last = A.lastCol();
  //
  for (Integer i=first_row; i<=last_row; i++)
  {
    for (Integer j=first_col; j<=last_col; j++)
    {
      Real sum = 0.0;
      for (Integer k = first; k<=last; k++)
        sum += A(i, k) * B(k, j);
      (*res)(i, j) = sum;
    }
  }
  return res;
}

Matrix * STK::mult	(	MatrixSquare const &	A,
		Matrix const &	B
	)

Matrix multiplication.

Perform the matrix product of the square matrix A with the matrix B. The matrix A and B must have the correct range.

The result is created on the stack.

Parameters:

[in]	A	Left operand
[in]	B	Right operand

Returns:: a pointer on the result as a Matrix

Definition at line 341 of file STK_LinAlgebra2D.cpp.

References STK::IContainer2D::firstCol(), STK::IContainer2D::firstRow(), STK::IContainer2D::lastCol(), STK::IContainer2D::lastRow(), STK::MatrixSquare::range(), STK::IContainer2D::rangeHo(), STK::IContainer2D::rangeVe(), and STK::sum().

{
#ifdef STK_DEBUG
  if (A.range() != B.rangeVe())
    throw std::runtime_error("In Matrix* mult(A, B) "
                             "A.range() != B.rangeHo()");
#endif
  // create result
  Matrix* res = new Matrix(A.range(), B.rangeHo());
  // indexes
  const Integer first_row = A.firstRow(), last_row = A.lastRow();
  const Integer first_col = B.firstCol(), last_col = B.lastCol();
  const Integer first = A.firstCol(), last = A.lastCol();
  for (Integer i=first_row; i<=last_row; i++)
  {
    for (Integer j=first_col; j<=last_col; j++)
    {
      Real sum = 0.0;
      for (Integer k = first; k<=last; k++)
        sum += A(i, k) * B(k, j);
      (*res)(i, j) = sum;
    }
  }
  return res;
}

template<class Container1D >

Vector* STK::mult	(	Matrix const &	A,
		ITContainer1D< Real, Container1D > const &	X
	)

Matrix multiplication by a Vector.

Perform the product of the matrix A by the Vector X The matrix A and X must have the correct range.

The result is created on the stack.

Parameters:

[in]	A	Left operand
[in]	X	Right operand

Returns:: a pointer on the result as a Vector

Definition at line 115 of file STK_LinAlgebra2D.h.

References STK::IContainer2D::firstCol(), STK::IContainer2D::firstRow(), STK::IContainer2D::lastCol(), STK::IContainer2D::lastRow(), STK::ITContainer1D< TYPE, TContainer1D >::range(), STK::IContainer2D::rangeHo(), STK::IContainer2D::rangeVe(), and STK::sum().

{
#ifdef STK_DEBUG
  if (A.rangeHo() != X.range())
    throw std::runtime_error("In Vector* mult(A, X) "
                             "A.rangeHo() != X.range()");
#endif
  // create result
  Vector* res = new Vector(A.rangeVe());
  // indexes
  const Integer first_row = A.firstRow(), last_row = A.lastRow();
  const Integer first_col = A.firstCol(), last_col = A.lastCol();
  for (Integer i=first_row; i<=last_row; i++)
  {
      Real sum = 0.0;
      for (Integer k = first_col; k<=last_col; k++)
        sum += A(i, k) * X[k];
      (*res)[i] = sum;
  }
  return res;
}

template<class TContainerHo , class TContainerVe , class TContainer2D , class TContainer1D >

void STK::mult	(	ITContainer2D< Real, TContainerHo, TContainerVe, TContainer2D > const &	A,
		ITContainer1D< Real, TContainer1D > const &	X,
		ITContainer1D< Real, TContainer1D > &	Y
	)

Matrix multiplication by a Vector.

Perform the product of the matrix A by the Vector X, Y= AX.

Parameters:

[in]	A	Left operand
[in]	X	Right operand
[out]	Y	the result

Returns:: a pointer on the result as a Vector

Definition at line 148 of file STK_LinAlgebra2D.h.

References STK::IContainer2D::firstRow(), STK::IContainer2D::lastRow(), STK::ITContainer1D< TYPE, TContainer1D >::range(), STK::IContainer2D::rangeHo(), STK::IContainer2D::rangeVe(), and STK::IContainer1D::resize().

{
  if (A.rangeHo() != X.range())
    throw std::runtime_error("In void mult(A, X, Y) "
                             "A.rangeHo() != X.range()");
  // create result
  Y.resize(A.rangeVe());
  // indexes
  const Integer first = A.firstRow(), last = A.lastRow();
  for (Integer i=first; i<=last; i++)
  {
    Y[i] = dot<TContainerHo, TContainer1D>(A(i), X);
  }
}

template<class Container1D >

void STK::mult	(	MatrixSquare const &	A,
		ITContainer1D< Real, Container1D > const &	X,
		ITContainer1D< Real, Container1D > &	Y
	)

Square Matrix multiplication by a Vector.

Perform the product of the square matrix A by the Vector X, Y= AX.

Parameters:

[in]	A	Left operand
[in]	X	Right operand
[out]	Y	the result

Returns:: a pointer on the result as a Vector

Definition at line 177 of file STK_LinAlgebra2D.h.

References STK::MatrixSquare::first(), STK::MatrixSquare::last(), STK::ITContainer1D< TYPE, TContainer1D >::range(), STK::MatrixSquare::range(), and STK::IContainer1D::resize().

{
  if (A.range() != X.range())
    throw std::runtime_error("In void mult(A, X, Y) "
                             "A.range() != X.range()");
  // create result
  Y.resize(X.range());
  // indexes
  const Integer first = A.first(), last = A.last();
  for (Integer i=first; i<=last; i++)
  {
      Y[i] = dot<Point, Container1D>(A(i), X);
  }
}

template<class TContainerHo , class TContainerVe , class TContainer2D , class TContainer1D >

void STK::multLefTranspose	(	ITContainer2D< Real, TContainerHo, TContainerVe, TContainer2D > const &	A,
		ITContainer1D< Real, TContainer1D > const &	X,
		ITContainer1D< Real, TContainer1D > &	Y
	)

Matrix multiplication by a Vector.

Perform the product of the transposed matrix A by the Vector X, Y= A'X.

Parameters:

[in]	A	Left operand
[in]	X	Right operand
[out]	Y	the result

Returns:: a pointer on the result as a Vector

Definition at line 206 of file STK_LinAlgebra2D.h.

References STK::IContainer2D::firstRow(), STK::IContainer2D::lastRow(), STK::ITContainer1D< TYPE, TContainer1D >::range(), STK::IContainer2D::rangeVe(), and STK::IContainer1D::resize().

{
  if (A.rangeVe() != X.range())
    throw std::runtime_error("In void multLeftTranspose(A, X, Y) "
                             "A.rangeVe() != X.range()");
  // create result
  Y.resize(X.range());
  // indexes
  const Integer first = A.firstRow(), last = A.lastRow();
  for (Integer i=first; i<=last; i++)
  {
    Y[i] = dot<TContainerVe, TContainer1D>(A[i], X);
  }
}

template<class Container1D >

void STK::multLeftTranspose	(	MatrixSquare const &	A,
		ITContainer1D< Real, Container1D > const &	X,
		ITContainer1D< Real, Container1D > &	Y
	)

transposed square Matrix multiplication by a Vector.

Perform the product of the transposed square matrix A by the Vector X, Y= AX.

Parameters:

[in]	A	Left operand
[in]	X	Right operand
[out]	Y	the result

Returns:: a pointer on the result as a Vector

Definition at line 235 of file STK_LinAlgebra2D.h.

References STK::MatrixSquare::first(), STK::MatrixSquare::last(), STK::ITContainer1D< TYPE, TContainer1D >::range(), STK::MatrixSquare::range(), and STK::IContainer1D::resize().

Referenced by STK::Law::MultivariateNormal< Container1D >::logLikelihood(), STK::Law::MultivariateNormal< Container1D >::lpdf(), STK::MultidimRegression::regression(), and STK::BSplineRegression::regression().

{
  if (A.range() != X.range())
    throw std::runtime_error("In void mult(A, X, Y) "
                             "A.range() != X.range()");
  // create result
  Y.resize(X.range());
  // indexes
  const Integer first = A.first(), last = A.last();
  for (Integer i=first; i<=last; i++)
  {
      Y[i] = dot<Vector, Container1D>(A[i], X);
  }
}

Matrix * STK::multLeftTranspose	(	Matrix const &	A,
		Matrix const &	B
	)

Matrix multiplication.

Perform the matrix product $ A'B $ .

The result is created on the stack.

Parameters:

[in]	A	Left operand
[in]	B	Right operand

Returns:: a pointer on the result as a Matrix

Definition at line 105 of file STK_LinAlgebra2D.cpp.

References STK::IContainer2D::firstCol(), STK::IContainer2D::firstRow(), STK::IContainer2D::lastCol(), STK::IContainer2D::lastRow(), STK::IContainer2D::rangeHo(), STK::IContainer2D::rangeVe(), and STK::sum().

{
#ifdef STK_DEBUG
  if (A.rangeVe() != B.rangeVe())
    throw std::runtime_error("In Matrix* multTranspose(A, B) "
                             "A.rangeVe() != B.rangeVe()");
#endif
  // create result
  Matrix* res = new Matrix(A.rangeHo(), B.rangeHo());
  // indexes
  const Integer first = A.firstRow(), last = A.lastRow();
  const Integer first_row = A.firstCol(), last_row = A.lastCol();
  const Integer first_col = B.firstCol(), last_col = B.lastCol();
  //
  for (Integer i=first_row; i<=last_row; i++)
  {
    for (Integer j=first_col; j<=last_col; j++)
    {
      Real sum = 0.0;
      for (Integer k = first; k<=last; k++)
        sum += A(k, i) * B(k, j);
      (*res)(i, j) = sum;
    }
  }
  return res;
}

template<class Container1D >

Vector* STK::multLeftTranspose	(	Matrix const &	A,
		ITContainer1D< Real, Container1D > const &	X
	)

Matrix by Vector multiplication.

Perform the product $ A'X $ .

The result is created on the stack.

Parameters:

[in]	A	Left operand
[in]	X	Right operand

Returns:: a pointer on the result as a Vector

Definition at line 278 of file STK_LinAlgebra2D.h.

References STK::IContainer2D::firstCol(), STK::IContainer2D::firstRow(), STK::IContainer2D::lastCol(), STK::IContainer2D::lastRow(), STK::ITContainer1D< TYPE, TContainer1D >::range(), STK::IContainer2D::rangeHo(), STK::IContainer2D::rangeVe(), and STK::sum().

{
#ifdef STK_DEBUG
  if (A.rangeVe() != X.range())
    throw std::runtime_error("In Vector* multLeftTranspose(A, X) "
                             "A.rangeHo() != X.range()");
#endif
  // create result
  Vector* res = new Vector(A.rangeHo());
  // indexes
  const Integer first_row = A.firstRow(), last_row = A.lastRow();
  const Integer first_col = A.firstCol(), last_col = A.lastCol();

  for (Integer i=first_col; i<=last_col; i++)
  {
      Real sum = 0.0;
      for (Integer k = first_row; k<=last_row; k++)
        sum += A(k, i) * X[k];
      (*res)[i] = sum;
  }
  return res;
}

Matrix * STK::weightedMultLeftTranspose	(	Matrix const &	A,
		Matrix const &	B,
		Vector const &	weights
	)

Weighted matrix multiplication.

Perform the weighted matrix product $ A'WB $ .

The result is created on the stack.

Parameters:

[in]	A	Left operand
[in]	B	Right operand
[in]	weights	the weights to apply

Returns:: a pointer on the result as a Matrix

Definition at line 145 of file STK_LinAlgebra2D.cpp.

References STK::IContainer2D::firstCol(), STK::IContainer2D::firstRow(), STK::IContainer2D::lastCol(), STK::IContainer2D::lastRow(), STK::IContainer2D::rangeHo(), STK::IContainer2D::rangeVe(), and STK::sum().

Referenced by STK::MultidimRegression::wregression(), and STK::BSplineRegression::wregression().

{
#ifdef STK_DEBUG
  if (A.rangeVe() != B.rangeVe())
    throw std::runtime_error("In Matrix* multTranspose(A, B) "
                             "A.rangeVe() != B.rangeVe()");
#endif
  // create result
  Matrix* res = new Matrix(A.rangeHo(), B.rangeHo());
  // indexes
  const Integer first = A.firstRow(), last = A.lastRow();
  const Integer first_row = A.firstCol(), last_row = A.lastCol();
  const Integer first_col = B.firstCol(), last_col = B.lastCol();
  //
  for (Integer i=first_row; i<=last_row; i++)
  {
    for (Integer j=first_col; j<=last_col; j++)
    {
      Real sum = 0.0;
      for (Integer k = first; k<=last; k++)
        sum += weights[k] * A(k, i) * B(k, j);
      (*res)(i, j) = sum;
    }
  }
  return res;
}

MatrixSquare * STK::multLeftTranspose ( Matrix const & A )

Matrix multiplication.

Perform the matrix product $ A'A $ .

The result is created on the stack.

Parameters:

[in] A the matrix to multiply by itself

Returns:: a pointer on the result as a MatrixSquare

Definition at line 223 of file STK_LinAlgebra2D.cpp.

References STK::IContainer2D::firstCol(), STK::IContainer2D::firstRow(), STK::IContainer2D::lastCol(), STK::IContainer2D::lastRow(), STK::IContainer2D::rangeHo(), and STK::sum().

{
  // create result
  MatrixSquare* res = new MatrixSquare(A.rangeHo());
  // indexes
  const Integer first = A.firstCol(), last = A.lastCol();
  const Integer first_row = A.firstRow(), last_row = A.lastRow();
  //
  for (Integer i=first; i<=last; i++)
  {
    // diagonal
    Real sum = 0.0;
    for (Integer k = first_row; k<=last_row; k++)
      sum += A(k, i) * A(k, i);
    (*res)(i, i) = sum;
    // outside diagonal
    for (Integer j=first; j<i; j++)
    {
      Real sum = 0.0;
      for (Integer k = first_row; k<=last_row; k++)
        sum += A(k, i) * A(k, j);
      (*res)(j, i) = ((*res)(i, j) = sum);
    }
  }
  return res;
}

MatrixSquare * STK::multRightTranspose ( Matrix const & A )

Matrix multiplication.

Perform the matrix product $ AA' $ .

The result is created on the stack.

Parameters:

[in] A the matrix to multiply by itself

Returns:: a pointer on the result as a MatrixSquare

Definition at line 260 of file STK_LinAlgebra2D.cpp.

References STK::IContainer2D::firstCol(), STK::IContainer2D::firstRow(), STK::IContainer2D::lastCol(), STK::IContainer2D::lastRow(), STK::IContainer2D::rangeVe(), and STK::sum().

Referenced by STK::LinearAAModel::simul().

{
  // create result
  MatrixSquare* res = new MatrixSquare(A.rangeVe());
  // get dimensions
  const Integer first = A.firstRow(), last = A.lastRow();
  const Integer first_col = A.firstCol(), last_col = A.lastCol();
  // compute AA'
  for (Integer i=first; i<=last; i++)
  {
    // compute diagonal
    Real sum = 0.0;
    for (Integer k = first_col; k<=last_col; k++)
      sum += A(i, k) * A(i, k);
    (*res)(i, i) = sum;
    // compute outside diagonal
    for (Integer j=first; j<i; j++)
    {
      Real sum = 0.0;
      for (Integer k = first_col; k<=last_col; k++)
        sum += A(i, k) * A(j, k);
      (*res)(j, i) = ((*res)(i, j) = sum);
    }
  }
  return res;
;
}

MatrixSquare * STK::weightedMultLeftTranspose	(	Matrix const &	A,
		Vector const &	weights
	)

Weighted matrix multiplication.

Perform the matrix product $ A'WA $ .

The result is created on the stack.

Parameters:

[in]	A	Left operand
[in]	weights	the weights to apply

Returns:: a pointer on the result as a MatrixSquare

Definition at line 299 of file STK_LinAlgebra2D.cpp.

References STK::IContainer2D::firstCol(), STK::IContainer2D::firstRow(), STK::IContainer2D::lastCol(), STK::IContainer2D::lastRow(), STK::IContainer2D::rangeHo(), and STK::sum().

{
  // create result
  MatrixSquare* res = new MatrixSquare(A.rangeHo());
  // indexes
  const Integer first = A.firstCol(), last = A.lastCol();
  const Integer first_row = A.firstRow(), last_row = A.lastRow();
  //
  for (Integer i=first; i<=last; i++)
  {
    // diagonal
    Real sum = 0.0;
    for (Integer k = first_row; k<=last_row; k++)
      sum += weights[k] * A(k, i) * A(k, i);
    (*res)(i, i) = sum;
    // outside diagonal
    for (Integer j=first; j<i; j++)
    {
      Real sum = 0.0;
      for (Integer k = first_row; k<=last_row; k++)
        sum += weights[k] * A(k, i) * A(k, j);
      (*res)(j, i) = ((*res)(i, j) = sum);
    }
  }
  return res;
}

template<class TContainerHo , class TContainerVe , class TContainer2D >

Real STK::normInf ( ITContainer2D< Real, TContainerHo, TContainerVe, TContainer2D > const & A )

Compute the infinity norm of a 2D container.

Return the maximal absolute value of the 2D container x

$s= \max_{i=1}^n |x_i|$

Parameters:

[in] A matrix to treat.

Returns:: the infinity norm of A

Definition at line 362 of file STK_LinAlgebra2D.h.

References STK::abs(), STK::IContainer2D::firstCol(), STK::IContainer2D::firstRow(), STK::IContainer2D::lastCol(), STK::IContainer2D::lastRow(), and STK::max().

{
  Real scale = 0.0;
  for (Integer j=A.firstRow(); j<=A.lastRow(); j++)
    for (Integer i=A.firstCol(); i<=A.lastCol(); i++)
      scale = STK::max(scale, STK::abs(A(i,j)));
  return (scale);
}

template<class TContainerHo_1 , class TContainerVe_1 , class TContainer2D_1 , class TContainerHo_2 , class TContainerVe_2 , class TContainer2D_2 >

void STK::transpose	(	ITContainer2D< Real, TContainerHo_1, TContainerVe_1, TContainer2D_1 > const &	A,
		ITContainer2D< Real, TContainerHo_2, TContainerVe_2, TContainer2D_2 > &	At
	)

transpose a matrix

Transpose the Matrix A and give the result in At.

Parameters:

[in]	A	the matrix to transpose
[out]	At	the transposed matrix

Definition at line 382 of file STK_LinAlgebra2D.h.

References STK::IContainer2D::firstCol(), STK::IContainer2D::lastCol(), STK::IContainer2D::rangeHo(), STK::IContainer2D::rangeVe(), and STK::IContainer2D::resize().

Referenced by STK::transpose().

{
  // Resize At.
  At.resize(A.rangeHo(), A.rangeVe());

  // copy each column of A in each row of At
  for (Integer j=A.firstCol(); j<=A.lastCol(); j++)
    At(j) = A[j];
}

template<class TContainerHo , class TContainerVe , class TContainer2D >

ITContainer2D< Real, TContainerHo, TContainerVe, TContainer2D >& STK::transpose ( ITContainer2D< Real, TContainerHo, TContainerVe, TContainer2D > & Q )

Transpose a Matrix.

Transpose a matrix and overload it with the result

Parameters:

[in,out] Q the matrix to transpose

Returns:: a reference on the matrix transposed

Definition at line 404 of file STK_LinAlgebra2D.h.

References STK::IContainer2D::firstCol(), STK::IContainer2D::firstRow(), STK::IContainer2D::lastCol(), STK::IContainer2D::lastRow(), STK::Funct::Q, STK::Funct::R, STK::IContainer2D::sizeHo(), STK::IContainer2D::sizeVe(), and STK::transpose().

{
  // check if we have to create an auxiliary cont
  if (Q.sizeVe() != Q.sizeHo())
  {
    // create temporary Matrix
    Matrix R;
    transpose(R, Q);
    Q = R;
    return Q;
  }
  // Square Matrix
  const Integer rbeg = Q.firstRow(), rend = Q.lastRow()
              , cbeg = Q.firstCol(), cend = Q.lastCol();
  // for every rows and cols
  for (Integer  i=rbeg, j=cbeg; j<cend; i++, j++)
  {
    // Create refs to current column and current row
    Vector Qcol(Q, Range(rbeg-cbeg+1 +j, rend), j);
    Point  Qrow(Q, Range(cbeg-rbeg+1 +i, cend), i);
    for ( Integer k1=rbeg-cbeg+1 +j, k2=cbeg-rbeg+1 +i
        ; k1<=rend
        ; k1++, k2++
        )
    {
      Real aux(Qcol[k1]);
      Qcol[k1]  = Qrow[k2];
      Qrow[k2] = aux;
    }
  }
  return Q;
}

template<class TContainerHo1 , class TContainerVe1 , class TContainer2D1 , class TContainerHo2 , class TContainerVe2 , class TContainer2D2 , class TContainerHo3 , class TContainerVe3 , class TContainer2D3 >

void STK::mult	(	ITContainer2D< Real, TContainerHo1, TContainerVe1, TContainer2D1 > &	C,
		const ITContainer2D< Real, TContainerHo2, TContainerVe2, TContainer2D2 > &	A,
		const ITContainer2D< Real, TContainerHo3, TContainerVe3, TContainer2D3 > &	B
	)

Matrix multiplication.

Perform the matricial product of the Matrix A with the Matrix B $ C = AB $ .

The range of C gives the range of the iterations. This method should be specialized for upper and lower triangular Matrix.

Parameters:

[out]	C	Matrix result
[in]	A	Matrix to multiply
[in]	B	Matrix which multiply

Definition at line 455 of file STK_LinAlgebra2D.h.

References STK::dot(), STK::IContainer2D::firstCol(), STK::IContainer2D::firstRow(), STK::IContainer2D::lastCol(), STK::IContainer2D::lastRow(), STK::IContainer2D::rangeHo(), and STK::IContainer2D::rangeVe().

{
#ifdef STK_DEBUG
  if (C.rangeVe() != A.rangeVe())
    throw std::runtime_error("In mult(C, A, B) "
                             "C.rangeVe() != A.rangeVe()");
  if (C.rangeHo() != B.rangeHo())
    throw std::runtime_error("In mult(C, A, B) "
                             "C.rangeHo() != B.rangeHo()");
#endif
  // indexes
  Integer ifirst_col = C.firstCol(), ilast_col = C.lastCol();
  Integer ifirst_row = C.firstRow(), ilast_row = C.lastRow();
  // for all cols
  for (Integer j=ifirst_col; j<=ilast_col; j++)
  {
    // for all lines
    for (Integer i=ifirst_row; i<=ilast_row; i++)
      C(i, j) = dot(A(i), B[j]);
  }
}

Vector* STK::multLeftTranspose	(	Matrix const &	A,
		Vector const &	X
	)

*

Definition at line 412 of file STK_LinAlgebra2D.cpp.

References STK::IContainer2D::firstCol(), STK::IContainer2D::firstRow(), STK::IContainer2D::lastCol(), STK::IContainer2D::lastRow(), STK::ITContainer1D< TYPE, TContainer1D >::range(), STK::IContainer2D::rangeHo(), STK::IContainer2D::rangeVe(), and STK::sum().

{
#ifdef STK_DEBUG
  if (A.rangeVe() != X.range())
    throw std::runtime_error("In Vector* multLeftTranspose(A, X) "
                             "A.rangeHo() != X.range()");
#endif
  // create result
  Vector* res = new Vector(A.rangeHo());
  // indexes
  const Integer first_row = A.firstRow(), last_row = A.lastRow();
  const Integer first_col = A.firstCol(), last_col = A.lastCol();

  for (Integer i=first_col; i<=last_col; i++)
  {
      Real sum = 0.0;
      for (Integer k = first_row; k<=last_row; k++)
        sum += A(k, i) * X[k];
      (*res)[i] = sum;
  }
  return res;
}

Algebra.

Classes

Typedefs

Functions

Detailed Description

Typedef Documentation

Function Documentation