Changeset 2f436b in git for Singular/sparsmat.cc


Ignore:
Timestamp:
Dec 31, 2000, 4:14:47 PM (23 years ago)
Author:
Olaf Bachmann <obachman@…>
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', '38dfc5131670d387a89455159ed1e071997eec94')
Children:
e609098c45a74ac91c002ffa7ece5eebe7f8c002
Parents:
33ec1145a109507ad3e3cf4a69a847b703358e93
Message:
* version 1-3-13: sparsemat improivements


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

Legend:

Unmodified
Added
Removed

  • Singular/sparsmat.cc

    r33ec11 r2f436b  
    22*  Computer Algebra System SINGULAR     *
    33****************************************/
    4 /* $Id: sparsmat.cc,v 1.46 2000-12-06 11:03:30 Singular Exp $ */
     4/* $Id: sparsmat.cc,v 1.47 2000-12-31 15:14:45 obachman Exp $ */
    55
    66/*
    … …  
    2222#include "sparsmat.h"
    2323#include "prCopy.h"
    24 
    25 /* ----------------- general definitions ------------------ */
    26 /* in structs.h
    27 typedef struct smprec sm_prec;
    28 typedef sm_prec * smpoly;
    29 struct smprec{
    30   smpoly n;            // the next element
    31   int pos;             // position
    32   int e;               // level
    33   poly m;              // the element
    34   float f;             // complexity of the element
    35 };
    36 */
    37 
    38 #if OLD > 0
    39 static poly smSelectCopy(poly, const poly);
     24#include "p_Procs.h"
     25#include "kbuckets.h"
     26#include "p_Mult_q.h"
     27
     28// define SM_NO_BUCKETS, if sparsemat stuff should not use buckets
     29// #define SM_NO_BUCKETS
     30
     31// this is also influenced by TEST_OPT_NOTBUCKETS
     32#ifndef SM_NO_BUCKETS
     33// buckets do carry a small additional overhead: only use them if
     34// min-length of polys is >= SM_MIN_LENGTH_BUCKET
     35#define SM_MIN_LENGTH_BUCKET MIN_LENGTH_BUCKET - 5
    4036#else
    41 #define smSelectCopy ppMult_Coeff_mm_DivSelect
     37#define SM_MIN_LENGTH_BUCKET    INT_MAX
    4238#endif
     39
    4340
    4441/* declare internal 'C' stuff */
    … …  
    5855static BOOLEAN smHaveDenom(poly);
    5956static number smCleardenom(ideal);
     57
     58static poly pp_Mult_Coeff_mm_DivSelect_MultDiv(poly p, int &lp, poly m,
     59                                               poly a, poly b)
     60{
     61  if (rOrd_is_c_dp(currRing))
     62  {
     63    int shorter;
     64    p = currRing->p_Procs->pp_Mult_Coeff_mm_DivSelectMult(p, m, a, b,
     65                                                          shorter, currRing);
     66    lp -= shorter;
     67  }
     68  else
     69  {
     70    p = pp_Mult_Coeff_mm_DivSelect(p, lp, m, currRing);
     71    smExpMultDiv(p, a, b);
     72  }
     73  return p;
     74}
     75
     76static poly smSelectCopy_ExpMultDiv(poly p, poly m, poly a, poly b)
     77{
     78  int lp = 0;
     79  return pp_Mult_Coeff_mm_DivSelect_MultDiv(p, lp, m, a, b);
     80}
     81
    6082
    6183/* class for sparse matrix:
    … …  
    134156};
    135157
     158Exponent_t smExpBound(ideal Id)
     159{
     160  Exponent_t max = 0;
     161  long ldeg;
     162  int dummy, i, n = IDELEMS(Id);
     163 
     164  for (i=0; i<n; i++)
     165  {
     166    if (Id->m[i] != NULL)
     167    {
     168      ldeg = pLDeg(Id->m[i], &dummy);
     169      if (ldeg > max) max = ldeg;
     170    }
     171  }
     172  return max*2*n;
     173}
     174
     175// #define HOMOG_LP
    136176/* ----------------- ops with rings ------------------ */
    137177ideal smRingCopy(ideal I, ring *ri, sip_sring &tmpR)
    … …  
    139179  ring origR =NULL;
    140180  ideal II;
    141   if (currRing->order[0]!=ringorder_c)
    142   {
    143     origR =currRing;
    144     tmpR=*origR;
    145     int *ord=(int*)omAlloc0(3*sizeof(int));
    146     int *block0=(int*)omAlloc(3*sizeof(int));
    147     int *block1=(int*)omAlloc(3*sizeof(int));
    148     ord[0]=ringorder_c;
     181  origR =currRing;
     182  tmpR=*origR;
     183  int *ord=(int*)omAlloc0(3*sizeof(int));
     184  int *block0=(int*)omAlloc(3*sizeof(int));
     185  int *block1=(int*)omAlloc(3*sizeof(int));
     186  ord[0]=ringorder_c;
     187#ifdef HOMOG_LP
     188  if (idHomIdeal(I))
     189    ord[1] = ringorder_lp;
     190  else
     191#endif
    149192    ord[1]=ringorder_dp;
    150     tmpR.order=ord;
    151     tmpR.OrdSgn=1;
    152     block0[1]=1;
    153     tmpR.block0=block0;
    154     block1[1]=tmpR.N;
    155     tmpR.block1=block1;
    156     rComplete(&tmpR,1);
    157     rChangeCurrRing(&tmpR);
    158     // fetch data from the old ring
    159     II=idInit(IDELEMS(I),I->rank);
    160     int k;
    161     for (k=0;k<IDELEMS(I);k++) II->m[k] = prCopyR( I->m[k], origR);
    162   }
    163   else
    164   {
    165     II=idCopy(I);
    166   }
     193  tmpR.order=ord;
     194  tmpR.OrdSgn=1;
     195  block0[1]=1;
     196  tmpR.block0=block0;
     197  block1[1]=tmpR.N;
     198  tmpR.block1=block1;
     199
     200  tmpR.bitmask = smExpBound(I);
     201
     202  // unfortunately, we can not work (yet) with r->N == 0
     203  if (tmpR.bitmask < 1) tmpR.bitmask = 1;
     204  if (tmpR.bitmask > currRing->bitmask) tmpR.bitmask = currRing->bitmask;
     205
     206  rComplete(&tmpR,1);
     207  rChangeCurrRing(&tmpR);
     208  if (TEST_OPT_PROT)
     209    Print("[%d:%d]", (long) tmpR.bitmask, tmpR.ExpL_Size);
     210  // fetch data from the old ring
     211  II = idrCopyR(I, origR);
     212  idTest(II);
    167213  *ri = origR;
    168214  return II;
    … …  
    951997        {
    952998          res = res->n = smElemCopy(b);
    953           res->m = smMult(b->m, w);
     999          res->m = ppMult_qq(b->m, w);
    9541000          res->e = 1;
    9551001          res->f = smPolyWeight(res);
    … …  
    9661012      {
    9671013        res = res->n = smElemCopy(b);
    968         res->m = smMult(b->m, w);
     1014        res->m = ppMult_qq(b->m, w);
    9691015        res->e = 1;
    9701016        res->f = smPolyWeight(res);
    … …  
    9731019      else
    9741020      {
    975         ha = smMult(a->m, p);
     1021        ha = ppMult_qq(a->m, p);
    9761022        pDelete(&a->m);
    977         hb = smMult(b->m, w);
     1023        hb = ppMult_qq(b->m, w);
    9781024        ha = pAdd(ha, hb);
    9791025        if (ha != NULL)
    … …  
    16801726
    16811727/* ----------------- arithmetic ------------------ */
    1682 
    1683 /*
    1684 *  returns a*b
    1685 *  a,b NOT destroyed
    1686 */
    1687 poly smMult(poly a, poly b)
    1688 {
    1689   poly pa, res, r;
    1690 
    1691   if (smSmaller(a, b))
    1692   {
    1693     r = a;
    1694     a = b;
    1695     b = r;
    1696   }
    1697   if (pNext(b) == NULL)
    1698   {
    1699     if (pLmIsConstantComp(b))
    1700       return ppMult_nn(a, pGetCoeff(b));
    1701     else
    1702       return ppMult_mm(a, b);
    1703   }
    1704   pa = res = ppMult_mm(a, b);
    1705   pIter(b);
    1706   do
    1707   {
    1708     r = ppMult_mm(a, b);
    1709     smCombineChain(&pa, r);
    1710     pIter(b);
    1711   } while (b != NULL);
    1712   return res;
    1713 }
    1714 
    17151728/*2
    17161729* exact division a/b
    … …  
    17701783}
    17711784
     1785#define X_MAS
     1786#ifdef X_MAS
     1787
    17721788/*
    17731789*  returns the part of (a*b)/exp(lead(c)) with nonegative exponents
    … …  
    17761792{
    17771793  poly pa, e, res, r;
    1778   BOOLEAN lead;
    1779 
    1780   if (smSmaller(a, b))
     1794  BOOLEAN lead;\
     1795  int la, lb, lr;
     1796
     1797  if ((c == NULL) || pLmIsConstantComp(c))
     1798  {
     1799    return ppMult_qq(a, b);
     1800  }
     1801
     1802  pqLength(a, b, la, lb, SM_MIN_LENGTH_BUCKET);
     1803
     1804  // we iter over b, so, make sure b is the shortest
     1805  // such that we minimize our iterations
     1806  if (lb > la)
    17811807  {
    17821808    r = a;
    17831809    a = b;
    17841810    b = r;
    1785   }
    1786   if ((c == NULL) || pLmIsConstantComp(c))
    1787   {
    1788     if (pNext(b) == NULL)
    1789     {
    1790       if (pLmIsConstantComp(b))
    1791         return ppMult_nn(a, pGetCoeff(b));
    1792       else
    1793         return ppMult_mm(a, b);
    1794     }
    1795     pa = res = ppMult_mm(a, b);
    1796     pIter(b);
    1797     do
    1798     {
    1799       r = ppMult_mm(a, b);
    1800       smCombineChain(&pa, r);
    1801       pIter(b);
    1802     } while (b != NULL);
    1803     return res;
     1811    lr = la;
     1812    la = lb;
     1813    lb = lr;
    18041814  }
    18051815  res = NULL;
    … …  
    18121822    {
    18131823      lead = pLmDivisibleByNoComp(e, a);
    1814       r = smSelectCopy(a, e);
    1815       smExpMultDiv(r, b, c);
     1824      r = smSelectCopy_ExpMultDiv(a, e, b, c);
    18161825    }
    18171826    else
    … …  
    18421851    }
    18431852  }
    1844   do
     1853 
     1854  if (!TEST_OPT_NOT_BUCKETS && lb >= SM_MIN_LENGTH_BUCKET)
     1855  {
     1856    // use buckets if minimum length is smaller than threshold
     1857    spolyrec rp;
     1858    poly append;
     1859    // find the last monomial before pa
     1860    if (res == pa)
     1861    {
     1862      append = &rp;
     1863      pNext(append) = res;
     1864    }
     1865    else
     1866    {
     1867      append = res;
     1868      while (pNext(append) != pa)
     1869      {
     1870        assume(pNext(append) != NULL);
     1871        pIter(append);
     1872      }
     1873    }
     1874    kBucket_pt bucket = kBucketCreate(currRing);
     1875    kBucketInit(bucket, pNext(append), 0);
     1876    do
     1877    {
     1878      pSetCoeff0(e,pGetCoeff(b));
     1879      if (smIsNegQuot(e, b, c))
     1880      {
     1881        lr = la;
     1882        r = pp_Mult_Coeff_mm_DivSelect_MultDiv(a, lr, e, b, c);
     1883        if (pLmDivisibleByNoComp(e, a))
     1884          append = kBucket_ExtractLarger_Add_q(bucket, append, r, &lr);
     1885        else
     1886          kBucket_Add_q(bucket, r, &lr);
     1887      }
     1888      else
     1889      {
     1890        r = ppMult_mm(a, e);
     1891        append = kBucket_ExtractLarger_Add_q(bucket, append, r, &la);
     1892      }
     1893      pIter(b);
     1894    } while (b != NULL);
     1895    pNext(append) = kBucketClear(bucket);
     1896    kBucketDestroy(&bucket);
     1897  }
     1898  else
     1899  {
     1900    // use old sm stuff
     1901    do
     1902    {
     1903      pSetCoeff0(e,pGetCoeff(b));
     1904      if (smIsNegQuot(e, b, c))
     1905      {
     1906        r = smSelectCopy_ExpMultDiv(a, e, b, c);
     1907        if (pLmDivisibleByNoComp(e, a))
     1908          smCombineChain(&pa, r);
     1909        else
     1910          pa = pAdd(pa,r);
     1911      }
     1912      else
     1913      {
     1914        r = ppMult_mm(a, e);
     1915        smCombineChain(&pa, r);
     1916      }
     1917      pIter(b);
     1918    } while (b != NULL);
     1919  }
     1920  pLmFree(e);
     1921  return res;
     1922}
     1923
     1924#else
     1925
     1926/*
     1927*  returns the part of (a*b)/exp(lead(c)) with nonegative exponents
     1928*/
     1929poly smMultDiv(poly a, poly b, const poly c)
     1930{
     1931  poly pa, e, res, r;
     1932  BOOLEAN lead;
     1933
     1934  if ((c == NULL) || pLmIsConstantComp(c))
     1935  {
     1936    return ppMult_qq(a, b);
     1937  }
     1938  if (smSmaller(a, b))
     1939  {
     1940    r = a;
     1941    a = b;
     1942    b = r;
     1943  }
     1944  res = NULL;
     1945  e = pInit();
     1946  lead = FALSE;
     1947  while (!lead)
    18451948  {
    18461949    pSetCoeff0(e,pGetCoeff(b));
    18471950    if (smIsNegQuot(e, b, c))
    18481951    {
    1849       r = smSelectCopy(a, e);
    1850       smExpMultDiv(r, b, c);
     1952      lead = pLmDivisibleByNoComp(e, a);
     1953      r = smSelectCopy_ExpMultDiv(a, e, b, c);
     1954    }
     1955    else
     1956    {
     1957      lead = TRUE;
     1958      r = ppMult_mm(a, e);
     1959    }
     1960    if (lead)
     1961    {
     1962      if (res != NULL)
     1963      {
     1964        smFindRef(&pa, &res, r);
     1965        if (pa == NULL)
     1966          lead = FALSE;
     1967      }
     1968      else
     1969      {
     1970        pa = res = r;
     1971      }
     1972    }
     1973    else
     1974      res = pAdd(res, r);
     1975    pIter(b);
     1976    if (b == NULL)
     1977    {
     1978      pLmFree(e);
     1979      return res;
     1980    }
     1981  }
     1982  do
     1983  {
     1984    pSetCoeff0(e,pGetCoeff(b));
     1985    if (smIsNegQuot(e, b, c))
     1986    {
     1987      r = smSelectCopy_ExpMultDiv(a, e, b, c);
    18511988      if (pLmDivisibleByNoComp(e, a))
    18521989        smCombineChain(&pa, r);
    … …  
    18642001  return res;
    18652002}
    1866 
     2003#endif
    18672004/*n
    18682005* exact division a/b
    … …  
    18802017}
    18812018
     2019
    18822020/* ------------ internals arithmetic ------------- */
    18832021static void smExactPolyDiv(poly a, poly b)
    … …  
    18872025  poly h;
    18882026  number y, yn;
    1889 
    1890   do
    1891   {
    1892     y = nDiv(pGetCoeff(a), x);
    1893     nNormalize(y);
    1894     pSetCoeff(a,y);
    1895     yn = nNeg(nCopy(y));
    1896     pSetCoeff0(e,yn);
    1897     if (smIsNegQuot(e, a, b))
    1898     {
    1899       h = smSelectCopy(tail, e);
    1900       smExpMultDiv(h, a, b);
    1901     }
    1902     else
    1903       h = ppMult_mm(tail, e);
    1904     nDelete(&yn);
    1905     a = pNext(a) = pAdd(pNext(a), h);
    1906   } while (a!=NULL);
     2027  int lt = pLength(tail);
     2028
     2029  if (lt + 1 >= SM_MIN_LENGTH_BUCKET &&  !TEST_OPT_NOT_BUCKETS)
     2030  {
     2031    kBucket_pt bucket = kBucketCreate(currRing);
     2032    kBucketInit(bucket, pNext(a), 0);
     2033    int lh = 0;
     2034    do
     2035    {
     2036      y = nDiv(pGetCoeff(a), x);
     2037      nNormalize(y);
     2038      pSetCoeff(a,y);
     2039      yn = nNeg(nCopy(y));
     2040      pSetCoeff0(e,yn);
     2041      lh = lt;
     2042      if (smIsNegQuot(e, a, b))
     2043      {
     2044        h = pp_Mult_Coeff_mm_DivSelect_MultDiv(tail, lh, e, a, b);
     2045      }
     2046      else
     2047        h = ppMult_mm(tail, e);
     2048      nDelete(&yn);
     2049      kBucket_Add_q(bucket, h, &lh);
     2050   
     2051      a = pNext(a) = kBucketExtractLm(bucket);
     2052    } while (a!=NULL);
     2053    kBucketDestroy(&bucket);
     2054  }
     2055  else
     2056  {
     2057    do
     2058    {
     2059      y = nDiv(pGetCoeff(a), x);
     2060      nNormalize(y);
     2061      pSetCoeff(a,y);
     2062      yn = nNeg(nCopy(y));
     2063      pSetCoeff0(e,yn);
     2064      if (smIsNegQuot(e, a, b))
     2065        h = smSelectCopy_ExpMultDiv(tail, e, a, b);
     2066      else
     2067        h = ppMult_mm(tail, e);
     2068      nDelete(&yn);
     2069      a = pNext(a) = pAdd(pNext(a), h);
     2070    } while (a!=NULL);
     2071  }
    19072072  pLmFree(e);
    19082073}
    … …  
    19402105  pLmTest(b);
    19412106  pLmTest(c);
     2107  poly bc = p_New(currRing);
     2108 
     2109  p_ExpVectorDiff(bc, b, c, currRing);
     2110 
    19422111  while(t!=NULL)
    19432112  {
    1944     pExpVectorAddSub(t, b, c);
     2113    pExpVectorAdd(t, bc);
    19452114    pIter(t);
    19462115  }
    1947 }
     2116  p_LmFree(bc, currRing);
     2117}
     2118
    19482119
    19492120static void smPolyDivN(poly a, const number x)
    … …  
    20162187  *px = pa;
    20172188}
     2189
    20182190
    20192191static void smFindRef(poly *ref, poly *px, poly r)
Note: See TracChangeset for help on using the changeset viewer.