Changeset 4c4a80 in git for kernel


Ignore:
Timestamp:
Jan 26, 2010, 7:07:58 PM (14 years ago)
Author:
Martin Monerjan
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', '4188d308699580d975efd0f6cca8dcb41c396f70')
Children:
9008029955f32c6d5ef414a2e5d9bf43a5c892f6
Parents:
5f6e18a13e164d6edb16c75b117c5a59fa98cbf1
Message:
gfanp::parallelButNotEqual
Homogeneity preprocessing now works
changed to intvec::compare in gcone::areEqual
Correct way to set gcone::dd_LinealitySpace



git-svn-id: file:///usr/local/Singular/svn/trunk@12459 2c84dea3-7e68-4137-9b89-c4e89433aadc
Location:
kernel
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • kernel/gfan.cc

    r5f6e18 r4c4a80  
    198198inline bool gcone::areEqual(facet *f, facet *s)
    199199{
     200#ifdef gfanp
     201        gcone::numberOfFacetChecks++;
     202#endif
    200203        bool res = TRUE;
    201204        int notParallelCtr=0;
     
    210213//      intvec *foo=ivNeg(sNormal);
    211214//      if(fNormal->compare(foo)!=0) //facet normals
    212 //              return TRUE;
     215//              return TRUE;   
    213216        else//parallelity, so we check the codim2-facets
    214217        {
     
    226229                                intvec* s2Normal;
    227230                                s2Normal = s2Act->getFacetNormal();
    228                                 bool foo=areEqual(f2Normal,s2Normal);
    229                                 if(foo)
     231//                              bool foo=areEqual(f2Normal,s2Normal);
     232                                int foo=f2Normal->compare(s2Normal);
     233                                if(foo==0)
    230234                                        ctr++;
    231235                                s2Act = s2Act->next;
     
    241245                res=TRUE;
    242246        else
    243                 res=FALSE;                             
     247        {
     248#ifdef gfanp
     249                gcone::parallelButNotEqual++;
     250#endif
     251                res=FALSE;
     252        }
    244253        return res;
    245254}       
     
    276285}
    277286               
    278 /** Return the flipped GB
     287/** Returns a pointer to the flipped GB
    279288Seems not be used
    280289Anyhow idCopy would make sense here.
     
    317326}
    318327               
     328/** Returns a pointer to this->interiorPoint
     329* MIND: ivCopy returns a new intvec
     330* @see facet::getFacetNormal
     331*/
    319332inline intvec *facet::getInteriorPoint()
    320333{
    321         return this->interiorPoint;
     334        return ivCopy(this->interiorPoint);
    322335}       
    323336               
     
    689702        */
    690703//      ideal initialForm=idInit(IDELEMS(I),1);
    691         intvec *gamma=new intvec(this->numVars);
     704//      intvec *gamma=new intvec(this->numVars);
    692705        int falseGammaCounter=0;
    693706        int *redRowsArray=NULL;
     
    698711                ideal initialForm=idInit(IDELEMS(I),1);
    699712                //read row ii into gamma
    700                 double tmp;     
     713                double tmp;
     714                intvec *gamma=new intvec(this->numVars);
    701715                for(int jj=1;jj<=this->numVars;jj++)
    702716                {
     
    781795                delete P;
    782796                idDelete(&initialForm);
    783                 //idDelete(L);         
     797                //idDelete(L);
     798                delete gamma;
    784799        }//for(ii<ddineq-rowsize
    785         delete gamma;
     800//      delete gamma;
    786801        int offset=0;//needed for correction of redRowsArray[ii]
    787802        for( int ii=0;ii<num_elts;ii++ )
     
    846861#endif
    847862
    848         dd_MatrixCanonicalize(&ddineq, &ddlinset, &ddredrows, &ddnewpos, &dderr);
    849         //time(&canonicalizeTac);
    850         //cout << "dd_MatrixCanonicalize time: " << difftime(canonicalizeTac,canonicalizeTic) << "sec" << endl;
     863        dd_MatrixCanonicalize(&ddineq, &ddlinset, &ddredrows, &ddnewpos, &dderr);       
    851864        ddrows = ddineq->rowsize;       //Size of the matrix with redundancies removed
    852865        ddcols = ddineq->colsize;
    853                        
    854         //ddCreateMatrix(ddrows,ddcols+1);
    855         ddFacets = dd_CopyMatrix(ddineq);
     866       
     867        this->ddFacets = dd_CopyMatrix(ddineq);
    856868                       
    857869        /*Write the normals into class facet*/ 
     
    860872        for (int kk = 0; kk<ddrows; kk++)
    861873        {
    862                 intvec *load = new intvec(this->numVars);       //intvec to store a single facet normal that will then be stored via setFacetNormal
     874                intvec *load = new intvec(this->numVars);//intvec to store a single facet normal that will then be stored via setFacetNormal
    863875                for (int jj = 1; jj <ddcols; jj++)
    864876                {
     
    868880                }//for (int jj = 1; jj <ddcols; jj++)
    869881                               
    870                 /*Quick'n'dirty hack for flippability*/
    871 //              bool isFlip=FALSE;
    872 //              if(gcone::hasHomInput==FALSE)
    873 //              {
     882                /*Quick'n'dirty hack for flippability. Executed only if gcone::hasHomInput==FALSE
     883                * Otherwise every facet intersects the positive orthant
     884                */     
     885                if(gcone::hasHomInput==FALSE)
     886                {
    874887                        //TODO: No dP needed
    875888                        bool isFlip=FALSE;
    876                         for (int jj = 0; jj<load->length(); jj++)
     889                        for(int jj = 0; jj<load->length(); jj++)
    877890                        {
    878891                                intvec *ivCanonical = new intvec(load->length());
    879892                                (*ivCanonical)[jj]=1;
    880         //                                      cout << "dotProd=" << dotProduct(*load,*ivCanonical) << endl;
    881893                                if (dotProduct(*load,*ivCanonical)<0)   
    882                                                 //if (ivMult(load,ivCanonical)<0)
    883894                                {
    884895                                        isFlip=TRUE;
     
    887898                                delete ivCanonical;
    888899                        }/*End of check for flippability*/
    889                         if (isFlip==FALSE)
     900                        if(isFlip==FALSE)
    890901                        {
    891902                                this->numFacets++;
     
    911922#endif
    912923                        }
    913 //              }
    914                 else
     924                        else
     925                        {
     926                                this->numFacets++;
     927                                if(this->numFacets==1)
     928                                {
     929                                        facet *fRoot = new facet();
     930                                        this->facetPtr = fRoot;
     931                                        fAct = fRoot;
     932                                }
     933                                else
     934                                {
     935                                        fAct->next = new facet();
     936                                        fAct = fAct->next;
     937                                }
     938                                fAct->isFlippable=TRUE;
     939                                fAct->setFacetNormal(load);
     940                                fAct->setUCN(this->getUCN());                                   
     941                        }
     942                }//hasHomInput==FALSE
     943                else    //Every facet is flippable
    915944                {       /*Now load should be full and we can call setFacetNormal*/                                     
    916945                        this->numFacets++;
     
    9711000        time_getConeNormals += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
    9721001#endif
    973        
    9741002
    9751003}//gcone::getConeNormals(ideal I)
     
    11421170void gcone::getExtremalRays(const gcone &gc)
    11431171{
     1172#ifdef gfanp
     1173        timeval start, end;
     1174        gettimeofday(&start, 0);
     1175#endif
    11441176        //Add lineality space - dd_LinealitySpace
    11451177        dd_MatrixPtr ddineq;
     
    12171249                        fAct->isFlippable=FALSE;
    12181250                fAct = fAct->next;
    1219         }       
     1251        }
     1252#ifdef gfanp
     1253        gettimeofday(&end, 0);
     1254        t_getExtremalRays += (end.tv_sec - start.tv_sec + 1e-6*(end.tv_usec - start.tv_usec));
     1255#endif 
    12201256}               
    12211257               
     
    21322168        }//for(k
    21332169        /*Output kernel, i.e. set gcone::dd_LinealitySpace here*/       
    2134         gcone::dd_LinealitySpace = dd_CreateMatrix(dimKer,this->numVars+1);
     2170//      gcone::dd_LinealitySpace = dd_CreateMatrix(dimKer,this->numVars+1);
     2171        res = dd_CreateMatrix(dimKer,this->numVars+1);
    21352172        int row=-1;
    21362173        for(int kk=1;kk<n;kk++)
     
    21422179                        {
    21432180                                if(d[ii]>0)
    2144                                         mpq_set(dd_LinealitySpace->matrix[row][ii],ddineq->matrix[d[ii]][kk]);
     2181                                        mpq_set(res->matrix[row][ii],ddineq->matrix[d[ii]][kk]);
    21452182                                else if(ii==kk)                         
    2146                                         mpq_set_str(dd_LinealitySpace->matrix[row][ii],"1",10);
    2147                                 mpq_canonicalize(dd_LinealitySpace->matrix[row][ii]);
     2183                                        mpq_set_str(res->matrix[row][ii],"1",10);
     2184                                mpq_canonicalize(res->matrix[row][ii]);
    21482185                        }
    21492186                }
    21502187        }
    21512188        dd_FreeMatrix(ddineq);
     2189        return(res);
    21522190        //Better safe than sorry:
    21532191        //NOTE dd_LinealitySpace contains RATIONAL ENTRIES
    21542192        //Therefore if you need integer ones: CALL gcone::makeInt(...) method
    21552193}       
     2194
     2195/** Compute the kernel of a given mxn-matrix */
     2196// dd_MatrixPtr gcone::kernel(const dd_MatrixPtr &M)
     2197// {
     2198// }
     2199
    21562200/** \brief The new method of Markwig and Jensen
    21572201 * Compute gcBasis and facets for the arbitrary starting cone. Store \f$(codim-1)\f$-facets as normals.
     
    22022246        * Afterwards static dd_MatrixPtr gcone::dd_LinealitySpace is set
    22032247        */
    2204         gcAct->computeLinealitySpace();
     2248        dd_LinealitySpace = gcAct->computeLinealitySpace();
    22052249        /*End of lineality space computation*/         
    22062250//      gcAct->getCodim2Normals(*gcAct);
     
    32223266float gcone::time_getConeNormals;
    32233267float gcone::time_getCodim2Normals;
     3268float gcone::t_getExtremalRays;
    32243269float gcone::time_flip;
    32253270float gcone::t_markings;
     
    32313276float gcone::t_mI;
    32323277float gcone::t_iP;
     3278unsigned gcone::parallelButNotEqual=0;
     3279unsigned gcone::numberOfFacetChecks=0;
    32333280#endif
    32343281bool gcone::hasHomInput=FALSE;
     
    32403287        if(rHasGlobalOrdering(currRing))
    32413288        {       
    3242         int numvar = pVariables;
    3243         gfanHeuristic = h;
     3289                int numvar = pVariables;
     3290                gfanHeuristic = h;
     3291               
     3292                enum searchMethod {
     3293                        reverseSearch,
     3294                        noRevS
     3295                };
    32443296       
    3245         enum searchMethod {
    3246                 reverseSearch,
    3247                 noRevS
    3248         };
     3297                searchMethod method;
     3298                method = noRevS;
    32493299       
    3250         searchMethod method;
    3251         method = noRevS;
    3252 
    3253         ring inputRing=currRing;        // The ring the user entered
    3254         ring rootRing;                  // The ring associated to the target ordering
    3255         ideal res;     
    3256         facet *fRoot;   
    3257         dd_set_global_constants();
    3258         if(method==noRevS)
    3259         {
    3260                 gcone *gcRoot = new gcone(currRing,inputIdeal);
    3261                 gcone *gcAct;
    3262                 gcAct = gcRoot;
    3263                 gcAct->numVars=pVariables;
    3264                 gcAct->getGB(inputIdeal);
    3265                 /*Check whether input is homogeneous
    3266                 if TRUE each facet intersects the positive orthant, so we don't need the
    3267                 flippability test in getConeNormals
    3268                 */
    3269                 if(idHomIdeal(gcAct->gcBasis,NULL))
    3270                         gcone::hasHomInput=TRUE;
    3271 #ifndef NDEBUG
    3272                 cout << "GB of input ideal is:" << endl;
    3273 //              idShow(gcAct->gcBasis);
    3274 #endif
    3275                 if(gcAct->isMonomial(gcAct->gcBasis))
    3276                 {
    3277                         WerrorS("Monomial input - terminating");
    3278                         exit(-1);
     3300                ring inputRing=currRing;        // The ring the user entered
     3301                ring rootRing;                  // The ring associated to the target ordering
     3302//              ideal res;
     3303                dd_set_global_constants();
     3304                if(method==noRevS)
     3305                {
     3306                        gcone *gcRoot = new gcone(currRing,inputIdeal);
     3307                        gcone *gcAct;
     3308                        gcAct = gcRoot;
     3309                        gcAct->numVars=pVariables;
     3310                        gcAct->getGB(inputIdeal);
     3311                        /*Check whether input is homogeneous
     3312                        if TRUE each facet intersects the positive orthant, so we don't need the
     3313                        flippability test in getConeNormals & getExtremalRays
     3314                        */
     3315                        if(idHomIdeal(gcAct->gcBasis,NULL))
     3316                                gcone::hasHomInput=TRUE;
     3317        #ifndef NDEBUG
     3318        //              cout << "GB of input ideal is:" << endl;
     3319        //              idShow(gcAct->gcBasis);
     3320        #endif
     3321                        if(gcAct->isMonomial(gcAct->gcBasis))
     3322                        {
     3323                                WerrorS("Monomial input - terminating");
     3324                                exit(-1);
     3325                                gcAct->getConeNormals(gcAct->gcBasis);
     3326                                lResList=lprepareResult(gcAct,1);
     3327                                dd_free_global_constants;
     3328                                //This is filthy
     3329                                return lResList;
     3330                        }                       
    32793331                        gcAct->getConeNormals(gcAct->gcBasis);
    3280                         lResList=lprepareResult(gcAct,1);
    3281                         dd_free_global_constants;
    3282                         //This is filthy
    3283                         return lResList;
    3284                 }
    3285                 cout << endl;
    3286                 gcAct->getConeNormals(gcAct->gcBasis);
    3287                 gcAct->noRevS(*gcAct);
    3288                 rChangeCurrRing(inputRing);
    3289                 //res=gcAct->gcBasis;
    3290                 //Below is a workaround, since gcAct->gcBasis gets deleted in noRevS
    3291                 res = inputIdeal;
    3292                 lResList=lprepareResult(gcRoot,gcRoot->getCounter());
    3293                 /*Cleanup*/
    3294                 gcone *gcDel;   
    3295                 gcDel = gcRoot;
    3296                 gcAct = gcRoot;
    3297                 while(gcAct!=NULL)
    3298                 {
    3299                         gcDel = gcAct;
    3300                         gcAct = gcAct->next;
    3301                         delete gcDel;
    3302                 }
    3303         }
    3304         dd_FreeMatrix(gcone::dd_LinealitySpace);
    3305         dd_free_global_constants();
     3332                        gcAct->noRevS(*gcAct);  //Here we go!
     3333                        rChangeCurrRing(inputRing);
     3334                        //res=gcAct->gcBasis;
     3335                        //Below is a workaround, since gcAct->gcBasis gets deleted in noRevS
     3336//                      res = inputIdeal;
     3337                        lResList=lprepareResult(gcRoot,gcRoot->getCounter());
     3338                        /*Cleanup*/
     3339                        gcone *gcDel;   
     3340                        gcDel = gcRoot;
     3341                        gcAct = gcRoot;
     3342                        while(gcAct!=NULL)
     3343                        {
     3344                                gcDel = gcAct;
     3345                                gcAct = gcAct->next;
     3346                                delete gcDel;
     3347                        }
     3348                }//method==noRevS
     3349                dd_FreeMatrix(gcone::dd_LinealitySpace);
     3350                dd_free_global_constants();
    33063351        }//rHasGlobalOrdering
    33073352        else
     
    33173362#ifdef gfanp
    33183363        cout << "t_getConeNormals:" << gcone::time_getConeNormals << endl;
    3319         cout << "t_getCodim2Normals:" << gcone::time_getCodim2Normals << endl;
    3320         cout << "  t_ddMC:" << gcone::t_ddMC << endl;
    3321         cout << "  t_mI:" << gcone::t_mI << endl;
    3322         cout << "  t_iP:" << gcone::t_iP << endl;
     3364//      cout << "t_getCodim2Normals:" << gcone::time_getCodim2Normals << endl;
     3365//      cout << "  t_ddMC:" << gcone::t_ddMC << endl;
     3366//      cout << "  t_mI:" << gcone::t_mI << endl;
     3367//      cout << "  t_iP:" << gcone::t_iP << endl;
     3368        cout << "t_getExtremalRays:" << gcone::t_getExtremalRays << endl;
    33233369        cout << "t_Flip:" << gcone::time_flip << endl;
    33243370        cout << "  t_markings:" << gcone::t_markings << endl;
     
    33273373        cout << "t_computeInv:" << gcone::time_computeInv << endl;
    33283374        cout << "t_enqueue:" << gcone::time_enqueue << endl;
     3375        cout << endl;
     3376        cout << "Checked " << gcone::numberOfFacetChecks << " Facets" << endl;
     3377        cout << " out of which there were " << gcone::parallelButNotEqual << " parallel but not equal." << endl;
    33293378#endif
    33303379        cout << "Maximum lenght of list of facets: " << gcone::maxSize << endl;
  • kernel/gfan.h

    r5f6e18 r4c4a80  
    141141                static float time_getConeNormals;
    142142                static float time_getCodim2Normals;
     143                static float t_getExtremalRays;
    143144                static float time_flip;
    144145                static float t_ffG;
     
    151152                static float t_mI;
    152153                static float t_iP;
     154                static unsigned parallelButNotEqual;
     155                static unsigned numberOfFacetChecks;
    153156#endif
    154157                /** Matrix to contain the homogeneity/lineality space */
     
    169172                /**
    170173                 * At least as a workaround we store the irredundant facets of a matrix here.
    171                  * Otherwise, since we throw away non-flippable facets, facets2Matrix will not
    172                  * yield all the necessary information
     174                 * This is needed to compute an interior points of a cone. Note that there
     175                 * will be non-flippable facets in it!           
    173176                 */
    174177                dd_MatrixPtr ddFacets;  //Matrix to store irredundant facets of the cone
Note: See TracChangeset for help on using the changeset viewer.