// Compute the singular value decomposition of A.  
// In VXL this is implemented as an class extending the basic matrix class.

vnl_svd<double> mySVD( A );

// The SVD class holds three matrices U, W, V such that the original 
// matrix A can be written A = U*W*V^T.  
//
// The matrices U and V are orthogonal (i.e. U*U^T = U^T*U = I), and 
// the matrix W stores the singular values in decreasing order along 
// its diagonal.

vnl_matrix<double>      myU = mySVD.U();
vnl_diag_matrix<double> myW = mySVD.W();  
vnl_matrix<double>      myV = mySVD.V();

// Actually vnl_matrix<> would have been fine for A.W() too.
// It's just defined as vnl_diag_matrix<> for efficiency ..

// Extract the largest singular value

double s_max = myW(1,1);  

// Extract the smallest singular value

unsigned rows = A.rows();
unsigned cols = A.cols();
double s_min = (rows > cols) ? myW(cols,cols) : myW(rows,rows);

// Note, we could have simply used the helper methods mySVD.sigma_max() 
// and mySVD.sigma_min() instead! 

// The SVD class has other helper methods too; see the docs
// for more details ....