source: git/Singular/dyn_modules/singmathic/singmathic.cc @ fcdeea

spielwiese
Last change on this file since fcdeea was fcdeea, checked in by Hans Schoenemann <hannes@…>, 5 years ago
fix: threadsupport
  • Property mode set to 100644
File size: 16.3 KB
Line 
1#include "kernel/mod2.h"
2
3#ifdef HAVE_MATHICGB
4
5#include "kernel/mod2.h"
6
7#include "misc/options.h"
8
9#include "kernel/ideals.h"
10#include "kernel/polys.h"
11
12#include "Singular/ipid.h"
13#include "Singular/feOpt.h"
14#include "Singular/mod_lib.h"
15
16#include <mathicgb.h>
17
18typedef mgb::GroebnerConfiguration::Coefficient Coefficient;
19typedef mgb::GroebnerConfiguration::VarIndex VarIndex;
20typedef mgb::GroebnerConfiguration::Exponent Exponent;
21typedef mgb::GroebnerConfiguration::BaseOrder BaseOrder;
22
23// Constructs a Singular ideal.
24class MathicToSingStream
25{
26public:
27  MathicToSingStream(Coefficient modulus, VarIndex varCount):
28    mModulus(modulus),
29    mVarCount(varCount),
30    mPolyCount(0),
31    mTerm(0),
32    mIdeal(0)
33  {}
34
35  ~MathicToSingStream() {deleteIdeal();}
36
37  // Mathic stream interface
38
39  Coefficient modulus() const {return mModulus;}
40  VarIndex varCount() const {return mModulus;}
41
42  void idealBegin(size_t polyCount)
43  {
44    deleteIdeal();
45    mIdeal = idInit(polyCount);
46    mPolyCount = 0;
47  }
48
49  void appendPolynomialBegin(size_t termCount) {}
50
51  void appendTermBegin(const mgb::GroebnerConfiguration::Component c)
52  {
53    if (mTerm == 0)
54      mTerm = mIdeal->m[mPolyCount] = pInit();
55    else
56      mTerm = mTerm->next = pInit();
57    pSetComp(mTerm,c);
58  }
59
60  void appendExponent(VarIndex index, Exponent exponent)
61  {
62    pSetExp(mTerm, index + 1, exponent);
63  }
64
65  void appendTermDone(Coefficient coefficient)
66  {
67    mTerm->coef = reinterpret_cast<number>(coefficient);
68    pSetm(mTerm);
69  }
70
71  void appendPolynomialDone()
72  {
73    ++mPolyCount;
74    mTerm = 0;
75  }
76
77  void idealDone() {}
78
79
80  // Singular interface
81
82  ::ideal takeIdeal()
83  {
84    ::ideal id = mIdeal;
85    mIdeal = 0;
86    return id;
87  }
88
89private:
90  void deleteIdeal()
91  {
92    if (mIdeal != 0)
93    {
94      idDelete(&mIdeal);
95      mIdeal = 0;
96    }
97  }
98
99  const Coefficient mModulus;
100  const VarIndex mVarCount;
101  size_t mPolyCount;
102  poly mTerm;
103  ::ideal mIdeal;
104};
105
106#include <iostream>
107
108bool setOrder(ring r, mgb::GroebnerConfiguration& conf)
109{
110  const VarIndex varCount = conf.varCount();
111
112  bool didSetComponentBefore = false;
113  mgb::GroebnerConfiguration::BaseOrder baseOrder =
114    mgb::GroebnerConfiguration::RevLexDescendingBaseOrder;
115
116  std::vector<Exponent> gradings;
117  for (int block = 0; r->order[block] != ringorder_no; ++block)
118  {
119    // *** ringorder_no
120
121    const rRingOrder_t type = static_cast<rRingOrder_t>(r->order[block]);
122    if (r->block0[block] < 0 || r->block1[block] < 0)
123    {
124      WerrorS("Unexpected negative block0/block1 in ring.");
125      return false;
126    }
127    const VarIndex block0 = static_cast<VarIndex>(r->block0[block]);
128    const VarIndex block1 = static_cast<VarIndex>(r->block1[block]);
129    const int* const weights = r->wvhdl[block];
130    if (block0 > block1)
131    {
132      WerrorS("Unexpected block0 > block1 in ring.");
133      return false;
134    }
135
136    // *** ringorder_c and ringorder_C
137    if (type == ringorder_c || type == ringorder_C)
138    {
139      if (block0 != 0 || block1 != 0 || weights != 0)
140      {
141        WerrorS("Unexpected non-zero fields on c/C block in ring.");
142        return false;
143      }
144      if (didSetComponentBefore)
145      {
146        WerrorS("Unexpected two c/C blocks in ring.");
147        return false;
148      }
149      didSetComponentBefore = true;
150      if (r->order[block + 1] == ringorder_no)
151      {
152        conf.setComponentBefore
153          (mgb::GroebnerConfiguration::ComponentAfterBaseOrder);
154      } else
155        conf.setComponentBefore(gradings.size() / varCount);
156      conf.setComponentsAscending(type == ringorder_C);
157      continue;
158    }
159    if (block0 == 0 || block1 == 0)
160    {
161      WerrorS("Expected block0 != 0 and block1 != 0 in ring.");
162      return false;
163    }
164    if (block1 > varCount)
165    {
166      // todo: first handle any block types where this is not true
167      WerrorS("Expected block1 <= #vars in ring.");
168      return false;
169    }
170
171    // dim is how many variables this block concerns.
172    const size_t dim = static_cast<size_t>(block1 - block0 + 1);
173
174    // *** single-graded/ungraded lex/revlex orders
175    // a(w): w-graded and that's it
176    // a64(w): w-graded with 64-bit weights (not supported here)
177    //    lp:               lex from  left (descending)
178    //    Dp:  1-graded,    lex from  left (descending)
179    //    Ds: -1-graded,    lex from  left (descending)
180    // Wp(w):  w-graded,    lex from  left (descending)
181    // Ws(w): -w-graded,    lex from  left (descending)
182    //    rp:               lex from right (ascending)
183    //    rs:            revlex from right (descending)
184    //    dp:  1-graded, revlex from right (descending)
185    //    ds: -1-graded, revlex from right (descending)
186    // wp(w):  w-graded, revlex from right (descending)
187    // ws(w): -w-graded, revlex from right (descending)
188    //    ls:            revlex from  left (ascending)
189
190    if (type == ringorder_a64)
191    {
192      WerrorS("Block type a64 not supported for MathicGB interface.");
193      return false;
194    }
195
196    // * handle the single-grading part
197    const bool oneGrading = (type == ringorder_Dp || type == ringorder_dp);
198    const bool minusOneGrading = (type == ringorder_Ds || type == ringorder_ds);
199    const bool wGrading =
200      (type == ringorder_a || type == ringorder_Wp || type == ringorder_wp);
201    const bool minusWGrading = (type == ringorder_ws || type == ringorder_Ws);
202    if (oneGrading || minusOneGrading || wGrading || minusWGrading)
203    {
204      const VarIndex begin = gradings.size();
205      gradings.resize(begin + varCount);
206      if (oneGrading || minusOneGrading)
207      {
208        if (weights != 0)
209        {
210          WerrorS("Expect wvhdl == 0 in Dp/dp/Ds/ds-block in ring.");
211          return false;
212        }
213        const Exponent value = oneGrading ? 1 : -1;
214        for (int var = block0 - 1; var < block1; ++var)
215          gradings[begin + var] = value;
216      }
217      else
218      {
219        if (weights == 0)
220        {
221          WerrorS("Expect wvhdl != 0 in a/Wp/wp/ws/Ws-block in ring.");
222          return false;
223        }
224        if (wGrading)
225        {
226          for (int var = 0; var < dim; ++var)
227            gradings[begin + (block0 - 1) + var] = weights[var];
228        }
229        else
230        {
231          for (int var = 0; var < dim; ++var)
232            gradings[begin + (block0 - 1) + var] = -weights[var];
233        }
234      }
235    }
236    if (type == ringorder_a)
237      continue; // a has only the grading, so we are done already
238
239    // * handle the lex/revlex part
240    const bool lexFromLeft =
241      type == ringorder_lp ||
242      type == ringorder_Dp ||
243      type == ringorder_Ds ||
244      type == ringorder_Wp ||
245      type == ringorder_Ws;
246    const bool lexFromRight = type == ringorder_rp;
247    const bool revlexFromLeft = type == ringorder_ls;
248    const bool revlexFromRight =
249      type == ringorder_rs ||
250      type == ringorder_dp ||
251      type == ringorder_ds ||
252      type == ringorder_wp ||
253      type == ringorder_ws;
254    if (lexFromLeft || lexFromRight || revlexFromLeft || revlexFromRight)
255    {
256      const int next = r->order[block + 1];
257      bool final = next == ringorder_no;
258      if (!final && r->order[block + 2] == ringorder_no)
259        final = next == ringorder_c || next == ringorder_C;
260      if (final)
261      {
262        if (lexFromRight)
263          baseOrder = mgb::GroebnerConfiguration::LexAscendingBaseOrder;
264        else if (revlexFromRight)
265          baseOrder = mgb::GroebnerConfiguration::RevLexDescendingBaseOrder;
266        else if (lexFromLeft)
267          baseOrder = mgb::GroebnerConfiguration::LexDescendingBaseOrder;
268        else
269          baseOrder = mgb::GroebnerConfiguration::RevLexAscendingBaseOrder;
270        continue;
271      }
272
273      const size_t begin = gradings.size();
274      gradings.resize(begin + dim * varCount);
275      const Exponent value = (lexFromLeft || lexFromRight) ? 1 : -1;
276      if (lexFromLeft || revlexFromLeft)
277      {
278        for (size_t row = 0; row < dim; ++row)
279          gradings[begin + row * varCount + (block0 - 1) + row] = value;
280      }
281      else
282      {
283        for (size_t row = 0; row < dim; ++row)
284          gradings[begin + row * varCount + (block1 - 1) - row] = value;
285      }
286      continue;
287    }
288
289    // *** ringorder_M: a square invertible matrix
290    if (type == ringorder_M)
291    {
292      if (weights == 0)
293      {
294        WerrorS("Expected wvhdl != 0 in M-block in ring.");
295        return false;
296      }
297      const size_t begin = gradings.size();
298      gradings.resize(begin + dim * varCount);
299      for (size_t row = 0; row < dim; ++row)
300        for (size_t col = block0 - 1; col < block1; ++col)
301          gradings[begin + row * varCount + col] = weights[row * dim + col];
302      continue;
303    }
304
305    // *** Miscellaneous unsupported or invalid block types
306    if (
307      type == ringorder_s ||
308      type == ringorder_S ||
309      type == ringorder_IS
310    )
311    {
312      // todo: Consider supporting this later.
313      WerrorS("Schreyer order s/S/IS not supported in MathicGB interface.");
314      return false;
315    }
316    if (type == ringorder_am)
317    {
318      // This block is a Schreyer-like ordering only used in Spielwiese.
319      // todo: Consider supporting it later.
320      WerrorS("Block type am not supported in MathicGB interface");
321      return false;
322    }
323    if (type == ringorder_L)
324    {
325      WerrorS("Invalid L-block found in order of ring.");
326      return false;
327    }
328    if (type == ringorder_aa)
329    {
330      // I don't know what an aa block is supposed to do.
331      WerrorS("aa ordering not supported by the MathicGB interface.");
332      return false;
333    }
334    if (type == ringorder_unspec)
335    {
336      WerrorS("Invalid unspec-block found in order of ring.");
337      return false;
338    }
339    WerrorS("Unknown block type found in order of ring.");
340    return false;
341  }
342
343  if (!didSetComponentBefore)
344  {
345    WerrorS("Expected to find a c/C block in ring.");
346    return false;
347  }
348
349  if (!conf.setMonomialOrder(baseOrder, gradings))
350  {
351    WerrorS("MathicGB does not support non-global orders.");
352    return false;
353  }
354  return true;
355}
356
357bool prOrderMatrix(ring r)
358{
359  const int varCount = r->N;
360  mgb::GroebnerConfiguration conf(101, varCount,0);
361  if (!setOrder(r, conf))
362    return false;
363  const std::vector<Exponent>& gradings = conf.monomialOrder().second;
364  if (gradings.size() % varCount != 0)
365  {
366    WerrorS("Expected matrix to be a multiple of varCount.");
367    return false;
368  }
369  const size_t rowCount = gradings.size() / varCount;
370  std::cout << "Order matrix:\n";
371  for (size_t row = 0; row < rowCount; ++row)
372  {
373    for (size_t col = 0; col < varCount; ++col)
374      std::cerr << ' ' << gradings[row * varCount + col];
375    std::cerr << '\n';
376  }
377  std::cerr
378    << "Base order: "
379    << mgb::GroebnerConfiguration::baseOrderName(conf.monomialOrder().first)
380    << '\n';
381  std::cerr << "Component before: " << conf.componentBefore() << '\n';
382  std::cerr << "Components ascending: " << conf.componentsAscending() << '\n';
383  std::cerr << "Schreyering: " << conf.schreyering() << '\n';
384  return true;
385}
386
387void prOrder(ring r)
388{
389  std::cout << "Printing order of ring.\n";
390  for (int block = 0; ; ++block)
391  {
392    switch (r->order[block])
393    {
394    case ringorder_no: // end of blocks
395      return;
396
397    case ringorder_a:
398      std::cout << "a";
399      break;
400
401    case ringorder_a64: ///< for int64 weights
402      std::cout << "a64";
403      break;
404
405    case ringorder_c:
406      std::cout << "c";
407      break;
408
409    case ringorder_C:
410      std::cout << "C";
411      break;
412
413    case ringorder_M:
414      std::cout << "M";
415      break;
416
417    case ringorder_S: ///< S?
418      std::cout << "S";
419      break;
420
421    case ringorder_s: ///< s?
422      std::cout << "s";
423      break;
424
425    case ringorder_lp:
426      std::cout << "lp";
427      break;
428
429    case ringorder_dp:
430      std::cout << "dp";
431      break;
432
433    case ringorder_rp:
434      std::cout << "rp";
435      break;
436
437    case ringorder_Dp:
438      std::cout << "Dp";
439      break;
440
441    case ringorder_wp:
442      std::cout << "wp";
443      break;
444
445    case ringorder_Wp:
446      std::cout << "Wp";
447      break;
448
449    case ringorder_ls:
450      std::cout << "ls"; // not global
451      break;
452
453    case ringorder_ds:
454      std::cout << "ds"; // not global
455      break;
456
457    case ringorder_Ds:
458      std::cout << "Ds"; // not global
459      break;
460
461    case ringorder_ws:
462      std::cout << "ws"; // not global
463      break;
464
465    case ringorder_Ws:
466      std::cout << "Ws"; // not global
467      break;
468
469    case ringorder_am:
470      std::cout << "am";
471      break;
472
473    case ringorder_L:
474      std::cout << "L";
475      break;
476
477    // the following are only used internally
478    case ringorder_aa: ///< for idElimination, like a, except pFDeg, pWeigths ignore it
479      std::cout << "aa";
480      break;
481
482    case ringorder_rs: ///< opposite of ls
483      std::cout << "rs";
484      break;
485
486    case ringorder_IS: ///< Induced (Schreyer) ordering
487      std::cout << "IS";
488      break;
489
490    case ringorder_unspec:
491      std::cout << "unspec";
492      break;
493    }
494    const int b0 = r->block0[block];
495    const int b1 = r->block1[block];
496    std::cout << ' ' << b0 << ':' << b1 << " (" << r->wvhdl[block] << ")" << std::flush;
497    if (r->wvhdl[block] != 0 && b0 != 0)
498    {
499      for (int v = 0; v <= b1 - b0; ++v)
500        std::cout << ' ' << r->wvhdl[block][v];
501    } else
502      std::cout << " null";
503    std::cout << '\n';
504  }
505}
506
507BOOLEAN prOrderX(leftv result, leftv arg)
508{
509  if (currRing == 0)
510  {
511    WerrorS("There is no current ring.");
512    return TRUE;
513  }
514  prOrder(currRing);
515  prOrderMatrix(currRing);
516  result->rtyp=NONE;
517  return FALSE;
518}
519
520BOOLEAN mathicgb(leftv result, leftv arg)
521{
522  result->rtyp=NONE;
523
524  if (arg == NULL || arg->next != NULL ||
525  ((arg->Typ() != IDEAL_CMD) &&(arg->Typ() != MODUL_CMD)))
526  {
527    WerrorS("Syntax: mathicgb(<ideal>/<module>)");
528    return TRUE;
529  }
530  if (!rField_is_Zp(currRing))
531  {
532    WerrorS("Polynomial ring must be over Zp.");
533    return TRUE;
534  }
535
536  const int characteristic = n_GetChar(currRing);
537  const int varCount = currRing->N;
538  const ideal I=(ideal) arg->Data();
539  mgb::GroebnerConfiguration conf(characteristic, varCount,I->rank);
540  feOptIndex fno=feGetOptIndex(FE_OPT_THREADS);
541  //conf.setMaxThreadCount(0); // default number of cores
542  conf.setMaxThreadCount((int)(long)feOptSpec[fno].value);
543  if (!setOrder(currRing, conf))
544    return TRUE;
545  if (TEST_OPT_PROT)
546    conf.setLogging("all");
547
548  mgb::GroebnerInputIdealStream toMathic(conf);
549
550  const ideal id = static_cast<const ideal>(arg->Data());
551  const int size = IDELEMS(id);
552  toMathic.idealBegin(size);
553  for (int i = 0; i < size; ++i)
554  {
555    const poly origP = id->m[i];
556    int termCount = 0;
557    for (poly p = origP; p != 0; p = pNext(p))
558      ++termCount;
559    toMathic.appendPolynomialBegin(termCount);
560
561    for (poly p = origP; p != 0; p = pNext(p))
562    {
563      toMathic.appendTermBegin(pGetComp(p));
564      for (int i = 1; i <= currRing->N; ++i)
565        toMathic.appendExponent(i - 1, pGetExp(p, i));
566      const long coefLong = reinterpret_cast<long>(pGetCoeff(p));
567      toMathic.appendTermDone(static_cast<int>(coefLong));
568    }
569    toMathic.appendPolynomialDone();
570  }
571  toMathic.idealDone();
572
573  MathicToSingStream fromMathic(characteristic, varCount);
574  mgb::computeGroebnerBasis(toMathic, fromMathic);
575
576  result->rtyp = arg->Typ();
577  result->data = fromMathic.takeIdeal();
578  return FALSE;
579}
580
581template class std::vector<Exponent>;
582template void mgb::computeGroebnerBasis<MathicToSingStream>
583  (mgb::GroebnerInputIdealStream&, MathicToSingStream&);
584
585extern "C" int SI_MOD_INIT(singmathic)(SModulFunctions* psModulFunctions)
586{
587  psModulFunctions->iiAddCproc(
588    (currPack->libname ? currPack->libname : ""),
589    "mathicgb",
590    FALSE,
591    mathicgb
592  );
593  psModulFunctions->iiAddCproc(
594    (currPack->libname ? currPack->libname : ""),
595    "mathicgb_prOrder",
596    FALSE,
597    prOrderX
598  );
599  return MAX_TOK;
600}
601
602/* #else
603
604int SI_MOD_INIT(singmathic)(SModulFunctions* psModulFunctions)
605{
606  WerrorS(
607    "Cannot initialize the Singular interface to MathicGB "
608    "as this Singular executable was built without support "
609    "for MathicGB."
610  );
611  return 1;
612}
613*/
614
615/* ressources: ------------------------------------------------------------
616
617http://stackoverflow.com/questions/3786408/number-of-threads-used-by-intel-tbb
618When you create the scheduler, you can specify the number of threads as
619tbb::task_scheduler_init init(nthread);
620
621    How do I know how many threads are available?
622
623    Do not ask!
624
625        Not even the scheduler knows how many threads really are available
626        There may be other processes running on the machine
627        Routine may be nested inside other parallel routines
628
629  conf.setMaxThreadCount(0); // default number of cores
630*/
631#endif /* HAVE_MATHICGB */
Note: See TracBrowser for help on using the repository browser.