Changeset 94aed9 in git for kernel/gfan.cc

Timestamp:

Jul 15, 2009, 11:49:09 AM (15 years ago)

Author:

Martin Monerjan

Branches:

(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')

Children:

d4f44f4f228fbb8258ce8e9927141d37d4470719

Parents:

7ae94bb7c3176ec1860c1447b4e220578a2affa3

Message:

Implemented codim2 check
Workaround for problems with intPointMatrix in flip


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

File:

• kernel/gfan.cc (modified) (18 diffs)

kernel/gfan.cc

@@ -2 +2 @@
 Compute the Groebner fan of an ideal
 $Author: monerjan $
-$Date: 2009-07-13 09:03:37 $
-$Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.76 2009-07-13 09:03:37 monerjan Exp $
-$Id: gfan.cc,v 1.76 2009-07-13 09:03:37 monerjan Exp $
+$Date: 2009-07-15 09:49:09 $
+$Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.77 2009-07-15 09:49:09 monerjan Exp $
+$Id: gfan.cc,v 1.77 2009-07-15 09:49:09 monerjan Exp $
 */
 
@@ -656 +656 @@
                         ddFacets = dd_CopyMatrix(ddineq);
 #ifdef gfan_DEBUG
-                        cout << "Having removed redundancies, the normals now read:" << endl;
-                        dd_WriteMatrix(stdout,ddineq);
+//                      cout << "Having removed redundancies, the normals now read:" << endl;
+//                      dd_WriteMatrix(stdout,ddineq);
 //                      cout << "Rows = " << ddrows << endl;
 //                      cout << "Cols = " << ddcols << endl;
@@ -664 +664 @@
                         /*Write the normals into class facet*/
 #ifdef gfan_DEBUG
-                        cout << "Creating list of normals" << endl;
+//                      cout << "Creating list of normals" << endl;
 #endif
                         /*The pointer *fRoot should be the return value of this function*/
@@ -673 +673 @@
                         //fAct = fRoot;                 //Seems to do the trick. fRoot and fAct have to point to the same adress!
                         //std::cout << "fRoot = " << fRoot << ", fAct = " << fAct << endl;
+                        int numNonFlip=0;
                         for (int kk = 0; kk<ddrows; kk++)
                         {
@@ -687 +688 @@
                                 
                                 /*Quick'n'dirty hack for flippability*/ 
-                                bool isFlip=FALSE;                              
+                                bool isFlip=FALSE;
+                                                                
                                 for (int jj = 0; jj<load->length(); jj++)
                                 {
@@ -703 +705 @@
                                 {
                                         this->numFacets++;
+                                        numNonFlip++;
                                         if(this->numFacets==1)
                                         {
@@ -742 +745 @@
                                         fAct->setUCN(this->getUCN());                                   
                                 }//if (isFlippable==FALSE)
-                                //delete load;
+                                //delete load;                          
                         }//for (int kk = 0; kk<ddrows; kk++)
+                        
+                        //In cases like I=<x-1,y-1> there are only non-flippable facets...
+                        if(numNonFlip==this->numFacets)
+                        {                                       
+                                cerr << "Only non-flippable facets. Terminating..." << endl;
+                        }
                         
                         /*
@@ -799 +808 @@
                         /*Now set appropriate linearity*/
                         dd_PolyhedraPtr ddpolyh;
-                        for (int ii=0; ii<this->numFacets; ii++)
-                        //while(fAct!=NULL)
+                        for (int ii=0; ii<this->numFacets; ii++)                        
                         {                               
                                 ddakt = dd_CopyMatrix(ddineq);
@@ -840 +848 @@
                                         codim2Act->setFacetNormal(n);
                                         delete n;                                       
-                                        /*intvec *n = new intvec(this->numVars);
-                                        makeInt(P,jj,*n);                                               
-                                        codim2Act->setFacetNormal(n);
-                                        this->facetPtr->numCodim2Facets++;      //Honor the creation of a codim-2-facet
-                                        codim2Act->next = new facet(2);
-                                        codim2Act = codim2Act->next;                                            
-                                        delete n;*/
-                                }                                                                               
+                                }               
+                                /*We check whether the facet spanned by the codim-2 facets
+                                * intersects with the positive orthant. Otherwise we define this
+                                * facet to be non-flippable
+                                */
+                                intvec *iv_intPoint = new intvec(this->numVars);
+                                dd_MatrixPtr shiftMatrix;
+                                dd_MatrixPtr intPointMatrix;
+                                shiftMatrix = dd_CreateMatrix(this->numVars,this->numVars+1);
+                                for(int kk=0;kk<this->numVars;kk++)
+                                {
+                                        dd_set_si(shiftMatrix->matrix[kk][0],1);
+                                        dd_set_si(shiftMatrix->matrix[kk][kk+1],1);
+                                }
+                                intPointMatrix=dd_MatrixAppend(ddakt,shiftMatrix);
+                                //dd_WriteMatrix(stdout,intPointMatrix);
+                                interiorPoint(intPointMatrix,*iv_intPoint);
+                                for(int ll=0;ll<this->numVars;ll++)
+                                {
+                                        if( (*iv_intPoint)[ll] < 0 )
+                                        {
+                                                fAct->isFlippable=FALSE;
+                                                break;
+                                        }
+                                }
+                                /*End of check*/                                        
                                 fAct = fAct->next;      
                                 dd_FreeMatrix(ddakt);
@@ -972 +998 @@
                         ina=idrCopyR(initialForm,srcRing);                      
 #ifdef gfan_DEBUG
-                        cout << "ina=";
-                        idShow(ina); cout << endl;
+//                      cout << "ina=";
+//                      idShow(ina); cout << endl;
 #endif
                         ideal H;
@@ -980 +1006 @@
                         idSkipZeroes(H);
 #ifdef gfan_DEBUG
-                        cout << "H="; idShow(H); cout << endl;
+//                      cout << "H="; idShow(H); cout << endl;
 #endif
                         /*Substep 2.2
@@ -990 +1016 @@
                         srcRing_H=idrCopyR(H,tmpRing);
 #ifdef gfan_DEBUG
-                        cout << "srcRing_H = ";
-                        idShow(srcRing_H); cout << endl;
+//                      cout << "srcRing_H = ";
+//                      idShow(srcRing_H); cout << endl;
 #endif
                         srcRing_HH=ffG(srcRing_H,this->gcBasis);                
 #ifdef gfan_DEBUG
-                        cout << "srcRing_HH = ";
-                        idShow(srcRing_HH); cout << endl;
+//                      cout << "srcRing_HH = ";
+//                      idShow(srcRing_HH); cout << endl;
 #endif
                         /*Substep 2.2.1
@@ -1019 +1045 @@
                         construction of the differences
                         */
-                        intPointMatrix = dd_CreateMatrix(iPMatrixRows+2,this->numVars+1);
+                        intPointMatrix = dd_CreateMatrix(iPMatrixRows+10,this->numVars+1);
                         intPointMatrix->numbtype=dd_Integer;    //NOTE: DO NOT REMOVE OR CHANGE TO dd_Rational
                         
@@ -1329 +1355 @@
                 
                 /** \brief Compute an interior point of a given cone
-                * Result will be written into intvec iv
+                * Result will be written into intvec iv. 
+                * Any rational point is automatically converted into an integer.
                 */
                 void interiorPoint(dd_MatrixPtr const &M, intvec &iv) //no const &M here since we want to remove redundant rows
@@ -1758 +1785 @@
                         
                         fAct = gcAct->facetPtr;                 
-                        gcAct->writeConeToFile(*gcAct);
+                        //gcAct->writeConeToFile(*gcAct);
                         
                         /*End of initialisation*/
@@ -1781 +1808 @@
                                         gcTmp->getCodim2Normals(*gcTmp);
                                         gcTmp->normalize();
-                                        gcTmp->writeConeToFile(*gcTmp);
+                                        //gcTmp->writeConeToFile(*gcTmp);
                                         /*add facets to SLA here*/
                                         SearchListRoot=gcTmp->enqueueNewFacets(*SearchListRoot);
-                                        gcTmp->showSLA(*SearchListRoot);
+                                        if(SearchListRoot!=NULL)
+                                                gcTmp->showSLA(*SearchListRoot);
                                         rChangeCurrRing(gcAct->baseRing);
                                         gcPtr->next=gcTmp;
@@ -1814 +1842 @@
                                 SearchListAct = SearchListRoot;
                         }
+                        cout << endl << "Found " << counter << " cones - terminating" << endl;
                 
                         //NOTE Hm, comment in and get a crash for free...
@@ -2038 +2067 @@
                                                                 {
                                                                         slHead = slAct->next;
-                                                                        slHead->prev = NULL;
+                                                                        if(slHead!=NULL)
+                                                                                slHead->prev = NULL;
                                                                         //set a bool flag to mark slAct as to be deleted
                                                                 }