Changeset 7f30e2 in git


Ignore:
Timestamp:
Nov 25, 2011, 5:22:38 PM (12 years ago)
Author:
Andreas Steenpass <steenpass@…>
Branches:
(u'spielwiese', '8d54773d6c9e2f1d2593a28bc68b7eeab54ed529')
Children:
f1cfef36a45967a6dc88672e6efc050bf9c592b1
Parents:
4093f96acfcc6419c56ee983c447099c3b87a36e
git-author:
Andreas Steenpass <steenpass@mathematik.uni-kl.de>2011-11-25 17:22:38+01:00
git-committer:
Andreas Steenpass <steenpass@mathematik.uni-kl.de>2012-08-02 18:26:02+02:00
Message:
improve modstd.lib and parallel.lib

modstd.lib:
  - change heuristics for variants
  - parallelize lifting
  - parallelize finding new primes
parallel.lib:
  fix behaviour for several arguments applied to a kernel command
Location:
Singular/LIB
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • Singular/LIB/modstd.lib

    r4093f9 r7f30e2  
    11////////////////////////////////////////////////////////////////////////////////
    2 version="$Id$";
     2version="$Id: modstd.lib 14375 2011-08-23 09:29:47Z steidel $";
    33category = "Commutative Algebra";
    44info="
    … …  
    88@*        G. Pfister      pfister@mathematik.uni-kl.de
    99@*        H. Schoenemann  hannes@mathematik.uni-kl.de
     10@*        A. Steenpass    steenpass@mathematik.uni-kl.de
    1011@*        S. Steidel      steidel@mathematik.uni-kl.de
    1112
    … …  
    2627LIB "poly.lib";
    2728LIB "ring.lib";
     29LIB "parallel.lib";
    2830
    2931////////////////////////////////////////////////////////////////////////////////
    … …  
    243245      if(printlevel >= 11)
    244246      {
    245          "isIncluded(K,J) takes "+string(timer - t)+" seconds";
     247         "isIncluded(J,K) takes "+string(timer - t)+" seconds";
    246248         "j = "+string(j);
    247249      }
    … …  
    373375////////////////////////////////////////////////////////////////////////////////
    374376
    375 proc primeTest(ideal I, bigint p)
     377proc primeTest(def II, bigint p)
    376378{
     379   if(typeof(II) == "string")
     380   {
     381      execute("ideal I = "+II+";");
     382   }
     383   else
     384   {
     385      ideal I = II;
     386   }
     387
    377388   int i,j;
     389   poly f;
     390   number cnt;
    378391   for(i = 1; i <= size(I); i++)
    379392   {
    380       for(j = 1; j <= size(I[i]); j++)
    381       {
    382          if((numerator(leadcoef(I[i][j])) mod p) == 0) { return(0); }
    383          if((denominator(leadcoef(I[i][j])) mod p) == 0) { return(0); }
     393      f = cleardenom(I[i]);
     394      if(f == 0) { return(0); }
     395      cnt = leadcoef(I[i])/leadcoef(f);
     396      if((numerator(cnt) mod p) == 0) { return(0); }
     397      if((denominator(cnt) mod p) == 0) { return(0); }
     398      for(j = size(f); j > 0; j--)
     399      {
     400         if((leadcoef(f[j]) mod p) == 0) { return(0); }
    384401      }
    385402   }
    … …  
    390407
    391408proc primeList(ideal I, int n, list #)
    392 "USAGE:  primeList(I,n); ( resp. primeList(I,n,L); ) I ideal, n integer
    393 RETURN:  the intvec of n greatest primes  <= 2147483647 (resp. n greatest primes
     409"USAGE:  primeList(I,n[,ncores]); ( resp. primeList(I,n[,L,ncores]); ) I ideal,
     410         n integer
     411RETURN:  the intvec of n greatest primes <= 2147483647 (resp. n greatest primes
    394412         < L[size(L)] union with L) such that none of these primes divides any
    395413         coefficient occuring in I
     414NOTE:    The number of cores to use can be defined by ncores, default is 1.
    396415EXAMPLE: example primList; shows an example
    397416"
    … …  
    399418   intvec L;
    400419   int i,p;
     420   int ncores = 1;
     421   
     422//-----------------  Initialize optional parameter ncores  ---------------------
     423   if(size(#) > 0)
     424   {
     425      if(size(#) == 1)
     426      {
     427         if(typeof(#[1]) == "int")
     428         {
     429            ncores = #[1];
     430            # = list();
     431         }
     432      }
     433      else
     434      {
     435         ncores = #[2];
     436      }
     437   }
     438   
    401439   if(size(#) == 0)
    402440   {
    … …  
    421459   }
    422460   if(p == 2) { ERROR("no more primes"); }
    423    for(i = 2; i <= n; i++)
    424    {
    425       p = prime(p-1);
    426       while(!primeTest(I,p))
     461   if(ncores == 1)
     462   {
     463      for(i = 2; i <= n; i++)
    427464      {
    428465         p = prime(p-1);
    429          if(p == 2) { ERROR("no more primes"); }
    430       }
    431       L[size(L)+1] = p;
     466         while(!primeTest(I,p))
     467         {
     468            p = prime(p-1);
     469            if(p == 2) { ERROR("no more primes"); }
     470         }
     471         L[size(L)+1] = p;
     472      }
     473   }
     474   else
     475   {
     476      int neededSize = size(L)+n-1;;
     477      list parallelResults;
     478      list arguments;
     479      int neededPrimes = neededSize-size(L);
     480      while(neededPrimes > 0)
     481      {
     482         arguments = list();
     483         for(i = ((neededPrimes div ncores)+1-(neededPrimes%ncores == 0))
     484            *ncores; i > 0; i--)
     485         {
     486            p = prime(p-1);
     487            if(p == 2) { ERROR("no more primes"); }
     488            arguments[i] = list("I", p);
     489         }
     490         parallelResults = parallelWaitAll("primeTest", arguments,
     491            list(list(list(ncores))));
     492         for(i = size(arguments); i > 0; i--)
     493         {
     494            if(parallelResults[i])
     495            {
     496               L[size(L)+1] = arguments[i][2];
     497            }
     498         }
     499         neededPrimes = neededSize-size(L);
     500      }
     501      if(size(L) > neededSize)
     502      {
     503         L = L[1..neededSize];
     504      }
    432505   }
    433506   return(L);
    … …  
    470543   ideal I = imap(R,I);
    471544
     545   intvec opt = option(get);
    472546   option(none);
    473547   option(prompt);
    … …  
    504578
    505579   setring R;
     580   option(set, opt);
    506581   return(imap(newR,sI),imap(newR,T));
    507582}
    … …  
    537612         The standard basis computation modulo p does also vary depending on the
    538613         integer variant, namely
    539 @*       - variant = 1/4: std(.,#[1]) resp. groebner,
    540 @*       - variant = 2/5: groebner,
    541 @*       - variant = 3/6: homog. - std(.,#[1]) resp. groebner - dehomog..
     614@*       - variant = 1: std(.,#[1]) resp. groebner,
     615@*       - variant = 2: groebner,
     616@*       - variant = 3: homog. - std(.,#[1]) resp. groebner - dehomog.,
     617@*       - variant = 4: fglm.
    542618EXAMPLE: example modpStd; shows an example
    543619"
    … …  
    552628   option(redSB);
    553629
    554    int t = timer;
    555    if((variant == 1) || (variant == 4))
     630   if(variant == 1)
    556631   {
    557632      if(size(#) > 0)
    … …  
    565640   }
    566641
    567    if((variant == 2) || (variant == 5))
     642   if(variant == 2)
    568643   {
    569644      i = groebner(i);
    570645   }
    571646
    572    if((variant == 3) || (variant == 6))
     647   if(variant == 3)
    573648   {
    574649      list rl = ringlist(@r);
    … …  
    607682      }
    608683
    609       t = timer;
    610684      i = subst(i, homvar, 1);
    611685      i = simplify(i, 34);
    … …  
    615689      i = interred(i);
    616690      kill HomR;
     691   }
     692   
     693   if(variant == 4)
     694   {
     695      def R1 = changeord("dp");
     696      setring R1;
     697      ideal i = fetch(@r,i);
     698      i = std(i);
     699      setring @r;
     700      i = fglm(R1,i);
    617701   }
    618702
    … …  
    646730RETURN:  a standard basis of I if no warning appears;
    647731NOTE:    The procedure computes a standard basis of I (over the rational
    648          numbers) by using  modular methods. If a warning appears then the
    649          result is a standard basis containing I and with high probability
    650          a standard basis of I.
     732         numbers) by using  modular methods.
    651733         By default the procedure computes a standard basis of I for sure, but
    652734         if the optional parameter #[2] = 0, it computes a standard basis of I
    653          with high probability and a warning appears at the end.
     735         with high probability.
    654736         The procedure distinguishes between different variants for the standard
    655737         basis computation in positive characteristic depending on the ordering
    656738         of the basering, the parameter #[2] and if the ideal I is homogeneous.
    657 @*       - variant = 1, if I is homogeneous and exactness = 0,
    658 @*       - variant = 2, if I is not homogeneous, 1-block-ordering and
    659                         exactness = 0,
     739@*       - variant = 1, if I is homogeneous,
     740@*       - variant = 2, if I is not homogeneous, 1-block-ordering,
    660741@*       - variant = 3, if I is not homogeneous, complicated ordering (lp or
    661                         > 1 block) and exactness = 0,
    662 @*       - variant = 4, if I is homogeneous and exactness = 1,
    663 @*       - variant = 5, if I is not homogeneous, 1-block-ordering and
    664                         exactness = 1,
    665 @*       - variant = 6, if I is not homogeneous, complicated ordering (lp or
    666                         > 1 block) and exactness = 1.
     742                        > 1 block),
     743@*       - variant = 4, if I is not homogeneous, ordering lp, dim(I) = 0.
    667744EXAMPLE: example modStd; shows an example
    668745"
    … …  
    752829   if(printlevel >= 10)
    753830   {
    754       "n1 = "+string(n1)+", n2 = "+string(n2)+", n3 = "+string(n3);
     831      "n1 = "+string(n1)+", n2 = "+string(n2)+", n3 = "+string(n3)
     832       +", exactness = "+string(exactness);
    755833   }
    756834
    … …  
    761839
    762840//--------------------  Initialize the list of primes  -------------------------
    763    intvec L = primeList(I,n2);
     841   int tt = timer;
     842   int rt = rtimer;
     843   intvec L = primeList(I,n2,n1);
     844   if(printlevel >= 10)
     845   {
     846      "CPU-time for primeList: "+string(timer-tt)+" seconds.";
     847      "Real-time for primeList: "+string(rtimer-rt)+" seconds.";
     848   }
    764849   L[5] = prime(random(an,en));
    765850
    … …  
    768853   int h = homog(I);
    769854
    770    int tt = timer;
    771    int rt = rtimer;
     855   tt = timer;
     856   rt = rtimer;
    772857
    773858   if(!mixedTest())
    … …  
    775860      if(h)
    776861      {
    777          if(exactness == 0)
    778          {
    779             variant = 1;
    780             if(printlevel >= 10) { "variant = 1"; }
    781          }
    782          if(exactness == 1)
    783          {
    784             variant = 4;
    785             if(printlevel >= 10) { "variant = 4"; }
    786          }
     862         variant = 1;
     863         if(printlevel >= 10) { "variant = 1"; }
     864         
    787865         rl[1] = L[5];
    788866         def @r = ring(rl);
    … …  
    802880         if((find(ordstr_R0, "M") > 0) || (find(ordstr_R0, "a") > 0) || neg)
    803881         {
    804             if(exactness == 0)
    805             {
    806                variant = 2;
    807                if(printlevel >= 10) { "variant = 2"; }
    808             }
    809             if(exactness == 1)
    810             {
    811                variant = 5;
    812                if(printlevel >= 10) { "variant = 5"; }
    813             }
     882            variant = 2;
     883            if(printlevel >= 10) { "variant = 2"; }
    814884         }
    815885         else
    … …  
    827897            if((order == "simple") || (size(rl) > 4))
    828898            {
    829                if(exactness == 0)
     899               variant = 2;
     900               if(printlevel >= 10) { "variant = 2"; }
     901            }
     902            else
     903            {
     904               rl[1] = L[5];
     905               def @r = ring(rl);
     906               setring @r;
     907               
     908               def @s = changeord("dp");
     909               setring @s;
     910               ideal I = std(fetch(R0,I));
     911               if(dim(I) == 0)
    830912               {
    831                   variant = 2;
    832                   if(printlevel >= 10) { "variant = 2"; }
     913                  variant = 4;
     914                  if(printlevel >= 10) { "variant = 4"; }
    833915               }
    834                if(exactness == 1)
    835                {
    836                   variant = 5;
    837                   if(printlevel >= 10) { "variant = 5"; }
    838                }
    839             }
    840             else
    841             {
    842                if(exactness == 0)
     916               else
    843917               {
    844918                  variant = 3;
    845919                  if(printlevel >= 10) { "variant = 3"; }
     920                 
     921                  int nvar@r = nvars(@r);
     922                  intvec w;
     923                  for(i = 1; i <= nvar@r; i++)
     924                  {
     925                     w[i] = deg(var(i));
     926                  }
     927                  w[nvar@r + 1] = 1;
     928
     929                  list hiRi = hilbRing(fetch(R0,I),w);
     930                  intvec W = hiRi[2];
     931                  @s = hiRi[1];
     932                  setring @s;
     933
     934                  Id(1) = std(Id(1));
     935                  intvec hi = hilb(Id(1), 1, W);
    846936               }
    847                if(exactness == 1)
    848                {
    849                   variant = 6;
    850                   if(printlevel >= 10) { "variant = 6"; }
    851                }
    852 
    853                rl[1] = L[5];
    854                def @r = ring(rl);
    855                setring @r;
    856                int nvar@r = nvars(@r);
    857                intvec w;
    858                for(i = 1; i <= nvar@r; i++)
    859                {
    860                   w[i] = deg(var(i));
    861                }
    862                w[nvar@r + 1] = 1;
    863 
    864                list hiRi = hilbRing(fetch(R0,I),w);
    865                intvec W = hiRi[2];
    866                def @s = hiRi[1];
    867                setring @s;
    868 
    869                Id(1) = std(Id(1));
    870                intvec hi = hilb(Id(1), 1, W);
    871 
     937               
    872938               setring R0;
    873939               kill @r,@s;
    … …  
    9711037         link l(i) = "ssi:fork";
    9721038         open(l(i));
    973          if((variant == 1) || (variant == 3) ||
    974             (variant == 4) || (variant == 6))
     1039         if((variant == 1) || (variant == 3))
    9751040         {
    9761041            write(l(i), quote(modpStd(I_for_fork, eval(L[i + 1]),
    9771042                                                  eval(variant), eval(hi))));
    9781043         }
    979          if((variant == 2) || (variant == 5))
     1044         if((variant == 2) || (variant == 4))
    9801045         {
    9811046            write(l(i), quote(modpStd(I_for_fork, eval(L[i + 1]),
    … …  
    9851050
    9861051      int t = timer;
    987       if((variant == 1) || (variant == 3) || (variant == 4) || (variant == 6))
     1052      if((variant == 1) || (variant == 3))
    9881053      {
    9891054         P = modpStd(I_for_fork, L[1], variant, hi);
    9901055      }
    991       if((variant == 2) || (variant == 5))
     1056      if((variant == 2) || (variant == 4))
    9921057      {
    9931058         P = modpStd(I_for_fork, L[1], variant);
    … …  
    10051070//--------------  Main standard basis computations in positive  ----------------
    10061071//----------------------  characteristic start here  ---------------------------
     1072
     1073   list arguments_farey, results_farey;
    10071074
    10081075   while(1)
    … …  
    10291096                  if(j <= size(L))
    10301097                  {
    1031                      if((variant == 1) || (variant == 3) ||
    1032                         (variant == 4) || (variant == 6))
     1098                     if((variant == 1) || (variant == 3))
    10331099                     {
    10341100                        write(l(i), quote(modpStd(I_for_fork, eval(L[j]),
    … …  
    10361102                        j++;
    10371103                     }
    1038                      if((variant == 2) || (variant == 5))
     1104                     if((variant == 2) || (variant == 4))
    10391105                     {
    10401106                        write(l(i), quote(modpStd(I_for_fork,
    … …  
    10621128         while(j <= size(L))
    10631129         {
    1064             if((variant == 1) || (variant == 3) ||
    1065                (variant == 4) || (variant == 6))
     1130            if((variant == 1) || (variant == 3))
    10661131            {
    10671132               P = modpStd(I, L[j], variant, hi);
    10681133            }
    1069             if((variant == 2) || (variant == 5))
     1134            if((variant == 2) || (variant == 4))
    10701135            {
    10711136               P = modpStd(I, L[j], variant);
    … …  
    10951160//-------------------  Lift results to basering via farey  ---------------------
    10961161
    1097       tt = timer;
     1162      tt = timer; rt = rtimer;
    10981163      N = T2[1];
    10991164      for(i = 2; i <= size(T2); i++) { N = N*T2[i]; }
    11001165      H = chinrem(T1,T2);
    1101       J = farey(H,N);
    1102       if(printlevel >= 10) { "Lifting-process takes "+string(timer - tt)
    1103                              +" seconds"; }
     1166      if(n1 == 1)
     1167      {
     1168         J = farey(H,N);
     1169      }
     1170      else
     1171      {
     1172         for(i = ncols(H); i > 0; i--)
     1173         {
     1174            arguments_farey[i] = list(ideal(H[i]), N);
     1175         }
     1176         results_farey = parallelWaitAll("farey", arguments_farey,
     1177                                         list(list(list(n1))));
     1178         for(i = ncols(H); i > 0; i--)
     1179         {
     1180            J[i] = results_farey[i][1];
     1181         }
     1182      }
     1183      if(printlevel >= 10)
     1184      {
     1185         "CPU-time for lifting-process is "+string(timer - tt)+" seconds.";
     1186         "Real-time for lifting-process is "+string(rtimer - rt)+" seconds.";
     1187      }
    11041188
    11051189//----------------  Test if we already have a standard basis of I --------------
    11061190
    11071191      tt = timer; rt = rtimer;
    1108       if((variant == 1) || (variant == 3) || (variant == 4) || (variant == 6))
     1192      if((variant == 1) || (variant == 3))
    11091193      {
    11101194         pTest = pTestSB(I,J,L,variant,hi);
    11111195      }
    1112       if((variant == 2) || (variant == 5))
     1196      if((variant == 2) || (variant == 4))
    11131197      {
    11141198         pTest = pTestSB(I,J,L,variant);
    … …  
    11321216
    11331217         attrib(J,"isSB",1);
    1134          tt = timer; rt = rtimer;
    1135          sizeTest = 1 - isIncluded(I,J,n1);
    1136 
    1137          if(printlevel >= 10)
    1138          {
    1139             "CPU-time for checking if I subset <G> is "
    1140             +string(timer - tt)+" seconds.";
    1141             "Real-time for checking if I subset <G> is "
    1142             +string(rtimer - rt)+" seconds.";
    1143          }
    1144 
    1145          if(sizeTest == 0)
    1146          {
    1147             if((variant == 1) || (variant == 2) || (variant == 3))
    1148             {
    1149                "===========================================================";
    1150                "WARNING: Output might not be a standard basis of the input.";
    1151                "===========================================================";
    1152                option(set, opt);
    1153                if(n1 > 1) { kill I_for_fork; }
    1154                return(J);
    1155             }
    1156             if((variant == 4) || (variant == 5) || (variant == 6))
     1218         
     1219         if(exactness == 0)
     1220         {
     1221            option(set, opt);
     1222            if(n1 > 1) { kill I_for_fork; }
     1223            return(J);
     1224         }
     1225         
     1226         if(exactness == 1)
     1227         {
     1228            tt = timer; rt = rtimer;
     1229            sizeTest = 1 - isIncluded(I,J,n1);
     1230
     1231            if(printlevel >= 10)
     1232            {
     1233               "CPU-time for checking if I subset <G> is "
     1234               +string(timer - tt)+" seconds.";
     1235               "Real-time for checking if I subset <G> is "
     1236               +string(rtimer - rt)+" seconds.";
     1237            }
     1238
     1239            if(sizeTest == 0)
    11571240            {
    11581241               tt = timer; rt = rtimer;
    … …  
    11851268
    11861269      j = size(L) + 1;
    1187       L = primeList(I,n3,L);
     1270      tt = timer; rt = rtimer;
     1271      L = primeList(I,n3,L,n1);
     1272      if(printlevel >= 10)
     1273      {
     1274         "CPU-time for primeList: "+string(timer-tt)+" seconds.";
     1275         "Real-time for primeList: "+string(rtimer-rt)+" seconds.";
     1276      }
    11881277
    11891278      if(n1 > 1)
    … …  
    11921281         {
    11931282            open(l(i));
    1194             if((variant == 1) || (variant == 3) ||
    1195                (variant == 4) || (variant == 6))
     1283            if((variant == 1) || (variant == 3))
    11961284            {
    11971285               write(l(i), quote(modpStd(I_for_fork, eval(L[j+i-1]),
    11981286                                                     eval(variant), eval(hi))));
    11991287            }
    1200             if((variant == 2) || (variant == 5))
     1288            if((variant == 2) || (variant == 4))
    12011289            {
    12021290               write(l(i), quote(modpStd(I_for_fork, eval(L[j+i-1]),
    … …  
    13181406   def S = ring(rl);
    13191407   setring S;
     1408   intvec opt = option(get);
    13201409   option(redSB);
    13211410   module Z,M,Z2;
    … …  
    13621451               if(size(reduce(G3,testG1)) == 0)
    13631452               {
     1453                  option(set, opt);
    13641454                  return(G3);
    13651455               }

  • Singular/LIB/parallel.lib

    r4093f9 r7f30e2  
    368368    }
    369369    write(l(i), quote(execute("result = "+eval(commands[k])
    370                                     +"(currentargs[1..size(currentargs)]);")));
     370      +"("+argsToString("currentargs", size(currentargs))+");")));
    371371    assignment[i] = k;
    372372    k++;
    … …  
    455455      }
    456456      write(l(wait), quote(execute("def result = "+eval(commands[k])
    457                                      +"(currentargs[1..size(currentargs)]);")));
     457        +"("+argsToString("currentargs", size(currentargs))+");")));
    458458      assignment[wait] = k;
    459459      k++;
    … …  
    717717  return(i);
    718718}
     719
     720static proc argsToString(string name, int length)
     721{
     722  string arglist;
     723  if(length > 0) {
     724    arglist = name+"[1]";
     725  }
     726  int i;
     727  for(i = 2; i <= length; i++) {
     728    arglist = arglist+", "+name+"["+string(i)+"]";
     729  }
     730  return(arglist);
     731}
Note: See TracChangeset for help on using the changeset viewer.