Changeset 7fee4f6 in git


Ignore:
Timestamp:
Jun 14, 2009, 12:00:14 PM (15 years ago)
Author:
Martin Monerjan
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', '3720ae8bfcff4a4649ee98a15552089151d2d59b')
Children:
e28cf1ea73e36a7c3768b1ea8eb0e279b22ecede
Parents:
9918fd329f6178bcb0152873e21d810764461b77
Message:
Fixed computation of interior points
Tried rCopyAndAddWeight for computation of dstRing. Needs more work


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

Legend:

Unmodified
Added
Removed
  • kernel/gfan.cc

    r9918fd r7fee4f6  
    22Compute the Groebner fan of an ideal
    33$Author: monerjan $
    4 $Date: 2009-06-12 16:32:00 $
    5 $Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.59 2009-06-12 16:32:00 monerjan Exp $
    6 $Id: gfan.cc,v 1.59 2009-06-12 16:32:00 monerjan Exp $
     4$Date: 2009-06-14 10:00:14 $
     5$Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.60 2009-06-14 10:00:14 monerjan Exp $
     6$Id: gfan.cc,v 1.60 2009-06-14 10:00:14 monerjan Exp $
    77*/
    88
     
    940940                        turn the minimal basis into a reduced one
    941941                        */
    942                         int i,j;
    943                         ring dstRing=rCopy0(srcRing);
    944                         i=rBlocks(srcRing);
    945                        
    946                         dstRing->order=(int *)omAlloc((i+1)*sizeof(int));
    947                         for(j=i;j>0;j--)
    948                         {
    949                                 dstRing->order[j]=srcRing->order[j-1];
    950                                 dstRing->block0[j]=srcRing->block0[j-1];
    951                                 dstRing->block1[j]=srcRing->block1[j-1];
    952                                 if (srcRing->wvhdl[j-1] != NULL)
    953                                 {
    954                                         dstRing->wvhdl[j] = (int*) omMemDup(srcRing->wvhdl[j-1]);
    955                                 }
    956                         }
    957                         dstRing->order[0]=ringorder_a;
    958                         dstRing->order[1]=ringorder_dp;
    959                         dstRing->order[2]=ringorder_C;                 
    960                         dstRing->wvhdl[0] =( int *)omAlloc((iv_weight->length())*sizeof(int));
    961                        
    962                         for (int ii=0;ii<this->numVars;ii++)
    963                         {                               
    964                                 dstRing->wvhdl[0][ii]=(*iv_weight)[ii];                         
    965                         }
    966                         rComplete(dstRing);
     942                        ring dstRing=rCopyAndAddWeight(tmpRing,iv_weight);     
     943                       
     944//                      int i,j;
     945//                      ring dstRing=rCopy0(srcRing);
     946//                      i=rBlocks(srcRing);
     947//                     
     948//                      dstRing->order=(int *)omAlloc((i+1)*sizeof(int));
     949//                      for(j=i;j>0;j--)
     950//                      {
     951//                              dstRing->order[j]=srcRing->order[j-1];
     952//                              dstRing->block0[j]=srcRing->block0[j-1];
     953//                              dstRing->block1[j]=srcRing->block1[j-1];
     954//                              if (srcRing->wvhdl[j-1] != NULL)
     955//                              {
     956//                                      dstRing->wvhdl[j] = (int*) omMemDup(srcRing->wvhdl[j-1]);
     957//                              }
     958//                      }
     959//                      dstRing->order[0]=ringorder_a;
     960//                      dstRing->order[1]=ringorder_dp;
     961//                      dstRing->order[2]=ringorder_C;                 
     962//                      dstRing->wvhdl[0] =( int *)omAlloc((iv_weight->length())*sizeof(int));
     963//                     
     964//                      for (int ii=0;ii<this->numVars;ii++)
     965//                      {                               
     966//                              dstRing->wvhdl[0][ii]=(*iv_weight)[ii];                         
     967//                      }
     968//                      rComplete(dstRing);
    967969                       
    968970                        // NOTE May assume that at this point srcRing already has 3 blocks of orderins, starting with a
     
    10001002                        f->printFlipGB();
    10011003#endif                 
    1002                         rChangeCurrRing(srcRing);       //return to the ring we started the computation of flipGB in
     1004                        //rChangeCurrRing(srcRing);     //return to the ring we started the computation of flipGB in
    10031005                }//void flip(ideal gb, facet *f)
    10041006                               
     
    12121214                        //dd_LPSolve(lpInt,solver,&err);        //This will not result in a point from the relative interior
    12131215                        if (err!=dd_NoError){cout << "Error during dd_LPSolve" << endl;}
    1214                         //cout << "Tick 5" << endl;
    12151216                                                                       
    12161217                        //lpSol=dd_CopyLPSolution(lpInt);
    12171218                        if (err!=dd_NoError){cout << "Error during dd_CopyLPSolution" << endl;}                 
    1218                         //cout << "Tick 6" << endl;
    12191219#ifdef gfan_DEBUG
    12201220                        cout << "Interior point: ";
    1221 #endif
    12221221                        for (int ii=1; ii<(lpSol->d)-1;ii++)
    12231222                        {
    1224 #ifdef gfan_DEBUG
    12251223                                dd_WriteNumber(stdout,lpSol->sol[ii]);
    1226 #endif
    1227                                 /* NOTE This works only as long as gmp returns fractions with the same denominator*/
    1228                                 (iv)[ii-1]=(int)mpz_get_d(mpq_numref(lpSol->sol[ii]));  //looks evil, but does the trick
    1229                         }
     1224                        }
     1225#endif
     1226                       
     1227                        //NOTE The following strongly resembles parts of makeInt.
     1228                        //Maybe merge sometimes
     1229                        mpz_t kgV; mpz_init(kgV);
     1230                        mpz_set_str(kgV,"1",10);
     1231                        mpz_t den; mpz_init(den);
     1232                        mpz_t tmp; mpz_init(tmp);
     1233                        mpq_get_den(tmp,lpSol->sol[1]);
     1234                        for (int ii=1;ii<(lpSol->d)-1;ii++)
     1235                        {
     1236                                mpq_get_den(den,lpSol->sol[ii+1]);
     1237                                mpz_lcm(kgV,tmp,den);
     1238                                mpz_set(tmp, kgV);
     1239                        }
     1240                        mpq_t qkgV;
     1241                        mpq_init(qkgV);
     1242                        mpq_set_z(qkgV,kgV);
     1243                        for (int ii=1;ii<(lpSol->d)-1;ii++)
     1244                        {
     1245                                mpq_t product;
     1246                                mpq_init(product);
     1247                                mpq_mul(product,qkgV,lpSol->sol[ii]);
     1248                                iv[ii-1]=(int)mpz_get_d(mpq_numref(product));
     1249                                mpq_clear(product);
     1250                        }
     1251                        mpq_clear(qkgV);
     1252                        mpz_clear(tmp);
     1253                        mpz_clear(den);
     1254                        mpz_clear(kgV);                 
     1255                       
    12301256                        dd_FreeLPSolution(lpSol);
    12311257                        dd_FreeLPData(lpInt);
Note: See TracChangeset for help on using the changeset viewer.