#include "vinci.h"

Include dependency graph for MP/Source/vinci_lass.c:

Macros
#define	MAXIMUM 1.0e150 /* define the maximum used for testing infinity */

#define	EPSILON_LASS EPSILON /* Numbers smaller than this are treated as zero in rhs*/

#define	EPS1 EPSILON /* Numbers smaller than this are treated as zero in coefficient */

#define	EPS_NORM EPSILON /* EPSILON used for constraint normalization*/

#define	LaShiftLevel 0 /* Shifting is possible if d>=LaShiftLevel */

#define	LaShift

#define	verboseFirstLevel /* output the intermediate volume of the first level */

Functions
void *	my_malloc (long int size)

void	my_free (void *pointer, long int size)

void	create_key (T_Key *key, int key_choice)

void	free_key (T_Key key, int key_choice)

void	add_hypervar (T_LassInt hyperplane, T_LassInt variable, T_Key *key)

void	delete_hypervar (T_LassInt hyperplane, T_LassInt variable, T_Key *key)

int *	create_int_vector (n)

void	free_int_vector (int *v, int n)

static T_LassInt	add_reduced_index (T_LassInt red, T_LassInt indices, T_LassInt ref_indices)

static void	del_original_indices (T_LassInt indices, T_LassInt org_indices)

static void	del_original (T_LassInt base, T_LassInt *indices)

static void	rm_original_inElAll_index (T_LassInt baserow)

static void	rm_constraint (rational A, int LastPlane_, int d, int rm_index)

static rational *	compact ()

static int	notInPivot (int *pivot, int col, int i)

static void	shift_P (rational *A, int LastPlane_, int d)

static int	norm_and_clean_constraints (rational A, int LastPlane_, int d, T_LassInt *Del_index, int Index_needed)

static rational	lass (rational *A, int LastPlane_, int d)

void	VOLUME_LASSERRE_FILE (double planes2[7][4], rational *volume)

void	volume_lasserre_file (double planes2[7][4], rational *volume)

void	VOLUME_LASSERRE_FILE_ (double planes2[7][4], rational *volume)

void	volume_lasserre_file_ (double planes2[7][4], rational *volume)

void	tree_out (T_Tree *ppr, int pi_balance, T_Key key, rational volume, T_Key keyfound, int key_choice)

T_VertexSet	duplicate_set (T_VertexSet s)

static int	compare_key (T_Key key1, T_Key key2, int key_choice)

Variables
int	G_d

int	G_m

real **	G_Hyperplanes = NULL

T_VertexSet *	G_Incidence = NULL

T_VertexSet	G_Vertices

int	G_Precomp = 0

int	G_Storage = -1

int	G_RandomSeed = 4

boolean	G_OutOfMem = TRUE

void *	G_MemRes = NULL

rational	G_Minus1 = -1

rational *	A

rational *	pivotrow

T_LassInt *	All_index

T_LassInt *	Pivot

int **	p2c

rational *	planescopy

static T_Key	key

static T_Key *	keyfound

static T_Tree *	tree_volumes

static T_VertexSet *	face

Macro Definition Documentation

#define EPS1 EPSILON /* Numbers smaller than this are treated as zero in coefficient */

Definition at line 27 of file MP/Source/vinci_lass.c.

Referenced by norm_and_clean_constraints().

#define EPS_NORM EPSILON /* EPSILON used for constraint normalization*/

Definition at line 28 of file MP/Source/vinci_lass.c.

Referenced by norm_and_clean_constraints().

#define EPSILON_LASS EPSILON /* Numbers smaller than this are treated as zero in rhs*/

Definition at line 26 of file MP/Source/vinci_lass.c.

Referenced by lass().

#define LaShift

Value:

1 /* shift polytope to make at least d components of

the rhs zero if there are less than d-LaShift zeros */

Definition at line 30 of file MP/Source/vinci_lass.c.

Referenced by lass().

#define LaShiftLevel 0 /* Shifting is possible if d>=LaShiftLevel */

Definition at line 29 of file MP/Source/vinci_lass.c.

Referenced by lass().

#define MAXIMUM 1.0e150 /* define the maximum used for testing infinity */

Definition at line 25 of file MP/Source/vinci_lass.c.

Referenced by lass().

#define verboseFirstLevel /* output the intermediate volume of the first level */

Definition at line 34 of file MP/Source/vinci_lass.c.

Function Documentation

void add_hypervar	(	T_LassInt	hyperplane,
		T_LassInt	variable,
		T_Key *	key
	)

Definition at line 168 of file MP/Source/vinci_lass.c.

References G_d, G_m, i, and j.

Referenced by lass().

 {  int i, j;
 
    if (hyperplane != G_m+1)
    {  
       for (i = 0; (*key).hypervar.hyperplanes [i] < hyperplane; i++);
       if ((*key).hypervar.hyperplanes [i] != hyperplane)
       {  /* insert index */
          for (j = G_d; j > i; j--)
             (*key).hypervar.hyperplanes [j] = (*key).hypervar.hyperplanes [j-1];
          (*key).hypervar.hyperplanes [i] = hyperplane;
       }
    }
 
    if (variable != G_d+1)
    {  
       for (i = 0; (*key).hypervar.variables [i] < variable; i++);
       if ((*key).hypervar.variables [i] != variable)
       {  /* insert index */
          for (j = G_d; j > i; j--)
             (*key).hypervar.variables [j] = (*key).hypervar.variables [j-1];
          (*key).hypervar.variables [i] = variable;
       }
    }
 }

Here is the caller graph for this function:

static T_LassInt add_reduced_index	(	T_LassInt	red,
		T_LassInt *	indices,
		T_LassInt *	ref_indices
	)

static

Definition at line 258 of file MP/Source/vinci_lass.c.

References G_m, and i.

Referenced by lass(), and norm_and_clean_constraints().

 { register int i;
   T_LassInt xch, base;
     
   for (i=0; red>=ref_indices[i]; i++) red++;  /* reduced index -> original index */
   base=red;  
   while (ref_indices[i]<=G_m) {
       xch=ref_indices[i];
       ref_indices[i]=red;
       red=xch;
       i++;
   };
   ref_indices[i]=red;
   ref_indices[i+1]=G_m+2;
   if (indices==NULL) return base;
   red=base;
   for (i=0; base>indices[i]; i++);
   while (indices[i]<=G_m) {
       xch=indices[i];
       indices[i]=red;
       red=xch;
       i++;
   };
   indices[i]=red;
   indices[i+1]=G_m+2;
   return base;
 }

Here is the caller graph for this function:

static rational* compact ( )

static

Definition at line 345 of file MP/Source/vinci_lass.c.

References G_d, G_Hyperplanes, G_m, i, j, and my_malloc().

Referenced by VOLUME_LASSERRE_FILE().

 {   register int i, j;
     register rational *po,*pc;
 
     if (!(pc = (rational *) my_malloc (G_m*(G_d+1)*sizeof(rational)))) {
         fprintf (stderr, "\n***** ERROR: Out of memory in 'compact.*pc'");
         exit(0); 
     }
     po=pc;
     for (i=0; i<G_m; i++) {
         for (j=0; j<=G_d; j++,pc++) *pc= G_Hyperplanes [i][j];
     };
     return po;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

static int compare_key	(	T_Key	key1,
		T_Key	key2,
		int	key_choice
	)

static

Definition at line 991 of file MP/Source/vinci_lass.c.

References T_Key::d, G_m, T_Key::hyperplanes, T_Key::hypervar, i, KEY_PLANES, KEY_PLANES_VAR, KEY_VERTICES, T_VertexSet::lastel, T_VertexSet::loe, p1, T_Key::set, and T_Key::vertices.

Referenced by tree_out().

 {  int       i;
    T_LassInt *p1, *p2;
    
    switch (key_choice)
    {
    case KEY_PLANES:
       for (i=0; ; i++)
          if      (key1.hyperplanes [i] < key2.hyperplanes [i]) return -1;
          else if (key1.hyperplanes [i] > key2.hyperplanes [i]) return 1;
          else if (key1.hyperplanes [i] > G_m)                  return 0;
       break;
    case KEY_VERTICES:
       if      (key1.vertices.d < key2.vertices.d) return -1;
       else if (key1.vertices.d > key2.vertices.d) return 1;
       else /* both volumes are for the same dimension */
       if      (key1.vertices.set.lastel < key2.vertices.set.lastel) return -1;
       else if (key1.vertices.set.lastel > key2.vertices.set.lastel) return 1;
       else
       {  /* both sets have the same cardinality */
          for (i=0; i <= key1.vertices.set.lastel; i++)
             if      (key1.vertices.set.loe [i] -> no < key2.vertices.set.loe [i] -> no)
                return -1;
             else if (key1.vertices.set.loe [i] -> no > key2.vertices.set.loe [i] -> no)
                return 1;
       }
       return 0;
       break;
    case KEY_PLANES_VAR:
 /*
       for (i=0; ; i++)
          if      (key1.hypervar.hyperplanes [i] < key2.hypervar.hyperplanes [i])
                                                                         return -1;
          else if (key1.hypervar.hyperplanes [i] > key2.hypervar.hyperplanes [i])
                                                                         return 1;
          else if (key1.hypervar.hyperplanes [i] > G_m)                  return 0;
 */
       for (p1 = key1.hypervar.hyperplanes, p2 = key2.hypervar.hyperplanes;;
            p1++, p2++)
          if      ((*p1) < (*p2)) return -1;
          else if ((*p1) > (*p2)) return  1;
          else if ((*p1) > G_m)   return  0;
       break;
    }
 }

Here is the caller graph for this function:

int* create_int_vector ( n )

Definition at line 236 of file MP/Source/vinci_lass.c.

References my_malloc(), n, and v.

 {  int *v;
    
    v = (int *) my_malloc (n * sizeof (int));
    
    return v;
 }

Here is the call graph for this function:

void create_key	(	T_Key *	key,
		int	key_choice
	)

Definition at line 131 of file MP/Source/vinci_lass.c.

References create_int_vector(), G_d, G_Storage, T_Key::hyperplanes, KEY_PLANES, KEY_PLANES_VAR, my_malloc(), and T_Key::variables.

Referenced by tree_out(), and VOLUME_LASSERRE_FILE().

 {
    switch (key_choice)
    {
    case KEY_PLANES:
       key -> hyperplanes = create_int_vector (G_d + 1);
       break;
    case KEY_PLANES_VAR:
       key -> hypervar.hyperplanes =
                                  (T_LassInt *) my_malloc ((G_Storage + 2) * sizeof (T_LassInt));
       key -> hypervar.variables =
                                  (T_LassInt *) my_malloc ((G_Storage + 2) * sizeof (T_LassInt));
       break;
    }
 }

Here is the call graph for this function:

Here is the caller graph for this function:

static void del_original	(	T_LassInt	base,
		T_LassInt *	indices
	)

static

Definition at line 308 of file MP/Source/vinci_lass.c.

References G_m, and i.

Referenced by lass(), and rm_original_inElAll_index().

 { int i;
 
   for (i=0; ((base!=indices[i])&&(indices[i]<=G_m)); i++);  /* search original index */
   if (base!=indices[i]) {
       fprintf(stderr, "ERROR: Deletion index not found!\n");
       exit(0);
   };
   for (;indices[i]<=G_m;i++) indices[i]=indices[i+1];
 }  

Here is the caller graph for this function:

static void del_original_indices	(	T_LassInt *	indices,
		T_LassInt *	org_indices
	)

static

Definition at line 292 of file MP/Source/vinci_lass.c.

References G_m, and i.

Referenced by lass().

 {   register int i, cnt;
 
     i=cnt=0;
     while (org_indices[i]<=G_m) {
         while ((org_indices[cnt+i]==indices[cnt])&&(indices[cnt]<=G_m)) cnt++;
         org_indices[i]=org_indices[i+cnt];
         i++;
     }
 }

Here is the caller graph for this function:

void delete_hypervar	(	T_LassInt	hyperplane,
		T_LassInt	variable,
		T_Key *	key
	)

Definition at line 204 of file MP/Source/vinci_lass.c.

References G_d, G_m, i, and j.

Referenced by lass().

 {  int i, j;
 
    if (hyperplane != G_m+1)
    {  
       for (i = 0; i <= G_d && (*key).hypervar.hyperplanes [i] != hyperplane; i++);
       if (i != G_d + 1)
       {  /* index found, delete it */
          for (j = i; (*key).hypervar.hyperplanes [j] != G_m + 1; j++)
             (*key).hypervar.hyperplanes [j] = (*key).hypervar.hyperplanes [j+1];
       }
    }
 
    if (variable != G_d+1)
    {  
       for (i = 0; i <= G_d && (*key).hypervar.variables [i] != variable; i++);
       if (i != G_d + 1)
       {  /* index found, delete it */
          for (j = i; (*key).hypervar.variables [j] != G_d + 1; j++)
             (*key).hypervar.variables [j] = (*key).hypervar.variables [j+1];
       }
    }
 }

Here is the caller graph for this function:

T_VertexSet duplicate_set ( T_VertexSet s )

Definition at line 976 of file MP/Source/vinci_lass.c.

References T_VertexSet::lastel, T_VertexSet::loe, T_VertexSet::maxel, and my_malloc().

Referenced by tree_out().

 {  T_VertexSet newset;
 
    newset.maxel = s.lastel;
    newset.lastel = s.lastel;
    newset.loe = (T_Vertex **) my_malloc ((s.lastel + 1) * sizeof (T_Vertex *));
    memcpy (newset.loe, s.loe, (s.lastel + 1) * sizeof (T_Vertex *));
    return newset;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void free_int_vector	(	int *	v,
		int	n
	)

Definition at line 248 of file MP/Source/vinci_lass.c.

References my_free().

Referenced by free_key().

 {  
    my_free (v, n * sizeof (int));
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void free_key	(	T_Key	key,
		int	key_choice
	)

Definition at line 151 of file MP/Source/vinci_lass.c.

References free_int_vector(), G_d, G_Storage, T_Key::hyperplanes, T_Key::hypervar, KEY_PLANES, KEY_PLANES_VAR, my_free(), and T_Key::variables.

Referenced by VOLUME_LASSERRE_FILE().

 {
    switch (key_choice)
    {
    case KEY_PLANES:
       free_int_vector (key.hyperplanes, G_d + 1);
       break;
    case KEY_PLANES_VAR:
       my_free (key.hypervar.hyperplanes, (G_Storage + 2) * sizeof (T_LassInt));
       my_free (key.hypervar.variables, (G_Storage + 2) * sizeof (T_LassInt));
       break;
    }
 }

Here is the call graph for this function:

Here is the caller graph for this function:

static rational lass	(	rational *	A,
		int	LastPlane_,
		int	d
	)

static

Definition at line 572 of file MP/Source/vinci_lass.c.

References A, add_hypervar(), add_reduced_index(), All_index, d, del_original(), del_original_indices(), delete_hypervar(), EPSILON_LASS, FALSE, G_d, G_m, G_MemRes, G_OutOfMem, G_Storage, i, j, KEY_PLANES_VAR, LaShift, LaShiftLevel, MAXIMUM, MIN_PIVOT_LASS, my_free(), my_malloc(), norm_and_clean_constraints(), p2c, Pivot, pivotrow, rm_original_inElAll_index(), shift_P(), tree_out(), TRUE, and volume().

Referenced by VOLUME_LASSERRE_FILE().

 {   rational * redA;            /* A reduced by one dimension and constraint */
     int i, j;
     T_LassInt baserow, basecol, col;
     int dimdiff, row;         /* dimension difference and boolean */
     int i_balance = FALSE;
     rational ma, mi, *volume, *realp1, *realp2;
     int Index_needed;         /* Boolean, if index operations are needed */
     T_LassInt * Del_index; /* contains the indices of the deleted planes */
 
     /* test if volume is already known and return it if so */
 
     dimdiff = G_d-d;
     if ((G_Storage > (dimdiff-2)) && (dimdiff >= 2)) {
         tree_out (&tree_volumes, &i_balance, key, &volume, &keyfound, KEY_PLANES_VAR);
         if ((*volume)>=0)  {  /* this volume has already been computed */       
           printf("\nNeed Scale routine!\n");
         }
         if (!G_OutOfMem) {
             (*volume)=0;      /* initialize */
             dimdiff=TRUE;
         }
         else dimdiff=FALSE;
     }
     else dimdiff=FALSE;
 
     /* if d==1 compute the volume and give it back */
 
     if (d == 1) {
         ma=-MAXIMUM;
         mi= MAXIMUM;
         for (i=0; i<=LastPlane_; i++,A+=2) { 
             if (*A>EPSILON_LASS) { if ((*(A+1)/ *A)<mi) mi=(*(A+1)/ *A); }
             else if (*A<-EPSILON_LASS) { if ((*(A+1)/ *A)>ma) ma=*(A+1)/ *A; } 
             else if ((*(A+1))<-(100000*EPSILON_LASS)) return 0; 
         }
         if ((ma<-.5*MAXIMUM)||(mi>.5*MAXIMUM)) {
             printf("\nVolume is unbounded!\n");
             exit(0);
         }
         if ((mi-ma)>EPSILON_LASS) {
             if (dimdiff) (*volume)=mi-ma;
             return mi-ma;
         }
         return 0;
     }
 
     /* if d>1 apply the recursive scheme by fixing constraints. */
 
     Index_needed = (G_Storage>(G_d-d-1));
     if (Index_needed){
         if (!(Del_index = (T_LassInt *) my_malloc ((LastPlane_ + 2) * sizeof (T_LassInt)))){
             fprintf (stderr, "\n***** ERROR/WARNING: Out of memory in 'lass'\n");
             if (G_OutOfMem) exit(0);
             G_OutOfMem = TRUE;
             free(G_MemRes);
             Del_index = (T_LassInt *) malloc ((LastPlane_ + 2) * sizeof (T_LassInt));
         };
         Del_index[0]=G_m+2;   /* initialize: mark end */
     }
     ma=0;                                         /* used to sum up the summands */
     if (norm_and_clean_constraints(A, &LastPlane_, d, Del_index, Index_needed)!=0)
         goto label2;
 
     /* if appropriate shift polytope */
 
     if (d>=LaShiftLevel) {
         realp1=A+d;
         realp2=realp1+LastPlane_*(d+1);
         j=0;
         while (realp1<=realp2) {
             if (fabs(*realp1)<EPSILON_LASS) j++;
             realp1+=d+1;
         }
         if (d-j>=LaShift) shift_P(A, LastPlane_, d);
     }
 
 
     redA = (rational *) my_malloc (LastPlane_* d*sizeof(rational));
     if (redA == NULL) {
         fprintf (stderr, "\n***** ERROR/WARNING: Out of memory in 'lass.*redA'\n");
         if (G_OutOfMem) exit(0);
         G_OutOfMem = TRUE;
         free(G_MemRes);
         redA = (rational *) malloc (LastPlane_* d*sizeof(rational));
     }
 #ifdef ReverseLass
     for (row=LastPlane_; row>=0; row--) {
 #else
     for (row=0; row<=LastPlane_; row++) {
 #endif
         if (fabs(*(A+row*(d+1)+d))<EPSILON_LASS) 
             continue;                        /* skip this constraint if b_row == 0 */
         if (Index_needed)
         {  baserow=add_reduced_index(row, NULL, All_index);
            p2c[G_d-d][1] = baserow;
            add_hypervar (baserow, G_d+1, &key);
         }       
         memcpy(&pivotrow[0], A+row*(d+1), sizeof(rational)*(d+1));
         col=0;                               /* search for pivot column */
         for (i=0; i<d; i++) {        
 #if PIVOTING_LASS == 0
             if (fabs(pivotrow[i])>=MIN_PIVOT_LASS) {col=i; break;};
 #endif
             if (fabs(pivotrow[i])>fabs(pivotrow[col])) col=i;
         };
         if (G_Storage>(G_d-d-1))
         {  basecol=add_reduced_index(col, NULL, Pivot);
            p2c[G_d-d][0] = basecol;
            add_hypervar (G_m+1, basecol, &key);
         }
 
         /* copy A onto redA and at the same time perform pivoting */
          
         mi=1.0/pivotrow[col];
         for (i=0; i<=d; i++) pivotrow[i]*=mi;
         realp1=A;
         realp2=redA;
         for (i=0; i<=LastPlane_; i++) {
             if (i==row) {
                 realp1+=d+1;
                 continue;
             };
             mi=*(A+(i*(d+1))+col);
             for (j=0; j<=d; j++) {
                 if (j==col) {
                     realp1++;
                     continue;
                 };
                 *realp2=(*realp1)-pivotrow[j]*mi;
                 realp1++;
                 realp2++;
             };
         };
         ma+= *(A+row*(d+1)+d)/(d*fabs(*(A+row*(d+1)+col)))
              *lass(redA, LastPlane_-1, d-1);
         if (Index_needed)
         {  rm_original_inElAll_index(baserow);
            delete_hypervar (baserow, G_d+1, &key);
         }
         if (G_Storage>(G_d-d-1))
         {  del_original(basecol, Pivot);
            delete_hypervar (G_m+1, basecol, &key);
         }
         /*$#ifdef verboseFirstLevel
             if (d==G_d) 
                 printf("\nVolume accumulated to iteration %i is %20.12f",row,ma );
         #endif$*/
     };
     my_free (redA, LastPlane_* d * sizeof (rational));
     label2: 
     if (Index_needed) {
         del_original_indices(Del_index, All_index);
         my_free (Del_index, (LastPlane_ + 2) * sizeof (T_LassInt));
     };
     if (dimdiff)(*volume)=ma;
     return ma;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void my_free	(	void *	pointer,
		long int	size
	)

Definition at line 121 of file MP/Source/vinci_lass.c.

Referenced by free_int_vector(), free_key(), and lass().

 {
    free (pointer);
 }

Here is the caller graph for this function:

void* my_malloc ( long int size )

Definition at line 90 of file MP/Source/vinci_lass.c.

References G_MemRes, G_OutOfMem, my_malloc(), and TRUE.

Referenced by compact(), create_int_vector(), create_key(), duplicate_set(), lass(), my_malloc(), tree_out(), and VOLUME_LASSERRE_FILE().

 {  void *pointer;
 
    pointer = malloc (size);
 
    if (pointer == NULL)
       if (!G_OutOfMem)
       {
          fprintf (stderr, "\n***** WARNING: Out of memory; freeing memory reserve\n");
          free (G_MemRes);
          G_OutOfMem = TRUE;
          pointer = my_malloc (size);
       }
       else
       {
          fprintf (stderr, "\n***** ERROR: Out of memory\n");
          exit (0);
       }
    else
    {
 
    }
 
    return pointer;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

static int norm_and_clean_constraints	(	rational *	A,
		int *	LastPlane_,
		int	d,
		T_LassInt *	Del_index,
		int	Index_needed
	)

static

Definition at line 483 of file MP/Source/vinci_lass.c.

References A, add_reduced_index(), All_index, d, EPS1, EPS_NORM, i, j, p1, rm_constraint(), and sqrt().

Referenced by lass().

 {   register int i, j, row = 0;
     register rational r0, *p1, *p2;
 
     /* find nonzero[][] and maximal elements and normalize */
   
     p1=A;                                  /* begin of first constraint */
     while (row<=(*LastPlane_)) {           /* remove zeros and normalize */
         r0=0.0;                            /* norm of vector */
         for (j=0; j<d; j++,p1++) 
             r0+=(*p1)*(*p1);               /* compute euclidean norm */
         r0=sqrt(r0);
         if (r0<EPS_NORM) {
             if ((*p1)<-100000*EPS1){      /* if negative rhs */
                 return 1;                  /* infeasible constraint */
             }
             rm_constraint(A, LastPlane_, d,row);
             if (Index_needed) add_reduced_index(row, Del_index, All_index);
             p1-=d;
         }
         else {
             r0=1.0/r0;
             p1-=d;
             for (j=0; j<=d; j++,p1++)
                 (*p1)*=r0;
             row++; 
         }
     }
 
     /* detect identical or reverse constraints */
     
     for (row=0; row<(*LastPlane_); row++) {
         i=row+1;
         while (i<=*LastPlane_) {        /* test all subsequent rows i if equal to row */
             r0=0.0;
             p1=A+row*(d+1);
             p2=A+i*(d+1);
             for (j=0;j<d;j++,p1++,p2++)
                 r0+=(*p1)*(*p2);        /* cosinus of arc among those two vectors */
             if (r0>0) {
                 /* NEW VERSION of removing constraints */ 
                 if (fabs(r0-1.0)<EPS_NORM) {
                     if ((*p1)>(*p2)){
                         if (Index_needed) add_reduced_index(row, Del_index, All_index);
                         rm_constraint(A, LastPlane_, d,row);
                         i=row+1;
                     }
                     else {
                         if (Index_needed) add_reduced_index(i, Del_index, All_index);
                         if (i<(*LastPlane_)) 
                             rm_constraint(A, LastPlane_, d,i);
                         else (*LastPlane_)--;
                     }
                 }
                 else i++;
 
                 /* OLD VERSION :
                 if ((fabs(r0-1.0)<EPS_NORM) && (fabs((*p1)-(*p2))<EPS1)){
                     if (Index_needed) add_reduced_index(i, Del_index, All_index);
                     if (i<(*LastPlane_)) 
                         rm_constraint(A, LastPlane_, d,i);
                     else (*LastPlane_)--;
                 }
                 else i++;
                 */
             }
             else {
                 if (fabs(r0+1.0)<EPS_NORM){
                     if ((*p1)>0){
                         if ((*p2)<(EPS1-(*p1))) return 1; 
                      }
                      else {
                          if ((*p1)<(EPS1-(*p2))) return 1; 
                      }
                 }
                 i++;
             }
         }
     }
     return 0;  /* elimination succesful */
 }

Here is the call graph for this function:

Here is the caller graph for this function:

static int notInPivot	(	int *	pivot,
		int	col,
		int	i
	)

static

Definition at line 366 of file MP/Source/vinci_lass.c.

References FALSE, and TRUE.

Referenced by shift_P().

 { register int h;
   for (h=0;h<col;h++)
    if (pivot[h]==i) return FALSE;
   return TRUE;
 }

Here is the caller graph for this function:

static void rm_constraint	(	rational *	A,
		int *	LastPlane_,
		int	d,
		int	rm_index
	)

static

Definition at line 328 of file MP/Source/vinci_lass.c.

References i, and p1.

Referenced by norm_and_clean_constraints().

 {   register rational *p1, *p2; 
     register int i;
 
     p1=A+rm_index*(d+1);
     p2=A+(rm_index+1)*(d+1);
     for (i=0; i<(((*LastPlane_)-rm_index)*(d+1)); i++) {
         *p1=*p2;
         p1++;
         p2++;
     };
     (*LastPlane_)--;
 }

Here is the caller graph for this function:

static void rm_original_inElAll_index ( T_LassInt baserow )

static

Definition at line 323 of file MP/Source/vinci_lass.c.

References All_index, and del_original().

Referenced by lass().

325 { del_original(baserow, All_index); }

All_index

T_LassInt * All_index

Definition: vinci_lass.c:69

del_original

static void del_original(T_LassInt base, T_LassInt *indices)

Definition: MP/Source/vinci_lass.c:308

Here is the call graph for this function:

Here is the caller graph for this function:

static void shift_P	(	rational *	A,
		int	LastPlane_,
		int	d
	)

static

Definition at line 374 of file MP/Source/vinci_lass.c.

References A, create_int_vector(), d, G_d, i, j, MIN_PIVOT_LASS, notInPivot(), p1, and planescopy.

Referenced by lass().

 {   register rational  *p1, *p2, *p3, d1, d2, d3;
     register int col, i, j;
     static int *pivot = NULL;
                  /* contains the pivot row of each column */
 
 
 
     if (pivot == NULL) pivot = create_int_vector (G_d + 1);
     
     p1=A;                         /* search pivot of first column */
     pivot[0]=0; 
     d3=fabs(d1=*p1);
     for (i=0; i<=LastPlane_; i++) {
         d2=fabs(*p1);
 #if PIVOTING_LASS == 0
         if (d2>=MIN_PIVOT_LASS) {pivot[0]=i; d1=*p1; break;};
 #endif
         if (d2>d3) { pivot[0]=i; d1=*p1; d3=d2; };
         p1+=(d+1);
     }
     /* copy pivot row into planescopy */
     p1=A+pivot[0]*(d+1)+1;   
     p2=planescopy+pivot[0]*(d+1)+1;
     for (i=1,d2=1.0/d1; i<=d; i++,p1++,p2++) *p2 = (*p1)*d2;
     /* complete first pivoting and copying */
     p1=A+1;                          
     p2=planescopy+1;
     for (i=0; i<=LastPlane_; i++, p1++, p2++) {
         if (i==pivot[0]) {
             p1+=d;
             p2+=d;
             continue;   /* pivot row already done */
         }
         d1=*(p1-1); 
         p3=planescopy+pivot[0]*(d+1)+1;
         for (j=1; j<=d; j++, p1++, p2++, p3++) (*p2)=(*p1)-d1*(*p3);
     }
     
     /* subsequent elimination below */
   
     for (col=1;col<d;col++) {
         for (i=0;i<=LastPlane_;i++)       /* search first row not already used as pivot row*/
             if (notInPivot(pivot,col,i)) {
                 pivot[col]=i; 
                 break; 
             }
         p1=planescopy+i*(d+1)+col;               /* search subsequent pivot row */
         d3=fabs(d1=*p1);
         for (; i<=LastPlane_; i++, p1+=(d+1))  
             if (notInPivot(pivot,col,i)) {
                 d2=fabs(*(p1));
 #if PIVOTING_LASS == 0
                 if (d2>=MIN_PIVOT_LASS) {
                     pivot[col]=i; 
                     d1=*p1;
                     break; 
                 }
 #endif
                 if (d2>d3) { 
                     pivot[col]=i;
                     d1=*p1;
                     d3=d2;
                 }
             };
         /* update pivot row */
         p1=planescopy+pivot[col]*(d+1)+col+1;
         d2=1.0/d1;
         for (j=col+1; j<=d; j++, p1++) (*p1) *= d2;
         if (col==(d-1)) break;   /* the rest is not needed in the last case */
         /* update rest of rows */
         p1=planescopy+col+1;
         p2=planescopy+pivot[col]*(d+1)+col+1;
         for (i=0; i<=LastPlane_; i++, p1+=(col+1)) {
             if (!notInPivot(pivot,col+1,i)) {
                 p1+=d-col;
                 continue;
             }
             d1=*(p1-1);
             for (j=col+1; j<=d; j++, p1++, p2++) *p1=(*p1)-d1*(*p2);
             p2-=d-col;
         }
     };
 
     /* compute x* by backward substitution; result goes into rhs of planescopy */
   
     for (i=d-2; 0<=i; i--){
         p1=planescopy+pivot[i]*(d+1)+d;
         p2=p1-d+i+1;
         for (j=i+1; j<d; j++, p2++)
             *(p1)-= (*p2)*(*(planescopy+pivot[j]*(d+1)+d));
     }
  
     /* compute shifted b  */
 
     for (i=0; i<=LastPlane_; i++) {
         p1=A+i*(d+1);
         p2=p1+d;
         if (notInPivot(pivot,d,i)) 
             for (j=0; j<d; j++,p1++) {
                 *p2 -= (*p1)*(*(planescopy+pivot[j]*(d+1)+d));
             }
         else *p2=0;
     }
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void tree_out	(	T_Tree **	ppr,
		int *	pi_balance,
		T_Key	key,
		rational **	volume,
		T_Key **	keyfound,
		int	key_choice
	)

Definition at line 818 of file MP/Source/vinci_lass.c.

References compare_key(), create_key(), T_Key::d, duplicate_set(), FALSE, G_d, G_Minus1, G_OutOfMem, G_Storage, T_Key::hyperplanes, T_Key::hypervar, if(), key, KEY_PLANES, KEY_PLANES_VAR, KEY_VERTICES, my_malloc(), p1, T_Key::set, T_Tree::tree_b, T_Tree::tree_l, tree_out(), T_Tree::tree_r, TRUE, T_Key::variables, T_Key::vertices, and T_Tree::vol.

Referenced by lass(), and tree_out().

 {  T_Tree *p1, *p2;
    int    cmp;
 
    /* Are we grounded? If so, add the node here and set the rebalance flag, then exit.  */
    if (!*ppr)
    {
       if (G_OutOfMem)
       {  /* don't allocate if out of memory */
          *volume = &G_Minus1;
       }
       else
       {  (*ppr) = (T_Tree *) my_malloc (sizeof (T_Tree));
          (*ppr) -> tree_l = NULL;
          (*ppr) -> tree_r = NULL;
          (*ppr) -> tree_b = 0;
          /* copy the key into the new node */
          create_key (&((*ppr) -> key), key_choice);
          switch (key_choice)
          {
          case KEY_PLANES:
             memcpy ((*ppr) -> key.hyperplanes, key.hyperplanes,
                     (G_d + 1) * sizeof (int));
             break;
          case KEY_VERTICES:
             (*ppr) -> key.vertices.set = duplicate_set (key.vertices.set);
             (*ppr) -> key.vertices.d = key.vertices.d;
             break;
          case KEY_PLANES_VAR:
             memcpy ((*ppr) -> key.hypervar.hyperplanes, key.hypervar.hyperplanes,
                     (G_Storage + 2) * sizeof (T_LassInt));
             memcpy ((*ppr) -> key.hypervar.variables, key.hypervar.variables,
                     (G_Storage + 2) * sizeof (T_LassInt));
             break;
          }
          (*ppr) -> vol = -1;       /* to recognise that element is newly created */
          *volume = &((*ppr) -> vol);
          *pi_balance = TRUE;
       }
       return;
    }
 
    cmp = compare_key ((*ppr) -> key, key, key_choice);
 
    /* if LESS, prepare to move to the left. */
    if (cmp < 0)
    {
       tree_out (&((*ppr) -> tree_l), pi_balance, key, volume, keyfound, key_choice);
       if (*pi_balance)
       {  /* left branch has grown longer */
          switch ((*ppr) -> tree_b)
          {
          case 1:  /* right branch WAS longer; balance is ok now */
                   /* LESS: case 1.. balance restored implicitly */
             (*ppr) -> tree_b = 0;
             *pi_balance = FALSE;
             break;
          case 0:  /* balance WAS okay; now left branch longer */
                   /* LESS: case 0.. balance bad but still ok */
             (*ppr) -> tree_b = -1;
             break;
          case -1: /* left branch was already too long. rebalance */
             p1 = (*ppr) -> tree_l;
             if (p1 -> tree_b == -1)
             {  /* LESS: single LL */
                (*ppr) -> tree_l = p1->tree_r;
                p1 -> tree_r = (*ppr);
                (*ppr) -> tree_b = 0;
                (*ppr) = p1;
             }
             else
             {  /* LESS: real LR */
                p2 = p1 -> tree_r;
                p1 -> tree_r = p2 -> tree_l;
                p2 -> tree_l = p1;
                (*ppr) -> tree_l = p2 -> tree_r;
                p2 -> tree_r = (*ppr);
                if (p2 -> tree_b == -1)
                   (*ppr) -> tree_b = 1;
                else (*ppr) -> tree_b = 0;
                if (p2->tree_b == 1)
                   p1 -> tree_b = -1;
                else p1 -> tree_b = 0;
                (*ppr) = p2;
             }
             (*ppr) -> tree_b = 0;
             *pi_balance = FALSE;
          } /* switch */
       } /* if */
    } /* cmp < 0 */
 
    /* if MORE, prepare to move to the right. */
    else if (cmp > 0)
    {
       tree_out (&((*ppr) -> tree_r), pi_balance, key, volume, keyfound, key_choice);
       if (*pi_balance)
       {  /* right branch has grown longer */
          switch ((*ppr) -> tree_b)
          {
          case -1: /* MORE: balance was off, fixed implicitly */
             (*ppr) -> tree_b = 0;
             *pi_balance = FALSE;
             break;
          case 0:  /* MORE: balance was okay, now off but ok */
             (*ppr)->tree_b = 1;
             break;
          case 1:  /* MORE: balance was off, need to rebalance */
             p1 = (*ppr) -> tree_r;
             if (p1 -> tree_b == 1)
             {  /* MORE: single RR */
                (*ppr) -> tree_r = p1 -> tree_l;
                p1 -> tree_l = (*ppr);
                (*ppr) -> tree_b = 0;
                (*ppr) = p1;
             }
             else
             {  /* MORE: real RL */
                p2 = p1 -> tree_l;
                p1 -> tree_l = p2 -> tree_r;
                p2 -> tree_r = p1;
                (*ppr) -> tree_r = p2 -> tree_l;
                p2 -> tree_l = (*ppr);
                if (p2 -> tree_b == 1)
                   (*ppr) -> tree_b = -1;
                else (*ppr) -> tree_b = 0;
                if (p2 -> tree_b == -1)
                   p1 -> tree_b = 1;
                else p1 -> tree_b = 0;
                (*ppr) = p2;
             }
             (*ppr) -> tree_b = 0;
             *pi_balance = FALSE;
          } /* switch */
       } /* if */
    } /* cmp > 0 */
 
    /* not less, not more: this is the same key! give volume back! */
    else
    {
       *pi_balance = FALSE;
       *volume = &((*ppr) -> vol);
       *keyfound = &((*ppr) -> key);
    }
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void VOLUME_LASSERRE_FILE	(	double	planes2[7][4],
		rational *	volume
	)

Definition at line 737 of file MP/Source/vinci_lass.c.

References All_index, compact(), create_key(), FALSE, free_key(), G_d, G_Hyperplanes, G_m, G_MemRes, G_OutOfMem, G_Storage, T_Key::hyperplanes, T_Key::hypervar, i, j, KEY_PLANES_VAR, lass(), my_malloc(), p2c, Pivot, pivotrow, planescopy, and T_Key::variables.

   {  int i, j;
 
   G_d = 3;
   G_m = 7;
   G_Storage = 0;
 
    /* reserves memory space for the global variable G_Hyperplanes; G_m and G_d must be  */
    /* set correctly                                                                     */
    
   G_Hyperplanes = (real **) my_malloc (G_m * sizeof (real *));
    
   for (i = 0; i < G_m; i++)
     G_Hyperplanes [i] = (real *) my_malloc ((G_d + 1) * sizeof (real));
   /* The last entry is needed for the right hand side. */
 
   /* Transfer input from Fortran to C array */
 
   for( i = 0; i <= G_d; i++)
     for( j= 0; j < G_m; j++){
       G_Hyperplanes [j][i] = planes2[j][i];
    /*     printf("\n%10.3e",planes2[j][i]); */
     }
   
 /*    fprintf (stderr, "\n****** Hyperplanes are: ******\n"); */
 /*    for (j = 0; j < G_m; j++) */
 /*      {  fprintf (stderr, "  Hyperplane [%i]:  ", j); */
 /*      for (i=0; i < G_d; i++) fprintf (stderr, "%10.3lf", G_Hyperplanes [j] [i]); */
 /*      fprintf (stderr,"\n"); */
 /*      fprintf (stderr, " :%10.3lf", G_Hyperplanes [j] [G_d]); */
 /*      fprintf (stderr,"\n"); */
 /*      } */
 
 /*     if (G_Storage > G_d - 3) */
 /*        G_Storage = G_d - 3; */
       /* necessary to prevent memory waste because in the tree arrays of length         */
       /* G_Storage + 2 are allocated                                                    */
       
    pivotrow = (rational *) my_malloc ((G_d + 1) * sizeof (rational));
    All_index = (T_LassInt *) my_malloc ((G_m + 1) * sizeof (T_LassInt));
    Pivot = (T_LassInt *) my_malloc ((G_d + 1) * sizeof (T_LassInt));
    p2c = (int **) my_malloc (G_d * sizeof (int *));
    for (i=0; i<G_d; i++){
        p2c[i] = (int *) my_malloc (2 * sizeof (int));
    }
    G_OutOfMem=FALSE;
    G_MemRes=malloc(G_d*G_d*G_m*sizeof(rational)); /* memory reserve; */
 
    A=compact();
    planescopy=compact();
    tree_volumes = NULL;
    create_key (&key, KEY_PLANES_VAR);
    key.hypervar.hyperplanes [0] = G_m + 1;
    key.hypervar.variables [0] = G_d + 1;
    All_index[0]=G_m+2;  /* initialization (end mark) */
    Pivot[0]=G_m+2;      /* initialization (end mark) */
    *volume = lass (A, G_m-1, G_d);
     
    if (!G_OutOfMem) free (G_MemRes);
    free_key (key, KEY_PLANES_VAR);
 }

Here is the call graph for this function:

void volume_lasserre_file	(	double	planes2[7][4],
		rational *	volume
	)

Definition at line 800 of file MP/Source/vinci_lass.c.

References VOLUME_LASSERRE_FILE().

Referenced by volratio().

                                                                      {
   VOLUME_LASSERRE_FILE( planes2, volume);
 }

Here is the call graph for this function:

Here is the caller graph for this function:

void VOLUME_LASSERRE_FILE_	(	double	planes2[7][4],
		rational *	volume
	)

Definition at line 804 of file MP/Source/vinci_lass.c.

References VOLUME_LASSERRE_FILE().

                                                                      {
   VOLUME_LASSERRE_FILE( planes2, volume);
 }

Here is the call graph for this function:

void volume_lasserre_file_	(	double	planes2[7][4],
		rational *	volume
	)

Definition at line 808 of file MP/Source/vinci_lass.c.

References VOLUME_LASSERRE_FILE().

                                                                      {
   VOLUME_LASSERRE_FILE( planes2, volume);
 }

Here is the call graph for this function:

Variable Documentation

rational* A

Definition at line 67 of file MP/Source/vinci_lass.c.

T_LassInt* All_index

Definition at line 69 of file MP/Source/vinci_lass.c.

T_VertexSet* face

static

Definition at line 79 of file MP/Source/vinci_lass.c.

int G_d

Definition at line 50 of file MP/Source/vinci_lass.c.

Referenced by add_hypervar(), compact(), create_key(), delete_hypervar(), free_key(), lass(), shift_P(), tree_out(), and VOLUME_LASSERRE_FILE().

real** G_Hyperplanes = NULL

Definition at line 53 of file MP/Source/vinci_lass.c.

Referenced by compact(), and VOLUME_LASSERRE_FILE().

T_VertexSet* G_Incidence = NULL

Definition at line 54 of file MP/Source/vinci_lass.c.

int G_m

Definition at line 51 of file MP/Source/vinci_lass.c.

Referenced by add_hypervar(), add_reduced_index(), compact(), compare_key(), del_original(), del_original_indices(), delete_hypervar(), lass(), and VOLUME_LASSERRE_FILE().

void* G_MemRes = NULL

Definition at line 62 of file MP/Source/vinci_lass.c.

Referenced by lass(), my_malloc(), and VOLUME_LASSERRE_FILE().

rational G_Minus1 = -1

Definition at line 63 of file MP/Source/vinci_lass.c.

Referenced by tree_out().

boolean G_OutOfMem = TRUE

Definition at line 61 of file MP/Source/vinci_lass.c.

Referenced by lass(), my_malloc(), tree_out(), and VOLUME_LASSERRE_FILE().

int G_Precomp = 0

Definition at line 57 of file MP/Source/vinci_lass.c.

int G_RandomSeed = 4

Definition at line 59 of file MP/Source/vinci_lass.c.

int G_Storage = -1

Definition at line 58 of file MP/Source/vinci_lass.c.

Referenced by create_key(), free_key(), lass(), tree_out(), and VOLUME_LASSERRE_FILE().

T_VertexSet G_Vertices

Definition at line 55 of file MP/Source/vinci_lass.c.

T_Key key

static

Definition at line 76 of file MP/Source/vinci_lass.c.

Referenced by tree_out().

T_Key * keyfound

static

Definition at line 76 of file MP/Source/vinci_lass.c.

int** p2c

Definition at line 71 of file MP/Source/vinci_lass.c.

T_LassInt* Pivot

Definition at line 70 of file MP/Source/vinci_lass.c.

rational* pivotrow

Definition at line 68 of file MP/Source/vinci_lass.c.

rational* planescopy

Definition at line 75 of file MP/Source/vinci_lass.c.

T_Tree* tree_volumes

static

Definition at line 78 of file MP/Source/vinci_lass.c.

Macros

Functions

Variables

Macro Definition Documentation

Function Documentation

Variable Documentation