Changeset 224d153 in git for kernel/gfan.cc


Ignore:
Timestamp:
May 29, 2009, 9:52:12 AM (15 years ago)
Author:
Martin Monerjan
Branches:
(u'spielwiese', '17f1d200f27c5bd38f5dfc6e8a0879242279d1d8')
Children:
3ea5cee1504eb5dced1aa095461b2eafa2c20a37
Parents:
69aadadc4549786f157e7b86a56775084cb15d54
Message:
More bugfixes in normalize
Implemented list of (codim-2)-facets using the facet structure.
NOTE: Need to save this when deleting cones in order to preserve the list
of comparable facets


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

Legend:

Unmodified
Added
Removed

  • kernel/gfan.cc

    r69aada r224d153  
    22Compute the Groebner fan of an ideal
    33$Author: monerjan $
    4 $Date: 2009-05-28 06:27:09 $
    5 $Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.56 2009-05-28 06:27:09 monerjan Exp $
    6 $Id: gfan.cc,v 1.56 2009-05-28 06:27:09 monerjan Exp $
     4$Date: 2009-05-29 07:52:12 $
     5$Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.57 2009-05-29 07:52:12 monerjan Exp $
     6$Id: gfan.cc,v 1.57 2009-05-29 07:52:12 monerjan Exp $
    77*/
    88
    … …  
    7777                int UCN;
    7878               
     79                /** \brief The codim of the facet
     80                */
     81                int codim;
     82               
    7983                /** \brief The Groebner basis on the other side of a shared facet
    8084                 *
    … …  
    9094                bool isIncoming;        //Is the facet incoming or outgoing?
    9195                facet *next;            //Pointer to next facet
     96                facet *codim2Ptr;       //Pointer to (codim-2)-facet. Bit of recursion here ;-)
    9297                //intvec **codim2Normals =(intvec**)omAlloc0(sizeof(intvec*));  //Integer matrix containing the (codim-2)-facets
    93                 struct c2N{
    94                         intvec normal;
    95                         intvec *next;
    96                 };
    97                                
     98                                               
    9899                /** The default constructor. Do I need a constructor of type facet(intvec)? */
    99100                facet()                 
    … …  
    103104                        this->next=NULL;
    104105                        this->UCN=NULL;
    105                         //this->codim2Normals=NULL;                     
     106                        this->codim2Ptr=NULL;
     107                        this->codim=1;          //default for (codim-1)-facets
     108                }
     109               
     110                /** \brief Constructor for facets of codim >= 2
     111                */
     112                facet(int const &n)
     113                {
     114                        this->next=NULL;
     115                        this->UCN=NULL;
     116                        this->codim2Ptr=NULL;
     117                        if(n>1)
     118                        {
     119                                this->codim=n;
     120                        }//NOTE Handle exception here!
    106121                }
    107122               
    … …  
    112127                void setFacetNormal(intvec *iv)
    113128                {
    114                         this->fNormal = ivCopy(iv);
    115                         //return;
     129                        this->fNormal = ivCopy(iv);                     
    116130                }
    117131               
    118132                /** Hopefully returns the facet normal */
    119133                intvec *getFacetNormal()
    120                 {       
    121                         //this->fNormal->show();        cout << endl;   
     134                {                               
    122135                        return this->fNormal;
    123136                }
    … …  
    165178                {
    166179                        return this->interiorPoint;
    167                 }
     180                }               
    168181               
    169182                /*bool isFlippable(intvec &load)
    … …  
    251264                {
    252265                        this->next=NULL;
    253                         this->facetPtr=NULL;
     266                        this->facetPtr=NULL;    //maybe this->facetPtr = new facet();
    254267                        this->baseRing=currRing;
    255268                        this->UCN=1;
     269                        this->numFacets=0;
    256270                }
    257271               
    … …  
    271285                        this->baseRing=currRing;
    272286                        this->UCN=1;
     287                        this->numFacets=0;
    273288                }
    274289               
    … …  
    279294                */                     
    280295                //NOTE Prolly need to specify the facet to flip over
    281                 gcone(const gcone& gc)         
     296                gcone(const gcone& gc, const facet &f)         
    282297                {
    283298                        this->next=NULL;
    284299                        this->numVars=gc.numVars;
    285                         this->UCN=(gc.UCN)+1;   //add 1 to the UCN of previous cone
     300                        this->UCN=(gc.UCN)+1;   //add 1 to the UCN of previous cone. This is NOT UNIQUE!
    286301                        facet *fAct= new facet();                       
    287302                        this->facetPtr=fAct;
    … …  
    305320                }
    306321               
    307                 /** \brief Default destructor */
    308                 ~gcone(){;}             //destructor
     322                /** \brief Default destructor
     323                */
     324                ~gcone()
     325                {
     326                        //NOTE SAVE THE FACET STRUCTURE!!!
     327                }
    309328                       
    310329               
    … …  
    340359                        //delete iv;
    341360                        //delete f;
    342                 }
     361                }               
    343362               
    344363                /** \brief Set gcone::numFacets */
    … …  
    531550                                        fAct->next = new facet();                                       
    532551                                        fAct=fAct->next;        //this should definitely not be called in the above while-loop :D
     552                                        this->numFacets++;
    533553                                }//if (isFlippable==FALSE)
    534554                                //delete load;
    … …  
    551571                       
    552572                        //Compute the number of facets
    553                         this->numFacets=ddineq->rowsize;
     573                        // wrong because ddineq->rowsize contains only irredundant facets. There can still be
     574                        // non-flippable ones!
     575                        //this->numFacets=ddineq->rowsize;
    554576                       
    555577                        //Clean up but don't delete the return value! (Whatever it will turn out to be)                 
    … …  
    14171439                                for (int ii=0; ii<this->numFacets; ii++)
    14181440                                {       
    1419                                         cout << "------------" << endl;
     1441                                        cout << endl << "------------" << endl;
    14201442                                        ddakt = dd_CopyMatrix(ddineq);
    14211443                                        set_addelem(ddakt->linset,ii+1);
    … …  
    14231445                                        dd_MatrixCanonicalize(&ddakt, &impl_linset, &redset, &newpos, &err);                   
    14241446                                       
    1425                                         dd_WriteMatrix(stdout,ddakt);
     1447                                        //dd_WriteMatrix(stdout,ddakt);
    14261448                                        ddpolyh=dd_DDMatrix2Poly(ddakt, &err);
    14271449                                        P=dd_CopyGenerators(ddpolyh);
    … …  
    14321454                                        * and add the resulting vector to the int matrix facet::codim2Facets
    14331455                                        */
     1456                                        this->facetPtr->codim2Ptr = new facet(2);       //construct a (codim-2)-facet
     1457                                        facet *codim2Act;
     1458                                        codim2Act = this->facetPtr->codim2Ptr;
    14341459                                        for (int jj=1;jj<=P->rowsize;jj++)
    14351460                                        {
    14361461                                                intvec *n = new intvec(this->numVars);
    1437                                                 normalize(P,jj,*n);
    1438                                                 //fAct->addCodim2Facet(n);                                             
     1462                                                normalize(P,jj,*n);                                             
     1463                                                codim2Act->setFacetNormal(n);
     1464                                                codim2Act->next = new facet(2);
     1465                                                codim2Act = codim2Act->next;
    14391466                                                n->show();
    14401467                                                delete n;                                                                       
    1441                                         }
    1442                                        
    1443                                         //intvec *load = new intvec(this->numVars);
     1468                                        }                                       
    14441469                                       
    14451470                                        dd_FreeMatrix(ddakt);
    14461471                                        dd_FreePolyhedra(ddpolyh);
    1447                                 }
    1448                                 //gc.showFacets();
    1449                                
     1472                                }                               
    14501473                        }                       
    14511474       
    … …  
    14631486                        //NOTE: gcTmp may be deleted, gcRoot from main proc should probably not!
    14641487                       
    1465                 }//void noRevS(gcone &gc)
    1466                
    1467                 void normalize(dd_MatrixPtr const &M, int line, intvec &n)
     1488                }//void noRevS(gcone &gc)       
     1489               
     1490               
     1491                /** \brief Make a set of rational vectors into integers
     1492                *
     1493                * We compute the lcm of the denominators and multiply with this to get integer values
     1494                * \param dd_MatrixPtr,intvec
     1495                */
     1496                void normalize(dd_MatrixPtr const &M, int const line, intvec &n)
    14681497                {                       
    14691498                        mpz_t denom[this->numVars];
    … …  
    14761505                        mpz_init(kgV);
    14771506                        mpz_init(tmp);
    1478                         //intvec *ivres = new intvec(this->numVars);
    1479 //                      intvec ivres(this->numVars);
    14801507                       
    14811508                        for (int ii=0;ii<(M->colsize)-1;ii++)
    … …  
    14841511                                mpz_init(z);
    14851512                                mpq_get_den(z,M->matrix[line-1][ii+1]);
    1486                                 //mpz_out_str(stdout,10,z);
    14871513                                mpz_set( denom[ii], z);
    14881514                                mpz_clear(z);                           
    14891515                        }
     1516                       
    14901517                        /*Compute lcm of the denominators*/
    14911518                        mpz_set(tmp,denom[0]);
    … …  
    14941521                                mpz_lcm(kgV,tmp,denom[ii]);                             
    14951522                        }
     1523                       
    14961524                        /*Multiply the nominators by kgV*/
    14971525                        mpq_t qkgV,res;
    14981526                        mpq_init(qkgV);
    1499 //                      mpq_canonicalize(qkgV);
    15001527                        mpq_set_str(qkgV,"1",10);
    1501 //                      mpq_canonicalize(qkgV);
     1528                        mpq_canonicalize(qkgV);
     1529                       
    15021530                        mpq_init(res);
    15031531                        mpq_set_str(res,"1",10);
    1504 //                      mpq_canonicalize(res);
    1505                         //mpq_set_z(qkgV,kgV);
     1532                        mpq_canonicalize(res);
     1533                       
    15061534                        mpq_set_num(qkgV,kgV);
    1507                         //mpq_set_den(qkgV,1);
     1535                       
    15081536//                      mpq_canonicalize(qkgV);
    15091537                        for (int ii=0;ii<(M->colsize)-1;ii++)
    15101538                        {
    15111539                                mpq_mul(res,qkgV,M->matrix[line-1][ii+1]);
    1512 //                              (*ivres)[ii]=(int)mpz_get_d(mpq_numref(res));
    15131540                                n[ii]=(int)mpz_get_d(mpq_numref(res));
    15141541                        }
    15151542                        //mpz_clear(denom[this->numVars]);
    15161543                        mpz_clear(kgV);
    1517                         mpq_clear(qkgV); mpq_clear(res);
    1518                        
    1519                         //return ivres;
     1544                        mpq_clear(qkgV); mpq_clear(res);                       
     1545                       
    15201546                }
    15211547               
    … …  
    15411567                       
    15421568                        for (int jj=0; jj<M->rowsize; jj++)
    1543                         {
     1569                        {                               
    15441570                                comp = fAct->getFacetNormal();
    15451571                                for (int ii=0; ii<this->numVars; ii++)
    … …  
    15471573                                        dd_set_si(M->matrix[jj][ii+1],(*comp)[ii]);
    15481574                                }
    1549                                 fAct = fAct->next;                             
     1575                                if(fAct->next!=NULL)
     1576                                {
     1577                                        fAct = fAct->next;
     1578                                }
    15501579                        }//jj
    15511580                       
    … …  
    16901719                gcAct->numVars=pVariables;
    16911720                gcAct->getGB(inputIdeal);
    1692                 gcAct->getConeNormals(gcAct->gcBasis);
    1693                 cout << "ding" << endl;         
     1721                gcAct->getConeNormals(gcAct->gcBasis);         
    16941722                gcAct->noRevS(*gcAct);         
    16951723                res=gcAct->gcBasis;     
Note: See TracChangeset for help on using the changeset viewer.