tgbgauss.cc File Reference

#include "kernel/mod2.h"
#include "misc/options.h"
#include "kernel/GBEngine/tgbgauss.h"
#include <stdlib.h>
#include "kernel/GBEngine/kutil.h"
#include "kernel/polys.h"

Go to the source code of this file.

Functions
mac_poly	mac_p_add_ff_qq (mac_poly a, number f, mac_poly b)
void	mac_mult_cons (mac_poly p, number c)
int	mac_length (mac_poly p)
void	mac_destroy (mac_poly p)
void	simple_gauss (tgb_sparse_matrix *mat, slimgb_alg *)
void	simple_gauss2 (tgb_matrix *mat)
static int	row_cmp_gen (const void *a, const void *b)

Variables
static const int	bundle_size =100

Function Documentation

mac_destroy()

void mac_destroy

(

mac_poly

p

)

Definition at line 113 of file tgbgauss.cc.

{
  mac_poly iter=p;
  while(iter)
  {
    mac_poly next=iter->next;
    nDelete(&iter->coef);
    delete iter;
    iter=next;
  }
}

mac_length()

int mac_length

(

mac_poly

p

)

Definition at line 102 of file tgbgauss.cc.

{
  int l=0;
  while(p){
    l++;
    p=p->next;
  }
  return l;
}

mac_mult_cons()

void mac_mult_cons	(	mac_poly	p,
		number	c
	)

Definition at line 91 of file tgbgauss.cc.

{
  while(p)
  {
    number m=nMult(p->coef,c);
    nDelete(&(p->coef));
    p->coef=m;
    p=p->next;
  }
}

mac_p_add_ff_qq()

mac_poly mac_p_add_ff_qq	(	mac_poly	a,
		number	f,
		mac_poly	b
	)

Definition at line 16 of file tgbgauss.cc.

{
  mac_poly erg;
  mac_poly* set_this;
  set_this=&erg;
  while((a!=NULL) &&(b!=NULL))
  {
    if (a->exp<b->exp)
    {
      (*set_this)=a;
      a=a->next;
      set_this= &((*set_this)->next);
    }
    else
    {
      if (a->exp>b->exp)
      {
        mac_poly in =new mac_poly_r();
        in->exp=b->exp;
        in->coef=nMult(b->coef,f);
        (*set_this)=in;
        b=b->next;
        set_this= &((*set_this)->next);
      }
      else
      {
        //a->exp==b->ecp
        number n=nMult(b->coef,f);
        number n2=nAdd(a->coef,n);
        nDelete(&n);
        nDelete(&(a->coef));
        if (nIsZero(n2))
        {
          nDelete(&n2);
          mac_poly ao=a;
          a=a->next;
          delete ao;
          b=b->next;
        }
        else
        {
          a->coef=n2;
          b=b->next;
          (*set_this)=a;
          a=a->next;
          set_this= &((*set_this)->next);
        }
      }
    }
  }
  if((a==NULL)&&(b==NULL))
  {
    (*set_this)=NULL;
    return erg;
  }
  if (b==NULL)
  {
    (*set_this=a);
    return erg;
  }
 
  //a==NULL
  while(b!=NULL)
  {
    mac_poly mp= new mac_poly_r();
    mp->exp=b->exp;
    mp->coef=nMult(f,b->coef);
    (*set_this)=mp;
    set_this=&(mp->next);
    b=b->next;
  }
  (*set_this)=NULL;
  return erg;
}

row_cmp_gen()

static int row_cmp_gen ( const void *  a, const void *  b  )

static

Definition at line 683 of file tgbgauss.cc.

{
  const mac_poly ap= *((mac_poly*) a);
  const mac_poly bp=*((mac_poly*) b);
  if (ap==NULL) return 1;
  if (bp==NULL) return -1;
  if (ap->exp<bp->exp) return -1;
  return 1;
}

simple_gauss()

void simple_gauss	(	tgb_sparse_matrix *	mat,
		slimgb_alg *	c
	)

Definition at line 125 of file tgbgauss.cc.

{
  int col, row;
  int* row_cache=(int*) omAlloc(mat->get_rows()*sizeof(int));
  col=0;
  row=0;
  int i;
  int pn=mat->get_rows();
  int matcol=mat->get_columns();
  int* area=(int*) omAlloc(sizeof(int)*((matcol-1)/bundle_size+1));
  const int max_area_index=(matcol-1)/bundle_size;
    //rows are divided in areas
  //if row begins with columns col, it is located in [area[col/bundle_size],area[col/bundle_size+1]-1]
  assume(pn>0);
  //first clear zeroes
  for(i=0;i<pn;i++)
  {
    if(mat->zero_row(i))
    {
      mat->perm_rows(i,pn-1);
      pn--;
      if(i!=pn){i--;}
    }
 
  }
  mat->sort_rows();
  for(i=0;i<pn;i++)
  {
      row_cache[i]=mat->min_col_not_zero_in_row(i);
      // Print("row_cache:%d\n",row_cache[i]);
  }
  int last_area=-1;
  for(i=0;i<pn;i++)
  {
    int this_area=row_cache[i]/bundle_size;
    assume(this_area>=last_area);
    if(this_area>last_area)
    {
      int j;
      for(j=last_area+1;j<=this_area;j++)
        area[j]=i;
      last_area=this_area;
    }
  }
  for(i=last_area+1;i<=max_area_index;i++)
  {
    area[i]=pn;
  }
  while(row<pn-1)
  {
    //row is the row where pivot should be
    // row== pn-1 means we have only to act on one row so no red nec.
    //we assume further all rows till the pn-1 row are non-zero
 
    //select column
 
    //col=mat->min_col_not_zero_in_row(row);
    int max_in_area;
    {
      int tai=row_cache[row]/bundle_size;
      assume(tai<=max_area_index);
      if(tai==max_area_index)
        max_in_area=pn-1;
      else
        max_in_area=area[tai+1]-1;
    }
    assume(row_cache[row]==mat->min_col_not_zero_in_row(row));
    col=row_cache[row];
 
    assume(col!=matcol);
    int found_in_row;
 
    found_in_row=row;
    BOOLEAN must_reduce=FALSE;
    assume(pn<=mat->get_rows());
    for(i=row+1;i<=max_in_area;i++)
    {
      int first;//=mat->min_col_not_zero_in_row(i);
      assume(row_cache[i]==mat->min_col_not_zero_in_row(i));
      first=row_cache[i];
      assume(first!=matcol);
      if(first<col)
      {
        col=first;
        found_in_row=i;
        must_reduce=FALSE;
      }
      else
      {
        if(first==col)
          must_reduce=TRUE;
      }
    }
    //select pivot
    int act_l=nSize(mat->get(found_in_row,col))*mat->non_zero_entries(found_in_row);
    if(must_reduce)
    {
      for(i=found_in_row+1;i<=max_in_area;i++)
      {
        assume(mat->min_col_not_zero_in_row(i)>=col);
        assume(row_cache[i]==mat->min_col_not_zero_in_row(i));
#ifndef SING_NDEBUG
        int first=row_cache[i];
        assume(first!=matcol);
#endif
        //      if((!(mat->is_zero_entry(i,col)))&&(mat->non_zero_entries(i)<act_l))
        int nz;
        if((row_cache[i]==col)&&((nz=nSize(mat->get(i,col))*mat->non_zero_entries(i))<act_l))
        {
          found_in_row=i;
          act_l=nz;
        }
 
      }
    }
    mat->perm_rows(row,found_in_row);
    int h=row_cache[row];
    row_cache[row]=row_cache[found_in_row];
    row_cache[found_in_row]=h;
 
    if(!must_reduce)
    {
      row++;
      continue;
    }
    //reduction
    //must extract content and normalize here
    mat->row_content(row);
    //mat->row_normalize(row); // row_content does normalize
 
    //for(i=row+1;i<pn;i++){
    for(i=max_in_area;i>row;i--)
    {
      int col_area_index=col/bundle_size;
      assume(col_area_index<=max_area_index);
      assume(mat->min_col_not_zero_in_row(i)>=col);
      assume(row_cache[i]==mat->min_col_not_zero_in_row(i));
#ifndef SING_NDEBUG
      int first=row_cache[i];
      assume(first!=matcol);
#endif
      if(row_cache[i]==col)
      {
 
        number c1=mat->get(i,col);
        number c2=mat->get(row,col);
        number n1=c1;
        number n2=c2;
 
        ksCheckCoeff(&n1,&n2,currRing->cf);
        //nDelete(&c1);
        n1=nInpNeg(n1);
        mat->mult_row(i,n2);
        mat->add_lambda_times_row(i,row,n1);
        nDelete(&n1);
        nDelete(&n2);
        assume(mat->is_zero_entry(i,col));
        row_cache[i]=mat->min_col_not_zero_in_row(i);
        assume(mat->min_col_not_zero_in_row(i)>col);
        if(row_cache[i]==matcol)
        {
          int index;
          index=i;
          int last_in_area;
          int this_cai=col_area_index;
          while(this_cai<max_area_index)
          {
            last_in_area=area[this_cai+1]-1;
            int h_c=row_cache[last_in_area];
            row_cache[last_in_area]=row_cache[index];
            row_cache[index]=h_c;
            mat->perm_rows(index,last_in_area);
            index=last_in_area;
            this_cai++;
            area[this_cai]--;
          }
          mat->perm_rows(index,pn-1);
          row_cache[index]=row_cache[pn-1];
          row_cache[pn-1]=matcol;
          pn--;
        }
        else
        {
          int index;
          index=i;
          int last_in_area;
          int this_cai=col_area_index;
          int final_cai=row_cache[index]/bundle_size;
          assume(final_cai<=max_area_index);
          while(this_cai<final_cai)
          {
            last_in_area=area[this_cai+1]-1;
            int h_c=row_cache[last_in_area];
            row_cache[last_in_area]=row_cache[index];
            row_cache[index]=h_c;
            mat->perm_rows(index,last_in_area);
            index=last_in_area;
            this_cai++;
            area[this_cai]--;
          }
        }
      }
      else
        assume(mat->min_col_not_zero_in_row(i)>col);
    }
//     for(i=row+1;i<pn;i++)
//     {
//       assume(mat->min_col_not_zero_in_row(i)==row_cache[i]);
//       // if(mat->zero_row(i))
//       assume(matcol==mat->get_columns());
//       if(row_cache[i]==matcol)
//       {
//         assume(mat->zero_row(i));
//         mat->perm_rows(i,pn-1);
//         row_cache[i]=row_cache[pn-1];
//         row_cache[pn-1]=matcol;
//         pn--;
//         if(i!=pn){i--;}
//       }
//     }
#ifdef TGB_DEBUG
  {
    int last=-1;
    for(i=0;i<pn;i++)
    {
      int act=mat->min_col_not_zero_in_row(i);
      assume(act>last);
    }
    for(i=pn;i<mat->get_rows();i++)
    {
      assume(mat->zero_row(i));
    }
  }
#endif
    row++;
  }
  omFree(area);
  omFree(row_cache);
}

simple_gauss2()

void simple_gauss2

(

tgb_matrix *

mat

)

Definition at line 365 of file tgbgauss.cc.

{
  int col, row;
  col=0;
  row=0;
  int i;
  int pn=mat->get_rows();
  assume(pn>0);
  //first clear zeroes
//   for(i=0;i<pn;i++)
//   {
//     if(mat->zero_row(i))
//     {
//       mat->perm_rows(i,pn-1);
//       pn--;
//       if(i!=pn){i--;}
//     }
//   }
  while((row<pn-1)&&(col<mat->get_columns())){
    //row is the row where pivot should be
    // row== pn-1 means we have only to act on one row so no red nec.
    //we assume further all rows till the pn-1 row are non-zero
 
    //select column
 
    //    col=mat->min_col_not_zero_in_row(row);
    assume(col!=mat->get_columns());
    int found_in_row=-1;
 
    //    found_in_row=row;
    assume(pn<=mat->get_rows());
    for(i=row;i<pn;i++)
    {
      //    int first=mat->min_col_not_zero_in_row(i);
      //  if(first<col)
      if(!(mat->is_zero_entry(i,col)))
      {
        found_in_row=i;
        break;
      }
    }
    if(found_in_row!=-1)
    {
    //select pivot
      int act_l=mat->non_zero_entries(found_in_row);
      for(i=i+1;i<pn;i++)
      {
        int vgl;
        assume(mat->min_col_not_zero_in_row(i)>=col);
        if((!(mat->is_zero_entry(i,col)))
        &&((vgl=mat->non_zero_entries(i))<act_l))
        {
          found_in_row=i;
          act_l=vgl;
        }
 
      }
      mat->perm_rows(row,found_in_row);
 
      //reduction
      for(i=row+1;i<pn;i++){
        assume(mat->min_col_not_zero_in_row(i)>=col);
        if(!(mat->is_zero_entry(i,col)))
        {
          number c1=nCopy(mat->get(i,col));
          c1=nInpNeg(c1);
          number c2=mat->get(row,col);
          number n1=c1;
          number n2=c2;
 
          ksCheckCoeff(&n1,&n2,currRing->cf);
          nDelete(&c1);
          mat->mult_row(i,n2);
          mat->add_lambda_times_row(i,row,n1);
          assume(mat->is_zero_entry(i,col));
        }
        assume(mat->min_col_not_zero_in_row(i)>col);
      }
      row++;
    }
    col++;
    // for(i=row+1;i<pn;i++)
//     {
//       if(mat->zero_row(i))
//       {
//         mat->perm_rows(i,pn-1);
//         pn--;
//         if(i!=pn){i--;}
//       }
//     }
  }
}

Variable Documentation

bundle_size

const int bundle_size =100

static

Definition at line 14 of file tgbgauss.cc.