Changeset 6ed8c4 in git


Ignore:
Timestamp:
Jul 19, 2011, 4:45:41 PM (13 years ago)
Author:
mlee <martinlee84@…>
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', '79dfb9a6d258bfeb991428bdb25b8f55e9e809b2')
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:
1 added
5 edited

Legend:

Unmodified
Added
Removed
  • kernel/linearAlgebra.cc

    rbb5c28 r6ed8c4  
    1717
    1818// include header files
    19 #include <kernel/mod2.h>
     19#include "mod2.h"
     20
     21#include <coeffs/coeffs.h>
     22#include <coeffs/numbers.h>
     23
     24#include <coeffs/mpr_complex.h>
     25#include <polys/monomials/p_polys.h>
     26#include <polys/simpleideals.h>
     27#include <polys/matpol.h>
     28
     29// #include <polys/polys.h>
     30
    2031#include <kernel/structs.h>
    21 #include <polys/polys.h>
    2232#include <kernel/ideals.h>
    23 #include <coeffs/numbers.h>
    24 #include <polys/matpol.h>
    25 #include <coeffs/mpr_complex.h>
    2633#include <kernel/linearAlgebra.h>
    2734
     
    3946 * the negative of nSize will be returned.
    4047 **/
    41 int pivotScore(number n)
    42 {
    43   int s = nSize(n);
    44   if (rField_is_long_C(currRing) ||
    45       rField_is_long_R(currRing) ||
    46       rField_is_R(currRing))
     48int pivotScore(number n, const ring r)
     49{
     50  int s = n_Size(n, r->cf);
     51  if (rField_is_long_C(r) ||
     52      rField_is_long_R(r) ||
     53      rField_is_R(r))
    4754    return -s;
    4855  else
     
    5865 **/
    5966bool pivot(const matrix aMat, const int r1, const int r2, const int c1,
    60            const int c2, int* bestR, int* bestC)
     67           const int c2, int* bestR, int* bestC, const ring R)
    6168{
    6269  int bestScore; int score; bool foundBestScore = false; poly matEntry;
     
    6976      if (matEntry != NULL)
    7077      {
    71         score = pivotScore(pGetCoeff(matEntry));
     78        score = pivotScore(pGetCoeff(matEntry), R);
    7279        if ((!foundBestScore) || (score < bestScore))
    7380        {
     
    8491}
    8592
    86 bool unitMatrix(const int n, matrix &unitMat)
     93bool unitMatrix(const int n, matrix &unitMat, const ring R)
    8794{
    8895  if (n < 1) return false;
    8996  unitMat = mpNew(n, n);
    90   for (int r = 1; r <= n; r++) MATELEM(unitMat, r, r) = pOne();
     97  for (int r = 1; r <= n; r++) MATELEM(unitMat, r, r) = p_One(R);
    9198  return true;
    9299}
    93100
    94 void luDecomp(const matrix aMat, matrix &pMat, matrix &lMat, matrix &uMat)
     101void luDecomp(const matrix aMat, matrix &pMat, matrix &lMat, matrix &uMat,
     102              const ring R)
    95103{
    96104  int rr = aMat->rows();
     
    101109  for (int r = 1; r <= rr; r++)
    102110    for (int c = 1; c <= cc; c++)
    103       MATELEM(uMat, r, c) = pCopy(aMat->m[c - 1 + (r - 1) * cc]);
     111      MATELEM(uMat, r, c) = p_Copy(aMat->m[c - 1 + (r - 1) * cc], R);
    104112
    105113  /* we use an int array to store all row permutations;
     
    109117
    110118  /* fill lMat with the (rr x rr) unit matrix */
    111   unitMatrix(rr, lMat);
     119  unitMatrix(rr, lMat,R);
    112120
    113121  int bestR; int bestC; int intSwap; poly pSwap; int cOffset = 0;
     
    116124    if (r > cc) break;
    117125    while ((r + cOffset <= cc) &&
    118            (!pivot(uMat, r, rr, r + cOffset, r + cOffset, &bestR, &bestC)))
     126           (!pivot(uMat, r, rr, r + cOffset, r + cOffset, &bestR, &bestC, R)))
    119127      cOffset++;
    120128    if (r + cOffset <= cc)
     
    155163        if (p != NULL)
    156164        {
    157           n = nDiv(pGetCoeff(p), pivotElement);
    158           nNormalize(n);
     165          n = n_Div(pGetCoeff(p), pivotElement, R->cf);
     166          n_Normalize(n,R->cf);
    159167
    160168          /* filling lMat;
    161169             old entry was zero, so no need to call pDelete(.) */
    162           MATELEM(lMat, rGauss, r) = pNSet(nCopy(n));
     170          MATELEM(lMat, rGauss, r) = p_NSet(n_Copy(n,R->cf),R);
    163171
    164172          /* adjusting uMat: */
    165           MATELEM(uMat, rGauss, r + cOffset) = NULL; pDelete(&p);
    166           n = nNeg(n);
     173          MATELEM(uMat, rGauss, r + cOffset) = NULL; p_Delete(&p,R);
     174          n = n_Neg(n,R->cf);
    167175          for (int cGauss = r + cOffset + 1; cGauss <= cc; cGauss++)
    168176          {
    169177            MATELEM(uMat, rGauss, cGauss)
    170               = pAdd(MATELEM(uMat, rGauss, cGauss),
    171                      ppMult_nn(MATELEM(uMat, r, cGauss), n));
    172             pNormalize(MATELEM(uMat, rGauss, cGauss));
     178              = p_Add_q(MATELEM(uMat, rGauss, cGauss),
     179                     pp_Mult_nn(MATELEM(uMat, r, cGauss), n, R), R);
     180            p_Normalize(MATELEM(uMat, rGauss, cGauss),R);
    173181          }
    174182
    175           nDelete(&n); /* clean up n */
     183          n_Delete(&n,R->cf); /* clean up n */
    176184        }
    177185      }
     
    181189  /* building the permutation matrix from 'permut' */
    182190  for (int r = 1; r <= rr; r++)
    183     MATELEM(pMat, r, permut[r]) = pOne();
     191    MATELEM(pMat, r, permut[r]) = p_One(R);
    184192  delete[] permut;
    185193
     
    192200 * is given by its LU-decomposition.
    193201 **/
    194 bool luInverse(const matrix aMat, matrix &iMat)
     202bool luInverse(const matrix aMat, matrix &iMat, const ring R)
    195203{ /* aMat is guaranteed to be an (n x n)-matrix */
    196204  matrix pMat;
    197205  matrix lMat;
    198206  matrix uMat;
    199   luDecomp(aMat, pMat, lMat, uMat);
    200   bool result = luInverseFromLUDecomp(pMat, lMat, uMat, iMat);
     207  luDecomp(aMat, pMat, lMat, uMat, R);
     208  bool result = luInverseFromLUDecomp(pMat, lMat, uMat, iMat, R);
    201209
    202210  /* clean-up */
    203   idDelete((ideal*)&pMat);
    204   idDelete((ideal*)&lMat);
    205   idDelete((ideal*)&uMat);
     211  id_Delete((ideal*)&pMat,R);
     212  id_Delete((ideal*)&lMat,R);
     213  id_Delete((ideal*)&uMat,R);
    206214
    207215  return result;
     
    222230}
    223231
    224 int luRank(const matrix aMat, const bool isRowEchelon)
     232int luRank(const matrix aMat, const bool isRowEchelon, const ring R)
    225233{
    226234  if (isRowEchelon) return rankFromRowEchelonForm(aMat);
     
    231239    matrix lMat;
    232240    matrix uMat;
    233     luDecomp(aMat, pMat, lMat, uMat);
     241    luDecomp(aMat, pMat, lMat, uMat,R);
    234242    int result = rankFromRowEchelonForm(uMat);
    235243
    236244    /* clean-up */
    237     idDelete((ideal*)&pMat);
    238     idDelete((ideal*)&lMat);
    239     idDelete((ideal*)&uMat);
     245    id_Delete((ideal*)&pMat,R);
     246    id_Delete((ideal*)&lMat,R);
     247    id_Delete((ideal*)&uMat,R);
    240248
    241249    return result;
     
    244252
    245253bool upperRightTriangleInverse(const matrix uMat, matrix &iMat,
    246                                bool diagonalIsOne)
     254                               bool diagonalIsOne, const ring R)
    247255{
    248256  int d = uMat->rows();
     
    270278    {
    271279      if (diagonalIsOne)
    272         MATELEM(iMat, r, r) = pOne();
     280        MATELEM(iMat, r, r) = p_One(R);
    273281      else
    274         MATELEM(iMat, r, r) = pNSet(nInvers(pGetCoeff(MATELEM(uMat, r, r))));
     282        MATELEM(iMat, r, r) = p_NSet(n_Invers(p_GetCoeff(MATELEM(uMat, r, r),R),R->cf),R);
    275283      for (int c = r + 1; c <= d; c++)
    276284      {
     
    278286        for (int k = r + 1; k <= c; k++)
    279287        {
    280           q = ppMult_qq(MATELEM(uMat, r, k), MATELEM(iMat, k, c));
    281           p = pAdd(p, q);
     288          q = pp_Mult_qq(MATELEM(uMat, r, k), MATELEM(iMat, k, c),R);
     289          p = p_Add_q(p, q,R);
    282290        }
    283         p = pNeg(p);
    284         p = pMult(p, pCopy(MATELEM(iMat, r, r)));
    285         pNormalize(p);
     291        p = p_Neg(p,R);
     292        p = p_Mult_q(p, p_Copy(MATELEM(iMat, r, r),R),R);
     293        p_Normalize(p,R);
    286294        MATELEM(iMat, r, c) = p;
    287295      }
     
    345353 **/
    346354bool luInverseFromLUDecomp(const matrix pMat, const matrix lMat,
    347                            const matrix uMat, matrix &iMat)
     355                           const matrix uMat, matrix &iMat, const ring R)
    348356{ /* uMat is guaranteed to be quadratic */
    349357  //int d = uMat->rows();
     
    359367       entries equal to 1 */
    360368    lowerLeftTriangleInverse(lMat, lMatInverse, true);
    361     iMat = mpMult(mpMult(uMatInverse, lMatInverse), pMat);
     369    iMat = mp_Mult(mp_Mult(uMatInverse, lMatInverse,R), pMat,R);
    362370
    363371    /* clean-up */
     
    892900
    893901void hessenberg(const matrix aMat, matrix &pMat, matrix &hessenbergMat,
    894                 const number tolerance)
     902                const number tolerance, const ring R)
    895903{
    896904  int n = MATROWS(aMat);
     
    934942        idDelete((ideal*)&u); idDelete((ideal*)&pTmp);
    935943        /* now include pTmpFull in pMat (by letting it act from the left) */
    936         pTmp = mpMult(pTmpFull, pMat); idDelete((ideal*)&pMat);
     944        pTmp = mp_Mult(pTmpFull, pMat,R); idDelete((ideal*)&pMat);
    937945        pMat = pTmp;
    938946        /* now let pTmpFull act on hessenbergMat from the left and from the
    939947           right (note that pTmpFull is self-inverse) */
    940         pTmp = mpMult(pTmpFull, hessenbergMat);
     948        pTmp = mp_Mult(pTmpFull, hessenbergMat,R);
    941949        idDelete((ideal*)&hessenbergMat);
    942         hessenbergMat = mpMult(pTmp, pTmpFull);
     950        hessenbergMat = mp_Mult(pTmp, pTmpFull, R);
    943951        idDelete((ideal*)&pTmp); idDelete((ideal*)&pTmpFull);
    944952        /* as there may be inaccuracy, we erase those entries of hessenbergMat
     
    964972      matrix &H,             /**< [in/out]  the matrix to be transformed */
    965973      int it,                /**< [in]      iteration index */
    966       const number tolerance /**< [in]      accuracy for square roots */
     974      const number tolerance,/**< [in]      accuracy for square roots */
     975      const ring R
    967976            )
    968977{
     
    10481057    tmp1 = hessenbergStep(c, uVec, hMat, tolerance); nDelete(&tmp1);
    10491058    /* now replace H by hMat * H * hMat: */
    1050     matrix wMat = mpMult(hMat, H); idDelete((ideal*)&H);
    1051     matrix H1 = mpMult(wMat, hMat);
     1059    matrix wMat = mp_Mult(hMat, H,R); idDelete((ideal*)&H);
     1060    matrix H1 = mp_Mult(wMat, hMat,R);
    10521061    idDelete((ideal*)&wMat); idDelete((ideal*)&hMat);
    10531062    /* now need to re-establish Hessenberg form of H1 and put it in H */
    1054     hessenberg(H1, wMat, H, tolerance);
     1063    hessenberg(H1, wMat, H, tolerance,R);
    10551064    idDelete((ideal*)&wMat); idDelete((ideal*)&H1);
    10561065  }
     
    10781087       int& eigenValuesL,
    10791088       const number tol1,
    1080        const number tol2
     1089       const number tol2,
     1090       const ring R
    10811091         )
    10821092{
     
    11141124      /* bring currentMat into Hessenberg form to fasten computations: */
    11151125      matrix mm1; matrix mm2;
    1116       hessenberg(currentMat, mm1, mm2, tol2);
     1126      hessenberg(currentMat, mm1, mm2, tol2,R);
    11171127      idDelete((ideal*)&currentMat); idDelete((ideal*)&mm1);
    11181128      currentMat = mm2;
     
    11431153        else   /* no deflation found yet */
    11441154        {
    1145           mpTrafo(currentMat, it, tol2);
     1155          mpTrafo(currentMat, it, tol2,R);
    11461156          it++;   /* try again */
    11471157        }
     
    11971207}
    11981208
    1199 lists qrDoubleShift(const matrix A, const number tol1, const number tol2,
    1200                     const number tol3)
    1201 {
    1202   int n = MATROWS(A);
    1203   matrix* queue = new matrix[n];
    1204   queue[0] = mpCopy(A); int queueL = 1;
    1205   number* eigenVs = new number[n]; int eigenL = 0;
    1206   /* here comes the main call: */
    1207   bool worked = qrDS(n, queue, queueL, eigenVs, eigenL, tol1, tol2);
    1208   lists result = (lists)omAlloc(sizeof(slists));
    1209   if (!worked)
    1210   {
    1211     for (int i = 0; i < eigenL; i++)
    1212       nDelete(&eigenVs[i]);
    1213     delete [] eigenVs;
    1214     for (int i = 0; i < queueL; i++)
    1215       idDelete((ideal*)&queue[i]);
    1216     delete [] queue;
    1217     result->Init(1);
    1218     result->m[0].rtyp = INT_CMD;
    1219     result->m[0].data = (void*)0;   /* a list with a single entry
    1220                                              which is the int zero */
    1221   }
    1222   else
    1223   {
    1224     /* now eigenVs[0..eigenL-1] contain all eigenvalues; among them, there
    1225        may be equal entries */
    1226     number* distinctEVs = new number[n]; int distinctC = 0;
    1227     int* mults = new int[n];
    1228     for (int i = 0; i < eigenL; i++)
    1229     {
    1230       int index = similar(distinctEVs, distinctC, eigenVs[i], tol3);
    1231       if (index == -1) /* a new eigenvalue */
    1232       {
    1233         distinctEVs[distinctC] = nCopy(eigenVs[i]);
    1234         mults[distinctC++] = 1;
    1235       }
    1236       else mults[index]++;
    1237       nDelete(&eigenVs[i]);
    1238     }
    1239     delete [] eigenVs;
    1240 
    1241     lists eigenvalues = (lists)omAlloc(sizeof(slists));
    1242     eigenvalues->Init(distinctC);
    1243     lists multiplicities = (lists)omAlloc(sizeof(slists));
    1244     multiplicities->Init(distinctC);
    1245     for (int i = 0; i < distinctC; i++)
    1246     {
    1247       eigenvalues->m[i].rtyp = NUMBER_CMD;
    1248       eigenvalues->m[i].data = (void*)nCopy(distinctEVs[i]);
    1249       multiplicities->m[i].rtyp = INT_CMD;
    1250       multiplicities->m[i].data = (void*)mults[i];
    1251       nDelete(&distinctEVs[i]);
    1252     }
    1253     delete [] distinctEVs; delete [] mults;
    1254 
    1255     result->Init(2);
    1256     result->m[0].rtyp = LIST_CMD;
    1257     result->m[0].data = (char*)eigenvalues;
    1258     result->m[1].rtyp = LIST_CMD;
    1259     result->m[1].data = (char*)multiplicities;
    1260   }
    1261   return result;
    1262 }
    1263 
    12641209/* This codes assumes that there are at least two variables in the current
    12651210   base ring. No assumption is made regarding the monomial ordering. */
     
    14401385        {
    14411386          t = gg;
    1442           gg = nGcd(t, pGetCoeff(MATELEM(uMat, r, col)), currRing);
     1387          gg = nGcd(t, pGetCoeff(MATELEM(uMat, r, col)));
    14431388          nDelete(&t);
    14441389        }
     
    14601405        {
    14611406          number g = nGcd(pGetCoeff(MATELEM(uMat, row, col)),
    1462                           pGetCoeff(MATELEM(uMat, r, col)), currRing);
     1407                          pGetCoeff(MATELEM(uMat, r, col)));
    14631408          number f1 = nDiv(pGetCoeff(MATELEM(uMat, r, col)), g);
    14641409          nNormalize(f1);   /* this division works without remainder */
  • kernel/linearAlgebra.h

    rbb5c28 r6ed8c4  
    5050       matrix &pMat,      /**< [out] the row permutation matrix P   */
    5151       matrix &lMat,      /**< [out] the lower triangular matrix L  */
    52        matrix &uMat       /**< [out] the upper row echelon matrix U */
     52       matrix &uMat,      /**< [out] the upper row echelon matrix U */
     53       const ring r= currRing       /**< [in] current ring */
    5354             );
    5455
     
    6263 **/
    6364int pivotScore(
    64                number n  /**< [in] a non-zero matrix entry */
     65               number n,    /**< [in] a non-zero matrix entry */
     66               const ring r= currRing /**< [in] current ring */
    6567              );
    6668
     
    8789           int* bestR,        /**< [out] address of row index of best
    8890                                         pivot element                  */
    89            int* bestC         /**< [out] address of column index of
     91           int* bestC,        /**< [out] address of column index of
    9092                                         best pivot element             */
     93           const ring r= currRing       /**< [in] current ring */
    9194          );
    9295
     
    108111bool luInverse(
    109112               const matrix aMat, /**< [in]  matrix to be inverted */
    110                matrix &iMat       /**< [out] inverse of aMat if
     113               matrix &iMat,      /**< [out] inverse of aMat if
    111114                                             invertible            */
     115               const ring r=currRing /**< [in] current ring */
    112116              );
    113117
     
    136140                                     triangular form               */
    137141       matrix &iMat,      /**< [out] inverse of uMat if invertible */
    138        bool diagonalIsOne /**< [in]  if true, then all diagonal
     142       bool diagonalIsOne,/**< [in]  if true, then all diagonal
    139143                                     entries of uMat are 1         */
     144       const ring r= currRing /**< [in] current ring */
    140145                              );
    141146
     
    197202       const matrix uMat, /**< [in]  upper right matrix of an LU-
    198203                                     decomposition                */
    199        matrix &iMat       /**< [out] inverse of A if invertible   */
     204       matrix &iMat,      /**< [out] inverse of A if invertible   */
     205       const ring r= currRing
    200206                          );
    201207
     
    215221int luRank(
    216222       const matrix aMat,      /**< [in]  input matrix */
    217        const bool isRowEchelon /**< [in]  if true then aMat is assumed to be
     223       const bool isRowEchelon,/**< [in]  if true then aMat is assumed to be
    218224                                          already in row echelon form */
     225       const ring r= currRing
    219226          );
    220227
     
    337344bool unitMatrix(
    338345       const int n,     /**< [in]  size of the matrix */
    339        matrix &unitMat  /**< [out] the new (nxn) unit matrix */
     346       matrix &unitMat,  /**< [out] the new (nxn) unit matrix */
     347       const ring r= currRing /** [in] current ring */
    340348               );
    341349
     
    457465       matrix &pMat,           /**< [out] the transformation matrix */
    458466       matrix &hessenbergMat,  /**< [out] the Hessenberg form of aMat */
    459        const number tolerance  /**< [in]  accuracy */
     467       const number tolerance, /**< [in]  accuracy */
     468       const ring r
    460469               );
    461470
     
    484493                  );
    485494
    486 /**
    487  * Computes all eigenvalues of a given real quadratic matrix with
    488  * multiplicites.
    489  *
    490  * The method assumes that the current ground field is the complex numbers.
    491  * Computations are based on the QR double shift algorithm involving
    492  * Hessenberg form and householder transformations.
    493  * If the algorithm works, then it returns a list with two entries which
    494  * are again lists of the same size:
    495  * _[1][i] is the i-th mutually distinct eigenvalue that was found,
    496  * _[2][i] is the (int) multiplicity of _[1][i].
    497  * If the algorithm does not work (due to an ill-posed matrix), a list with
    498  * the single entry (int)0 is returned.
    499  * 'tol1' is used for detection of deflation in the actual qr double shift
    500  * algorithm.
    501  * 'tol2' is used for ending Heron's iteration whenever square roots
    502  * are being computed.
    503  * 'tol3' is used to distinguish between distinct eigenvalues: When
    504  * the Euclidean distance between two computed eigenvalues is less then
    505  * tol3, then they will be regarded equal, resulting in a higher
    506  * multiplicity of the corresponding eigenvalue.
    507  *
    508  * @return a list with one entry (int)0, or two entries which are again lists
    509  **/
    510 lists qrDoubleShift(
    511        const matrix A,     /**< [in]  the quadratic matrix         */
    512        const number tol1,  /**< [in]  tolerance for deflation      */
    513        const number tol2,  /**< [in]  tolerance for square roots   */
    514        const number tol3   /**< [in]  tolerance for distinguishing
    515                                       eigenvalues                  */
    516                    );
    517495
    518496/**
  • libpolys/coeffs/numbers.h

    rbb5c28 r6ed8c4  
    1010#include <coeffs/coeffs.h>
    1111
    12 /*
    1312// the access methods
    1413//
     
    2827#define nIsMOne(n)        n_IsMOne(n, currRing->cf)
    2928#define nGreaterZero(n)   n_GreaterZero(n, currRing->cf)
     29#define nGreater(a, b)   n_Greater (a,b,currRing->cf)
    3030#define nWrite(n)         n_Write(n,currRing->cf)
    3131#define nNormalize(n)     n_Normalize(n,currRing->cf)
     
    4444
    4545#define nSetMap(R)        n_SetMap(R,currRing->cf)
    46 */
    4746
    4847
  • libpolys/polys/monomials/ring.h

    rbb5c28 r6ed8c4  
    3030struct p_Procs_s;
    3131typedef struct p_Procs_s p_Procs_s;
    32 class slists;
    33 typedef slists *           lists;
     32//class slists;
     33//typedef slists *           lists;
    3434class kBucket;
    3535typedef kBucket*           kBucket_pt;
     
    716716void rSetWeightVec(ring r, int64 *wv);
    717717
    718 lists rDecompose(const ring r);
    719 ring rCompose(const lists  L);
     718//lists rDecompose(const ring r);
     719//ring rCompose(const lists  L);
    720720/////////////////////////////
    721721// Auxillary functions
  • libpolys/polys/polys.h

    rbb5c28 r6ed8c4  
    1818 *
    1919 ***************************************************************/
    20 /*
     20
    2121// deletes old coeff before setting the new one
    2222#define pSetCoeff(p,n)      p_SetCoeff(p,n,currRing)
     
    3939#define pGetExpSum(p1, p2, i)    p_GetExpSum(p1, p2, i, currRing)
    4040#define pGetExpDiff(p1, p2, i)   p_GetExpDiff(p1, p2, i, currRing)
    41 */
     41
    4242
    4343/***************************************************************
Note: See TracChangeset for help on using the changeset viewer.