Changeset b2e4bd in git


Ignore:
Timestamp:
Jun 16, 2011, 12:01:59 PM (12 years ago)
Author:
Stefan Steidel <steidel@…>
Branches:
(u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'a800fe4b3e9d37a38c5a10cc0ae9dfa0c15a4ee6')
Children:
ccd1b0f4ccf9a7be41b783b8f77aef7cc476ffe8
Parents:
9189e93eec836e1bd1e7b9961d718cb4ecce03c9
Message:
Some parallelization features and required input for procedure assPrimes changed. New procedure pTestPoly included.

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

Legend:

Unmodified
Added
Removed

  • Singular/LIB/assprimeszerodim.lib

    r9189e93 rb2e4bd  
    3131         parallel computing)
    3232ASSUME:  I is zero-dimensional in Q[variables]
    33 NOTE:    Parallelization is just applicable using 32-bit Singular version since
    34          MP-links are not compatible with 64-bit Singular version.
    3533RETURN:  the radical of I
    3634EXAMPLE: example zeroRadical; shows an example
    … …  
    6664   if(size(#) > 0)
    6765   {
    68       int n = #[1];
    69       if((n > 1) && (1 - system("with","MP")))
    70       {
    71      "========================================================================";
    72      "There is no MP available on your system. Since this is necessary to     ";
    73      "parallelize the algorithm, the computation will be done without forking.";
    74      "========================================================================";
    75          n = 1;
     66      int n1 = #[1];
     67      if(n1 >= 10)
     68      {
     69         int n2 = n1 + 1;
     70         int n3 = n1;
     71      }
     72      else
     73      {
     74         int n2 = 10;
     75         int n3 = 10;
    7676      }
    7777   }
    7878   else
    7979   {
    80       int n = 1;
     80      int n1 = 1;
     81      int n2 = 10;
     82      int n3 = 10;
    8183   }
    8284
    8385//--------------------  Initialize the list of primes  -------------------------
    84    intvec L = primeList(I,10);
     86   intvec L = primeList(I,n2);
    8587   L[5] = prime(random(100000000,536870912));
    8688
    87    if(n > 1)
    88    {
    89 
    90 //-----  Create n links l(1),...,l(n), open all of them and compute  -----------
    91 //-----  polynomial F for the primes L[2],...,L[n + 1].              -----------
    92 
    93       for(i = 1; i <= n; i++)
    94       {
    95          link l(i) = "MPtcp:fork";
     89   if(n1 > 1)
     90   {
     91
     92//-----  Create n links l(1),...,l(n1), open all of them and compute  ----------
     93//-----  polynomial F for the primes L[2],...,L[n1 + 1].              ----------
     94
     95      for(i = 1; i <= n1; i++)
     96      {
     97         //link l(i) = "MPtcp:fork";
     98         link l(i) = "ssi:fork";
    9699         open(l(i));
    97100         write(l(i), quote(zeroRadP(eval(I), eval(L[i + 1]))));
    … …  
    107110      index++;
    108111
    109       j = j + n + 1;
     112      j = j + n1 + 1;
    110113   }
    111114
    … …  
    114117   while(!crit)
    115118   {
    116       if(n > 1)
     119      if(n1 > 1)
    117120      {
    118121         while(j <= size(L) + 1)
    119122         {
    120             for(i = 1; i <= n; i++)
     123            for(i = 1; i <= n1; i++)
    121124            {
    122125               if(status(l(i), "read", "ready"))
    … …  
    141144            }
    142145            //--- k describes the number of closed links ---
    143             if(k == n)
     146            if(k == n1)
    144147            {
    145148               j++;
    … …  
    151154      else
    152155      {
    153          while(j<=size(L))
     156         while(j <= size(L))
    154157         {
    155158            P = zeroRadP(I, L[j]);
    … …  
    164167      G = chinrem(CO1,CO2);
    165168      N = CO2[1];
    166       for(i = 2; i <= size(CO2); i++){N = N*CO2[i];}
     169      for(i = 2; i <= size(CO2); i++){ N = N*CO2[i]; }
    167170      F = farey(G,N);
    168171
    … …  
    180183
    181184         j = size(L) + 1;
    182          L = primeList(I,10,L);
    183 
    184          if(n > 1)
     185         L = primeList(I,n3,L);
     186
     187         if(n1 > 1)
    185188         {
    186             for(i = 1; i <= n; i++)
     189            for(i = 1; i <= n1; i++)
    187190            {
    188191               open(l(i));
    189192               write(l(i), quote(zeroRadP(eval(I), eval(L[j+i-1]))));
    190193            }
    191             j = j + n;
     194            j = j + n1;
    192195            k = 0;
    193196         }
    … …  
    204207
    205208   if(k == 0) { return(I); }
    206    else { return(modStd(I + F)); }
     209   else { return(modStd(I + F, n1)); }
    207210}
    208211
    … …  
    210213
    211214proc assPrimes(list #)
    212 "USAGE:  assPrimes(I,[a],[n]); I ideal or module,
     215"USAGE:  assPrimes(I,[n],[a]); I ideal or module,
    213216         optional: int n: number of processors (for parallel computing), int a:
    214217         - a = 1: method of Eisenbud/Hunecke/Vasconcelos
    … …  
    217220         assPrimes(I) chooses n = a = 1
    218221ASSUME:  I is zero-dimensional over Q[variables]
    219 NOTE:    Parallelization is just applicable using 32-bit Singular version since
    220          MP-links are not compatible with 64-bit Singular version.
    221222RETURN:  a list Re of associated primes of I:
    222223EXAMPLE: example assPrimes; shows an example
    … …  
    239240      if(size(#) == 2)
    240241      {
    241          int alg = #[2];
    242          int n = 1;
     242         int n1 = #[2];
     243         int alg = 1;
     244         if(n1 >= 10)
     245         {
     246            int n2 = n1 + 1;
     247            int n3 = n1;
     248         }
     249         else
     250         {
     251            int n2 = 10;
     252            int n3 = 10;
     253         }
    243254      }
    244255      if(size(#) == 3)
    245256      {
    246          int alg = #[2];
    247          int n = #[3];
    248          if((n > 1) && (1 - system("with","MP")))
     257         int n1 = #[2];
     258         int alg = #[3];
     259         if(n1 >= 10)
     260         {
     261            int n2 = n1 + 1;
     262            int n3 = n1;
     263         }
     264         else
    249265         {
    250      "========================================================================";
    251      "There is no MP available on your system. Since this is necessary to     ";
    252      "parallelize the algorithm, the computation will be done without forking.";
    253      "========================================================================";
    254             n = 1;
     266            int n2 = 10;
     267            int n3 = 10;
    255268         }
    256269      }
    … …  
    258271   else
    259272   {
     273      int n1 = 1;
    260274      int alg = 1;
    261       int n = 1;
     275      int n2 = 10;
     276      int n3 = 10;
     277   }
     278   
     279   if(printlevel >= 10)
     280   {
     281      "n1 = "+string(n1)+", n2 = "+string(n2)+", n3 = "+string(n3);
    262282   }
    263283
    … …  
    265285   int RT = rtimer;
    266286   int TT;
     287   int t_sleep;
    267288
    268289   def SPR = basering;
     290   map phi;
     291   list H = ideal(0);
     292   ideal F;
     293   poly F1;
     294   
    269295   if(printlevel >= 10) { "========== Start modStd =========="; }
    270    I = modStd(I,n);
     296   I = modStd(I,n1);
    271297   if(printlevel >= 10) { "=========== End modStd ==========="; }
    272298   if(printlevel >= 9) { "modStd takes "+string(rtimer-RT)+" seconds."; }
    273299   int d = vdim(I);
    274300   if(d == -1) { ERROR("Ideal is not zero-dimensional."); }
    275    if(homog(I) == 1){ return(list(maxideal(1))); }
     301   if(homog(I) == 1) { return(list(maxideal(1))); }
    276302   poly f = findGen(I);
    277303   if(printlevel >= 9) { "Coordinate change: "+string(f); }
    … …  
    293319            {
    294320               TT = timer;
    295                I = zeroR(I,n);
     321               I = zeroR(I,n1);
    296322               if(printlevel >= 9)
    297323               {
    298                   "zeroR(I,n) takes "+string(timer-TT)+" seconds.";
     324                  "zeroR(I,n1) takes "+string(timer-TT)+" seconds.";
    299325               }
    300326               TT = timer;
    … …  
    311337         {
    312338            TT = timer;
    313             I = zeroR(I,n);
     339            I = zeroR(I,n1);
    314340            if(printlevel >= 9)
    315341            {
    316                "zeroR(I,n) takes "+string(timer-TT)+" seconds.";
     342               "zeroR(I,n1) takes "+string(timer-TT)+" seconds.";
    317343            }
    318344            TT = timer;
    … …  
    340366
    341367//--------------------  Initialize the list of primes  -------------------------
    342    intvec L = primeList(I,10);
     368   intvec L = primeList(I,n2);
    343369   L[5] = prime(random(1000000000,2134567879));
    344370
    … …  
    358384//-----  main standard basis computations in positive characteristic        ----
    359385
    360    if(n > 1)
     386   if(n1 > 1)
    361387   {
    362388
    363389//-----  Create n1 links l(1),...,l(n1), open all of them and compute  ---------
    364 //-----  standard basis for the primes L[2],...,L[n + 1].              ---------
    365 
    366       for(i = 1; i <= n; i++)
    367       {
    368          link l(i) = "MPtcp:fork";
     390//-----  standard basis for the primes L[2],...,L[n1 + 1].             ---------
     391
     392      for(i = 1; i <= n1; i++)
     393      {
     394         //link l(i) = "MPtcp:fork";
     395         link l(i) = "ssi:fork";
    369396         open(l(i));
    370397         write(l(i), quote(modpSpecialAlgDep(eval(ringL), eval(L[i + 1]),
    … …  
    381408      index++;
    382409
    383       j = j + n + 1;
     410      j = j + n1 + 1;
    384411   }
    385412
    … …  
    388415   int tt = timer;
    389416   int rt = rtimer;
    390 
     417   
    391418   while(1)
    392419   {
    393420      tt = timer;
    394421      rt = rtimer;
    395 
    396       if(n > 1)
     422     
     423      if(printlevel >= 9) { "size(L) = "+string(size(L)); }
     424     
     425      if(n1 > 1)
    397426      {
    398427         while(j <= size(L) + 1)
    399428         {
    400             for(i = 1; i <= n; i++)
     429            for(i = 1; i <= n1; i++)
    401430            {
    402431               //--- ask if link l(i) is ready otherwise sleep for t seconds ---
    … …  
    424453            }
    425454            //--- k describes the number of closed links ---
    426             if(k == n)
     455            if(k == n1)
    427456            {
    428457               j++;
    … …  
    433462      else
    434463      {
    435          while(j<=size(L))
     464         while(j <= size(L))
    436465         {
    437466            P = modpSpecialAlgDep(ringL, L[j], alg);
    … …  
    442471         }
    443472      }
     473     
    444474      if(printlevel >= 9)
    445475      {
    … …  
    450480      }
    451481
    452       // insert deleteUnluckyPrimes
     482//-------------------  Lift results to basering via farey ----------------------
     483
     484      tt = timer;
    453485      G = chinrem(CO1,CO2);
    454486      N = CO2[1];
    455       for(j = 2; j <= size(CO2); j++){N = N*CO2[j];}
     487      for(j = 2; j <= size(CO2); j++){ N = N*CO2[j]; }
    456488      F = farey(G,N);
    457       if(F[1]-testF[1]==0)
    458       {
    459          if(printlevel >= 9) { "size(L) = "+string(size(L)); }
    460 
     489      if(printlevel >= 10) { "Lifting-process takes "+string(timer - tt)
     490                             +" seconds"; }
     491
     492      if(pTestPoly(F[1], ringL, alg, L))                             
     493      {
    461494         F = cleardenom(F[1]);
    462495
    … …  
    478511            {
    479512               setring SPR;
    480                map phi = rHelp,var(nvars(SPR));
    481                list H = phi(H);
     513               phi = rHelp,var(nvars(SPR));
     514               H = phi(H);
     515               
    482516               if(printlevel >= 9)
    483517               {
    … …  
    485519                  "CPU-time without test is "+string(timer - T)+" seconds.";
    486520               }
     521               
    487522               T = timer;
    488523               RT = rtimer;
    489524
    490                ideal F = phi(F);
    491                poly F1;
    492 
    493                if(n > 1)
     525               F = phi(F);
     526
     527               if(n1 > 1)
    494528               {
    495529                  open(l(1));
    496530                  write(l(1), quote(quickSubst(eval(F[1]), eval(f), eval(I))));
    497                   int t_sleep = timer;
     531                  t_sleep = timer;
    498532               }
    499533               else
    … …  
    511545                  }
    512546
    513                   if(n > 1)
     547                  if(n1 > 1)
    514548                  {
    515549                     t_sleep = timer - t_sleep;
    … …  
    546580
    547581      int_break = 0;
     582      setring rHelp;
    548583      testF = F;
    549584      CO1 = G;
    … …  
    553588      j = size(L) + 1;
    554589
    555       setring(SPR);
    556       L = primeList(I,10,L);
     590      setring SPR;
     591      L = primeList(I,n3,L);
    557592      setring rHelp;
    558593
    559       if(n > 1)
    560       {
    561          for(i = 1; i <= n; i++)
     594      if(n1 > 1)
     595      {
     596         for(i = 1; i <= n1; i++)
    562597         {
    563598            open(l(i));
    … …  
    565600                                                             eval(alg))));
    566601         }
    567          j = j + n;
     602         j = j + n1;
    568603         k = 0;
    569604      }
    … …  
    572607example
    573608{ "EXAMPLE:";  echo = 2;
    574    ring R=0,(a,b,c,d,e,f),dp;
    575    ideal I=
     609   ring R = 0,(a,b,c,d,e,f),dp;
     610   ideal I =
    576611   2fb+2ec+d2+a2+a,
    577612   2fc+2ed+2ba+b,
    … …  
    880915////////////////////////////////////////////////////////////////////////////////
    881916
     917static proc pTestPoly(poly testF, list ringL, int alg, list L)
     918{
     919   int i,j,p;
     920   def R0 = basering;
     921   
     922   while(!i)
     923   {
     924      i = 1;
     925      p = prime(random(1000000000,2134567879));
     926      for(j = 1; j <= size(L); j++)
     927      {
     928         if(p == L[j]) { i = 0; break; }
     929      }
     930   }
     931   
     932   ringL[1] = p;
     933   def @R = ring(ringL);
     934   setring @R;
     935   poly f = fetch(SPR,f_for_fork);
     936   ideal I = fetch(SPR,I_for_fork);
     937   if(alg == 1) { list M = specialAlgDepEHV(f,I); }
     938   if(alg == 2) { list M = specialAlgDepGTZ(f,I); }
     939   if(alg == 3) { list M = specialAlgDepMonico(f,I); }
     940   def @S = M[1];
     941   setring @S;
     942   poly testF = fetch(R0,testF);
     943   int k = (testF == F);
     944   
     945   setring R0;
     946   return(k);
     947}
     948
     949////////////////////////////////////////////////////////////////////////////////
     950
    882951/*
    883952Examples:
Note: See TracChangeset for help on using the changeset viewer.