STK::Svd Class Reference

The class Svd compute the Singular Value Decomposition of a Matrix with the Golub-Reinsch Algorithm. More...

#include <STK_Svd.h>

List of all members.

Public Member Functions
	Svd (Matrix const &A=Matrix(), bool ref=false, bool withU=true, bool withV=true)
	Default constructor.
	Svd (const Svd &S)
	Copy Constructor.
virtual	~Svd ()
	destructor.
Svd &	operator= (const Svd &S)
	Operator = : overwrite the Svd with S.
void	clearU ()
	clear U_.
void	clearV ()
	clear V_.
void	clear ()
	clear U_, V_ and D_.
void	newSvd (Matrix const &A=Matrix(), bool withU=true, bool withV=true)
	Compute the svd of the Matrix A and copy the data see the corresponding constructor Take care that if U_ was previously a reference, it cannot be modified.
Integer	nrowU () const
	Number of rows of U_.
Integer	ncolU () const
	Number of cols of U_.
Integer	ncolD () const
	Number of rows of D_.
Integer	nrowV () const
	Number of rows of V_.
Integer	ncolV () const
	Number of cols of V_.
Real	normSup () const
	Norm of the matrix.
Integer	rank () const
	rank of the matrix
bool	error () const
	if an error occur during svd()
Matrix const &	getU () const
	get U (const)
MatrixSquare const &	getV () const
	get V (const)
const Point &	getD () const
	get D (const)
Matrix &	getU ()
	get U
MatrixSquare &	getV ()
	get V
Point &	getD ()
	get D
Static Public Member Functions
static Real	bidiag (const Matrix &M, Point &D, Vector &F)
	Computing the bidiagonalisation of M.
static void	rightEliminate (Point &D, Vector &F, Integer const &nrow, MatrixSquare &V, bool withV=true, Real const &tol=Arithmetic< Real >::epsilon())
	right eliminate the element on the subdiagonal of the row nrow
static void	leftEliminate (Point &D, Vector &F, Integer const &nrow, Matrix &U, bool withU=true, Real const &tol=Arithmetic< Real >::epsilon())
	left eliminate the element on the subdiagonal of the row nrow
static bool	diag (Point &D, Vector &F, Matrix &U, MatrixSquare &V, bool withU=true, bool withV=true, Real const &tol=Arithmetic< Real >::epsilon())
	Computing the diagonalisation of a bidiagnal matrix.
Protected Attributes
Matrix	U_
	U_ matrix.
MatrixSquare	V_
	V_ square matrix.
Point	D_
	Array of the singular values.
Integer	ncolV_
	Number of cols (and rows) of V.
Integer	ncolD_
	Number of cols of D_.
Integer	ncolU_
	Number of cols of U.
Integer	nrowU_
	Number of rows of U_.
bool	withU_
	Compute U_ ?
bool	withV_
	Compute V_ ?
bool	ref_
	Is this structure just a pointer on U_ ?
Real	norm_
	norm of the matrix (largest singular value)
Integer	rank_
	rank of the matrix
bool	error_
	Everything OK during computation ?
Private Member Functions
void	init ()
	Initialize the containers.
void	compSvd ()
	Svd main steps.
void	compU ()
	Compute U (if withU_ is true)
void	compV ()
	Compute V (if withV_ is true)
Private Attributes
Vector	F_
	Values of the Surdiagonal.

---

Detailed Description

The class Svd compute the Singular Value Decomposition of a Matrix with the Golub-Reinsch Algorithm.

The method take as:

• input: A matrix A(nrow,ncol)
• output:
  1. U Matrix (nrow,ncolU).
  2. D Vector (ncol)
  3. V Matrix (ncol,ncol). and perform the decomposition:
• A = UDV' (transpose V). U can have more cols than A, and it is possible to ompute some (all) vectors of Ker(A).

Definition at line 63 of file STK_Svd.h.

---

Constructor & Destructor Documentation

home iovleff Developpement workspace stkpp projects Algebra src STK_Svd cpp STK::Svd::Svd	(	Matrix const &	A = Matrix(),
		bool	ref = false,
		bool	withU = true,
		bool	withV = true
	)

Default constructor.

Parameters:

A	the matrix to decompose.
ref	if true, U_ is a reference of A.
withU	if true, we save the left housolder transforms in U_.
withV	if true, we save the right housolder transforms in V_.

Definition at line 81 of file STK_Svd.cpp.

References compSvd(), and init().

        : U_(A, ref)
        , withU_(withU)
        , withV_(withV)
        , ref_(ref)
{
  // init containers
  init();
  // compute svd
  compSvd();
}

STK::Svd::Svd

(

const Svd &

S

)

Copy Constructor.

Parameters:

S	the Svd to copy

Definition at line 98 of file STK_Svd.cpp.

        : U_(S.U_, S.ref_)
        , V_(S.V_)
        , D_(S.D_)
        , ncolV_(S.ncolV_)
        , ncolD_(S.ncolD_)
        , ncolU_(S.ncolU_)
        , nrowU_(S.nrowU_)
        , withU_(S.withU_)
        , withV_(S.withV_)
        , ref_(S.ref_)
        , norm_(S.norm_)
        , rank_(S.rank_)
        , error_(S.error_)
{ ; }

STK::Svd::~Svd

(

)

[virtual]

destructor.

Definition at line 135 of file STK_Svd.cpp.

{ ;}

---

Member Function Documentation

Svd & STK::Svd::operator=

(

const Svd &

S

)

Operator = : overwrite the Svd with S.

Parameters:

S	the Svd to copy

Definition at line 115 of file STK_Svd.cpp.

References D_, error_, ncolD_, ncolU_, ncolV_, norm_, nrowU_, rank_, ref_, U_, V_, withU_, and withV_.

{
  U_ = S.U_;
  V_ = S.V_;
  D_ = S.D_;
  ncolV_ = S.ncolV_;
  ncolD_ = S.ncolD_;
  ncolU_ = S.ncolU_;
  nrowU_ = S.nrowU_;
  withU_ = S.withU_;
  withV_ = S.withV_;
  ref_ = false;        // no reference possible with operator =
  norm_ = S.norm_;
  rank_ = S.rank_;
  error_ = S.error_;

  return *this;
}

void STK::Svd::clearU

(

)

clear U_.

Definition at line 221 of file STK_Svd.cpp.

References STK::IArray2D< TYPE, TArray2D >::clear(), ncolU_, nrowU_, U_, and withU_.

Referenced by clear().

{ U_.clear();
  nrowU_ = 0;
  ncolU_ = 0;
  withU_ = false;
}

void STK::Svd::clearV

(

)

clear V_.

Definition at line 230 of file STK_Svd.cpp.

References STK::IArray2D< TYPE, TArray2D >::clear(), ncolV_, V_, and withV_.

Referenced by clear().

{ V_.clear(); withV_ = false;
  ncolV_ = 0;
}

void STK::Svd::clear

(

)

clear U_, V_ and D_.

Definition at line 237 of file STK_Svd.cpp.

References STK::ArrayHo< Real >::clear(), clearU(), clearV(), and D_.

Referenced by newSvd().

{ clearU();
  clearV();
  D_.clear();
}

void STK::Svd::newSvd	(	Matrix const &	A = Matrix(),
		bool	withU = true,
		bool	withV = true
	)

Compute the svd of the Matrix A and copy the data see the corresponding constructor Take care that if U_ was previously a reference, it cannot be modified.

Parameters:

A	is the matrix to decompose.
withU	if true, we save the left housolder transforms in U_.
withV	if true, we save the right housolder transforms in V_.

Definition at line 201 of file STK_Svd.cpp.

References clear(), compSvd(), init(), ref_, U_, withU_, and withV_.

{ // clear old results
  clear();

  // create U_
  U_   = A;           // Copy A in U_
  ref_ = false;       // this is not a ref_

  withU_   = withU;   // copy withU_ value
  withV_   = withV;   // copy withV_ value

  init();             // initialize (U_) and dimensions
  compSvd();          // compute the svd
}

Real STK::Svd::bidiag	(	const Matrix &	M,
		Point &	D,
		Vector &	F
	)		[static]

Computing the bidiagonalisation of M.

The diagonal and the subdiagonal are stored in D and F

Parameters:

M	the matrix to bidiagonalize, the matrix is overwritten with the left and right Householder vectors. The method return the estimate of the inf norm of M.
D	the element of the diagonal
F	the element of the surdiagnal

Definition at line 245 of file STK_Svd.cpp.

References STK::abs(), STK::IContainer2D::firstCol(), STK::IContainer2D::firstRow(), STK::ITContainer1D< TYPE, TContainer1D >::front(), STK::house(), STK::Range::incFirst(), STK::IContainer2D::lastCol(), STK::IContainer2D::lastRow(), STK::leftHouseholder(), STK::max(), STK::min(), STK::norm(), STK::Funct::P, STK::IContainer2D::rangeVe(), STK::IContainer1D::resize(), STK::rightHouseholder(), STK::IContainer2D::sizeHo(), and STK::IContainer2D::sizeVe().

Referenced by compSvd().

{
  // norm of the matrix M
  Real norm  = 0.0;
  // compute the number of iteration
  Integer first_iter = M.firstCol();
  Integer last_iter  = M.firstCol() + min(M.sizeHo(), M.sizeVe()) -1;
  // Diagonal values
  D.resize(Range(first_iter, last_iter));
  // Upper diagonal values
  F.resize(Range(first_iter-1, last_iter));
  F.front() = 0.0;
  // Bidiagonalisation of M
  // loop on the cols and rows
  Range rowRange0(M.rangeVe())
    , rowRange1(Range(M.firstRow()+1, M.lastRow()))
    , colRange1(Range(M.firstCol()+1, M.lastCol()));
  for ( Integer iter=first_iter ; iter<=last_iter
      ; iter++
      , rowRange0.incFirst()
      , rowRange1.incFirst()
      , colRange1.incFirst()
      )
  {
    // reference on the current column iter
    Vector X( M, rowRange0, iter);
    // Left Householder
    D[iter] = house(X);
    // apply Householder to next cols
    leftHouseholder(M[colRange1], X);
    // check if there is a row
    if ((iter < last_iter)||(M.sizeHo()>M.sizeVe()))
    {
      // ref on the current row iter
      Point P(M, colRange1, iter);
      // Right Householder
      F[iter] = house(P);
      // apply Householder to next rows
      rightHouseholder(M(rowRange1), P);
    }
    else
      F[iter] = 0.0;
    // Estimation of the norm of M
    norm = max(abs(D[iter])+abs(F[iter]), norm);
  }
  // return estimated norm
  return norm;
}

void STK::Svd::rightEliminate	(	Point &	D,
		Vector &	F,
		Integer const &	nrow,
		MatrixSquare &	V,
		bool	withV = true,
		Real const &	tol = Arithmetic<Real>::epsilon()
	)		[static]

right eliminate the element on the subdiagonal of the row nrow

Parameters:

D	the diagonal of the matrix
F	the surdiagonal of the matrix
nrow	the number of the row were we want to rightEliminate
V	a right orthogonal Square Matrix
withV	true if we want to update V
tol	the tolerance to use

Definition at line 401 of file STK_Svd.cpp.

References STK::abs(), STK::norm(), and STK::rightGivens().

Referenced by compSvd(), and diag().

{
  // the element to eliminate
  Real z = F[nrow];
  // if the element is not 0.0
  if (abs(z)+tol != tol)
  {
    // column of the element to eliminate
    Integer ncol1 = nrow+1;
    // begin the Givens rotations
    for (Integer k=nrow, k1=nrow-1; k>=1 ; k--, k1--)
    {
      // compute and apply Givens rotation to the rows (k, k+1)
      Real aux, sinus, cosinus;
      Real y = D[k];
      D[k]   = (aux     = norm(y,z));
      z      = (sinus   = -z/aux) * F[k1];
      F[k1] *= (cosinus =  y/aux);
      // Update V_
      if (withV)
        rightGivens(V, ncol1, k, cosinus, sinus);
      // if 0.0 we can break now
      if (abs(z)+tol == tol) break;
    }
  }
  // the element is now 0
  F[nrow] = 0.0;        // is 0.0
}

void STK::Svd::leftEliminate	(	Point &	D,
		Vector &	F,
		Integer const &	nrow,
		Matrix &	U,
		bool	withU = true,
		Real const &	tol = Arithmetic<Real>::epsilon()
	)		[static]

left eliminate the element on the subdiagonal of the row nrow

Parameters:

D	the diagonal of the matrix
F	the surdiagonal of the matrix
nrow	the number of the row were we want to rightEliminate
U	a left orthogonal Matrix
withU	true if we want to update U
tol	the tolerance to use

Definition at line 438 of file STK_Svd.cpp.

References STK::abs(), STK::norm(), and STK::rightGivens().

Referenced by diag().

{
  //the element to eliminate
  Real z = F[nrow];
  // if the element is not 0.0
  if (abs(z)+tol != tol)
  {
    // begin the Givens rotations
    for (Integer k=nrow+1; k <=D.last(); k++)
    {
      // compute and apply Givens rotation to the rows (nrow, k)
      Real y = D[k];
      Real aux, cosinus, sinus;
      D[k]  = (aux     = norm(y,z));
      z     = (sinus   = -z/aux) * F[k];
      F[k] *= (cosinus = y/aux);
      // Update U_
      if (withU)
        rightGivens(U, nrow, k, cosinus, sinus);
      if (abs(z)+tol == tol) break;
    }
  }
  F[nrow] = 0.0;
}

bool STK::Svd::diag	(	Point &	D,
		Vector &	F,
		Matrix &	U,
		MatrixSquare &	V,
		bool	withU = true,
		bool	withV = true,
		Real const &	tol = Arithmetic<Real>::epsilon()
	)		[static]

Computing the diagonalisation of a bidiagnal matrix.

Parameters:

D	the diagoanl of the matrix
F	the subdiagonal of the matrix
U	a left orthogonal Matrix
withU	true if we want to update U
V	a right orthogonal Square Matrix
withV	true if we want to update V
tol	the tolerance to use

Definition at line 471 of file STK_Svd.cpp.

References STK::abs(), d1, d2, error(), leftEliminate(), STK::norm(), rightEliminate(), STK::rightGivens(), STK::sign(), STK::IArray2DBase< TYPE, PTRCOL, TArrayHo, TArrayVe, TArray2D >::swap(), and STK::IArray2DBase< TYPE, PTRCOL, TArrayHo, TArrayVe, TArray2D >::swapCols().

Referenced by compSvd().

{
  // result of the diag process
  bool error = false;
  // Diagonalization of A : Reduction of la matrice bidiagonale
  for (Integer end=D.last(); end>=D.first(); --end)
  { // 30 iter max
    Integer iter;
    for (iter=1; iter<=30; iter++)
    { // if  the last element of the surdiagonale is 0.0
      // stop the iterations
      Integer beg;
      if (abs(F[end-1])+tol == tol)  { F[end-1] = 0.0; break;}
      // now F[end-1] !=0
      // if  D[end] == 0, we can annulate F[end-1]
      // with rotations of the columns.
      if (abs(D[end])+tol == tol)
      {
        D[end]   = 0.0;
        rightEliminate(D, F, end-1, V, withV, tol);
        break; // Last element of the surdiagonal is 0
      }
      // now D[end] != 0 and F[end-1] != 0
      // look for the greatest matrix such that all the elements
      // of the diagonal and surdiagonal are not zeros
      for (beg = end-1; beg>D.first(); --beg)
      {
        if ((abs(D[beg])+tol == tol)||(abs(F[beg])+tol == tol))
        break;
      }
      // now F[beg-1]!=0
      // if D[beg] == 0 and F[beg] != 0,
      // we can eliminate the element F[beg]
      // with rotations of the rows
      if ((abs(D[beg])+tol == tol) && (abs(F[beg])+tol != tol))
      {
        D[beg] = 0.0;
        leftEliminate(D, F, beg, U, withU, tol);
      }

      // Si F[beg]==0, on augmente beg
      if (abs(F[beg])+tol == tol) { F[beg] = 0.0; beg++;}

      // On peut commencer les rotations QR entre les lignes beg et end
      // Shift computation
      // easy shift : commented
      // Real aux = norm(D[end],F[end-1]);
      // Real y   = (D[beg]+aux)*(D[beg]-aux);
      // Real z   = D[beg]*F[beg];
      // Wilkinson shift : look at the doc
      Real dd1 = D[end-1];      // d_1
      Real dd2 = D[end];        // d_2
      Real ff1 = F[end-2];      // f_1
      Real ff2 = F[end-1];      // f_2
      // g
      Real aux = ( (dd1+dd2)*(dd1-dd2)
                 + (ff1-ff2)*(ff1+ff2))/(2*dd1*ff2);
      // A - mu
      Real d1 = D[beg];
      Real f1 = F[beg];
      Real y  = (d1+dd2)*(d1-dd2)
              + ff2*(dd1/(aux+sign(aux,norm(1.0,aux)))- ff2);
      Real z  = d1*f1;
      // chase no null element
      Integer k, k1;
      for (k=beg, k1 = beg+1; k<end; ++k, ++k1)
      { // Rotation colonnes (k,k+1)
        // Input :  d1 contient D[k],
        //          f1 contient F[k],
        //          y  contient F[k-1]
        // Output : d1 contient F[k],
        //          d2 contient D[k+1],
        //          y  contient D[k]
        Real cosinus, sinus;
        Real d2 = D[k1];
        F[k-1] = (aux = norm(y,z));                            // F[k-1]
       // arbitrary rotation if y = z = 0.0
        if (aux)
          y  = (cosinus = y/aux) * d1 - (sinus = -z/aux) * f1; // D[k]
        else
          y  = cosinus * d1 - sinus * f1;                      // D[k]

        z      = -sinus * d2;                                  // z
        d1     =  sinus * d1 + cosinus * f1;                   // F[k]
        d2    *=  cosinus;                                     // D[k+1]
        // Update V
        if (withV)
          rightGivens(V, k1, k, cosinus, sinus);
        // avoid underflow
        // Rotation lignes (k,k+1)
        // Input :  d1 contient F[k],
        //          d2 contient D[k+1],
        //          y  contient D[k]
        // Output : d1 contient D[k+1],
        //          f1 contient F[k+1],
        //          y  contient F[k]
        f1   = F[k1];
        D[k] = (aux = norm(y,z));                              // D[k]
        // arbitrary rotation if y = z = 0.0
        if (aux)
          y  = (cosinus = y/aux) * d1 - (sinus = -z/aux) * d2; // F[k]
        else
          y   = cosinus * d1 - sinus * d2;                     // F[k]

        z   = -sinus * f1;                                    // z
        d1  = sinus *d1 + cosinus * d2;                       // D[k+1]
        f1 *= cosinus;                                        // F[k+1]
        // Update U
        if (withU)
          rightGivens(U, k1, k, cosinus, sinus);
      } // end of the QR updating iteration
      D[end]   = d1;
      F[end-1] = y;
      F[beg-1] = 0.0;  // F[beg-1] is overwritten, we have to set 0.0
    } // iter
    // too many iterations
    if (iter >= 30) { error = true;}
    // positive singular value only
    if (D[end]< 0.0)
    {
      // change sign of the singular value
      D[end]= -D[end];
      // change sign of the column end of V
      if (withV) V[end] = -V[end];
    }

    // We have to sort the singular value : we use a basic strategy
    Real z = D[end];        // current value
    for (Integer i=end+1; i<=D.last(); i++)
    { if (D[i]> z)                // if the ith singular value is greater
      { D.swap(i-1, i);       // swap the cols
        if (withU) U.swapCols(i-1, i);
        if (withV) V.swapCols(i-1, i);
      }
      else break;
    } // end sort
  } // boucle end
  return error;
}

Integer STK::Svd::nrowU

(

)

const [inline]

Number of rows of U_.

Definition at line 200 of file STK_Svd.h.

References STK::IContainer2D::sizeVe(), and U_.

{ return U_.sizeVe();}

Integer STK::Svd::ncolU

(

)

const [inline]

Number of cols of U_.

Definition at line 202 of file STK_Svd.h.

References STK::IContainer2D::sizeHo(), and U_.

{ return U_.sizeHo();}

Integer STK::Svd::ncolD

(

)

const [inline]

Number of rows of D_.

Definition at line 204 of file STK_Svd.h.

References ncolD_.

{ return ncolD_;}

Integer STK::Svd::nrowV

(

)

const [inline]

Number of rows of V_.

Definition at line 206 of file STK_Svd.h.

References STK::IContainer2D::sizeVe(), and V_.

{ return V_.sizeVe();}

Integer STK::Svd::ncolV

(

)

const [inline]

Number of cols of V_.

Definition at line 208 of file STK_Svd.h.

References STK::IContainer2D::sizeHo(), and V_.

{ return V_.sizeHo();}

Real STK::Svd::normSup

(

)

const [inline]

Norm of the matrix.

Definition at line 210 of file STK_Svd.h.

References norm_.

{ return norm_;}

Integer STK::Svd::rank

(

)

const [inline]

rank of the matrix

Definition at line 212 of file STK_Svd.h.

References rank_.

{ return rank_;}

bool STK::Svd::error

(

)

const [inline]

if an error occur during svd()

Definition at line 214 of file STK_Svd.h.

References error_.

Referenced by diag().

{ return error_;}

Matrix const& STK::Svd::getU

(

)

const [inline]

get U (const)

Definition at line 216 of file STK_Svd.h.

References U_.

{ return U_;}

MatrixSquare const& STK::Svd::getV

(

)

const [inline]

get V (const)

Definition at line 218 of file STK_Svd.h.

References V_.

{ return V_;}

const Point& STK::Svd::getD

(

)

const [inline]

get D (const)

Definition at line 220 of file STK_Svd.h.

References D_.

{ return D_;}

Matrix& STK::Svd::getU

(

)

[inline]

get U

Definition at line 223 of file STK_Svd.h.

References U_.

{ return U_;}

MatrixSquare& STK::Svd::getV

(

)

[inline]

get V

Definition at line 225 of file STK_Svd.h.

References V_.

{ return V_;}

Point& STK::Svd::getD

(

)

[inline]

get D

Definition at line 227 of file STK_Svd.h.

References D_.

{ return D_;}

void STK::Svd::init

(

)

[private]

Initialize the containers.

Definition at line 141 of file STK_Svd.cpp.

References STK::IContainer2D::empty(), error_, STK::min(), ncolD_, ncolU_, ncolV_, nrowU_, STK::IArray2D< TYPE, TArray2D >::shift(), STK::IContainer2D::sizeHo(), STK::IContainer2D::sizeVe(), and U_.

Referenced by newSvd(), and Svd().

{
  // U_ is just a copy of A, translate begin to 1
  U_.shift(1,1);   // if U_ is a ref on A, this can generate an error

  // If the container is empty, set default
  if (U_.empty())
  { ncolV_ = 0;
    ncolD_ = 0;
    ncolU_ = U_.sizeHo();
    nrowU_ = 0;
    return;
  }

  // The container is not empty, thus we have to compute the dimension
  // and eventually to adjust the container (U_)
  ncolU_ = U_.sizeHo();         // Number of cols of (U_)
  nrowU_ = U_.sizeVe();         // Number of rows of (U_)
  ncolV_ = U_.sizeHo();         // Number of cols of V_
  ncolD_ = min(nrowU_, ncolV_); // Maximal number of singular value
  error_ = false;               // no problems...
}

void STK::Svd::compSvd

(

)

[private]

Svd main steps.

Definition at line 166 of file STK_Svd.cpp.

References bidiag(), STK::RecursiveArray1D< TYPE, Container1D >::clear(), compU(), compV(), D_, diag(), STK::IContainer2D::empty(), error_, F_, ncolD_, ncolV_, norm_, nrowU_, rank_, STK::MatrixSquare::resize(), rightEliminate(), U_, V_, withU_, and withV_.

Referenced by newSvd(), and Svd().

{
  // if the container is empty, there is nothing to do
  if (U_.empty())
  { rank_ = 0;
    norm_ = 0.0;
    return;
  }
  // Bidiagonalize (U_)
  norm_ = bidiag(U_, D_, F_);
  // We need to create V_ before rightEliminate
  if (withV_) {
    V_.resize(ncolV_);
    compV();
  }
  // rightEliminate last element of F_ if any
  if (nrowU_ < ncolV_)
  { rightEliminate(D_, F_, nrowU_, V_, withV_, norm_);}
  // If (U_) is not needed, we can destroy the storage
  if (withU_) { compU();}
  // Diagonalize
  error_ = diag(D_, F_, U_, V_, withU_, withV_, norm_);
  // Compute the true inf norm
  norm_ = D_[1];
  // Compute the rank
  rank_ = 0;
  for (Integer iter=1; iter<=ncolD_; iter++)
    if (norm_+D_[iter]!=norm_) { rank_++;}
    else break;
  // The sub diagonal is now zero
  F_.clear();
}

void STK::Svd::compU

(

)

[private]

Compute U (if withU_ is true)

Definition at line 355 of file STK_Svd.cpp.

References beta(), D_, STK::dot(), STK::min(), ncolU_, nrowU_, STK::Funct::P, STK::ITContainer1D< TYPE, TContainer1D >::size(), and U_.

Referenced by compSvd().

{
  Integer niter = D_.size();            // Number of iterations
  Integer ncol  = min(nrowU_, ncolU_); // number of non zero cols of U_

  // initialization of the remaining cols of U_ to 0.0
  // put 0 to unused cols
  U_[Range(ncol+1, ncolU_)] = 0.0;
  // Computation of U_
  for (Integer iter=niter, iter1=niter+1; iter>=1; iter--, iter1--)
  {
    // ref of the column iter
    Vector X(U_, Range(iter1,nrowU_), iter);
    // ref of the row iter
    Point P(U_, Range(iter,ncolU_), iter);
    // Get beta and test
    Real beta = P[iter];
    if (beta)
    {
      // update the column iter
      P[iter] = 1.0 + beta;
      // Updating the cols iter+1 to ncolU_
      for (Integer j=iter1; j<=niter; j++)
      { // product of U_iter by U_j
        Real aux;
        Vector Y(U_, Range(iter1, nrowU_), j); // ref on the column j
        // U_j = aux = beta * X'Y
        P[j] = (aux = dot( X, Y) *beta);
        // U^j += aux * U^iter
        Y += X * aux;
      }
      // compute the vector v
      X *= beta;
    }
    else // U^iter = identity
    {
      P[iter] = 1.0;
      X = 0.0;
    }
    // update the column iter
    U_(Range(1,iter-1), iter) = 0.0;
  }
}

void STK::Svd::compV

(

)

[private]

Compute V (if withV_ is true)

Definition at line 296 of file STK_Svd.cpp.

References beta(), STK::Range::decFirst(), STK::dot(), ncolV_, nrowU_, STK::IContainer2D::rangeHo(), U_, and V_.

Referenced by compSvd().

{ // Construction of V_
  // Number of right Householder rotations
  Integer  niter = (ncolV_>nrowU_) ? (nrowU_) : (ncolV_-1);

  // initialization of the remaining rows and cols of V_ to Identity
  for (Integer iter=niter+2; iter<=ncolV_; iter++)
  {
    Vector W(V_, V_.rangeHo(), iter);
    W       = 0.0;
    W[iter] = 1.0;
  }

  Range range1(niter+1, ncolV_), range2(niter+2, ncolV_);
  for ( Integer iter0=niter, iter1=niter+1, iter2=niter+2; iter0>=1
      ; iter0--, iter1--, iter2--
      , range1.decFirst(), range2.decFirst()
      )
  {
    // get beta and test
    Real beta = U_(iter0, iter1);
    if (beta)
    {
      // ref on the row iter1 of V_
      Point  Vrow1(V_, range1, iter1);
      // diagonal element
      Vrow1[iter1] = 1.0+beta;
      // ref on the column iter1
      Vector Vcol1(V_, range2, iter1);
      // get the Householder vector
      Vcol1 = Point(U_, range2, iter0);
      // Apply housholder to next cols
      for (Integer j=iter2; j<=ncolV_; j++)
      {
        Real aux;
        // ref on the column j
        Vector Vcolj( V_, range2, j);
        // update column j
        Vrow1[j] = (aux = dot(Vcol1, Vcolj) * beta);
        Vcolj   += Vcol1 * aux;
      }
      // compute the Householder vector
      Vcol1 *= beta;
    }
    else // nothing to do
    {
      V_(range2, iter1) = 0.0;
      V_(iter1, iter1)  = 1.0;
      V_(iter1, range2) = 0.0;
    }
  }
  // First column and rows
  V_(1,1) =1.0;
  V_(Range(2,ncolV_),1) =0.0;
  V_(1,Range(2,ncolV_)) =0.0;
}

---

Member Data Documentation

Matrix STK::Svd::U_ [protected]

U_ matrix.

Definition at line 67 of file STK_Svd.h.

Referenced by clearU(), compSvd(), compU(), compV(), getU(), init(), ncolU(), newSvd(), nrowU(), and operator=().

MatrixSquare STK::Svd::V_ [protected]

V_ square matrix.

Definition at line 68 of file STK_Svd.h.

Referenced by clearV(), compSvd(), compV(), getV(), ncolV(), nrowV(), and operator=().

Point STK::Svd::D_ [protected]

Array of the singular values.

Definition at line 69 of file STK_Svd.h.

Referenced by clear(), compSvd(), compU(), getD(), and operator=().

Integer STK::Svd::ncolV_ [protected]

Number of cols (and rows) of V.

Definition at line 72 of file STK_Svd.h.

Referenced by clearV(), compSvd(), compV(), init(), and operator=().

Integer STK::Svd::ncolD_ [protected]

Number of cols of D_.

Definition at line 73 of file STK_Svd.h.

Referenced by compSvd(), init(), ncolD(), and operator=().

Integer STK::Svd::ncolU_ [protected]

Number of cols of U.

Definition at line 74 of file STK_Svd.h.

Referenced by clearU(), compU(), init(), and operator=().

Integer STK::Svd::nrowU_ [protected]

Number of rows of U_.

Definition at line 75 of file STK_Svd.h.

Referenced by clearU(), compSvd(), compU(), compV(), init(), and operator=().

bool STK::Svd::withU_ [protected]

Compute U_ ?

Definition at line 78 of file STK_Svd.h.

Referenced by clearU(), compSvd(), newSvd(), and operator=().

bool STK::Svd::withV_ [protected]

Compute V_ ?

Definition at line 79 of file STK_Svd.h.

Referenced by clearV(), compSvd(), newSvd(), and operator=().

bool STK::Svd::ref_ [protected]

Is this structure just a pointer on U_ ?

Definition at line 80 of file STK_Svd.h.

Referenced by newSvd(), and operator=().

Real STK::Svd::norm_ [protected]

norm of the matrix (largest singular value)

Definition at line 83 of file STK_Svd.h.

Referenced by compSvd(), normSup(), and operator=().

Integer STK::Svd::rank_ [protected]

rank of the matrix

Definition at line 84 of file STK_Svd.h.

Referenced by compSvd(), operator=(), and rank().

bool STK::Svd::error_ [protected]

Everything OK during computation ?

Definition at line 85 of file STK_Svd.h.

Referenced by compSvd(), error(), init(), and operator=().

Vector STK::Svd::F_ [private]

Values of the Surdiagonal.

Definition at line 231 of file STK_Svd.h.

Referenced by compSvd().

---

The documentation for this class was generated from the following files:

• STK_Svd.h
• STK_Svd.cpp