vinci_lass.c File Reference

#include "vinci.h"

Include dependency graph for vinci_lass.c:

Go to the source code of this file.

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)
static int	compare_key (T_Key key1, T_Key key2, int key_choice)
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)

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 vinci_lass.c.

Referenced by norm_and_clean_constraints().

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

Definition at line 28 of file 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 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 vinci_lass.c.

Referenced by lass().

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

Definition at line 29 of file vinci_lass.c.

Referenced by lass().

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

Definition at line 25 of file vinci_lass.c.

Referenced by QualityMetric::average_metric_and_weights(), QualityAssessor::get_QAFunction_name(), lass(), and QualityAssessor::loop_over_mesh().

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

Definition at line 34 of file vinci_lass.c.

Function Documentation

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

Definition at line 168 of file vinci_lass.c.

{  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;
      }
   }
}

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

static

Definition at line 258 of file 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 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 818 of file 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 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 vinci_lass.c.

{
   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;
   }
}

static void del_original ( T_LassInt  base, T_LassInt *  indices  )

static

Definition at line 308 of file 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 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 vinci_lass.c.

{  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];
      }
   }
}

T_VertexSet duplicate_set

(

T_VertexSet

s

)

Definition at line 1028 of file vinci_lass.c.

{  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;
}

void free_int_vector	(	int *	v,
		int	n
	)

Definition at line 248 of file vinci_lass.c.

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

void free_key	(	T_Key	key,
		int	key_choice
	)

Definition at line 151 of file vinci_lass.c.

{
   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;
   }
}

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

static

Definition at line 572 of file 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 vinci_lass.c.

{
   free (pointer);
}

void* my_malloc

(

long int

size

)

Definition at line 90 of file vinci_lass.c.

{  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;
}

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 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 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 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 vinci_lass.c.

References All_index, and del_original().

Referenced by lass().

{   del_original(baserow, All_index); }

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 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 870 of file vinci_lass.c.

{  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);
   }
}

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

Definition at line 737 of file 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.

Referenced by volume_lasserre_file(), VOLUME_LASSERRE_FILE_(), and volume_lasserre_file_().

  {  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:

Here is the caller graph for this function:

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

Definition at line 800 of file vinci_lass.c.

                                                                     {
  VOLUME_LASSERRE_FILE( planes2, volume);
}

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

Definition at line 804 of file 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 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 vinci_lass.c.

Referenced by CImg< uintT >::_solve(), CImgDisplay::_update_iskey(), TargetCalculator::compute_reference_corner_matrices(), PatchData::get_average_Lambda_3d(), PatchData::get_barrier_delta(), Mesquite::h_fcn_2e(), Mesquite::h_fcn_2i(), Mesquite::h_fcn_3e(), Mesquite::h_fcn_3i(), h_fcn_ridft2(), h_fcn_ridft3(), Mesquite::h_gdft_2(), Mesquite::h_gdft_2_v0(), Mesquite::h_gdft_2_v1(), Mesquite::h_gdft_2_v2(), Mesquite::h_gdft_3(), Mesquite::h_gdft_3_v0(), Mesquite::h_gdft_3_v1(), Mesquite::h_gdft_3_v2(), Alg_Metric_Base_2::initialize(), Alg_Metric_Base_3::initialize(), CImg< uintT >::invert(), lass(), norm_and_clean_constraints(), Mesquite::QR(), serialbs95_(), serialbs95mat_(), shift_P(), CImg< uintT >::solve(), and write_attr_CGNS().

T_LassInt* All_index

Definition at line 69 of file vinci_lass.c.

Referenced by lass(), norm_and_clean_constraints(), rm_original_inElAll_index(), and VOLUME_LASSERRE_FILE().

T_VertexSet* face

static

Definition at line 79 of file vinci_lass.c.

Referenced by Mesh::addFace(), Node::addFace(), block::buildFaces(), Face::create(), TopologyInfo::face_vertices(), mblockData::getPatchListWith(), Element::setFace(), Remesher_Simmetrix::simmesh_to_window(), and FluRegion::WriteGSPASCII().

int G_d

Definition at line 50 of file vinci_lass.c.

real** G_Hyperplanes = NULL

Definition at line 53 of file vinci_lass.c.

T_VertexSet* G_Incidence = NULL

Definition at line 54 of file vinci_lass.c.

int G_m

Definition at line 51 of file vinci_lass.c.

void* G_MemRes = NULL

Definition at line 62 of file vinci_lass.c.

rational G_Minus1 = -1

Definition at line 63 of file vinci_lass.c.

boolean G_OutOfMem = TRUE

Definition at line 61 of file vinci_lass.c.

int G_Precomp = 0

Definition at line 57 of file vinci_lass.c.

int G_RandomSeed = 4

Definition at line 59 of file vinci_lass.c.

int G_Storage = -1

Definition at line 58 of file vinci_lass.c.

T_VertexSet G_Vertices

Definition at line 55 of file vinci_lass.c.

T_Key key

static

Definition at line 76 of file vinci_lass.c.

Referenced by CImg< uintT >::_display(), CImg< uintT >::_display_object3d(), CImg< uintT >::_get_select(), StopWatchCollection::add(), conchr(), dchar(), CImg< uintT >::display_graph(), StopWatchCollection::get_key(), CImg< uintT >::get_select_graph(), locchr(), matmodel_hyperelastic(), rchar(), rflu_modboundlists::rflu_buildbvertexlists(), rflu_modtopologyutils::rflu_buildcellvertlist(), rflu_modtopologyutils::rflu_buildconnvertlist(), rflu_modedgelist::rflu_buildedgelist(), rflu_modfacelist::rflu_buildfacelist(), rflu_modtopologyutils::rflu_buildfacevertlist(), rflu_computedcuhreinfo(), rflu_computeintegral1oles(), rflu_computeintegral2oles(), rflu_computeintegral4oles(), rflu_computeintegral5oles(), rflu_modcobalt::rflu_convcobalt2rocflu(), rflu_modrenumberlist::rflu_denumberlist(), rflu_modhashtable::rflu_hashbuildkey(), rflu_modhashtable::rflu_hashbuildkey1(), rflu_modhashtable::rflu_hashedge(), rflu_modhashtable::rflu_hashface(), rflu_modhashtable::rflu_hashfuncprimary(), rflu_modhashtable::rflu_hashfuncsecondary(), rflu_modhashtable::rflu_hashvertex(), rflu_modhashtable::rflu_hashvertexfancy(), rflu_modpartitionregion::rflu_part_buildvertexlists(), rflu_modrenumberlist::rflu_renumberlist(), rflu_modstencilsutils::rflu_sortbfaces(), rflu_modhashtable::rflu_unhashbface(), rvav_readcomparisonssection(), tree_out(), and updateHash().

T_Key * keyfound

static

Definition at line 76 of file vinci_lass.c.

int** p2c

Definition at line 71 of file vinci_lass.c.

Referenced by lass(), and VOLUME_LASSERRE_FILE().

T_LassInt* Pivot

Definition at line 70 of file vinci_lass.c.

Referenced by lass(), and VOLUME_LASSERRE_FILE().

rational* pivotrow

Definition at line 68 of file vinci_lass.c.

Referenced by lass(), and VOLUME_LASSERRE_FILE().

rational* planescopy

Definition at line 75 of file vinci_lass.c.

Referenced by shift_P(), and VOLUME_LASSERRE_FILE().

T_Tree* tree_volumes

static

Definition at line 78 of file vinci_lass.c.