STK::Svd Class Reference
[Subproject STKAlgebra::Algebra]

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 (const Matrix &A=Matrix(), const bool &ref=false, const bool &withU=true, const bool &withV=true)

Svd (const Svd &S)

virtual ~Svd ()

Svd & operator= (const Svd &S)

void clearU ()

void clearV ()

void clear ()

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

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()

const Matrix & getU () const

get U (const)

const MatrixSquare & 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)

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

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

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

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 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 77 of file STK_Svd.h.

Constructor & Destructor Documentation

STK::Svd::Svd	(	const Matrix &	A = `Matrix()`,
		const bool &	ref = `false`,
		const bool &	withU = `true`,
		const bool &	withV = `true`
	)

Default Ctor.

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 100 of file STK_Svd.cpp.

References compSvd(), and init().

00105         : U_(A, ref)
00106         , withU_(withU)
00107         , withV_(withV)
00108         , ref_(ref)
00109 {
00110   // init containers
00111   init();
00112   // compute svd
00113   compSvd();
00114 }

STK::Svd::Svd ( const Svd & S )

Copy Ctor

Parameters:

S

the Svd to copy

Definition at line 117 of file STK_Svd.cpp.

00118         : U_(S.U_, S.ref_)
00119         , V_(S.V_)
00120         , D_(S.D_)
00121         , ncolV_(S.ncolV_)
00122         , ncolD_(S.ncolD_)
00123         , ncolU_(S.ncolU_)
00124         , nrowU_(S.nrowU_)
00125         , withU_(S.withU_)
00126         , withV_(S.withV_)
00127         , ref_(S.ref_)
00128         , norm_(S.norm_)
00129         , rank_(S.rank_)
00130         , error_(S.error_)
00131 { ; }

STK::Svd::~Svd ( ) [virtual]

Dtor.

Definition at line 154 of file STK_Svd.cpp.

00155 { ;}

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 134 of file STK_Svd.cpp.

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

00135 {
00136   U_ = S.U_;
00137   V_ = S.V_;
00138   D_ = S.D_;
00139   ncolV_ = S.ncolV_;
00140   ncolD_ = S.ncolD_;
00141   ncolU_ = S.ncolU_;
00142   nrowU_ = S.nrowU_;
00143   withU_ = S.withU_;
00144   withV_ = S.withV_;
00145   ref_ = false;        // no reference possible with operator =
00146   norm_ = S.norm_;
00147   rank_ = S.rank_;
00148   error_ = S.error_;
00149 
00150   return *this;
00151 }

void STK::Svd::clearU ( )

clear U_.

Definition at line 241 of file STK_Svd.cpp.

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

Referenced by clear().

00242 { U_.clear();
00243   nrowU_ = 0;
00244   ncolU_ = 0;
00245   withU_ = false;
00246 }

void STK::Svd::clearV ( )

clear V_.

Definition at line 250 of file STK_Svd.cpp.

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

Referenced by clear().

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

void STK::Svd::clear ( )

clear U_, V_ and D_.

Definition at line 257 of file STK_Svd.cpp.

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

Referenced by newSvd().

00258 { clearU();
00259   clearV();
00260   D_.clear();
00261 }

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

Compute the svd of the Matrix A and copy the data see the corresponding Ctor. 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 221 of file STK_Svd.cpp.

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

00225 { // clear old results
00226   clear();
00227 
00228   // create U_
00229   U_   = A;           // Copy A in U_
00230   ref_ = false;       // this is not a ref_
00231 
00232   withU_   = withU;   // copy withU_ value
00233   withV_   = withV;   // copy withV_ value
00234 
00235   init();             // initialize (U_) and dimensions
00236   compSvd();          // compute the svd
00237 }

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 265 of file STK_Svd.cpp.

References STK::abs(), STK::ITContainer2D< TYPE, TContainerHo, TContainerVe, TContainer2D >::firstCol(), STK::ITContainer2D< TYPE, TContainerHo, TContainerVe, TContainer2D >::firstRow(), STK::ITContainer1D< TYPE, TContainer1D >::front(), STK::ITContainer2D< TYPE, TContainerHo, TContainerVe, TContainer2D >::getRangeVe(), STK::house(), STK::Inx::incFirst(), STK::ITContainer2D< TYPE, TContainerHo, TContainerVe, TContainer2D >::lastCol(), STK::ITContainer2D< TYPE, TContainerHo, TContainerVe, TContainer2D >::lastRow(), STK::leftHouseholder(), STK::max(), STK::min(), STK::norm(), STK::ITContainer1D< TYPE, TContainer1D >::resize(), STK::rightHouseholder(), STK::ITContainer2D< TYPE, TContainerHo, TContainerVe, TContainer2D >::sizeHo(), and STK::ITContainer2D< TYPE, TContainerHo, TContainerVe, TContainer2D >::sizeVe().

Referenced by compSvd().

00266 {
00267   // norm of the matrix M
00268   Real norm  = 0.0;
00269   // compute the number of iteration
00270   Integer first_iter = M.firstCol();
00271   Integer last_iter  = M.firstCol() + min(M.sizeHo(), M.sizeVe()) -1;
00272   // Diagonal values
00273   D.resize(Inx(first_iter, last_iter));
00274   // Upper diagonal values
00275   F.resize(Inx(first_iter-1, last_iter));
00276   F.front() = 0.0;
00277   // Bidiagonalisation of M
00278   // loop on the cols and rows
00279   Inx rowRange0(M.getRangeVe())
00280     , rowRange1(Inx(M.firstRow()+1, M.lastRow()))
00281     , colRange1(Inx(M.firstCol()+1, M.lastCol()));
00282   for ( Integer iter=first_iter ; iter<=last_iter
00283       ; iter++
00284       , rowRange0.incFirst()
00285       , rowRange1.incFirst()
00286       , colRange1.incFirst()
00287       )
00288   {
00289     // reference on the current col iter
00290     Vector X( M, rowRange0, iter);
00291     // Left Householder
00292     D[iter] = house(X);
00293     // apply Householder to next cols
00294     leftHouseholder(M[colRange1], X);
00295     // check if there is a row
00296     if ((iter < last_iter)||(M.sizeHo()>M.sizeVe()))
00297     {
00298       // ref on the current row iter
00299       Point P(M, colRange1, iter);
00300       // Right Householder
00301       F[iter] = house(P);
00302       // apply Householder to next rows
00303       rightHouseholder(M(rowRange1), P);
00304     }
00305     else
00306       F[iter] = 0.0;
00307     // Estimation of the norm of M
00308     norm = max(abs(D[iter])+abs(F[iter]), norm);
00309   }
00310   return norm;
00311 }

void STK::Svd::rightEliminate	(	Point &	D,
		Vector &	F,
		const Integer &	nrow,
		MatrixSquare &	V,
		const bool &	withV = `true`,
		const Real &	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 420 of file STK_Svd.cpp.

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

Referenced by compSvd(), and diag().

00427 {
00428   // the element to eliminate
00429   Real z = F[nrow];
00430   // if the element is not 0.0
00431   if (abs(z)+tol != tol)
00432   {
00433     // column of the element to eliminate
00434     Integer ncol1 = nrow+1;
00435     // begin the Givens rotations
00436     for (Integer k=nrow, k1=nrow-1; k>=1 ; k--, k1--)
00437     {
00438       // compute and apply Givens rotation to the rows (k, k+1)
00439       Real aux, sinus, cosinus;
00440       Real y = D[k];
00441       D[k]   = (aux     = norm(y,z));
00442       z      = (sinus   = -z/aux) * F[k1];
00443       F[k1] *= (cosinus =  y/aux);
00444       // Update V_
00445       if (withV)
00446         rightGivens(V, ncol1, k, cosinus, sinus);
00447       // if 0.0 we can break now
00448       if (abs(z)+tol == tol) break;
00449     }
00450   }
00451   // the element is now 0
00452   F[nrow] = 0.0;        // is 0.0
00453 }

void STK::Svd::leftEliminate	(	Point &	D,
		Vector &	F,
		const Integer &	nrow,
		Matrix &	U,
		const bool &	withU = `true`,
		const Real &	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 457 of file STK_Svd.cpp.

References STK::abs(), STK::ITContainer1D< TYPE, TContainer1D >::last(), STK::norm(), and STK::rightGivens().

Referenced by diag().

00464 {
00465   //the element to eliminate
00466   Real z = F[nrow];
00467   // if the element is not 0.0
00468   if (abs(z)+tol != tol)
00469   {
00470     // begin the Givens rotations
00471     for (Integer k=nrow+1; k <=D.last(); k++)
00472     {
00473       // compute and apply Givens rotation to the rows (nrow, k)
00474       Real y = D[k];
00475       Real aux, cosinus, sinus;
00476       D[k]  = (aux     = norm(y,z));
00477       z     = (sinus   = -z/aux) * F[k];
00478       F[k] *= (cosinus = y/aux);
00479       // Update U_
00480       if (withU)
00481         rightGivens(U, nrow, k, cosinus, sinus);
00482       if (abs(z)+tol == tol) break;
00483     }
00484   }
00485   F[nrow] = 0.0;
00486 }

bool STK::Svd::diag	(	Point &	D,
		Vector &	F,
		Matrix &	U,
		MatrixSquare &	V,
		const bool &	withU = `true`,
		const bool &	withV = `true`,
		const Real &	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 490 of file STK_Svd.cpp.

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

Referenced by compSvd().

00498 {
00499   // result of the diag process
00500   bool error = false;
00501   // Diagonalization of A : Reduction of la matrice bidiagonale
00502   for (Integer end=D.last(); end>=D.first(); --end)
00503   { // 30 iter max
00504     Integer iter;
00505     for (iter=1; iter<=30; iter++)
00506     { // if  the last element of the surdiagonale is 0.0
00507       // stop the iterations
00508       Integer beg;
00509       if (abs(F[end-1])+tol == tol)  { F[end-1] = 0.0; break;}
00510       // now F[end-1] !=0
00511       // if  D[end] == 0, we can annulate F[end-1]
00512       // with rotations of the cols.
00513       if (abs(D[end])+tol == tol)
00514       {
00515         D[end]   = 0.0;
00516         rightEliminate(D, F, end-1, V, withV, tol);
00517         break; // Last element of the surdiagonal is 0
00518       }
00519       // now D[end] != 0 and F[end-1] != 0
00520       // look for the greatest matrix such that all the elements
00521       // of the diagonal and surdiagonal are not zeros
00522       for (beg = end-1; beg>D.first(); --end)
00523       {
00524         if ((abs(D[beg])+tol == tol)||(abs(F[beg])+tol == tol))
00525         break;
00526       }
00527       // now F[beg-1]!=0
00528       // if D[beg] == 0 and F[beg] != 0,
00529       // we can eliminate the element F[beg]
00530       // with rotations of the rows
00531       if ((abs(D[beg])+tol == tol) && (abs(F[beg])+tol != tol))
00532       {
00533         D[beg] = 0.0;
00534         leftEliminate(D, F, beg, U, withU, tol);
00535       }
00536 
00537       // Si F[beg]==0, on augmente beg
00538       if (abs(F[beg])+tol == tol) { F[beg] = 0.0; beg++;}
00539 
00540       // On peut commencer les rotations QR entre les lignes beg et end
00541       // Shift computation
00542       // easy shift : commented
00543       // Real aux = norm(D[end],F[end-1]);
00544       // Real y   = (D[beg]+aux)*(D[beg]-aux);
00545       // Real z   = D[beg]*F[beg];
00546       // Wilkinson shift : look at the doc
00547       Real dd1 = D[end-1];      // d_1
00548       Real dd2 = D[end];        // d_2
00549       Real ff1 = F[end-2];      // f_1
00550       Real ff2 = F[end-1];      // f_2
00551       // g
00552       Real aux = ( (dd1+dd2)*(dd1-dd2)
00553                  + (ff1-ff2)*(ff1+ff2))/(2*dd1*ff2);
00554       // A - mu
00555       Real d1 = D[beg];
00556       Real f1 = F[beg];
00557       Real y  = (d1+dd2)*(d1-dd2)
00558               + ff2*(dd1/(aux+sign(aux,norm(1.0,aux)))- ff2);
00559       Real z  = d1*f1;
00560       // chase no null element
00561       Integer k, k1;
00562       for (k=beg, k1 = beg+1; k<end; ++k, ++k1)
00563       { // Rotation colonnes (k,k+1)
00564         // Input :  d1 contient D[k],
00565         //          f1 contient F[k],
00566         //          y  contient F[k-1]
00567         // Output : d1 contient F[k],
00568         //          d2 contient D[k+1],
00569         //          y  contient D[k]
00570         Real cosinus, sinus;
00571         Real d2 = D[k1];
00572         F[k-1] = (aux = norm(y,z));                            // F[k-1]
00573        // arbitrary rotation if y = z = 0.0
00574         if (aux)
00575           y  = (cosinus = y/aux) * d1 - (sinus = -z/aux) * f1; // D[k]
00576         else
00577           y  = cosinus * d1 - sinus * f1;                      // D[k]
00578 
00579         z      = -sinus * d2;                                  // z
00580         d1     =  sinus * d1 + cosinus * f1;                   // F[k]
00581         d2    *=  cosinus;                                     // D[k+1]
00582         // Update V
00583         if (withV)
00584           rightGivens(V, k1, k, cosinus, sinus);
00585         // avoid underflow
00586         // Rotation lignes (k,k+1)
00587         // Input :  d1 contient F[k],
00588         //          d2 contient D[k+1],
00589         //          y  contient D[k]
00590         // Output : d1 contient D[k+1],
00591         //          f1 contient F[k+1],
00592         //          y  contient F[k]
00593         f1   = F[k1];
00594         D[k] = (aux = norm(y,z));                              // D[k]
00595         // arbitrary rotation if y = z = 0.0
00596         if (aux)
00597           y  = (cosinus = y/aux) * d1 - (sinus = -z/aux) * d2; // F[k]
00598         else
00599           y   = cosinus * d1 - sinus * d2;                     // F[k]
00600 
00601         z   = -sinus * f1;                                    // z
00602         d1  = sinus *d1 + cosinus * d2;                       // D[k+1]
00603         f1 *= cosinus;                                        // F[k+1]
00604         // Update U
00605         if (withU)
00606           rightGivens(U, k1, k, cosinus, sinus);
00607       } // end of the QR updating iteration
00608       D[end]   = d1;
00609       F[end-1] = y;
00610       F[beg-1] = 0.0;  // F[beg-1] is overwritten, we have to set 0.0
00611     } // iter
00612     // too many iterations
00613     if (iter >= 30) { error = true;}
00614     // positive singular value only
00615     if (D[end]< 0.0)
00616     {
00617       // change sign of the singular value
00618       D[end]= -D[end];
00619       // change sign of the col end of V
00620       if (withV) V[end] = -V[end];
00621     }
00622 
00623     // We have to sort the singular value : we use a basic strategy
00624     Real z = D[end];        // current value
00625     for (Integer i=end+1; i<=D.last(); i++)
00626     { if (D[i]> z)                // if the ith singular value is greater
00627       { D.swapElts(i-1, i);       // swap the cols
00628         if (withU) U.swapCols(i-1, i);
00629         if (withV) V.swapCols(i-1, i);
00630       }
00631       else break;
00632     } // end sort
00633   } // boucle end
00634   return error;
00635 }

Integer STK::Svd::nrowU ( ) const [inline]

Definition at line 218 of file STK_Svd.h.

References STK::ITContainer2D< TYPE, TContainerHo, TContainerVe, TContainer2D >::sizeVe(), and U_.

00218 { return U_.sizeVe();}

Integer STK::Svd::ncolU ( ) const [inline]

Definition at line 220 of file STK_Svd.h.

References STK::ITContainer2D< TYPE, TContainerHo, TContainerVe, TContainer2D >::sizeHo(), and U_.

00220 { return U_.sizeHo();}

Integer STK::Svd::ncolD ( ) const [inline]

Definition at line 222 of file STK_Svd.h.

References ncolD_.

00222 { return ncolD_;}

Integer STK::Svd::nrowV ( ) const [inline]

Definition at line 224 of file STK_Svd.h.

References STK::ITContainer2D< TYPE, TContainerHo, TContainerVe, TContainer2D >::sizeVe(), and V_.

00224 { return V_.sizeVe();}

Integer STK::Svd::ncolV ( ) const [inline]

Definition at line 226 of file STK_Svd.h.

References STK::ITContainer2D< TYPE, TContainerHo, TContainerVe, TContainer2D >::sizeHo(), and V_.

00226 { return V_.sizeHo();}

Real STK::Svd::normSup ( ) const [inline]

Definition at line 228 of file STK_Svd.h.

References norm_.

00228 { return norm_;}

Integer STK::Svd::rank ( ) const [inline]

Definition at line 230 of file STK_Svd.h.

References rank_.

00230 { return rank_;}

bool STK::Svd::error ( ) const [inline]

Definition at line 232 of file STK_Svd.h.

References error_.

Referenced by diag().

00232 { return error_;}

const Matrix& STK::Svd::getU ( ) const [inline]

Definition at line 234 of file STK_Svd.h.

References U_.

00234 { return U_;}

const MatrixSquare& STK::Svd::getV ( ) const [inline]

Definition at line 236 of file STK_Svd.h.

References V_.

00236 { return V_;}

const Point& STK::Svd::getD ( ) const [inline]

Definition at line 238 of file STK_Svd.h.

References D_.

00238 { return D_;}

Matrix& STK::Svd::getU ( ) [inline]

Definition at line 243 of file STK_Svd.h.

References U_.

00243 { return U_;}

MatrixSquare& STK::Svd::getV ( ) [inline]

Definition at line 245 of file STK_Svd.h.

References V_.

00245 { return V_;}

Point& STK::Svd::getD ( ) [inline]

Definition at line 247 of file STK_Svd.h.

References D_.

00247 { return D_;}

void STK::Svd::init ( ) [private]

Definition at line 160 of file STK_Svd.cpp.

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

Referenced by newSvd(), and Svd().

00161 {
00162   // U_ is just a copy of A, translate begin to 1
00163   U_.shift(1,1);   // if U_ is a ref on A, this can generate an error
00164 
00165   // If the container is empty, set default
00166   if (U_.empty())
00167   { ncolV_ = 0;
00168     ncolD_ = 0;
00169     ncolU_ = U_.sizeHo();
00170     nrowU_ = 0;
00171     return;
00172   }
00173 
00174   // The container is not empty, thus we have to compute the dimension
00175   // and eventually to adjust the container (U_)
00176   ncolU_ = U_.sizeHo();         // Number of cols of (U_)
00177   nrowU_ = U_.sizeVe();         // Number of rows of (U_)
00178   ncolV_ = U_.sizeHo();         // Number of cols of V_
00179   ncolD_ = min(nrowU_, ncolV_); // Maximal number of singular value
00180   error_ = false;               // no problems...
00181 }

void STK::Svd::compSvd ( ) [private]

Definition at line 185 of file STK_Svd.cpp.

References bidiag(), STK::Array1D< Real >::clear(), STK::IArray2D< TYPE, TArray2D >::clear(), compU(), compV(), D_, diag(), STK::ITContainer2D< TYPE, TContainerHo, TContainerVe, TContainer2D >::empty(), error_, F_, ncolD_, ncolV_, norm_, nrowU_, rank_, STK::MatrixSquare::resize(), rightEliminate(), U_, V_, withU_, and withV_.

Referenced by newSvd(), and Svd().

00186 {
00187   // if the container is empty, there is nothing to do
00188   if (U_.empty())
00189   { rank_ = 0;
00190     norm_ = 0.0;
00191     return;
00192   }
00193   // Bidiagonalize (U_)
00194   norm_ = bidiag(U_, D_, F_);
00195   // We need to create V_ before rightEliminate
00196   if (withV_) {
00197     V_.resize(ncolV_);
00198     compV();
00199   }
00200   // rightEliminate last element of F_ if any
00201   if (nrowU_ < ncolV_)
00202   { rightEliminate(D_, F_, nrowU_, V_, withV_, norm_);}
00203   // If (U_) is not needed, we can destroy the storage
00204   if (withU_) { compU();}
00205   else        { U_.clear();}
00206   // Diagonalize
00207   error_ = diag(D_, F_, U_, V_, withU_, withV_, norm_);
00208   // Compute the true inf norm
00209   norm_ = D_[1];
00210   // Compute the rank
00211   rank_ = 0;
00212   for (Integer iter=1; iter<=ncolD_; iter++)
00213     if (norm_+D_[iter]!=norm_) { rank_++;}
00214     else break;
00215   // The sub diagonal is now zero
00216   F_.clear();
00217 }

void STK::Svd::compU ( ) [private]

Definition at line 374 of file STK_Svd.cpp.

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

Referenced by compSvd().

00375 {
00376   Integer niter = D_.size();            // Number of iterations
00377   Integer ncol  = min(nrowU_, ncolU_); // number of non zero cols of U_
00378 
00379   // initialization of the remaining cols of U_ to 0.0
00380   // put 0 to unused cols
00381   U_[Inx(ncol+1, ncolU_)] = 0.0;
00382   // Computation of U_
00383   for (Integer iter=niter, iter1=niter+1; iter>=1; iter--, iter1--)
00384   {
00385     // ref of the col iter
00386     Vector X(U_, Inx(iter1,nrowU_), iter);
00387     // ref of the row iter
00388     Point P(U_, Inx(iter,ncolU_), iter);
00389     // Get beta and test
00390     Real beta = P[iter];
00391     if (beta)
00392     {
00393       // update the col iter
00394       P[iter] = 1.0 + beta;
00395       // Updating the cols iter+1 to ncolU_
00396       for (Integer j=iter1; j<=niter; j++)
00397       { // product of U_iter by U_j
00398         Real aux;
00399         Vector Y(U_, Inx(iter1, nrowU_), j); // ref on the col j
00400         // U_j = aux = beta * X'Y
00401         P[j] = (aux = dot( X, Y) *beta);
00402         // U^j += aux * U^iter
00403         Y += X * aux;
00404       }
00405       // compute the vector v
00406       X *= beta;
00407     }
00408     else // U^iter = identity
00409     {
00410       P[iter] = 1.0;
00411       X = 0.0;
00412     }
00413     // update the col iter
00414     U_(Inx(1,iter-1), iter) = 0.0;
00415   }
00416 }

void STK::Svd::compV ( ) [private]

Definition at line 315 of file STK_Svd.cpp.

References beta(), STK::Inx::decFirst(), STK::dot(), STK::ITContainer2D< TYPE, TContainerHo, TContainerVe, TContainer2D >::getRangeHo(), ncolV_, nrowU_, U_, and V_.

Referenced by compSvd().

00316 { // Construction of V_
00317   // Number of right Householder rotations
00318   Integer  niter = (ncolV_>nrowU_) ? (nrowU_) : (ncolV_-1);
00319 
00320   // initialization of the remaining rows and cols of V_ to Identity
00321   for (Integer iter=niter+2; iter<=ncolV_; iter++)
00322   {
00323     Vector W(V_, V_.getRangeHo(), iter);
00324     W       = 0.0;
00325     W[iter] = 1.0;
00326   }
00327 
00328   Inx range1(niter+1, ncolV_), range2(niter+2, ncolV_);
00329   for ( Integer iter0=niter, iter1=niter+1, iter2=niter+2; iter0>=1
00330       ; iter0--, iter1--, iter2--
00331       , range1.decFirst(), range2.decFirst()
00332       )
00333   {
00334     // get beta and test
00335     Real beta = U_(iter0, iter1);
00336     if (beta)
00337     {
00338       // ref on the row iter1 of V_
00339       Point  Vrow1(V_, range1, iter1);
00340       // diagonal element
00341       Vrow1[iter1] = 1.0+beta;
00342       // ref on the col iter1
00343       Vector Vcol1(V_, range2, iter1);
00344       // get the Householder vector
00345       Vcol1 = Point(U_, range2, iter0);
00346       // Apply housholder to next cols
00347       for (Integer j=iter2; j<=ncolV_; j++)
00348       {
00349         Real aux;
00350         // ref on the col j
00351         Vector Vcolj( V_, range2, j);
00352         // update col j
00353         Vrow1[j] = (aux = dot(Vcol1, Vcolj) * beta);
00354         Vcolj   += Vcol1 * aux;
00355       }
00356       // compute the Householder vector
00357       Vcol1 *= beta;
00358     }
00359     else // nothing to do
00360     {
00361       V_(range2, iter1) = 0.0;
00362       V_(iter1, iter1)  = 1.0;
00363       V_(iter1, range2) = 0.0;
00364     }
00365   }
00366   // First col and rows
00367   V_(1,1) =1.0;
00368   V_(Inx(2,ncolV_),1) =0.0;
00369   V_(1,Inx(2,ncolV_)) =0.0;
00370 }

Member Data Documentation

Matrix STK::Svd::U_ [protected]

Definition at line 81 of file STK_Svd.h.

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

MatrixSquare STK::Svd::V_ [protected]

Definition at line 82 of file STK_Svd.h.

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

Point STK::Svd::D_ [protected]

Definition at line 83 of file STK_Svd.h.

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

Integer STK::Svd::ncolV_ [protected]

Definition at line 86 of file STK_Svd.h.

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

Integer STK::Svd::ncolD_ [protected]

Definition at line 87 of file STK_Svd.h.

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

Integer STK::Svd::ncolU_ [protected]

Definition at line 88 of file STK_Svd.h.

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

Integer STK::Svd::nrowU_ [protected]

Definition at line 89 of file STK_Svd.h.

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

bool STK::Svd::withU_ [protected]

Definition at line 92 of file STK_Svd.h.

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

bool STK::Svd::withV_ [protected]

Definition at line 93 of file STK_Svd.h.

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

bool STK::Svd::ref_ [protected]

Definition at line 94 of file STK_Svd.h.

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

Real STK::Svd::norm_ [protected]

Definition at line 97 of file STK_Svd.h.

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

Integer STK::Svd::rank_ [protected]

Definition at line 98 of file STK_Svd.h.

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

bool STK::Svd::error_ [protected]

Definition at line 99 of file STK_Svd.h.

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

Vector STK::Svd::F_ [private]

Definition at line 251 of file STK_Svd.h.

Referenced by compSvd().

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


Public Member Functions
	Svd (const Matrix &A=Matrix(), const bool &ref=false, const bool &withU=true, const bool &withV=true)
	Svd (const Svd &S)
virtual	~Svd ()
Svd &	operator= (const Svd &S)
void	clearU ()
void	clearV ()
void	clear ()
void	newSvd (const Matrix &A=Matrix(), const bool &withU=true, const bool &withV=true)
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()
const Matrix &	getU () const
	get U (const)
const MatrixSquare &	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)
static void	rightEliminate (Point &D, Vector &F, const Integer &nrow, MatrixSquare &V, const bool &withV=true, const Real &tol=Arithmetic< Real >::epsilon())
static void	leftEliminate (Point &D, Vector &F, const Integer &nrow, Matrix &U, const bool &withU=true, const Real &tol=Arithmetic< Real >::epsilon())
static bool	diag (Point &D, Vector &F, Matrix &U, MatrixSquare &V, const bool &withU=true, const bool &withV=true, const Real &tol=Arithmetic< Real >::epsilon())
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.

STK::Svd Class Reference [Subproject STKAlgebra::Algebra]

Public Member Functions

Static Public Member Functions

Protected Attributes

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation

STK::Svd Class Reference
[Subproject STKAlgebra::Algebra]