Changeset f91fddc in git

Timestamp:

Apr 9, 2010, 9:20:11 AM (13 years ago)

Author:

Martin Monerjan

Branches:

(u'spielwiese', '0d6b7fcd9813a1ca1ed4220cfa2b104b97a0a003')

Children:

cacfb67cf6bf7cd84cd8e725d9acf88690092d61

Parents:

807ee2919f6e854d2cfa4def4e5f3c99b8883a2d

Message:

Back to 64bit - #ifdef SIZEOF_LONG where required
Expanded pCopy(pHead so as to not to yield a memleak in preprocessing
step
Marked position of free-bug at the end of getConeNormals
Tried another way of computing inv, much slower
int64gcd
delete [] ringString


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

Location:

kernel

Files:

• gfan.cc (modified) (23 diffs)
• gfan.h (modified) (1 diff)

kernel/gfan.cc

@@ -76 +76 @@
 //NOTE Defining this will slow things down!
 //Only good for very coarse profiling
-//#define gfanp
+// #define gfanp
 #ifdef gfanp
 #include <sys/time.h>
@@ -157 +157 @@
         while(f2Copy!=NULL)
         {
-                if(f2Act==NULL || f2Act==(facet*)0xfefefefefefefefe)
+                if(f2Act==NULL
+#ifndef NDEBUG
+  #if SIZEOF_LONG==8
+                                 || f2Act==(facet*)0xfefefefefefefefe
+  #elif SIZEOF_LONG==4
+                                 || f2Act==(facet*)0xfefefefe
+  #endif
+#endif
+                  )
                 {
                         f2Act=new facet(2);
@@ -273 +281 @@
                 }
         }
+//      if( fIntP->compare(gIntP)!=0) res=FALSE;
 #ifdef gfanp
         gettimeofday(&end, 0);
@@ -840 +849 @@
         dd_rowset ddredrows;            // # of redundant rows in ddineq
         dd_rowset ddlinset;             // the opposite
-        dd_rowindex ddnewpos;             // all to make dd_Canonicalize happy
+        dd_rowindex ddnewpos=NULL;                // all to make dd_Canonicalize happy
         dd_NumberType ddnumb=dd_Integer; //Number type
         dd_ErrorType dderr=dd_NoError;          
@@ -894 +903 @@
         for(int ii=0;ii<ddineq->rowsize;ii++)
         {
-                ideal initialForm=idInit(IDELEMS(I),1);
+                ideal initialForm=idInit(IDELEMS(I),I->rank);
                 //read row ii into gamma
-                int64 tmp;
+//              int64 tmp;
                 int64vec *gamma=new int64vec(this->numVars);
                 for(int jj=1;jj<=this->numVars;jj++)
                 {
+                        int64 tmp;
                         tmp=(int64)mpq_get_d(ddineq->matrix[ii][jj]);
                         (*gamma)[jj-1]=(int64)tmp;
@@ -909 +919 @@
                 for(int jj=0;jj<IDELEMS(initialForm);jj++)
                 {
-                        L->m[jj]=pCopy(pHead(initialForm->m[jj]));
+                        poly p=pHead(initialForm->m[jj]);
+                        L->m[jj]=pCopy(/*pHead(initialForm->m[jj]))*/p);
+                        pDelete(&p);
                 }               
                 
@@ -946 +958 @@
                                                         }
                                                 }
+                                                pDelete(&q);
                                                 if(isMaybeFacet==TRUE)
-                                                {                                                       
+                                                {
                                                         break;//while(p!=NULL)
                                                 }
                                                 p=pNext(p);
-                                                pDelete(&q);                                            
                                         }//while
 //                                      pDelete(&p);//NOTE Better to use pDel and qDel. Commenting in this line will not work!
-                                        pDelete(&q);
+                                        if(q!=NULL) pDelete(&q);
                                         pDelete(&pDel);
                                         pDelete(&qDel);
@@ -960 +972 @@
                                         {
                                                 dd_set_si(ddineq->matrix[ii][0],1);                                             
-                                                if(num_alloc==0)
-                                                        num_alloc += 1;
-                                                else
-                                                        num_alloc += 1;
+//                                              if(num_alloc==0)
+//                                                      num_alloc += 1;
+//                                              else                                            
+//                                                      num_alloc += 1;
+                                                if(num_alloc==num_elts) num_alloc==0 ? num_alloc=1 : num_alloc*=2;
+                                                
                                                 void *tmp = realloc(redRowsArray,(num_alloc*sizeof(int)));
                                                 if(!tmp)
@@ -987 +1001 @@
 //      delete gamma;
         int offset=0;//needed for correction of redRowsArray[ii]
+#ifdef gfan_DEBUG
+        printf("Removed %i of %i in preprocessing step\n",num_elts,ddineq->rowsize);
+#endif
         for( int ii=0;ii<num_elts;ii++ )
         {
                 dd_MatrixRowRemove(&ddineq,redRowsArray[ii]+1-offset);//cddlib sucks at enumeration
                 offset++;
-        }
-        free(redRowsArray);
+        }       
+        free(redRowsArray);//NOTE May crash on some machines.
         /*And now for the strictly positive rows
         * Doesn't gain significant speedup
@@ -1063 +1080 @@
                 for (int jj = 1; jj <ddcols; jj++)
                 {
-                        double foo;
-                        foo = mpq_get_d(ddineq->matrix[kk][jj]);
-                        (*load)[jj-1] = (int64)foo;     //store typecasted entry at pos jj-1 of load
-                        ggT = intgcd(ggT,(int64&)foo);
+                        int64 val;
+                        val = (int64)mpq_get_d(ddineq->matrix[kk][jj]);
+                        (*load)[jj-1] = val;    //store typecasted entry at pos jj-1 of load
+                        ggT = int64gcd(ggT,/*(int64&)foo*/val);
                 }//for (int jj = 1; jj <ddcols; jj++)
                 if(ggT>1)
@@ -1190 +1207 @@
         }       
         //Clean up but don't delete the return value!   
-        dd_FreeMatrix(ddineq);
-        set_free(ddredrows);
-        set_free(ddlinset);
-        free(ddnewpos);
+        dd_FreeMatrix(ddineq);  
+        set_free(ddredrows);//check
+        set_free(ddlinset);//check
+        free(ddnewpos);//<-- NOTE Here the crash occurs omAlloc issue?
 #ifdef gfanp
         gettimeofday(&end, 0);
@@ -1478 +1495 @@
 //              mpq_clear(product);
                 //Compute intgcd
-                ggT=intgcd(ggT,(*ivIntPointOfCone)[ii]);
+                ggT=int64gcd(ggT,(*ivIntPointOfCone)[ii]);
         }
         
@@ -1520 +1537 @@
                 int64 ggT=(*ivIntPointOfCone)[0];
                 for(int ii=0;ii<this->numVars;ii++)
-                        ggT=intgcd( ggT, (*ivIntPointOfCone)[ii]);
+                        ggT=int64gcd( ggT, (*ivIntPointOfCone)[ii]);
                 if(ggT>1)
                 {
@@ -1602 +1619 @@
                 int64 ggT=(*ivIntPointOfFacet)[0];
                 for(int ii=0;ii<this->numVars;ii++)
-                        ggT=intgcd(ggT,(*ivIntPointOfFacet)[ii]);
+                        ggT=int64gcd(ggT,(*ivIntPointOfFacet)[ii]);
                 if(ggT>1)
                 {
@@ -1855 +1872 @@
                         pGetExpV(aktpoly,src_ExpV);
                         rChangeCurrRing(tmpRing);       //this ring change is crucial!
-                        pGetExpV(pCopy(H->m[ii]),dst_ExpV);
+                        poly p=pCopy(H->m[ii]);
+                        pGetExpV(p/*pCopy(H->m[ii])*/,dst_ExpV);
+                        pDelete(&p);
                         rChangeCurrRing(srcRing);
                         bool expVAreEqual=TRUE;
@@ -2252 +2271 @@
  * and in gcone::flip for obvious reasons.
 */
-inline void gcone::computeInv(const ideal &gb, ideal &initialForm, const int64vec &fNormal)
+/*inline*/ void gcone::computeInv(const ideal &gb, ideal &initialForm, const int64vec &fNormal)
 {
 #ifdef gfanp
@@ -2265 +2284 @@
                 pGetExpV(aktpoly,leadExpV);     //find the leading exponent in leadExpV[1],...,leadExpV[n], use pNext(p)
                 initialFormElement=pHead(aktpoly);
+//              int *v=(int *)omAlloc((this->numVars+1)*sizeof(int));
                 while(pNext(aktpoly)!=NULL)     /*loop trough terms and check for parallelity*/
                 {
                         int64vec *check = new int64vec(this->numVars);
                         aktpoly=pNext(aktpoly); //next term
-//                      pSetm(aktpoly);
                         int *v=(int *)omAlloc((this->numVars+1)*sizeof(int));
                         pGetExpV(aktpoly,v);            
-                        /* Convert (int)v into (int64vec)check */                       
+                        /* Convert (int)v into (int64vec)check */
+//                      bool notPar=FALSE;
                         for (int jj=0;jj<this->numVars;jj++)
                         {
                                 (*check)[jj]=v[jj+1]-leadExpV[jj+1];
+//                              register int64 foo=(fNormal)[jj];
+//                              if( ( (*check)[jj]  == /*fNormal[jj]*/foo ) 
+//                               || ( (/*fNormal[jj]*/foo!=0) && ( ( (*check)[jj] % /*fNormal[jj]*/foo ) !=0 ) ) ) 
+//                              {
+//                                      notPar=TRUE;
+//                                      break;
+//                              }
                         }
-                        if (isParallel(*check,fNormal)) //pass *check when 
-//                      if(isParallel((const int64vec*)&check,(const int64vec*)&fNormal))
-//                      if(fNormal.compare(check)==0)
-                        {
-                                //Found a parallel vector. Add it
+                        omFree(v);
+                        if (isParallel(*check,fNormal))//Found a parallel vector. Add it
+//                      if(notPar==FALSE)
+                        {
                                 initialFormElement = pAdd((initialFormElement),(poly)pHead(aktpoly));//pAdd = p_Add_q destroys args
                         }
-                        omFree(v);
                         delete check;
                 }//while
+//              omFree(v);
 #ifdef gfan_DEBUG
 //              cout << "Initial Form=";                                
@@ -2478 +2504 @@
         gettimeofday(&start, 0);
 #endif*/                
-        int lhs,rhs;
         bool res;
-        lhs=dotProduct(a,b)*dotProduct(a,b);
-        rhs=dotProduct(a,a)*dotProduct(b,b);
-                        //cout << "LHS="<<lhs<<", RHS="<<rhs<<endl;
+        int lhs=dotProduct(a,b)*dotProduct(a,b);
+        int rhs=dotProduct(a,a)*dotProduct(b,b);        
         if (lhs==rhs)
         {
@@ -3372 +3396 @@
                 mpq_mul(res,qkgV,M->matrix[line-1][ii+1]);
                 n[ii]=(int64)mpz_get_d(mpq_numref(res));
-//              ggT=intgcd(ggT,n[ii]);
+//              ggT=int64gcd(ggT,n[ii]);
         }
         int64 ggT=n[0];
         for(int ii=0;ii<this->numVars;ii++)
-                ggT=intgcd(ggT,n[ii]);  
+                ggT=int64gcd(ggT,n[ii]);        
         //Normalization
         if(ggT>1)
@@ -3859 +3883 @@
 /** \brief Compute the gcd of two ints
  */
-static int intgcd(const int64 &a, const int64 &b)
+static int64 int64gcd(const int64 &a, const int64 &b)
+{
+        int64 r, p=a, q=b;
+        if(p < 0)
+                p = -p;
+        if(q < 0)
+                q = -q;
+        while(q != 0)
+        {
+                r = p % q;
+                p = q;
+                q = r;
+        }
+        return p;
+}
+        
+static int intgcd(const int &a, const int &b)
 {
         int r, p=a, q=b;
@@ -3873 +3913 @@
         }
         return p;
-}               
+}
                 
 /** \brief Construct a dd_MatrixPtr from a cone's list of facets
@@ -3971 +4011 @@
                 while(fAct!=NULL)
                 {       
-                        const int64vec *iv;
-                        iv=fAct->getRef2FacetNormal();//->getFacetNormal();
+                        const int64vec *iv=fAct->getRef2FacetNormal();
+//                      iv=fAct->getRef2FacetNormal();//->getFacetNormal();
                         f2Act=fAct->codim2Ptr;
                         for (int ii=0;ii<iv->length();ii++)
@@ -4009 +4049 @@
         gcOutputFile.close();
         }
+        delete [] ringString;
                         
 }//writeConeToFile(gcone const &gc)

kernel/gfan.h

@@ -251 +251 @@
 };
 lists lprepareResult(gcone *gc, const int n);
-static int intgcd(const int64 &a, const int64 &b);
+static int64 int64gcd(const int64 &a, const int64 &b);
+static int intgcd(const int &a, const int &b);
 static int dotProduct(const int64vec &iva, const int64vec &ivb);
 static bool isParallel(const int64vec &a, const int64vec &b);