Changeset 6ed8c4 in git

Timestamp:

Jul 19, 2011, 4:45:41 PM (12 years ago)

Author:

mlee <martinlee84@…>

Branches:

(u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'a800fe4b3e9d37a38c5a10cc0ae9dfa0c15a4ee6')

Children:

1f637e0cf55bac728426f222f498b81de378e159

Parents:

bb5c2849ef4e9115c1f69238f0aaf689e2be7e2d

git-author:

mlee <martinlee84@web.de>2011-07-19 16:45:41+02:00

git-committer:

Mohamed Barakat <mohamed.barakat@rwth-aachen.de>2011-11-09 12:52:40+01:00

Message:

added const ring to some functions in linearAlgebra.* (TODO complete this)
switched on compatibility layers in numbers.h and polys.h
added nGreater macro to numbers.h
commented out lists in ring.h
moved lists dependend functions from linearAlgebra.* to linearAlgebra.cut

Files:

• Singular/linearAlgebra.cut (added)
• kernel/linearAlgebra.cc (modified) (28 diffs)
• kernel/linearAlgebra.h (modified) (10 diffs)
• libpolys/coeffs/numbers.h (modified) (3 diffs)
• libpolys/polys/monomials/ring.h (modified) (2 diffs)
• libpolys/polys/polys.h (modified) (2 diffs)

kernel/linearAlgebra.cc

@@ -17 +17 @@
 
 // include header files
-#include <kernel/mod2.h>
+#include "mod2.h"
+
+#include <coeffs/coeffs.h>
+#include <coeffs/numbers.h>
+
+#include <coeffs/mpr_complex.h>
+#include <polys/monomials/p_polys.h>
+#include <polys/simpleideals.h>
+#include <polys/matpol.h>
+
+// #include <polys/polys.h>
+
 #include <kernel/structs.h>
-#include <polys/polys.h>
 #include <kernel/ideals.h>
-#include <coeffs/numbers.h>
-#include <polys/matpol.h>
-#include <coeffs/mpr_complex.h>
 #include <kernel/linearAlgebra.h>
 
@@ -39 +46 @@
  * the negative of nSize will be returned.
  **/
-int pivotScore(number n)
-{
-  int s = nSize(n);
-  if (rField_is_long_C(currRing) ||
-      rField_is_long_R(currRing) ||
-      rField_is_R(currRing))
+int pivotScore(number n, const ring r)
+{
+  int s = n_Size(n, r->cf);
+  if (rField_is_long_C(r) ||
+      rField_is_long_R(r) ||
+      rField_is_R(r))
     return -s;
   else
@@ -58 +65 @@
  **/
 bool pivot(const matrix aMat, const int r1, const int r2, const int c1,
-           const int c2, int* bestR, int* bestC)
+           const int c2, int* bestR, int* bestC, const ring R)
 {
   int bestScore; int score; bool foundBestScore = false; poly matEntry;
@@ -69 +76 @@
       if (matEntry != NULL)
       {
-        score = pivotScore(pGetCoeff(matEntry));
+        score = pivotScore(pGetCoeff(matEntry), R);
         if ((!foundBestScore) || (score < bestScore))
         {
@@ -84 +91 @@
 }
 
-bool unitMatrix(const int n, matrix &unitMat)
+bool unitMatrix(const int n, matrix &unitMat, const ring R)
 {
   if (n < 1) return false;
   unitMat = mpNew(n, n);
-  for (int r = 1; r <= n; r++) MATELEM(unitMat, r, r) = pOne();
+  for (int r = 1; r <= n; r++) MATELEM(unitMat, r, r) = p_One(R);
   return true;
 }
 
-void luDecomp(const matrix aMat, matrix &pMat, matrix &lMat, matrix &uMat)
+void luDecomp(const matrix aMat, matrix &pMat, matrix &lMat, matrix &uMat,
+              const ring R)
 {
   int rr = aMat->rows();
@@ -101 +109 @@
   for (int r = 1; r <= rr; r++)
     for (int c = 1; c <= cc; c++)
-      MATELEM(uMat, r, c) = pCopy(aMat->m[c - 1 + (r - 1) * cc]);
+      MATELEM(uMat, r, c) = p_Copy(aMat->m[c - 1 + (r - 1) * cc], R);
 
   /* we use an int array to store all row permutations;
@@ -109 +117 @@
 
   /* fill lMat with the (rr x rr) unit matrix */
-  unitMatrix(rr, lMat);
+  unitMatrix(rr, lMat,R);
 
   int bestR; int bestC; int intSwap; poly pSwap; int cOffset = 0;
@@ -116 +124 @@
     if (r > cc) break;
     while ((r + cOffset <= cc) &&
-           (!pivot(uMat, r, rr, r + cOffset, r + cOffset, &bestR, &bestC)))
+           (!pivot(uMat, r, rr, r + cOffset, r + cOffset, &bestR, &bestC, R)))
       cOffset++;
     if (r + cOffset <= cc)
@@ -155 +163 @@
         if (p != NULL)
         {
-          n = nDiv(pGetCoeff(p), pivotElement);
-          nNormalize(n);
+          n = n_Div(pGetCoeff(p), pivotElement, R->cf);
+          n_Normalize(n,R->cf);
 
           /* filling lMat;
              old entry was zero, so no need to call pDelete(.) */
-          MATELEM(lMat, rGauss, r) = pNSet(nCopy(n));
+          MATELEM(lMat, rGauss, r) = p_NSet(n_Copy(n,R->cf),R);
 
           /* adjusting uMat: */
-          MATELEM(uMat, rGauss, r + cOffset) = NULL; pDelete(&p);
-          n = nNeg(n);
+          MATELEM(uMat, rGauss, r + cOffset) = NULL; p_Delete(&p,R);
+          n = n_Neg(n,R->cf);
           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));
+              = p_Add_q(MATELEM(uMat, rGauss, cGauss),
+                     pp_Mult_nn(MATELEM(uMat, r, cGauss), n, R), R);
+            p_Normalize(MATELEM(uMat, rGauss, cGauss),R);
           }
 
-          nDelete(&n); /* clean up n */
+          n_Delete(&n,R->cf); /* clean up n */
         }
       }
@@ -181 +189 @@
   /* building the permutation matrix from 'permut' */
   for (int r = 1; r <= rr; r++)
-    MATELEM(pMat, r, permut[r]) = pOne();
+    MATELEM(pMat, r, permut[r]) = p_One(R);
   delete[] permut;
 
@@ -192 +200 @@
  * is given by its LU-decomposition.
  **/
-bool luInverse(const matrix aMat, matrix &iMat)
+bool luInverse(const matrix aMat, matrix &iMat, const ring R)
 { /* aMat is guaranteed to be an (n x n)-matrix */
   matrix pMat;
   matrix lMat;
   matrix uMat;
-  luDecomp(aMat, pMat, lMat, uMat);
-  bool result = luInverseFromLUDecomp(pMat, lMat, uMat, iMat);
+  luDecomp(aMat, pMat, lMat, uMat, R);
+  bool result = luInverseFromLUDecomp(pMat, lMat, uMat, iMat, R);
 
   /* clean-up */
-  idDelete((ideal*)&pMat);
-  idDelete((ideal*)&lMat);
-  idDelete((ideal*)&uMat);
+  id_Delete((ideal*)&pMat,R);
+  id_Delete((ideal*)&lMat,R);
+  id_Delete((ideal*)&uMat,R);
 
   return result;
@@ -222 +230 @@
 }
 
-int luRank(const matrix aMat, const bool isRowEchelon)
+int luRank(const matrix aMat, const bool isRowEchelon, const ring R)
 {
   if (isRowEchelon) return rankFromRowEchelonForm(aMat);
@@ -231 +239 @@
     matrix lMat;
     matrix uMat;
-    luDecomp(aMat, pMat, lMat, uMat);
+    luDecomp(aMat, pMat, lMat, uMat,R);
     int result = rankFromRowEchelonForm(uMat);
 
     /* clean-up */
-    idDelete((ideal*)&pMat);
-    idDelete((ideal*)&lMat);
-    idDelete((ideal*)&uMat);
+    id_Delete((ideal*)&pMat,R);
+    id_Delete((ideal*)&lMat,R);
+    id_Delete((ideal*)&uMat,R);
 
     return result;
@@ -244 +252 @@
 
 bool upperRightTriangleInverse(const matrix uMat, matrix &iMat,
-                               bool diagonalIsOne)
+                               bool diagonalIsOne, const ring R)
 {
   int d = uMat->rows();
@@ -270 +278 @@
     {
       if (diagonalIsOne)
-        MATELEM(iMat, r, r) = pOne();
+        MATELEM(iMat, r, r) = p_One(R);
       else
-        MATELEM(iMat, r, r) = pNSet(nInvers(pGetCoeff(MATELEM(uMat, r, r))));
+        MATELEM(iMat, r, r) = p_NSet(n_Invers(p_GetCoeff(MATELEM(uMat, r, r),R),R->cf),R);
       for (int c = r + 1; c <= d; c++)
       {
@@ -278 +286 @@
         for (int k = r + 1; k <= c; k++)
         {
-          q = ppMult_qq(MATELEM(uMat, r, k), MATELEM(iMat, k, c));
-          p = pAdd(p, q);
+          q = pp_Mult_qq(MATELEM(uMat, r, k), MATELEM(iMat, k, c),R);
+          p = p_Add_q(p, q,R);
         }
-        p = pNeg(p);
-        p = pMult(p, pCopy(MATELEM(iMat, r, r)));
-        pNormalize(p);
+        p = p_Neg(p,R);
+        p = p_Mult_q(p, p_Copy(MATELEM(iMat, r, r),R),R);
+        p_Normalize(p,R);
         MATELEM(iMat, r, c) = p;
       }
@@ -345 +353 @@
  **/
 bool luInverseFromLUDecomp(const matrix pMat, const matrix lMat,
-                           const matrix uMat, matrix &iMat)
+                           const matrix uMat, matrix &iMat, const ring R)
 { /* uMat is guaranteed to be quadratic */
   //int d = uMat->rows();
@@ -359 +367 @@
        entries equal to 1 */
     lowerLeftTriangleInverse(lMat, lMatInverse, true);
-    iMat = mpMult(mpMult(uMatInverse, lMatInverse), pMat);
+    iMat = mp_Mult(mp_Mult(uMatInverse, lMatInverse,R), pMat,R);
 
     /* clean-up */
@@ -892 +900 @@
 
 void hessenberg(const matrix aMat, matrix &pMat, matrix &hessenbergMat,
-                const number tolerance)
+                const number tolerance, const ring R)
 {
   int n = MATROWS(aMat);
@@ -934 +942 @@
         idDelete((ideal*)&u); idDelete((ideal*)&pTmp);
         /* now include pTmpFull in pMat (by letting it act from the left) */
-        pTmp = mpMult(pTmpFull, pMat); idDelete((ideal*)&pMat);
+        pTmp = mp_Mult(pTmpFull, pMat,R); idDelete((ideal*)&pMat);
         pMat = pTmp;
         /* now let pTmpFull act on hessenbergMat from the left and from the
            right (note that pTmpFull is self-inverse) */
-        pTmp = mpMult(pTmpFull, hessenbergMat);
+        pTmp = mp_Mult(pTmpFull, hessenbergMat,R);
         idDelete((ideal*)&hessenbergMat);
-        hessenbergMat = mpMult(pTmp, pTmpFull);
+        hessenbergMat = mp_Mult(pTmp, pTmpFull, R);
         idDelete((ideal*)&pTmp); idDelete((ideal*)&pTmpFull);
         /* as there may be inaccuracy, we erase those entries of hessenbergMat
@@ -964 +972 @@
       matrix &H,             /**< [in/out]  the matrix to be transformed */
       int it,                /**< [in]      iteration index */
-      const number tolerance /**< [in]      accuracy for square roots */
+      const number tolerance,/**< [in]      accuracy for square roots */
+      const ring R
             )
 {
@@ -1048 +1057 @@
     tmp1 = hessenbergStep(c, uVec, hMat, tolerance); nDelete(&tmp1);
     /* now replace H by hMat * H * hMat: */
-    matrix wMat = mpMult(hMat, H); idDelete((ideal*)&H);
-    matrix H1 = mpMult(wMat, hMat);
+    matrix wMat = mp_Mult(hMat, H,R); idDelete((ideal*)&H);
+    matrix H1 = mp_Mult(wMat, hMat,R);
     idDelete((ideal*)&wMat); idDelete((ideal*)&hMat);
     /* now need to re-establish Hessenberg form of H1 and put it in H */
-    hessenberg(H1, wMat, H, tolerance);
+    hessenberg(H1, wMat, H, tolerance,R);
     idDelete((ideal*)&wMat); idDelete((ideal*)&H1);
   }
@@ -1078 +1087 @@
        int& eigenValuesL,
        const number tol1,
-       const number tol2
+       const number tol2,
+       const ring R
          )
 {
@@ -1114 +1124 @@
       /* bring currentMat into Hessenberg form to fasten computations: */
       matrix mm1; matrix mm2;
-      hessenberg(currentMat, mm1, mm2, tol2);
+      hessenberg(currentMat, mm1, mm2, tol2,R);
       idDelete((ideal*)&currentMat); idDelete((ideal*)&mm1);
       currentMat = mm2;
@@ -1143 +1153 @@
         else   /* no deflation found yet */
         {
-          mpTrafo(currentMat, it, tol2);
+          mpTrafo(currentMat, it, tol2,R);
           it++;   /* try again */
         }
@@ -1197 +1207 @@
 }
 
-lists qrDoubleShift(const matrix A, const number tol1, const number tol2,
-                    const number tol3)
-{
-  int n = MATROWS(A);
-  matrix* queue = new matrix[n];
-  queue[0] = mpCopy(A); int queueL = 1;
-  number* eigenVs = new number[n]; int eigenL = 0;
-  /* here comes the main call: */
-  bool worked = qrDS(n, queue, queueL, eigenVs, eigenL, tol1, tol2);
-  lists result = (lists)omAlloc(sizeof(slists));
-  if (!worked)
-  {
-    for (int i = 0; i < eigenL; i++)
-      nDelete(&eigenVs[i]);
-    delete [] eigenVs;
-    for (int i = 0; i < queueL; i++)
-      idDelete((ideal*)&queue[i]);
-    delete [] queue;
-    result->Init(1);
-    result->m[0].rtyp = INT_CMD;
-    result->m[0].data = (void*)0;   /* a list with a single entry
-                                             which is the int zero */
-  }
-  else
-  {
-    /* now eigenVs[0..eigenL-1] contain all eigenvalues; among them, there
-       may be equal entries */
-    number* distinctEVs = new number[n]; int distinctC = 0;
-    int* mults = new int[n];
-    for (int i = 0; i < eigenL; i++)
-    {
-      int index = similar(distinctEVs, distinctC, eigenVs[i], tol3);
-      if (index == -1) /* a new eigenvalue */
-      {
-        distinctEVs[distinctC] = nCopy(eigenVs[i]);
-        mults[distinctC++] = 1;
-      }
-      else mults[index]++;
-      nDelete(&eigenVs[i]);
-    }
-    delete [] eigenVs;
-
-    lists eigenvalues = (lists)omAlloc(sizeof(slists));
-    eigenvalues->Init(distinctC);
-    lists multiplicities = (lists)omAlloc(sizeof(slists));
-    multiplicities->Init(distinctC);
-    for (int i = 0; i < distinctC; i++)
-    {
-      eigenvalues->m[i].rtyp = NUMBER_CMD;
-      eigenvalues->m[i].data = (void*)nCopy(distinctEVs[i]);
-      multiplicities->m[i].rtyp = INT_CMD;
-      multiplicities->m[i].data = (void*)mults[i];
-      nDelete(&distinctEVs[i]);
-    }
-    delete [] distinctEVs; delete [] mults;
-
-    result->Init(2);
-    result->m[0].rtyp = LIST_CMD;
-    result->m[0].data = (char*)eigenvalues;
-    result->m[1].rtyp = LIST_CMD;
-    result->m[1].data = (char*)multiplicities;
-  }
-  return result;
-}
-
 /* This codes assumes that there are at least two variables in the current
    base ring. No assumption is made regarding the monomial ordering. */
@@ -1440 +1385 @@
         {
           t = gg;
-          gg = nGcd(t, pGetCoeff(MATELEM(uMat, r, col)), currRing);
+          gg = nGcd(t, pGetCoeff(MATELEM(uMat, r, col)));
           nDelete(&t);
         }
@@ -1460 +1405 @@
         {
           number g = nGcd(pGetCoeff(MATELEM(uMat, row, col)),
-                          pGetCoeff(MATELEM(uMat, r, col)), currRing);
+                          pGetCoeff(MATELEM(uMat, r, col)));
           number f1 = nDiv(pGetCoeff(MATELEM(uMat, r, col)), g);
           nNormalize(f1);   /* this division works without remainder */

kernel/linearAlgebra.h

@@ -50 +50 @@
        matrix &pMat,      /**< [out] the row permutation matrix P   */
        matrix &lMat,      /**< [out] the lower triangular matrix L  */
-       matrix &uMat       /**< [out] the upper row echelon matrix U */
+       matrix &uMat,      /**< [out] the upper row echelon matrix U */
+       const ring r= currRing       /**< [in] current ring */
              );
 
@@ -62 +63 @@
  **/
 int pivotScore(
-               number n  /**< [in] a non-zero matrix entry */
+               number n,    /**< [in] a non-zero matrix entry */
+               const ring r= currRing /**< [in] current ring */
               );
 
@@ -87 +89 @@
            int* bestR,        /**< [out] address of row index of best
                                          pivot element                  */
-           int* bestC         /**< [out] address of column index of
+           int* bestC,        /**< [out] address of column index of
                                          best pivot element             */
+           const ring r= currRing       /**< [in] current ring */
           );
 
@@ -108 +111 @@
 bool luInverse(
                const matrix aMat, /**< [in]  matrix to be inverted */
-               matrix &iMat       /**< [out] inverse of aMat if
+               matrix &iMat,      /**< [out] inverse of aMat if
                                              invertible            */
+               const ring r=currRing /**< [in] current ring */
               );
 
@@ -136 +140 @@
                                      triangular form               */
        matrix &iMat,      /**< [out] inverse of uMat if invertible */
-       bool diagonalIsOne /**< [in]  if true, then all diagonal
+       bool diagonalIsOne,/**< [in]  if true, then all diagonal
                                      entries of uMat are 1         */
+       const ring r= currRing /**< [in] current ring */
                               );
 
@@ -197 +202 @@
        const matrix uMat, /**< [in]  upper right matrix of an LU-
                                      decomposition                */
-       matrix &iMat       /**< [out] inverse of A if invertible   */
+       matrix &iMat,      /**< [out] inverse of A if invertible   */
+       const ring r= currRing
                           );
 
@@ -215 +221 @@
 int luRank(
        const matrix aMat,      /**< [in]  input matrix */
-       const bool isRowEchelon /**< [in]  if true then aMat is assumed to be
+       const bool isRowEchelon,/**< [in]  if true then aMat is assumed to be
                                           already in row echelon form */
+       const ring r= currRing
           );
 
@@ -337 +344 @@
 bool unitMatrix(
        const int n,     /**< [in]  size of the matrix */
-       matrix &unitMat  /**< [out] the new (nxn) unit matrix */
+       matrix &unitMat,  /**< [out] the new (nxn) unit matrix */
+       const ring r= currRing /** [in] current ring */
                );
 
@@ -457 +465 @@
        matrix &pMat,           /**< [out] the transformation matrix */
        matrix &hessenbergMat,  /**< [out] the Hessenberg form of aMat */
-       const number tolerance  /**< [in]  accuracy */
+       const number tolerance, /**< [in]  accuracy */
+       const ring r
                );
 
@@ -484 +493 @@
                   );
 
-/**
- * Computes all eigenvalues of a given real quadratic matrix with
- * multiplicites.
- *
- * The method assumes that the current ground field is the complex numbers.
- * Computations are based on the QR double shift algorithm involving
- * Hessenberg form and householder transformations.
- * If the algorithm works, then it returns a list with two entries which
- * are again lists of the same size:
- * _[1][i] is the i-th mutually distinct eigenvalue that was found,
- * _[2][i] is the (int) multiplicity of _[1][i].
- * If the algorithm does not work (due to an ill-posed matrix), a list with
- * the single entry (int)0 is returned.
- * 'tol1' is used for detection of deflation in the actual qr double shift
- * algorithm.
- * 'tol2' is used for ending Heron's iteration whenever square roots
- * are being computed.
- * 'tol3' is used to distinguish between distinct eigenvalues: When
- * the Euclidean distance between two computed eigenvalues is less then
- * tol3, then they will be regarded equal, resulting in a higher
- * multiplicity of the corresponding eigenvalue.
- *
- * @return a list with one entry (int)0, or two entries which are again lists
- **/
-lists qrDoubleShift(
-       const matrix A,     /**< [in]  the quadratic matrix         */
-       const number tol1,  /**< [in]  tolerance for deflation      */
-       const number tol2,  /**< [in]  tolerance for square roots   */
-       const number tol3   /**< [in]  tolerance for distinguishing
-                                      eigenvalues                  */
-                   );
 
 /**

libpolys/coeffs/numbers.h

@@ -10 +10 @@
 #include <coeffs/coeffs.h>
 
-/*
 // the access methods
 //
@@ -28 +27 @@
 #define nIsMOne(n)        n_IsMOne(n, currRing->cf)
 #define nGreaterZero(n)   n_GreaterZero(n, currRing->cf)
+#define nGreater(a, b)   n_Greater (a,b,currRing->cf)
 #define nWrite(n)         n_Write(n,currRing->cf)
 #define nNormalize(n)     n_Normalize(n,currRing->cf)
@@ -44 +44 @@
 
 #define nSetMap(R)        n_SetMap(R,currRing->cf)
-*/
 
 

libpolys/polys/monomials/ring.h

@@ -30 +30 @@
 struct p_Procs_s;
 typedef struct p_Procs_s p_Procs_s;
-class slists;
-typedef slists *           lists;
+//class slists;
+//typedef slists *           lists;
 class kBucket;
 typedef kBucket*           kBucket_pt;
@@ -716 +716 @@
 void rSetWeightVec(ring r, int64 *wv);
 
-lists rDecompose(const ring r);
-ring rCompose(const lists  L);
+//lists rDecompose(const ring r);
+//ring rCompose(const lists  L);
 /////////////////////////////
 // Auxillary functions

libpolys/polys/polys.h

@@ -18 +18 @@
  *
  ***************************************************************/
-/*
+
 // deletes old coeff before setting the new one
 #define pSetCoeff(p,n)      p_SetCoeff(p,n,currRing)
@@ -39 +39 @@
 #define pGetExpSum(p1, p2, i)    p_GetExpSum(p1, p2, i, currRing)
 #define pGetExpDiff(p1, p2, i)   p_GetExpDiff(p1, p2, i, currRing)
-*/
+
 
 /***************************************************************