Changeset 348034 in git


Ignore:
Timestamp:
May 6, 2009, 11:54:23 AM (14 years ago)
Author:
Martin Monerjan
Branches:
(u'spielwiese', '0d6b7fcd9813a1ca1ed4220cfa2b104b97a0a003')
Children:
c0180ff713d540e3f2e22d691024ebd3cf49c0a9
Parents:
946904abc570055ed6750cb3ecf70fd53d22275b
Message:
Partially fixed ring construction.
However now H=kStd(ina,...) does no longer work as expected and returns
the GB of the srcRing.
So apparently there is something going wrong with the ring change


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

Legend:

Unmodified
Added
Removed

  • kernel/gfan.cc

    r946904 r348034  
    22Compute the Groebner fan of an ideal
    33$Author: monerjan $
    4 $Date: 2009-05-05 15:15:02 $
    5 $Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.45 2009-05-05 15:15:02 monerjan Exp $
    6 $Id: gfan.cc,v 1.45 2009-05-05 15:15:02 monerjan Exp $
     4$Date: 2009-05-06 09:54:23 $
     5$Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.46 2009-05-06 09:54:23 monerjan Exp $
     6$Id: gfan.cc,v 1.46 2009-05-06 09:54:23 monerjan Exp $
    77*/
    88
    … …  
    1414#include "kutil.h"
    1515#include "intvec.h"
     16#include "int64vec.h"
    1617#include "polys.h"
    1718#include "ideals.h"
    … …  
    165166                                return TRUE;
    166167                }*/                                             
    167                 }//bool isFlippable(facet &f)
     168                //}//bool isFlippable(facet &f)
    168169               
    169170               
    … …  
    599600                        */
    600601                        ring srcRing=currRing;
    601                        
    602                         /* copied and modified from ring.cc::rAssureSyzComp */
    603                         //ring tmpRing=rCopyAndChangeWeight(srcRing,fNormal);
    604                         ring tmpRing=rCopy0(srcRing);
    605                         int i=rBlocks(srcRing);
    606                         int j;
    607                         tmpRing->order=(int *)omAlloc((i+1)*sizeof(int));
    608                         /*NOTE This should probably be set, but as of now crashes Singular*/
    609                         /*tmpRing->block0=(int *)omAlloc0((i+1)*sizeof(int));
    610                         tmpRing->block1=(int *)omAlloc0((i+1)*sizeof(int));*/
    611                         tmpRing->wvhdl[0] =( int *)omAlloc((fNormal->length())*sizeof(int)); //found in Singular/ipshell.cc
    612                         for(j=i;j>0;j--)
    613                         {
    614                                 tmpRing->order[j]=srcRing->order[j-1];
    615                                 tmpRing->block0[j]=srcRing->block0[j-1];
    616                                 tmpRing->block1[j]=srcRing->block1[j-1];
    617                                 if (srcRing->wvhdl[j-1] != NULL)
    618                                 {
    619                                         tmpRing->wvhdl[j] = (int*) omMemDup(srcRing->wvhdl[j-1]);
    620                                 }
    621                         }
    622                         tmpRing->order[0]=ringorder_a;
    623                         tmpRing->order[1]=ringorder_dp;
    624                         tmpRing->order[2]=ringorder_C;
    625                         //tmpRing->wvhdl[0] =( int *)omAlloc((fNormal->length())*sizeof(int)); //found in Singular/ipshell.cc
    626                        
    627                         for (int ii=0;ii<this->numVars;ii++)
    628                         {                               
    629                                 tmpRing->wvhdl[0][ii]=-(*fNormal)[ii];  //What exactly am I doing here?
    630                                 //cout << tmpring->wvhdl[0][ii] << endl;
    631                         }
    632                         rComplete(tmpRing);
     602
     603                        ring tmpRing=rCopyAndAddWeight(srcRing,fNormal);
    633604                        rChangeCurrRing(tmpRing);
    634605                       
    635606                        rWrite(currRing); cout << endl;
     607                       
    636608                        ideal ina;                     
    637609                        ina=idrCopyR(initialForm,srcRing);                     
    … …  
    640612                        idShow(ina); cout << endl;
    641613#endif
    642                         ideal H;
     614                        ideal H;  //NOTE WTF? rCopyAndAddWeight does not seem to be compatible with kStd... => wrong result!
    643615                        H=kStd(ina,NULL,isHomog,NULL);  //we know it is homogeneous
    644616                        idSkipZeroes(H);
    … …  
    702674                                        for (int kk=1;kk<=this->numVars;kk++)
    703675                                        {
    704                                                 cout << src_ExpV[kk] << "," << dst_ExpV[kk] << endl;
     676#ifdef gfan_DEBUG
     677                                                //cout << src_ExpV[kk] << "," << dst_ExpV[kk] << endl;
     678#endif
    705679                                                if (src_ExpV[kk]!=dst_ExpV[kk])
    706680                                                {
    … …  
    712686                                        {
    713687                                                markingsAreCorrect=TRUE; //everything is fine
    714                                                 cout << "correct markings" << endl;
     688#ifdef gfan_DEBUG
     689//                                              cout << "correct markings" << endl;
     690#endif
    715691                                        }//if (pHead(aktpoly)==pHead(H->m[jj])
    716692                                        delete src_ExpV;
    … …  
    772748                        turn the minimal basis into a reduced one
    773749                        */
     750                        int i,j;
    774751                        ring dstRing=rCopy0(srcRing);
    775752                        i=rBlocks(srcRing);
    … …  
    796773                        }
    797774                        rComplete(dstRing);
     775                       
     776                        // NOTE May assume that at this point srcRing already has 3 blocks of orderins, starting with a
     777                        // Thus:
     778                        //ring dstRing=rCopyAndChangeWeight(srcRing,iv_weight);
     779                        //cout << "PLING" << endl;
     780                        /*ring dstRing=rCopy0(srcRing);
     781                        for (int ii=0;ii<this->numVars;ii++)
     782                        {                               
     783                                dstRing->wvhdl[0][ii]=(*iv_weight)[ii];                         
     784                        }*/
    798785                        rChangeCurrRing(dstRing);
    799786#ifdef gfan_DEBUG
    … …  
    848835                                {                                       
    849836                                        if (pDivisibleBy(I->m[ii-1],p)) //does LM(I->m[ii]) divide LM(p) ?
    850                                         {
    851                                                 //cout << "TICK 3" << endl;
     837                                        {                                               
    852838                                                poly step1,step2,step3;
    853839                                                //cout << "dividing "; pWrite(pHead(p));cout << "by ";pWrite(pHead(I->m[ii-1])); cout << endl;
    854840                                                step1 = pDivideM(pHead(p),pHead(I->m[ii-1]));
    855                                                 //cout << "LT(p)/LT(f_i)="; pWrite(step1); cout << endl;
    856                                                 //cout << "TICK 3.1" << endl;
    857                                                 step2 = ppMult_qq(step1, I->m[ii-1]);
    858                                                 //cout << "TICK 3.2" << endl;
     841                                                //cout << "LT(p)/LT(f_i)="; pWrite(step1); cout << endl;                               
     842                                                step2 = ppMult_qq(step1, I->m[ii-1]);                                           
    859843                                                step3 = pSub(pCopy(p), step2);
    860                                                 //p=pSub(p,pMult( pDivide(pHead(p),pHead(I->m[ii])), I->m[ii]));
    861                                                 //cout << "TICK 4" << endl;
     844                                                //p=pSub(p,pMult( pDivide(pHead(p),pHead(I->m[ii])), I->m[ii]));                       
    862845                                                //pSetm(p);
    863846                                                pSort(step3); //must be here, otherwise strange behaviour with many +o+o+o+o+ terms
    … …  
    878861                                        pSort(r);
    879862                                        //cout << "r="; pWrite(r); cout << endl;
    880                                         //cout << "TICK 6" << endl;
     863                                       
    881864                                        if (pLength(p)!=1)
    882865                                        {
    … …  
    886869                                        {
    887870                                                p=NULL; //Hack to correct this situation                                               
    888                                         }
    889                                         //cout << "TICK 7" << endl;
     871                                        }                                       
    890872                                        //cout << "p="; pWrite(p);
    891873                                }//if (divOccured==FALSE)
    … …  
    1000982                        if (err!=dd_NoError){cout << "Error during dd_Matrix2LP in gcone::interiorPoint" << endl;}                     
    1001983                        if (lp==NULL){cout << "LP is NULL" << endl;}
    1002                         dd_WriteLP(stdout,lp);
    1003                         //cout << "Tick 3" << endl;
    1004                                                
     984#ifdef gfan_DEBUG
     985//                      dd_WriteLP(stdout,lp);
     986#endif 
     987                                       
    1005988                        lpInt=dd_MakeLPforInteriorFinding(lp);
    1006989                        if (err!=dd_NoError){cout << "Error during dd_MakeLPForInteriorFinding in gcone::interiorPoint" << endl;}
    1007                         dd_WriteLP(stdout,lpInt);
    1008                         //cout << "Tick 4" << endl;
    1009                        
     990#ifdef gfan_DEBUG
     991//                      dd_WriteLP(stdout,lpInt);
     992#endif                 
     993
    1010994                        dd_FindRelativeInterior(M,&ddlinset,&ddredrows,&lpSol,&err);
    1011995                        if (err!=dd_NoError)
    … …  
    10411025                }//void interiorPoint(dd_MatrixPtr const &M)
    10421026               
    1043                 ring rCopyAndChangeWeight(ring const r, intvec const *ivw)                             
    1044                 {
    1045                         ring res=rCopy0(r, FALSE, FALSE);
     1027                /** \brief Copy a ring and add a weighted ordering in first place
     1028                * Kudos to walkSupport.cc
     1029                */
     1030                ring rCopyAndAddWeight(ring const &r, intvec const *ivw)                               
     1031                {
     1032                        ring res=(ring)omAllocBin(ip_sring_bin);
     1033                        memcpy4(res,r,sizeof(ip_sring));
     1034                        res->VarOffset = NULL;
     1035                        res->ref=0;
     1036                       
     1037                        if (r->algring!=NULL)
     1038                                r->algring->ref++;
     1039                        if (r->parameter!=NULL)
     1040                        {
     1041                                res->minpoly=nCopy(r->minpoly);
     1042                                int l=rPar(r);
     1043                                res->parameter=(char **)omAlloc(l*sizeof(char_ptr));
     1044                                int i;
     1045                                for(i=0;i<rPar(r);i++)
     1046                                {
     1047                                        res->parameter[i]=omStrDup(r->parameter[i]);
     1048                                }
     1049                        }
     1050                       
    10461051                        int i=rBlocks(r);
    1047                         int j;
    1048 
    1049                         res->order=(int *)omAlloc((i+1)*sizeof(int));
    1050                         /*res->block0=(int *)omAlloc0((i+1)*sizeof(int));
    1051                         res->block1=(int *)omAlloc0((i+1)*sizeof(int));
    1052                         int ** wvhdl =(int **)omAlloc0((i+1)*sizeof(int**));*/
    1053                         for(j=i;j>0;j--)
    1054                         {
    1055                                 res->order[j]=r->order[j-1];
    1056                                 res->block0[j]=r->block0[j-1];
    1057                                 res->block1[j]=r->block1[j-1];
    1058                                 if (r->wvhdl[j-1] != NULL)
    1059                                 {
    1060                                         res->wvhdl[j] = (int*) omMemDup(r->wvhdl[j-1]);
    1061                                 }
    1062                         }
     1052                        int jj;
     1053                       
     1054                        res->order =(int *)omAlloc((i+1)*sizeof(int));
     1055                        res->block0=(int *)omAlloc((i+1)*sizeof(int));
     1056                        res->block1=(int *)omAlloc((i+1)*sizeof(int));
     1057                        res->wvhdl =(int **)omAlloc((i+1)*sizeof(int**));
     1058                        for(jj=0;jj<i;jj++)
     1059                        {                               
     1060                                if (r->wvhdl[jj] != NULL)
     1061                                {
     1062                                        res->wvhdl[jj] = (int*) omMemDup(r->wvhdl[jj-1]);
     1063                                }
     1064                                else
     1065                                {
     1066                                        res->wvhdl[jj+1]=NULL;
     1067                                }
     1068                        }
     1069                       
     1070                        for (jj=0;jj<i;jj++)
     1071                        {
     1072                                res->order[jj+1]=r->order[jj];
     1073                                res->block0[jj+1]=r->block0[jj];
     1074                                res->block1[jj+1]=r->block1[jj];
     1075                        }                                               
     1076                       
    10631077                        res->order[0]=ringorder_a;
    1064                         res->order[1]=ringorder_dp;
    1065                         res->order[2]=ringorder_C;
    1066                         //res->wvhdl = wvhdl;
    1067                        
    1068                         res->wvhdl[0] =( int *)omAlloc((ivw->length())*sizeof(int));
    1069                         for (int ii=0;ii<this->numVars;ii++)
     1078                        res->order[1]=ringorder_dp;     //basically useless, since that should never be used                   
     1079                        int length=ivw->length();
     1080                        int *A=(int *)omAlloc(length*sizeof(int));
     1081                        for (jj=0;jj<length;jj++)
    10701082                        {                               
    1071                                 res->wvhdl[0][ii]=(*ivw)[ii];                           
     1083                                A[jj]=(*ivw)[jj];                               
    10721084                        }                       
    1073                        
     1085                        res->wvhdl[0]=(int *)A;
     1086                        res->block0[0]=1;
     1087                        res->block1[0]=length;
     1088                       
     1089                        res->names = (char **)omAlloc0(rVar(r) * sizeof(char_ptr));
     1090                        for (i=rVar(res)-1;i>=0; i--)
     1091                        {
     1092                                res->names[i] = omStrDup(r->names[i]);
     1093                        }                       
    10741094                        rComplete(res);
    10751095                        return res;
    1076                 }//rCopyAndChange
     1096                }//rCopyAndAdd
     1097               
     1098                ring rCopyAndChangeWeight(ring const &r, intvec *ivw)
     1099                {
     1100                        ring res=rCopy0(currRing);
     1101                        rComplete(res);
     1102                        rSetWeightVec(res,(int64*)ivw);
     1103                        //rChangeCurrRing(rnew);
     1104                        return res;
     1105                }
    10771106               
    10781107                /**
    … …  
    10941123                        facet *fCmp;    //Ptr to alpha_j
    10951124                        fAct = fMin;
    1096                         fCmp = fMin->next;
    1097                        
    1098                         //cout << endl<< fMin->next << endl;
    1099                         //cout << gcTmp.facetPtr->next << endl;
    1100                         //gcTmp.facetPtr->printNormal();
    1101                         //fAct=gcTmp.facetPtr->next;
    1102                         //fAct->printNormal();                 
     1125                        fCmp = fMin->next;                             
    11031126                       
    11041127                        rChangeCurrRing(this->rootRing);        //because we compare the monomials in the rootring                     
    … …  
    11671190                                //fAct=fAct->next;
    11681191                        }//while(fAct.facetPtr->next!=NULL)
    1169                        
     1192                        delete alpha_i,alpha_j,sigma;
    11701193                        /*If testfacet was minimal then fMin should still point there */
    1171                         intvec *alpha_min = new intvec(this->numVars);
     1194                        //NOTE BUG: Comment in and -> out of memory error
     1195                        /*intvec *alpha_min = new intvec(this->numVars);
    11721196                        alpha_min=fMin->getFacetNormal();
     1197                        delete fCmp,fAct,fMin;
     1198                       
    11731199                        intvec *test = new intvec(this->numVars);
    1174                         test=testfacet.getFacetNormal();
    1175                         if (isParallel(alpha_min,test))
    1176                         //if (fMin==gcTmp.facetPtr)
     1200                        test=testfacet.getFacetNormal();                       
     1201                       
     1202                        if (isParallel(alpha_min,test))*/
     1203                        if (fMin==gcTmp.facetPtr)
    11771204                        {
    11781205                                rChangeCurrRing(actRing);
    … …  
    12801307        gcAct->getGB(rootIdeal);        //sets gcone::gcBasis
    12811308        idShow(gcAct->gcBasis);
    1282         gcAct->getConeNormals(gcAct->gcBasis);  //hopefully compute the normals
    1283         //gcAct->flip(gcAct->gcBasis,gcAct->facetPtr);
     1309        gcAct->getConeNormals(gcAct->gcBasis);  //hopefully compute the normals 
     1310        gcAct->flip(gcAct->gcBasis,gcAct->facetPtr);   
    12841311        /*Now it is time to compute the search facets, respectively start the reverse search.
    12851312        But since we are in the root all facets should be search facets. IS THIS TRUE?
    … …  
    12971324                fAct = fAct->next;             
    12981325        }*/
    1299         gcAct->reverseSearch(gcAct);
     1326        //gcAct->reverseSearch(gcAct);
    13001327        /*As of now extra.cc expects gfan to return type ideal. Probably this will change in near future.
    13011328        The return type will then be of type LIST_CMD
Note: See TracChangeset for help on using the changeset viewer.