Changeset f91fddc in git


Ignore:
Timestamp:
Apr 9, 2010, 9:20:11 AM (14 years ago)
Author:
Martin Monerjan
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', '648d28f488f6ff08f5607ff229b9ad9e4a5b93c2')
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:
2 edited

Legend:

Unmodified
Added
Removed
  • kernel/gfan.cc

    r807ee2 rf91fddc  
    7676//NOTE Defining this will slow things down!
    7777//Only good for very coarse profiling
    78 //#define gfanp
     78// #define gfanp
    7979#ifdef gfanp
    8080#include <sys/time.h>
     
    157157        while(f2Copy!=NULL)
    158158        {
    159                 if(f2Act==NULL || f2Act==(facet*)0xfefefefefefefefe)
     159                if(f2Act==NULL
     160#ifndef NDEBUG
     161  #if SIZEOF_LONG==8
     162                                 || f2Act==(facet*)0xfefefefefefefefe
     163  #elif SIZEOF_LONG==4
     164                                 || f2Act==(facet*)0xfefefefe
     165  #endif
     166#endif
     167                  )
    160168                {
    161169                        f2Act=new facet(2);
     
    273281                }
    274282        }
     283//      if( fIntP->compare(gIntP)!=0) res=FALSE;
    275284#ifdef gfanp
    276285        gettimeofday(&end, 0);
     
    840849        dd_rowset ddredrows;            // # of redundant rows in ddineq
    841850        dd_rowset ddlinset;             // the opposite
    842         dd_rowindex ddnewpos;             // all to make dd_Canonicalize happy
     851        dd_rowindex ddnewpos=NULL;                // all to make dd_Canonicalize happy
    843852        dd_NumberType ddnumb=dd_Integer; //Number type
    844853        dd_ErrorType dderr=dd_NoError;         
     
    894903        for(int ii=0;ii<ddineq->rowsize;ii++)
    895904        {
    896                 ideal initialForm=idInit(IDELEMS(I),1);
     905                ideal initialForm=idInit(IDELEMS(I),I->rank);
    897906                //read row ii into gamma
    898                 int64 tmp;
     907//              int64 tmp;
    899908                int64vec *gamma=new int64vec(this->numVars);
    900909                for(int jj=1;jj<=this->numVars;jj++)
    901910                {
     911                        int64 tmp;
    902912                        tmp=(int64)mpq_get_d(ddineq->matrix[ii][jj]);
    903913                        (*gamma)[jj-1]=(int64)tmp;
     
    909919                for(int jj=0;jj<IDELEMS(initialForm);jj++)
    910920                {
    911                         L->m[jj]=pCopy(pHead(initialForm->m[jj]));
     921                        poly p=pHead(initialForm->m[jj]);
     922                        L->m[jj]=pCopy(/*pHead(initialForm->m[jj]))*/p);
     923                        pDelete(&p);
    912924                }               
    913925               
     
    946958                                                        }
    947959                                                }
     960                                                pDelete(&q);
    948961                                                if(isMaybeFacet==TRUE)
    949                                                 {                                                       
     962                                                {
    950963                                                        break;//while(p!=NULL)
    951964                                                }
    952965                                                p=pNext(p);
    953                                                 pDelete(&q);                                           
    954966                                        }//while
    955967//                                      pDelete(&p);//NOTE Better to use pDel and qDel. Commenting in this line will not work!
    956                                         pDelete(&q);
     968                                        if(q!=NULL) pDelete(&q);
    957969                                        pDelete(&pDel);
    958970                                        pDelete(&qDel);
     
    960972                                        {
    961973                                                dd_set_si(ddineq->matrix[ii][0],1);                                             
    962                                                 if(num_alloc==0)
    963                                                         num_alloc += 1;
    964                                                 else
    965                                                         num_alloc += 1;
     974//                                              if(num_alloc==0)
     975//                                                      num_alloc += 1;
     976//                                              else                                           
     977//                                                      num_alloc += 1;
     978                                                if(num_alloc==num_elts) num_alloc==0 ? num_alloc=1 : num_alloc*=2;
     979                                               
    966980                                                void *tmp = realloc(redRowsArray,(num_alloc*sizeof(int)));
    967981                                                if(!tmp)
     
    9871001//      delete gamma;
    9881002        int offset=0;//needed for correction of redRowsArray[ii]
     1003#ifdef gfan_DEBUG
     1004        printf("Removed %i of %i in preprocessing step\n",num_elts,ddineq->rowsize);
     1005#endif
    9891006        for( int ii=0;ii<num_elts;ii++ )
    9901007        {
    9911008                dd_MatrixRowRemove(&ddineq,redRowsArray[ii]+1-offset);//cddlib sucks at enumeration
    9921009                offset++;
    993         }
    994         free(redRowsArray);
     1010        }       
     1011        free(redRowsArray);//NOTE May crash on some machines.
    9951012        /*And now for the strictly positive rows
    9961013        * Doesn't gain significant speedup
     
    10631080                for (int jj = 1; jj <ddcols; jj++)
    10641081                {
    1065                         double foo;
    1066                         foo = mpq_get_d(ddineq->matrix[kk][jj]);
    1067                         (*load)[jj-1] = (int64)foo;     //store typecasted entry at pos jj-1 of load
    1068                         ggT = intgcd(ggT,(int64&)foo);
     1082                        int64 val;
     1083                        val = (int64)mpq_get_d(ddineq->matrix[kk][jj]);
     1084                        (*load)[jj-1] = val;    //store typecasted entry at pos jj-1 of load
     1085                        ggT = int64gcd(ggT,/*(int64&)foo*/val);
    10691086                }//for (int jj = 1; jj <ddcols; jj++)
    10701087                if(ggT>1)
     
    11901207        }       
    11911208        //Clean up but don't delete the return value!   
    1192         dd_FreeMatrix(ddineq);
    1193         set_free(ddredrows);
    1194         set_free(ddlinset);
    1195         free(ddnewpos);
     1209        dd_FreeMatrix(ddineq); 
     1210        set_free(ddredrows);//check
     1211        set_free(ddlinset);//check
     1212        free(ddnewpos);//<-- NOTE Here the crash occurs omAlloc issue?
    11961213#ifdef gfanp
    11971214        gettimeofday(&end, 0);
     
    14781495//              mpq_clear(product);
    14791496                //Compute intgcd
    1480                 ggT=intgcd(ggT,(*ivIntPointOfCone)[ii]);
     1497                ggT=int64gcd(ggT,(*ivIntPointOfCone)[ii]);
    14811498        }
    14821499       
     
    15201537                int64 ggT=(*ivIntPointOfCone)[0];
    15211538                for(int ii=0;ii<this->numVars;ii++)
    1522                         ggT=intgcd( ggT, (*ivIntPointOfCone)[ii]);
     1539                        ggT=int64gcd( ggT, (*ivIntPointOfCone)[ii]);
    15231540                if(ggT>1)
    15241541                {
     
    16021619                int64 ggT=(*ivIntPointOfFacet)[0];
    16031620                for(int ii=0;ii<this->numVars;ii++)
    1604                         ggT=intgcd(ggT,(*ivIntPointOfFacet)[ii]);
     1621                        ggT=int64gcd(ggT,(*ivIntPointOfFacet)[ii]);
    16051622                if(ggT>1)
    16061623                {
     
    18551872                        pGetExpV(aktpoly,src_ExpV);
    18561873                        rChangeCurrRing(tmpRing);       //this ring change is crucial!
    1857                         pGetExpV(pCopy(H->m[ii]),dst_ExpV);
     1874                        poly p=pCopy(H->m[ii]);
     1875                        pGetExpV(p/*pCopy(H->m[ii])*/,dst_ExpV);
     1876                        pDelete(&p);
    18581877                        rChangeCurrRing(srcRing);
    18591878                        bool expVAreEqual=TRUE;
     
    22522271 * and in gcone::flip for obvious reasons.
    22532272*/
    2254 inline void gcone::computeInv(const ideal &gb, ideal &initialForm, const int64vec &fNormal)
     2273/*inline*/ void gcone::computeInv(const ideal &gb, ideal &initialForm, const int64vec &fNormal)
    22552274{
    22562275#ifdef gfanp
     
    22652284                pGetExpV(aktpoly,leadExpV);     //find the leading exponent in leadExpV[1],...,leadExpV[n], use pNext(p)
    22662285                initialFormElement=pHead(aktpoly);
     2286//              int *v=(int *)omAlloc((this->numVars+1)*sizeof(int));
    22672287                while(pNext(aktpoly)!=NULL)     /*loop trough terms and check for parallelity*/
    22682288                {
    22692289                        int64vec *check = new int64vec(this->numVars);
    22702290                        aktpoly=pNext(aktpoly); //next term
    2271 //                      pSetm(aktpoly);
    22722291                        int *v=(int *)omAlloc((this->numVars+1)*sizeof(int));
    22732292                        pGetExpV(aktpoly,v);           
    2274                         /* Convert (int)v into (int64vec)check */                       
     2293                        /* Convert (int)v into (int64vec)check */
     2294//                      bool notPar=FALSE;
    22752295                        for (int jj=0;jj<this->numVars;jj++)
    22762296                        {
    22772297                                (*check)[jj]=v[jj+1]-leadExpV[jj+1];
     2298//                              register int64 foo=(fNormal)[jj];
     2299//                              if( ( (*check)[jj]  == /*fNormal[jj]*/foo )
     2300//                               || ( (/*fNormal[jj]*/foo!=0) && ( ( (*check)[jj] % /*fNormal[jj]*/foo ) !=0 ) ) )
     2301//                              {
     2302//                                      notPar=TRUE;
     2303//                                      break;
     2304//                              }
    22782305                        }
    2279                         if (isParallel(*check,fNormal)) //pass *check when
    2280 //                      if(isParallel((const int64vec*)&check,(const int64vec*)&fNormal))
    2281 //                      if(fNormal.compare(check)==0)
    2282                         {
    2283                                 //Found a parallel vector. Add it
     2306                        omFree(v);
     2307                        if (isParallel(*check,fNormal))//Found a parallel vector. Add it
     2308//                      if(notPar==FALSE)
     2309                        {
    22842310                                initialFormElement = pAdd((initialFormElement),(poly)pHead(aktpoly));//pAdd = p_Add_q destroys args
    22852311                        }
    2286                         omFree(v);
    22872312                        delete check;
    22882313                }//while
     2314//              omFree(v);
    22892315#ifdef gfan_DEBUG
    22902316//              cout << "Initial Form=";                               
     
    24782504        gettimeofday(&start, 0);
    24792505#endif*/               
    2480         int lhs,rhs;
    24812506        bool res;
    2482         lhs=dotProduct(a,b)*dotProduct(a,b);
    2483         rhs=dotProduct(a,a)*dotProduct(b,b);
    2484                         //cout << "LHS="<<lhs<<", RHS="<<rhs<<endl;
     2507        int lhs=dotProduct(a,b)*dotProduct(a,b);
     2508        int rhs=dotProduct(a,a)*dotProduct(b,b);       
    24852509        if (lhs==rhs)
    24862510        {
     
    33723396                mpq_mul(res,qkgV,M->matrix[line-1][ii+1]);
    33733397                n[ii]=(int64)mpz_get_d(mpq_numref(res));
    3374 //              ggT=intgcd(ggT,n[ii]);
     3398//              ggT=int64gcd(ggT,n[ii]);
    33753399        }
    33763400        int64 ggT=n[0];
    33773401        for(int ii=0;ii<this->numVars;ii++)
    3378                 ggT=intgcd(ggT,n[ii]); 
     3402                ggT=int64gcd(ggT,n[ii]);       
    33793403        //Normalization
    33803404        if(ggT>1)
     
    38593883/** \brief Compute the gcd of two ints
    38603884 */
    3861 static int intgcd(const int64 &a, const int64 &b)
     3885static int64 int64gcd(const int64 &a, const int64 &b)
     3886{
     3887        int64 r, p=a, q=b;
     3888        if(p < 0)
     3889                p = -p;
     3890        if(q < 0)
     3891                q = -q;
     3892        while(q != 0)
     3893        {
     3894                r = p % q;
     3895                p = q;
     3896                q = r;
     3897        }
     3898        return p;
     3899}
     3900       
     3901static int intgcd(const int &a, const int &b)
    38623902{
    38633903        int r, p=a, q=b;
     
    38733913        }
    38743914        return p;
    3875 }               
     3915}
    38763916               
    38773917/** \brief Construct a dd_MatrixPtr from a cone's list of facets
     
    39714011                while(fAct!=NULL)
    39724012                {       
    3973                         const int64vec *iv;
    3974                         iv=fAct->getRef2FacetNormal();//->getFacetNormal();
     4013                        const int64vec *iv=fAct->getRef2FacetNormal();
     4014//                      iv=fAct->getRef2FacetNormal();//->getFacetNormal();
    39754015                        f2Act=fAct->codim2Ptr;
    39764016                        for (int ii=0;ii<iv->length();ii++)
     
    40094049        gcOutputFile.close();
    40104050        }
     4051        delete [] ringString;
    40114052                       
    40124053}//writeConeToFile(gcone const &gc)
  • kernel/gfan.h

    r807ee2 rf91fddc  
    251251};
    252252lists lprepareResult(gcone *gc, const int n);
    253 static int intgcd(const int64 &a, const int64 &b);
     253static int64 int64gcd(const int64 &a, const int64 &b);
     254static int intgcd(const int &a, const int &b);
    254255static int dotProduct(const int64vec &iva, const int64vec &ivb);
    255256static bool isParallel(const int64vec &a, const int64vec &b);
Note: See TracChangeset for help on using the changeset viewer.