added rank(Matrix)

miscmaths.cc

@@ -519,6 +519,27 @@ namespace MISCMATHS {
       return pinv;
     }
 
+  int rank(const Matrix& X)
+    {
+      // calculates the rank of matrix X
+      Tracer tr("rank");
+
+      DiagonalMatrix eigenvals;
+      SVD(X,eigenvals);
+
+      double tolerance = Max(X.Nrows(),X.Ncols()) * eigenvals.Maximum() * 1e-16;
+
+      int therank=0;
+
+      for(int i=0; i<eigenvals.Nrows(); i++)
+	if (eigenvals(i+1)>tolerance)
+	  therank++;
+
+      // cout << "tolerance = " << tolerance << "\n" << "eigenvalues = " << eigenvals << "\n" << "rank = " << therank << endl;
+
+      return therank;
+    }
+
 
   ReturnMatrix sqrtaff(const Matrix& mat)
     {

miscmaths.h

@@ -120,6 +120,7 @@ namespace MISCMATHS {
   int diag(DiagonalMatrix& m, const ColumnVector& diagvals);
   ReturnMatrix diag(const Matrix& mat);
   ReturnMatrix pinv(const Matrix& mat);
+  int rank(const Matrix& X);
   ReturnMatrix sqrtaff(const Matrix& mat);
 
   void reshape(Matrix& r, const Matrix& m, int nrows, int ncols);