Changeset fed143 in git

Timestamp:

Sep 20, 2010, 10:11:03 AM (13 years ago)

Author:

Frank Seelisch <seelisch@…>

Branches:

(u'spielwiese', '8d54773d6c9e2f1d2593a28bc68b7eeab54ed529')

Children:

cfd2e6ba5e3f99f46e106336cafaef280361dd61

Parents:

73f78da968d3bf51d9cad24fb4ecf386e3c13276

Message:

rank of a matrix with ground field entries; much faster than linalg.lib's mat_rk

git-svn-id: file:///usr/local/Singular/svn/trunk@13221 2c84dea3-7e68-4137-9b89-c4e89433aadc

Files:

• Singular/iparith.cc (modified) (5 diffs)
• Singular/tok.h (modified) (1 diff)
• kernel/linearAlgebra.cc (modified) (8 diffs)
• kernel/linearAlgebra.h (modified) (1 diff)

Singular/iparith.cc

@@ -412 +412 @@
   { "quotient",    0, QUOTIENT_CMD ,      CMD_2},
   { "random",      0, RANDOM_CMD ,        CMD_23},
+  { "rank",        0, RANK_CMD ,          CMD_12},
   { "read",        0, READ_CMD ,          CMD_12},
   { "reduce",      0, REDUCE_CMD ,        CMD_M},
@@ -2960 +2961 @@
   res->data =(char *)(long)((i > j) ? i : (rand() % (j-i+1)) + i);
 #endif /* buildin_rand */
+  return FALSE;
+}
+static BOOLEAN jjRANK2(leftv res, leftv u, leftv v)
+{
+  matrix m =(matrix)u->Data();
+  int isRowEchelon = (int)(long)v->Data();
+  if (isRowEchelon != 1) isRowEchelon = 0;
+  int rank = luRank(m, isRowEchelon);
+  res->data =(char *)(long)rank;
   return FALSE;
 }
@@ -3977 +3987 @@
 ,{jjQUOT,      QUOTIENT_CMD,   IDEAL_CMD,      MODUL_CMD,  MODUL_CMD, ALLOW_PLURAL |ALLOW_RING}
 ,{jjRANDOM,    RANDOM_CMD,     INT_CMD,        INT_CMD,    INT_CMD, ALLOW_PLURAL |ALLOW_RING}
+,{jjRANK2,     RANK_CMD,       INT_CMD,        MATRIX_CMD, INT_CMD, ALLOW_PLURAL |ALLOW_RING}
 ,{jjREAD2,     READ_CMD,       STRING_CMD,     LINK_CMD,   STRING_CMD, ALLOW_PLURAL |ALLOW_RING}
 ,{jjREDUCE_P,  REDUCE_CMD,     POLY_CMD,       POLY_CMD,   IDEAL_CMD, ALLOW_PLURAL |ALLOW_RING}
@@ -5043 +5054 @@
   }
   //res->data = (char *)0;
+  return FALSE;
+}
+static BOOLEAN jjRANK1(leftv res, leftv v)
+{
+  matrix m =(matrix)v->Data();
+  int rank = luRank(m, 0);
+  res->data =(char *)(long)rank;
   return FALSE;
 }
@@ -5971 +5989 @@
 ,{kQHWeight,    QHWEIGHT_CMD,    INTVEC_CMD,     MODUL_CMD     , ALLOW_PLURAL |ALLOW_RING}
 ,{jjWRONG,      QRING_CMD,       0,              ANY_TYPE      , ALLOW_PLURAL |ALLOW_RING}
+,{jjRANK1,      RANK_CMD,        INT_CMD,        MATRIX_CMD    , ALLOW_PLURAL |ALLOW_RING}
 ,{jjREAD,       READ_CMD,        STRING_CMD,     LINK_CMD      , ALLOW_PLURAL |ALLOW_RING}
 ,{jjREGULARITY, REGULARITY_CMD,  INT_CMD,        LIST_CMD      , NO_PLURAL |ALLOW_RING}

Singular/tok.h

@@ -148 +148 @@
   QRING_CMD,
   RANDOM_CMD,
+  RANK_CMD,
   READ_CMD,
   REPART_CMD,

kernel/linearAlgebra.cc

@@ -148 +148 @@
         {
           n = nDiv(pGetCoeff(p), pivotElement);
+          nNormalize(n);
 
           /* filling lMat;
@@ -157 +158 @@
           n = nNeg(n);
           for (int cGauss = r + cOffset + 1; cGauss <= cc; cGauss++)
+          {
             MATELEM(uMat, rGauss, cGauss)
               = pAdd(MATELEM(uMat, rGauss, cGauss),
                      ppMult_nn(MATELEM(uMat, r, cGauss), n));
+            pNormalize(MATELEM(uMat, rGauss, cGauss));
+          }
 
           nDelete(&n); /* clean up n */
@@ -182 +186 @@
 bool luInverse(const matrix aMat, matrix &iMat)
 { /* aMat is guaranteed to be an (n x n)-matrix */
-  int d = aMat->rows();
-
   matrix pMat;
   matrix lMat;
@@ -196 +198 @@
 
   return result;
+}
+
+/* Assumes that aMat is already in row echelon form */
+int rankFromRowEchelonForm(const matrix aMat)
+{
+  int rank = 0;
+  int rr = aMat->rows(); int cc = aMat->cols();
+  int r = 1; int c = 1;
+  while ((r <= rr) && (c <= cc))
+  {
+    if (MATELEM(aMat, r, c) == NULL) c++;
+    else { rank++; r++; }
+  }
+  return rank;
+}
+
+int luRank(const matrix aMat, const bool isRowEchelon)
+{
+  if (isRowEchelon) return rankFromRowEchelonForm(aMat);
+  else
+  { /* compute the LU-decomposition and read off the rank from
+       the upper triangular matrix of that decomposition */
+    matrix pMat;
+    matrix lMat;
+    matrix uMat;
+    luDecomp(aMat, pMat, lMat, uMat);
+    int result = rankFromRowEchelonForm(uMat);
+  
+    /* clean-up */
+    idDelete((ideal*)&pMat);
+    idDelete((ideal*)&lMat);
+    idDelete((ideal*)&uMat);
+  
+    return result;
+  }
 }
 
@@ -238 +275 @@
         p = pNeg(p);
         p = pMult(p, pCopy(MATELEM(iMat, r, r)));
+        pNormalize(p);
         MATELEM(iMat, r, c) = p;
       }
@@ -285 +323 @@
         p = pNeg(p);
         p = pMult(p, pCopy(MATELEM(iMat, c, c)));
+        pNormalize(p);
         MATELEM(iMat, r, c) = p;
       }
@@ -349 +388 @@
       p = pAdd(p, ppMult_qq(MATELEM(lMat, r, c), MATELEM(yVec, c, 1) ));
     MATELEM(yVec, r, 1) = pNeg(p);
+    pNormalize(MATELEM(yVec, r, 1));
   }
   
@@ -383 +423 @@
       poly q = pNSet(nInvers(pGetCoeff(MATELEM(uMat, r, nonZeroC))));
       MATELEM(xVec, nonZeroC, 1) = pMult(pNeg(p), q);
+      pNormalize(MATELEM(xVec, nonZeroC, 1));
       lastNonZeroC = nonZeroC;
     }

kernel/linearAlgebra.h

@@ -198 +198 @@
        matrix &iMat       /**< [out] inverse of A if invertible   */
                           );
-                          
+
+/**
+ * Computes the rank of a given (m x n)-matrix.
+ *
+ * The matrix may already be given in row echelon form, e.g., when
+ * the user has before called luDecomp and passes the upper triangular
+ * matrix U to luRank. In this case, the second argument can be set to
+ * true in order to pass this piece of information to the method.
+ * Otherwise, luDecomp will be called first to compute the matrix U.
+ * The rank is then read off the matrix U.
+ *
+ * @return the rank as an int
+ * @sa luDecomp(const matrix aMat, matrix &pMat, matrix &lMat, matrix &uMat)
+ **/
+int luRank(
+       const matrix aMat,      /**< [in]  input matrix */
+       const bool isRowEchelon /**< [in]  if true then aMat is assumed to be
+                                          already in row echelon form */
+          );
+
 /**
  * Solves the linear system A*x = b, where A is an (n x n)-matrix