Changeset 68eb0c in git


Ignore:
Timestamp:
Mar 6, 2010, 6:09:52 PM (13 years ago)
Author:
Martin Monerjan
Branches:
(u'spielwiese', '0d6b7fcd9813a1ca1ed4220cfa2b104b97a0a003')
Children:
5151ac95d6e9778dc832b31c28c9649d358d1b0c
Parents:
51346c4ca4b934458da6e3ad84176f094bfd5caa
Message:
getExtremalRays now stores the rays in intvec** gcone::gcRays, the
facets only point their fNormals below codim2Ptr there.
gcone::numRays
Cleanup of no longer used stuff: dotProduct, interiorPoint2


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

Legend:

Unmodified
Added
Removed

  • kernel/gfan.cc

    r51346c4 r68eb0c  
    7777//NOTE Defining this will slow things down!
    7878//Only good for very coarse profiling
    79 #define gfanp
     79// #define gfanp
    8080#ifdef gfanp
    8181#include <sys/time.h>
    … …  
    206206        this->next=NULL;
    207207        this->flipGB=NULL;
    208 //      delete this;
     208//      delete(this);
    209209}
    210210               
    … …  
    226226                        if(codim2Ptr->fNormal!=NULL)
    227227                        {
    228                                 delete codim2Ptr->fNormal;
     228                                delete codim2Ptr->fNormal;//NOTE Do not want this anymore since the rays are now in gcone!
    229229                                codim2Ptr = codim2Ptr->next;
    230230                        }
    … …  
    232232//              delete this->codim2Ptr;
    233233        }
    234         if(this->codim2Ptr!=NULL)
    235                 delete this->codim2Ptr;
     234        //The rays are stored in the cone!
     235//      if(this->codim2Ptr!=NULL)
     236//              delete this->codim2Ptr;
    236237        if(this->flipGB!=NULL)
    237238                idDelete((ideal *)&this->flipGB);
    … …  
    562563        this->numFacets=0;
    563564        this->ivIntPt=NULL;
     565        this->gcRays=NULL;
    564566}
    565567               
    … …  
    587589        this->numFacets=0;
    588590        this->ivIntPt=NULL;
     591        this->gcRays=NULL;
    589592}
    590593               
    … …  
    609612        this->numFacets=0;
    610613        this->ivIntPt=NULL;
     614        this->gcRays=NULL;
    611615}
    612616               
    … …  
    649653//      dd_FreeMatrix(this->ddFacets);
    650654        //dd_FreeMatrix(this->ddFacets);
     655        for(int ii=0;ii<this->numRays;ii++)
     656                delete(gcRays[ii]);
     657        omFree(gcRays);
    651658}                       
    652659
    … …  
    14351442        /* Compute interior point on the fly*/
    14361443        intvec *ivIntPointOfCone = new intvec(this->numVars);
    1437 //      for(int ii=0;ii<P->rowsize;ii++)
    1438 //      {}
    14391444        mpq_t *colSum = new mpq_t[this->numVars];
    14401445        int denom[this->numVars];//denominators of colSum
    … …  
    15281533        int rows=P->rowsize;
    15291534        facet *fAct=gc.facetPtr;
     1535        //Construct an array to hold the extremal rays of the cone
     1536        this->gcRays = (intvec**)omAlloc0(sizeof(intvec*)*P->rowsize);
     1537        for(int ii=0;ii<P->rowsize;ii++)
     1538        {
     1539                intvec *rowvec = new intvec(this->numVars);
     1540                makeInt(P,ii+1,*rowvec);//get an integer entry instead of rational, rowvec is primitve
     1541                this->gcRays[ii] = ivCopy(rowvec);
     1542                delete rowvec;
     1543        }
     1544        this->numRays=P->rowsize;
     1545        //Check which rays belong to which facet
    15301546        while(fAct!=NULL)
    15311547        {
    … …  
    15351551                for(int ii=0;ii<rows;ii++)
    15361552                {                       
    1537                         intvec *rowvec = new intvec(this->numVars);
    1538                         makeInt(P,ii+1,*rowvec);//get an integer entry instead of rational, rowvec is primitve                 
    1539                         if(dotProduct(*fNormal,*rowvec)==0)
     1553                        if(dotProduct(*fNormal,this->gcRays[ii])==0)
    15401554                        {
    15411555                                intvec *tmp = ivIntPointOfFacet;//Prevent memleak
    … …  
    15521566                                        codim2Act = codim2Act->next;
    15531567                                }
    1554                                 codim2Act->setFacetNormal(rowvec);
    1555                                 fAct->numRays++;
    1556                                 //TODO Add up to interior point here!
    1557                                 //Memleak sinve ivAdd returns a new intvec
    1558                                 ivIntPointOfFacet=ivAdd(ivIntPointOfFacet,rowvec);
    1559                                 //Now tmp still points to the OLD address of ivIntPt
     1568                                //codim2Act->setFacetNormal(rowvec);
     1569                                //Rather just let codim2Act point to the corresponding intvec of gcRays
     1570                                codim2Act->fNormal=this->gcRays[ii];
     1571                                fAct->numRays++;                                 
     1572                                //Memleak avoided via tmp
     1573                                ivIntPointOfFacet=ivAdd(ivIntPointOfFacet,this->gcRays[ii]);
     1574                                //Now tmp still points to the OLD address of ivIntPointOfFacet
    15601575                                delete(tmp);
    15611576                                       
    15621577                        }
    1563                         delete(rowvec);
    15641578                }//For non-homog input ivIntPointOfFacet should already be >0 here
    15651579                if(!hasHomInput) {assert(iv64isStrictlyPositive(ivIntPointOfFacet));}
    … …  
    15711585                        for(int ii=0;ii<this->numVars;ii++)
    15721586                                (*ivOne)[ii]=1;
    1573 //                      intvec *diff = new intvec(this->numVars);
    1574 //                      diff=ivSub(ivIntPointOfFacet,ivOne);
    15751587                        while( !iv64isStrictlyPositive(ivIntPointOfFacet) )
    15761588                        {
    … …  
    15791591                                {
    15801592                                        (*ivOne)[jj] = (*ivOne)[jj] << 1; //times 2
    1581 //                                      (*diff)[jj]= (*diff)[jj] << 1;
    15821593                                }
    15831594                                ivIntPointOfFacet = ivAdd(ivIntPointOfFacet/*diff*/,ivOne);
    15841595                                delete tmp;                             
    15851596                        }
    1586                         //fAct->setInteriorPoint(ivIntPointOfFacet);
    15871597                        delete ivOne;
    15881598                }
    1589                 //else
    15901599                int ggT=(*ivIntPointOfFacet)[0];
    15911600                for(int ii=0;ii<this->numVars;ii++)
    … …  
    16391648//              delete fNormal;         
    16401649                fAct = fAct->next;
    1641         }
     1650        }//end of facet checking
    16421651        dd_FreeMatrix(P);
    16431652        //Now all extremal rays should be set w.r.t their respective fNormal
    … …  
    24492458*
    24502459*/
    2451 // inline int gcone::dotProduct(intvec &iva, intvec &ivb)                               
    2452 // {                   
    2453 //      int res=0;
    2454 //      for (int i=0;i<this->numVars;i++)
    2455 //      {
    2456 //              res = res+(iva[i]*ivb[i]);
    2457 //      }
    2458 //      return res;
    2459 // }//int dotProduct
    24602460inline int gcone::dotProduct(const intvec &iva, const intvec &ivb)                             
    24612461{                       
    … …  
    24992499        return res;
    25002500}//bool isParallel
    2501 // inline int gcone::dotProduct(const intvec *a, const intvec *b)                               
    2502 // {                   
    2503 //      int res=0;
    2504 //      for (int i=0;i<this->numVars;i++)
    2505 //      {
    2506 //              res = res+((*a)[i]*(*b)[i]);
    2507 //      }
    2508 //      return res;
    2509 // }//int dotProduct
    2510 // inline bool gcone::isParallel(const intvec* a, const intvec* b)
    2511 // {
    2512 //      int lhs,rhs;
    2513 //      bool res;
    2514 //      lhs=dotProduct(a,b)*dotProduct(a,b);
    2515 //      rhs=dotProduct(a,a)*dotProduct(b,b);
    2516 //                      //cout << "LHS="<<lhs<<", RHS="<<rhs<<endl;
    2517 //      if (lhs==rhs)
    2518 //      {
    2519 //              res = TRUE;
    2520 //      }
    2521 //      else
    2522 //      {
    2523 //              res = FALSE;
    2524 //      }
    2525 //      return res;
    2526 // }
     2501
    25272502/** \brief Compute an interior point of a given cone
    25282503 * Result will be written into intvec iv.
    … …  
    27922767        mpz_clear(kgV);
    27932768        mpq_clear(qkgV); mpq_clear(res);
    2794        
    2795 /***************************************************/
    2796 //      //We need to find out if there are at least two pos rays
    2797 //      facet *f1 = this->facetPtr;
    2798 //      facet *f2 = NULL;
    2799 //      intvec *intF1=NULL;
    2800 //      intvec *intF2=NULL;
    2801 //      while(f1!=NULL)
    2802 //      {
    2803 //              if(f1->isFlippable)
    2804 //              {
    2805 //                      facet *f1Ray=f1->codim2Ptr;
    2806 //                      while(f1Ray!=NULL)
    2807 //                      {
    2808 //                              intvec *check = f1Ray->getFacetNormal();
    2809 //                              if(iv64isStrictlyPositive(check))
    2810 //                              {
    2811 //                                      intF1=ivCopy(check);
    2812 //                                      delete check;
    2813 //                                      break;
    2814 //                              }
    2815 //                              delete check;
    2816 //                              f1Ray = f1Ray->next;
    2817 //                      }       
    2818 //              }
    2819 //              if(intF1!=NULL)//We have found the first strictly positive ray
    2820 //                      break;
    2821 //              f1 = f1->next;
    2822 //      }
    2823 //      if(f1->next!=NULL)
    2824 //              f2=f1->next;
    2825 //      while(f2!=NULL)
    2826 //      {
    2827 //              if(f2->isFlippable)
    2828 //              {
    2829 //                      facet *f2Ray=f2->codim2Ptr;
    2830 //                      while(f2Ray!=NULL)
    2831 //                      {
    2832 //                              intvec *check = f2Ray->getFacetNormal();
    2833 //                              if(iv64isStrictlyPositive(check))
    2834 //                              {
    2835 //                                      //and not equal to intF1
    2836 //                                      if(check->compare(intF1)!=0)//we found a distinct 2nd ray
    2837 //                                      {
    2838 //                                              intF2=ivCopy(check);
    2839 //                                              delete check;
    2840 //                                      }
    2841 //                                      break;
    2842 //                              }
    2843 //                              delete check;
    2844 //                              f2Ray = f2Ray->next;
    2845 //                      }
    2846 //              }
    2847 //              if(intF2!=NULL)
    2848 //                      break;
    2849 //              f2=f2->next;
    2850 //      }
    2851 //      if(intF2==NULL)//do some linear algebra
    2852 //      {
    2853 //              f2=f1->next;
    2854 //              intvec *arbitraryRay = f2->codim2Ptr->getFacetNormal();
    2855 //              //Now add
    2856 //              intvec *ivSum = ivAdd(intF1, arbitraryRay);
    2857 //              while(iv64isStrictlyPositive(ivSum)==FALSE)
    2858 //              {
    2859 //                      intvec *tmp = intF1;
    2860 //                      for(int ii=0;ii<this->numVars;ii++)
    2861 //                              (*intF1)[ii] = (*intF1)[ii] << 1;       //times 2
    2862 //                      intF1 = ivAdd(intF1,ivSum);
    2863 //                      delete tmp;
    2864 //                     
    2865 //              }
    2866 //              this->setIntPoint(ivSum);
    2867 //
    2868 //      }
    2869 //      else    //just add intF1+intF2
    2870 //      {
    2871 //              intvec *sum = ivAdd(intF1, intF2);
    2872 //              this->setIntPoint(sum);
    2873 //              delete sum;
    2874 //      }
    2875 /***************************************************/
    2876 //      //Find 1st flippable facet
    2877 //      facet *f1 = this->facetPtr;
    2878 //      facet *f2 = NULL; //idPrint(this->gcBasis);
    2879 //      intvec *intF1=NULL;
    2880 //      intvec *intF2=NULL;
    2881 //      while(f1->isFlippable==FALSE && f1!=NULL)
    2882 //              f1=f1->next;
    2883 //      if(f1!=NULL)
    2884 //      {
    2885 //              intF1 = f1->getInteriorPoint();
    2886 //              f2 = f1->next;
    2887 //      }
    2888 //      while(f2!=NULL && f2->isFlippable==FALSE)
    2889 //              f2=f2->next;
    2890 //      if(f1==NULL || f2==NULL) //Is there only one flippable facet?
    2891 //      {       //f1==NULL would mean there ain't any flippable facet...
    2892 //              /*If f2==NULL we ignore f2 and try to find an int point
    2893 //              * using f1 and linear algebra
    2894 //              */
    2895 //              bool foundIntPoint=FALSE;
    2896 //              intvec *f1Ray=f1->codim2Ptr->getFacetNormal();
    2897 //              const intvec *fNormal=f1->getRef2FacetNormal();
    2898 //              while(foundIntPoint==FALSE)
    2899 //              {                       
    2900 //                      intvec *res;// = new intvec(this->numVars);
    2901 //                      res = ivAdd(const_cast<intvec*>(fNormal), f1Ray);
    2902 //                      int posCtr=0;   //Count the pos entries of res
    2903 //                      for(int ii=0;ii<this->numVars;ii++)
    2904 //                      {
    2905 //                              if( (*res)[ii]>0)
    2906 //                                      posCtr++;                       
    2907 //                      }                       
    2908 //                      if(posCtr==this->numVars)
    2909 //                      {
    2910 //                              foundIntPoint=TRUE;
    2911 //                              this->setIntPoint(res);
    2912 //                      }
    2913 //                      else
    2914 //                      {
    2915 //                              for(int ii=0;ii<this->numVars;ii++)
    2916 //                                      (*f1Ray)[ii] = ((*f1Ray)[ii]) << 1;
    2917 //                      }
    2918 //                      delete res;
    2919 //              }
    2920 //              //Now we should have an int point
    2921 //              delete f1Ray;
    2922 //      }
    2923 //      else    //ok, we have two flippable facets
    2924 //      {               
    2925 //              intF2=f2->getInteriorPoint();
    2926 //              mpq_t *ratIntPt = new mpq_t[this->numVars];
    2927 //              for(int ii=0;ii<this->numVars;ii++)
    2928 //                      mpq_init(ratIntPt[ii]);
    2929 //              for(int ii=0;ii<this->numVars;ii++)
    2930 //              {
    2931 //                      mpq_t a,b;
    2932 //                      mpq_init(a); mpq_init(b);
    2933 //                      mpq_set_si(a,(*intF1)[ii],1);
    2934 //                      mpq_set_si(b,(*intF2)[ii],1);
    2935 //                      mpq_t diff;
    2936 //                      mpq_init(diff);
    2937 //                      mpq_sub(diff,a,b);      //diff=a-b
    2938 //                      mpq_t quot;
    2939 //                      mpq_init(quot);
    2940 //                      mpq_div_2exp(quot,diff,1);      //quot=diff/2=(a-b)/2
    2941 //                      mpq_clear(diff);
    2942 //                      //Don't be clever and reuse diff here
    2943 //                      mpq_t sum; mpq_init(sum);
    2944 //                      mpq_add(sum,b,quot);    //sum=b+quot=b+(a-b)/2
    2945 //                      mpq_set(ratIntPt[ii],sum);
    2946 //                      mpq_clear(sum);
    2947 //                      mpq_clear(quot);
    2948 //                      mpq_clear(a); mpq_clear(b);
    2949 //              }       
    2950 //              /*Compute lcm of the denominators*/
    2951 //              mpz_t *denom = new mpz_t[this->numVars];
    2952 //              mpz_t tmp,kgV;
    2953 //              mpz_init(tmp); mpz_init(kgV);
    2954 //              for (int ii=0;ii<this->numVars;ii++)
    2955 //              {
    2956 //                      mpz_t z;
    2957 //                      mpz_init(z);
    2958 //                      mpq_get_den(z,ratIntPt[ii]);
    2959 //                      mpz_init(denom[ii]);
    2960 //                      mpz_set( denom[ii], z);
    2961 //                      mpz_clear(z);                           
    2962 //              }
    2963 //             
    2964 //              mpz_set(tmp,denom[0]);
    2965 //              for (int ii=0;ii<this->numVars;ii++)
    2966 //              {
    2967 //                      mpz_lcm(kgV,tmp,denom[ii]);
    2968 //                      mpz_set(tmp,kgV);                               
    2969 //              }
    2970 //              mpz_clear(tmp);
    2971 //              /*Multiply the nominators by kgV*/
    2972 //              mpq_t qkgV,res;
    2973 //              mpq_init(qkgV);
    2974 //              /*mpq_set_str(qkgV,"1",10);*/
    2975 //              mpq_canonicalize(qkgV);
    2976 //                             
    2977 //              mpq_init(res);
    2978 //              /*mpq_set_str(res,"1",10);*/
    2979 //              mpq_canonicalize(res);
    2980 //                             
    2981 //              mpq_set_num(qkgV,kgV);
    2982 //              intvec *n=new intvec(this->numVars);
    2983 //              for (int ii=0;ii<this->numVars;ii++)
    2984 //              {
    2985 //                      mpq_canonicalize(ratIntPt[ii]);mpq_mul(res,qkgV,ratIntPt[ii]);
    2986 //                      (*n)[ii]=(int)mpz_get_d(mpq_numref(res));
    2987 //              }
    2988 //              this->setIntPoint(n);
    2989 //              delete n;
    2990 //              delete [] ratIntPt;
    2991 //              delete [] denom;
    2992 //              mpz_clear(kgV);
    2993 //              mpq_clear(qkgV); mpq_clear(res);
    2994 //      }
    29952769}
    29962770       
    … …  
    39623736                                        cout << "Removing (";fAct->fNormal->show(1,1);cout << ") from list" << endl;
    39633737#endif
    3964                                         marker->shallowDelete();
     3738//                                      marker->shallowDelete();
    39653739//                                      delete(marker);
    39663740                                        break;

  • kernel/gfan.h

    r51346c4 r68eb0c  
    7575                facet *codim2Ptr;       //Pointer to (codim-2)-facet. Bit of recursion here ;-)
    7676                int numCodim2Facets;    //#of (codim-2)-facets of this facet. Set in getCodim2Normals()
    77                 unsigned numRays;       //Number of spanning rays of the cone
     77                unsigned numRays;       //Number of spanning rays of the facet
    7878                ring flipRing;          //the ring on the other side of the facet
    7979                                       
    … …  
    8484                /**  The copy constructor */
    8585                facet(const facet& f);
     86                /** A shallow copy of facets*/
    8687                facet* shallowCopy(const facet& f);
    8788                void shallowDelete();
    … …  
    9495                /** Stores the facet normal \param intvec*/
    9596                inline void setFacetNormal(intvec *iv);
    96                 /** Hopefully returns the facet normal */
     97                /** Returns the facet normal */
    9798                inline intvec *getFacetNormal();
    9899                /** Return a reference to the facet normal*/
    … …  
    189190                dd_MatrixPtr ddFacets;  //Matrix to store irredundant facets of the cone
    190191               
     192                /** Array of intvecs representing the rays of the cone*/
     193                intvec **gcRays;
     194                unsigned numRays;       //#rays of the cone
    191195                /** Contains the Groebner basis of the cone. Is set by gcone::getGB(ideal I)*/
    192196                ideal gcBasis;          //GB of the cone, set by gcone::getGB();
Note: See TracChangeset for help on using the changeset viewer.