source: git/kernel/gfan.cc @ 2ae96e

spielwiese
Last change on this file since 2ae96e was 0ad81f, checked in by Martin Monerjan, 15 years ago
Ring construction & change works for the lifting step in gcone::flip; git-svn-id: file:///usr/local/Singular/svn/trunk@11636 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 18.6 KB
Line 
1/*
2Compute the Groebner fan of an ideal
3$Author: monerjan $
4$Date: 2009-04-07 09:44:20 $
5$Header: /exports/cvsroot-2/cvsroot/kernel/gfan.cc,v 1.30 2009-04-07 09:44:20 monerjan Exp $
6$Id: gfan.cc,v 1.30 2009-04-07 09:44:20 monerjan Exp $
7*/
8
9#include "mod2.h"
10
11#ifdef HAVE_GFAN
12
13#include "kstd1.h"
14#include "intvec.h"
15#include "polys.h"
16#include "ideals.h"
17#include "kmatrix.h"
18#include "fast_maps.h"  //Mapping of ideals
19#include "maps.h"
20#include "iostream.h"   //deprecated
21
22//Hacks for different working places
23#define ITWM
24
25#ifdef UNI
26#include "/users/urmel/alggeom/monerjan/cddlib/include/setoper.h" //Support for cddlib. Dirty hack
27#include "/users/urmel/alggeom/monerjan/cddlib/include/cdd.h"
28#endif
29
30#ifdef HOME
31#include "/home/momo/studium/diplomarbeit/cddlib/include/setoper.h"
32#include "/home/momo/studium/diplomarbeit/cddlib/include/cdd.h"
33#endif
34
35#ifdef ITWM
36#include "/u/slg/monerjan/cddlib/include/setoper.h"
37#include "/u/slg/monerjan/cddlib/include/cdd.h"
38#endif
39
40#ifndef gfan_DEBUG
41#define gfan_DEBUG
42#endif
43
44//#include gcone.h
45
46/**
47*\brief Class facet
48*       Implements the facet structure as a linked list
49*
50*/
51class facet
52{
53        private:
54                /** \brief Inner normal of the facet, describing it uniquely up to isomorphism */
55                intvec *fNormal;
56               
57                /** \brief The Groebner basis on the other side of a shared facet
58                 *
59                 * In order not to have to compute the flipped GB twice we store the basis we already get
60                 * when identifying search facets. Thus in the next step of the reverse search we can
61                 * just copy the old cone and update the facet and the gcBasis.
62                 * facet::flibGB is set via facet::setFlipGB() and printed via facet::printFlipGB
63                 */
64                ideal flipGB;           //The Groebner Basis on the other side, computed via gcone::flip
65
66                       
67        public:
68                /** The default constructor. Do I need a constructor of type facet(intvec)? */
69                facet()                 
70                {
71                        // Pointer to next facet.  */
72                        /* Defaults to NULL. This way there is no need to check explicitly */
73                        this->next=NULL; 
74                }
75               
76                /** The default destructor */
77                ~facet(){;}
78               
79                /** Stores the facet normal \param intvec*/
80                void setFacetNormal(intvec *iv){
81                        this->fNormal = ivCopy(iv);
82                        //return;
83                }
84               
85                /** Hopefully returns the facet normal */
86                intvec *getFacetNormal()
87                {       
88                        //this->fNormal->show();        cout << endl;   
89                        return this->fNormal;
90                }
91               
92                /** Method to print the facet normal*/
93                void printNormal()
94                {
95                        fNormal->show();
96                }
97               
98                /** Store the flipped GB*/
99                void setFlipGB(ideal I)
100                {
101                        this->flipGB=I;
102                }
103               
104                /** Print the flipped GB*/
105                void printFlipGB()
106                {
107                        idShow(this->flipGB);
108                }
109               
110                bool isFlippable;       //flippable facet? Want to have cone->isflippable.facet[i]
111                bool isIncoming;        //Is the facet incoming or outgoing?
112                facet *next;            //Pointer to next facet
113};
114
115/**
116*\brief Implements the cone structure
117*
118* A cone is represented by a linked list of facet normals
119* @see facet
120*/
121/*class gcone
122finally this should become s.th. like gconelib.{h,cc} to provide an API
123*/
124class gcone
125{
126        private:
127                int numFacets;          //#of facets of the cone               
128
129        public:
130                /** \brief Default constructor.
131                *
132                * Initialises this->next=NULL and this->facetPtr=NULL
133                */
134                gcone()
135                {
136                        this->next=NULL;
137                        this->facetPtr=NULL;
138                }
139                ~gcone();               //destructor
140               
141                /** Pointer to the first facet */
142                facet *facetPtr;        //Will hold the adress of the first facet; set by gcone::getConeNormals
143                poly gcMarkedTerm;      //marked terms of the cone's Groebner basis
144                int numVars;            //#of variables in the ring
145               
146                /** Contains the Groebner basis of the cone. Is set by gcone::getGB(ideal I)*/
147                ideal gcBasis;          //GB of the cone, set by gcone::getGB();
148                gcone *next;            //Pointer to *previous* cone in search tree
149               
150                /** \brief Compute the normals of the cone
151                *
152                * This method computes a representation of the cone in terms of facet normals. It takes an ideal
153                * as its input. Redundancies are automatically removed using cddlib's dd_MatrixCanonicalize.
154                * Other methods for redundancy checkings might be implemented later. See Anders' diss p.44.
155                * Note that in order to use cddlib a 0-th column has to be added to the matrix since cddlib expects
156                * each row to represent an inequality of type const+x1+...+xn <= 0. While computing the normals we come across
157                * the set \f$ \partial\mathcal{G} \f$ which we might store for later use. C.f p71 of journal
158                * As a result of this procedure the pointer facetPtr points to the first facet of the cone.
159                */
160                void getConeNormals(const ideal &I)
161                {
162#ifdef gfan_DEBUG
163                        cout << "*** Computing Inequalities... ***" << endl;
164#endif         
165                        //All variables go here - except ineq matrix and *v, see below
166                        int lengthGB=IDELEMS(I);        // # of polys in the groebner basis
167                        int pCompCount;                 // # of terms in a poly
168                        poly aktpoly;
169                        int numvar = pVariables;        // # of variables in a polynomial (or ring?)
170                        int leadexp[numvar];            // dirty hack of exp.vects
171                        int aktexp[numvar];
172                        int cols,rows;                  // will contain the dimensions of the ineq matrix - deprecated by
173                        dd_rowrange ddrows;
174                        dd_colrange ddcols;
175                        dd_rowset ddredrows;            // # of redundant rows in ddineq
176                        dd_rowset ddlinset;             // the opposite
177                        dd_rowindex ddnewpos;           // all to make dd_Canonicalize happy
178                        dd_NumberType ddnumb=dd_Real;   //Number type
179                        dd_ErrorType dderr=dd_NoError;  //
180                        // End of var declaration
181#ifdef gfan_DEBUG
182                        cout << "The Groebner basis has " << lengthGB << " elements" << endl;
183                        cout << "The current ring has " << numvar << " variables" << endl;
184#endif
185                        cols = numvar;
186               
187                        //Compute the # inequalities i.e. rows of the matrix
188                        rows=0; //Initialization
189                        for (int ii=0;ii<IDELEMS(I);ii++)
190                        {
191                                aktpoly=(poly)I->m[ii];
192                                rows=rows+pLength(aktpoly)-1;
193                        }
194#ifdef gfan_DEBUG
195                        cout << "rows=" << rows << endl;
196                        cout << "Will create a " << rows << " x " << cols << " matrix to store inequalities" << endl;
197#endif
198                        dd_rowrange aktmatrixrow=0;     // needed to store the diffs of the expvects in the rows of ddineq
199                        dd_set_global_constants();
200                        ddrows=rows;
201                        ddcols=cols;
202                        dd_MatrixPtr ddineq;            //Matrix to store the inequalities
203                        ddineq=dd_CreateMatrix(ddrows,ddcols+1); //The first col has to be 0 since cddlib checks for additive consts there
204               
205                        // We loop through each g\in GB and compute the resulting inequalities
206                        for (int i=0; i<IDELEMS(I); i++)
207                        {
208                                aktpoly=(poly)I->m[i];          //get aktpoly as i-th component of I
209                                pCompCount=pLength(aktpoly);    //How many terms does aktpoly consist of?
210                                cout << "Poly No. " << i << " has " << pCompCount << " components" << endl;
211               
212                                int *v=(int *)omAlloc((numvar+1)*sizeof(int));
213                                pGetExpV(aktpoly,v);    //find the exp.vect in v[1],...,v[n], use pNext(p)
214                               
215                                //Store leadexp for aktpoly
216                                for (int kk=0;kk<numvar;kk++)
217                                {
218                                        leadexp[kk]=v[kk+1];
219                                        //Since we need to know the difference of leadexp with the other expvects we do nothing here
220                                        //but compute the diff below
221                                }
222               
223                               
224                                while (pNext(aktpoly)!=NULL) //move to next term until NULL
225                                {
226                                        aktpoly=pNext(aktpoly);
227                                        pSetm(aktpoly);         //doesn't seem to help anything
228                                        pGetExpV(aktpoly,v);
229                                        for (int kk=0;kk<numvar;kk++)
230                                        {
231                                                aktexp[kk]=v[kk+1];
232                                                //ineq[aktmatrixrow][kk]=leadexp[kk]-aktexp[kk];        //dito
233                                                dd_set_si(ddineq->matrix[(dd_rowrange)aktmatrixrow][kk+1],leadexp[kk]-aktexp[kk]); //because of the 1st col being const 0
234                                        }
235                                        aktmatrixrow=aktmatrixrow+1;
236                                }//while
237               
238                        } //for
239               
240                        //Maybe add another row to contain the constraints of the standard simplex?
241
242#ifdef gfan_DEBUG
243                        cout << "The inequality matrix is" << endl;
244                        dd_WriteMatrix(stdout, ddineq);
245#endif
246
247                        // The inequalities are now stored in ddineq
248                        // Next we check for superflous rows
249                        ddredrows = dd_RedundantRows(ddineq, &dderr);
250                        if (dderr!=dd_NoError)                  // did an error occur?
251                        {
252                                dd_WriteErrorMessages(stderr,dderr);    //if so tell us
253                        } else
254                        {
255                                cout << "Redundant rows: ";
256                                set_fwrite(stdout, ddredrows);          //otherwise print the redundant rows
257                        }//if dd_Error
258               
259                        //Remove reduntant rows here!
260                        dd_MatrixCanonicalize(&ddineq, &ddlinset, &ddredrows, &ddnewpos, &dderr);
261                        ddrows = ddineq->rowsize;       //Size of the matrix with redundancies removed
262                        ddcols = ddineq->colsize;
263#ifdef gfan_DEBUG
264                        cout << "Having removed redundancies, the normals now read:" << endl;
265                        dd_WriteMatrix(stdout,ddineq);
266                        cout << "Rows = " << ddrows << endl;
267                        cout << "Cols = " << ddcols << endl;
268#endif
269                       
270                        /*Write the normals into class facet*/
271#ifdef gfan_DEBUG
272                        cout << "Creating list of normals" << endl;
273#endif
274                        /*The pointer *fRoot should be the return value of this function*/
275                        facet *fRoot = new facet();     //instantiate new facet
276                        this->facetPtr = fRoot;         //set variable facetPtr of class gcone to first facet
277                        facet *fAct;                    //instantiate pointer to active facet
278                        fAct = fRoot;                   //Seems to do the trick. fRoot and fAct have to point to the same adress!
279                        std::cout << "fRoot = " << fRoot << ", fAct = " << fAct << endl;
280                        for (int kk = 0; kk<ddrows; kk++)
281                        {
282                                intvec *load = new intvec(numvar);      //intvec to store a single facet normal that will then be stored via setFacetNormal
283                                for (int jj = 1; jj <ddcols; jj++)
284                                {
285                                        double *foo;
286                                        foo = (double*)ddineq->matrix[kk][jj];  //get entry from actual position
287#ifdef gfan_DEBUG
288                                        std::cout << "fAct is " << *foo << " at " << fAct << std::endl;
289#endif 
290                                        (*load)[jj-1] = (int)*foo;              //store typecasted entry at pos jj-1 of load
291                                        //check for flipability here
292                                        if (jj<ddcols)                          //Is this facet NOT the last facet? Writing while instead of if is a really bad idea :)
293                                        {
294                                                //fAct->next = new facet();     //If so: instantiate new facet. Otherwise this->next=NULL due to the constructor                                               
295                                        }
296                                }//for (int jj = 1; jj <ddcols; jj++)
297                                /*Quick'n'dirty hack for flippability*/ 
298                                bool isFlippable;                       
299                                for (int jj = 0; jj<this->numVars; jj++)
300                                {                                       
301                                        intvec *ivCanonical = new intvec(this->numVars);
302                                        (*ivCanonical)[jj]=1;                                   
303                                        if (dotProduct(load,ivCanonical)>=0)
304                                        {
305                                                isFlippable=FALSE;                                             
306                                        }
307                                        else
308                                        {
309                                                isFlippable=TRUE;
310                                        }                                       
311                                        delete ivCanonical;
312                                }//for (int jj = 0; jj<this->numVars; jj++)
313                                if (isFlippable==FALSE)
314                                {
315                                        cout << "Ignoring facet";
316                                        load->show();
317                                        //fAct->next=NULL;
318                                }
319                                else
320                                {       /*Now load should be full and we can call setFacetNormal*/
321                                        fAct->setFacetNormal(load);
322                                        fAct->next = new facet();
323                                        //fAct->printNormal();
324                                        fAct=fAct->next;        //this should definitely not be called in the above while-loop :D
325                                }//if (isFlippable==FALSE)
326                                delete load;
327                        }//for (int kk = 0; kk<ddrows; kk++)
328                        /*
329                        Now we should have a linked list containing the facet normals of those facets that are
330                        -irredundant
331                        -flipable
332                        Adressing is done via *facetPtr
333                        */
334                       
335                        //Clean up but don't delete the return value! (Whatever it will turn out to be)
336                        dd_FreeMatrix(ddineq);
337                        set_free(ddredrows);
338                        free(ddnewpos);
339                        set_free(ddlinset);
340                        dd_free_global_constants();
341
342                }//method getConeNormals(ideal I)       
343               
344               
345                /** \brief Compute the Groebner Basis on the other side of a shared facet
346                *
347                * Implements algorithm 4.3.2 from Anders' thesis.
348                * As shown there it is not necessary to compute an interior point. The knowledge of the facet normal
349                * suffices. A term \f$ x^\gamma \f$ of \f$ g \f$ is in \f$  in_\omega(g) \f$ iff \f$ \gamma - leadexp(g)\f$
350                * is parallel to \f$ leadexp(g) \f$
351                * Parallelity is checked using basic linear algebra. See gcone::isParallel.
352                * Other possibilities includes  computing the rank of the matrix consisting of the vectors in question and
353                * computing an interior point of the facet and taking all terms having the same weight with respect
354                * to this interior point.
355                *\param ideal, facet
356                * Input: a marked,reduced Groebner basis and a facet
357                */
358                void flip(ideal gb, facet *f)           //Compute "the other side"
359                {                       
360                        intvec *fNormal = new intvec(this->numVars);    //facet normal, check for parallelity                   
361                        fNormal = f->getFacetNormal();  //read this->fNormal;
362#ifdef gfan_DEBUG
363                        cout << "===" << endl;
364                        cout << "running gcone::flip" << endl;
365                        cout << "fNormal=";
366                        fNormal->show();
367                        cout << endl;
368#endif                         
369                        /*1st step: Compute the initial ideal*/
370                        poly initialFormElement[IDELEMS(gb)];   //array of #polys in GB to store initial form
371                        ideal initialForm=idInit(IDELEMS(gb),this->gcBasis->rank);
372                        poly aktpoly;
373                        intvec *check = new intvec(this->numVars);      //array to store the difference of LE and v
374                       
375                        for (int ii=0;ii<IDELEMS(gb);ii++)
376                        {
377                                aktpoly = (poly)gb->m[ii];                                                             
378                                int *v=(int *)omAlloc((this->numVars+1)*sizeof(int));
379                                int *leadExpV=(int *)omAlloc((this->numVars+1)*sizeof(int));
380                                pGetExpV(aktpoly,leadExpV);     //find the leading exponent in leadExpV[1],...,leadExpV[n], use pNext(p)
381                                initialFormElement[ii]=pHead(aktpoly);
382                               
383                                while(pNext(aktpoly)!=NULL)     /*loop trough terms and check for parallelity*/
384                                {
385                                        aktpoly=pNext(aktpoly); //next term
386                                        pSetm(aktpoly);
387                                        pGetExpV(aktpoly,v);           
388                                        /* Convert (int)v into (intvec)check */                 
389                                        for (int jj=0;jj<this->numVars;jj++)
390                                        {
391                                                //cout << "v["<<jj+1<<"]="<<v[jj+1]<<endl;
392                                                //cout << "leadExpV["<<jj+1<<"]="<<leadExpV[jj+1]<<endl;
393                                                (*check)[jj]=v[jj+1]-leadExpV[jj+1];
394                                        }
395#ifdef gfan_DEBUG
396                                        cout << "check=";                       
397                                        check->show();
398                                        cout << endl;
399#endif
400                                        if (isParallel(check,fNormal)) //pass *check when
401                                        {
402                                                cout << "Parallel vector found, adding to initialFormElement" << endl;                         
403                                                initialFormElement[ii] = pAdd(initialFormElement[ii],(poly)pHead(aktpoly));
404                                        }                                               
405                                }//while
406                                cout << "Initial Form=";                               
407                                pWrite(initialFormElement[ii]);
408                                cout << "---" << endl;
409                                /*Now initialFormElement must be added to (ideal)initialForm */
410                                initialForm->m[ii]=initialFormElement[ii];
411                        }//for
412                        f->setFlipGB(initialForm);                     
413#ifdef gfan_DEBUG
414                        cout << "Initial ideal is: " << endl;
415                        idShow(initialForm);
416                        f->printFlipGB();
417                        cout << "===" << endl;
418#endif
419                        delete check;
420                       
421                        /*2nd step: lift initial ideal to a GB of the neighbouring cone using minus alpha as weight*/
422                        /*Substep 2.1
423                        compute $G_{-\alpha}(in_v(I))
424                        see journal p. 66
425                        */
426                        ring srcRing=currRing;
427                        ring dstRing=rCopy0(srcRing);
428                        dstRing->order[0]=ringorder_a;
429                        //tmpring->order[1]=ringorder_dp;
430                        //tmpring->order[2]=ringorder_C;
431                        dstRing->wvhdl[0] =( int *)omAlloc((fNormal->length())*sizeof(int)); //found in Singular/ipshell.cc
432                       
433                        for (int ii=0;ii<this->numVars;ii++)
434                        {                               
435                                dstRing->wvhdl[0][ii]=-(*fNormal)[ii];  //What exactly am I doing here?
436                                //cout << tmpring->wvhdl[0][ii] << endl;
437                        }
438                        rComplete(dstRing);
439                        rChangeCurrRing(dstRing);
440                        map theMap=(map)idMaxIdeal(1);
441                        rWrite(currRing); cout << endl;
442                        ideal ina;
443                        ina=fast_map(initialForm,srcRing,(ideal)theMap,dstRing);                       
444                        cout << "ina=";
445                        idShow(ina); cout << endl;
446                        ideal H;
447                        H=kStd(ina,NULL,isHomog,NULL);  //we know it is homogeneous
448                        idSkipZeroes(H);
449                        cout << "H="; idShow(H); cout << endl;
450                        /*Substep 2.2
451                        do the lifting
452                        */
453                        rChangeCurrRing(srcRing);
454                        ideal srcRing_H;
455                        ideal srcRing_HH;
456                        //map theMap = (map)idMaxIdeal(1);
457                        srcRing_H=fast_map(H,dstRing,(ideal)theMap,srcRing);
458                        srcRing_HH=srcRing_H-stdred(srcRing_H,this->gcBasis);
459                        /*Substep 2.3
460                        turn the minimal basis into a reduced one
461                        */
462                       
463                }//void flip(ideal gb, facet *f)
464                               
465                /** \brief Compute a Groebner Basis
466                *
467                * Computes the Groebner basis and stores the result in gcone::gcBasis
468                *\param ideal
469                *\return void
470                */
471                void getGB(ideal const &inputIdeal)
472                {
473                        ideal gb;
474                        gb=kStd(inputIdeal,NULL,testHomog,NULL);
475                        idSkipZeroes(gb);
476                        this->gcBasis=gb;       //write the GB into gcBasis
477                }//void getGB
478               
479                ideal GenGrbWlk(ideal, ideal);  //Implementation of the Generic Groebner Walk. Needed for a) Computing the sink and b) finding search facets
480
481
482                /** \brief Compute the dot product of two intvecs
483                *
484                */
485                int dotProduct(intvec const &iva, intvec const &ivb)                           
486                {
487                        //intvec iva=a;
488                        //intvec ivb=b;
489                        int res=0;
490                        for (int i=0;i<this->numVars;i++)
491                        {
492                                res = res+(iva[i]*ivb[i]);
493                        }
494                        return res;
495                }//int dotProduct
496
497                /** \brief Check whether two intvecs are parallel
498                *
499                * \f$ \alpha\parallel\beta\Leftrightarrow\langle\alpha,\beta\rangle^2=\langle\alpha,\alpha\rangle\langle\beta,\beta\rangle \f$
500                */
501                bool isParallel(intvec const &a, intvec const &b)
502                {                       
503                        int lhs,rhs;
504                        lhs=dotProduct(a,b)*dotProduct(a,b);
505                        rhs=dotProduct(a,a)*dotProduct(b,b);
506                        cout << "LHS="<<lhs<<", RHS="<<rhs<<endl;
507                        if (lhs==rhs)
508                        {
509                                return TRUE;
510                        }
511                        else
512                        {
513                                return FALSE;
514                        }
515                }//bool isParallel
516};//class gcone
517
518ideal gfan(ideal inputIdeal)
519{
520        int numvar = pVariables; 
521       
522        #ifdef gfan_DEBUG
523        cout << "Now in subroutine gfan" << endl;
524        #endif
525        ring inputRing=currRing;        // The ring the user entered
526        ring rootRing;                  // The ring associated to the target ordering
527        ideal res;
528        //matrix ineq;                  //Matrix containing the boundary inequalities
529        facet *fRoot;
530       
531        /*
532        1. Select target order, say dp.
533        2. Compute GB of inputIdeal wrt target order -> newRing, setCurrRing etc...
534        3. getConeNormals
535        */
536       
537        /* Construct a new ring which will serve as our root
538        Does not yet work as expected. Will work fine with order dp,Dp but otherwise hangs in getGB
539        */
540        rootRing=rCopy0(currRing);
541        rootRing->order[0]=ringorder_lp;
542        rComplete(rootRing);
543        rChangeCurrRing(rootRing);
544       
545        /* Fetch the inputIdeal into our rootRing */
546        map theMap=(map)idMaxIdeal(1);  //evil hack!
547        //idShow(idMaxIdeal(1));
548        /*for (int ii=0;ii<pVariables;ii++)
549        {
550                theMap->m[ii]=inputIdeal->m[ii];
551        }*/
552        theMap->preimage=NULL;
553        ideal rootIdeal;
554        rootIdeal=fast_map(inputIdeal,inputRing,(ideal)theMap, currRing);
555#ifdef gfan_DEBUG
556        cout << "Root ideal is " << endl;
557        idShow(rootIdeal);
558        cout << "The root ring is " << endl;
559        rWrite(rootRing);
560        cout << endl;
561#endif 
562       
563        gcone *gcRoot = new gcone();    //Instantiate the sink
564        gcone *gcAct;
565        gcAct = gcRoot;
566        gcAct->numVars=pVariables;
567        gcAct->getGB(rootIdeal);        //sets gcone::gcBasis
568        idShow(gcAct->gcBasis);
569        gcAct->getConeNormals(gcAct->gcBasis);  //hopefully compute the normals
570        gcAct->flip(gcAct->gcBasis,gcAct->facetPtr);
571        /*Now it is time to compute the search facets, respectively start the reverse search.
572        But since we are in the root all facets should be search facets. IS THIS TRUE?
573        MIND: AS OF NOW, THE LIST OF FACETS IS NOT PURGED OF NON-FLIPPAPLE FACETS
574        */
575       
576        /*As of now extra.cc expects gfan to return type ideal. Probably this will change in near future.
577        The return type will then be of type LIST_CMD
578        Assume gfan has finished, thus we have enumerated all the cones
579        Create an array of size of #cones. Let each entry in the array contain a pointer to the respective
580        Groebner Basis and merge this somehow with LIST_CMD
581        => Count the cones!
582        */
583        rChangeCurrRing(inputRing);
584        res=gcAct->gcBasis;
585        //cout << fRoot << endl;
586        return res;
587        //return GBlist;
588}
589/*
590Since gfan.cc is #included from extra.cc there must not be a int main(){} here
591*/
592#endif
Note: See TracBrowser for help on using the repository browser.