Changeset a9fa92 in git


Ignore:
Timestamp:
Apr 2, 2009, 11:44:23 AM (15 years ago)
Author:
Martin Monerjan
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'd25190065115c859833252500a64cfb7b11e3a50')
Children:
f6c6b01d649bc867417af16dfdfc71c7a9785e8d
Parents:
95a5d329460059206ad840fd798ce5fe36656f1f
Message:
Removed bug in gcone::getConeNormals due to which the actual data of facet::fNormal was lost.
Introduced void facet::*getFacetNormal() to return the facet normal.
More experimental code on parallelity. Will be removed by a method that is not inherently numerical instable like brute force division over all components.


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

Legend:

Unmodified
Added
Removed
  • kernel/gfan.cc

    r95a5d3 ra9fa92  
    22Compute the Groebner fan of an ideal
    33$Author: monerjan $
    4 $Date: 2009-04-01 14:07:20 $
    5 $Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.26 2009-04-01 14:07:20 monerjan Exp $
    6 $Id: gfan.cc,v 1.26 2009-04-01 14:07:20 monerjan Exp $
     4$Date: 2009-04-02 09:44:23 $
     5$Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.27 2009-04-02 09:44:23 monerjan Exp $
     6$Id: gfan.cc,v 1.27 2009-04-02 09:44:23 monerjan Exp $
    77*/
    88
     
    6868                /** Stores the facet normal \param intvec*/
    6969                void setFacetNormal(intvec *iv){
    70                         fNormal = iv;
     70                        this->fNormal = ivCopy(iv);
    7171                        //return;
     72                }
     73               
     74                /** Hopefully returns the facet normal */
     75                intvec *getFacetNormal()
     76                {       
     77                        //this->fNormal->show();        cout << endl;   
     78                        return this->fNormal;
    7279                }
    7380               
     
    250257#endif
    251258                        /*The pointer *fRoot should be the return value of this function*/
    252                         facet *fRoot = new facet();             //instantiate new facet with intvec with numvar rows, one column and initial values all 0
    253                         facetPtr = fRoot;               //set variable facetPtr of class gcone to first facet
     259                        facet *fRoot = new facet();     //instantiate new facet
     260                        this->facetPtr = fRoot;         //set variable facetPtr of class gcone to first facet
    254261                        facet *fAct;                    //instantiate pointer to active facet
    255262                        fAct = fRoot;           //This does not seem to do the trick. fRoot and fAct have to point to the same adress!
     
    269276                                        if (jj<ddcols)                          //Is this facet NOT the last facet? Writing while instead of if is a really bad idea :)
    270277                                        {
    271                                                 fAct->next = new facet();       //If so: instantiate new facet. Otherwise this->next=NULL due to the constructor
    272                                                 fAct = fAct->next;              //scary :)
     278                                                fAct->next = new facet();       //If so: instantiate new facet. Otherwise this->next=NULL due to the constructor                                               
    273279                                        }
    274280                                }//for jj<ddcols
     
    276282                                fAct->setFacetNormal(load);
    277283                                fAct->printNormal();
     284                                fAct=fAct->next;        //this should definitely not be called in the above while-loop :D
     285                                delete load;
    278286                        }
    279287                        /*
     
    312320                void flip(ideal gb, facet *f)           //Compute "the other side"
    313321                {                       
    314                         intvec fNormal; //facet normal, check for parallelity
    315                         /* EXTREMELY EXPERIMENTAL CODE AHEAD*/
     322                       
     323                        intvec *fNormal = new intvec(this->numVars);    //facet normal, check for parallelity                   
     324                        fNormal = f->getFacetNormal();  //read this->fNormal;                   
     325                        fNormal->show();
     326                       
     327                        bool isParallel=TRUE;
     328                        double lambda=0;
     329                        double lambda_temp=lambda;
     330                        /* EXTREMELY EXPERIMENTAL CODE AHEAD - NUMERICAL INSTABILITIES GUARANTEED DUE TO DIVISONS*/
    316331                        /*1st step: Compute the initial ideal*/
    317332                        poly initialForm[IDELEMS(gb)];  //array of #polys in GB to store initial form
    318333                        poly aktpoly;
    319334                        int check[this->numVars];       //array to store the difference of LE and v
     335                       
    320336                        for (int ii=0;ii<IDELEMS(gb);ii++)
    321337                        {
    322                                 aktpoly = (poly)gb->m[ii];
     338                                aktpoly = (poly)gb->m[ii];                                                             
    323339                                int *v=(int *)omAlloc((this->numVars+1)*sizeof(int));
    324340                                int *leadExpV=(int *)omAlloc((this->numVars+1)*sizeof(int));                           
     341                                intvec *iv=(intvec *)omAlloc((this->numVars+1)*sizeof(intvec)); //let's try this
     342                                intvec *ivLeadExpV=(intvec *)omAlloc((this->numVars+1)*sizeof(intvec));
    325343                                pGetExpV(aktpoly,leadExpV);     //find the leading exponent in leadExpV[1],...,leadExpV[n], use pNext(p)
     344                                //pGetExpV(aktpoly,ivLeadExpV);
    326345                                initialForm[ii]=pHead(aktpoly);
     346                               
    327347                                while(pNext(aktpoly)!=NULL)     /*loop trough terms and check for parallelity*/
    328348                                {
    329349                                        aktpoly=pNext(aktpoly); //next term
    330350                                        pGetExpV(aktpoly,v);
     351                                       
    331352                                        for (int i=0;i<pVariables;i++)
    332353                                        {
    333354                                                check[i] = v[i]-leadExpV[i];
     355                                                cout << "check["<<i<<"]="<<check[i]<<endl;
    334356                                        }
    335                                         if (isParallel(*check,fNormal)) //pass *check when
     357                                        lambda=check[0]/(*fNormal)[0];
     358                                        lambda_temp=lambda;
     359                                        cout << "lambda="<<lambda<<endl;
     360                                        for (int i=0;i<pVariables;i++)
    336361                                        {
     362                                                lambda_temp=check[i]/(*fNormal)[i];
     363                                                cout<<"lambda_temp="<<lambda_temp<<endl;
     364                                                if (lambda_temp!=lambda)
     365                                                {
     366                                                        isParallel=FALSE;
     367                                                }
     368                                        }
     369                                        if (isParallel==TRUE) //pass *check when
     370                                        {
     371                                                cout << "PARALLEL" << endl;
    337372                                                //initialForm[ii] += pHead(aktpoly);    //This should add the respective term to
    338373                                        }                               
     
    417452        gcAct->getGB(inputIdeal);
    418453        gcAct->getConeNormals(gcAct->gcBasis);  //hopefully compute the normals
    419        
     454        gcAct->flip(gcAct->gcBasis,gcAct->facetPtr);
    420455        /*Now it is time to compute the search facets, respectively start the reverse search.
    421456        But since we are in the root all facets should be search facets. IS THIS TRUE?
Note: See TracChangeset for help on using the changeset viewer.