zpspecial.hpp

Go to the documentation of this file.

#ifndef ALREADY_zpspecial
#define ALREADY_zpspecial

#ifdef MON_PREMIER

#ifdef tableop
Z<MON_PREMIER,short int> tableplus[2*MON_PREMIER];
Z<MON_PREMIER,short int> tablemoins[2*MON_PREMIER];

void inittableop()
{
  for(int i=0;i<MON_PREMIER;i++)
    {
      tableplus[i].rep=i;
      tableplus[MON_PREMIER+i].rep=i;
    }
  for(int i=0;i<MON_PREMIER;i++)
    {
      tablemoins[i].rep=i;
      tablemoins[MON_PREMIER+i].rep=MON_PREMIER-i;
    }
}
template<>
Z<MON_PREMIER,short int>&  
Z<MON_PREMIER,short int>::operator+=(const Z<MON_PREMIER,short int>  z) 
{
  rep=tableplus[rep+z.rep].rep;
  return *this;
}
template<>
Z<MON_PREMIER,short int > 
Z<MON_PREMIER,short int >::operator+(const  Z<MON_PREMIER,short int> &b) const
  {
    return tableplus[rep+b.rep];
  }
template<>
Z<MON_PREMIER,short int>&  
Z<MON_PREMIER,short int>::operator-=(const Z<MON_PREMIER,short int>  z) 
{
  rep=tablemoins[rep+z.rep].rep;
  return *this;
}

#else


template<>
Z<MON_PREMIER,short int>&  
Z<MON_PREMIER,short int>::operator+=(const Z<MON_PREMIER,short int>  z) 
  { 
    int tmp=rep;
    if(tmp) cout<<"";
    (tmp += z.rep); 
    if(tmp>=MON_PREMIER) tmp-=MON_PREMIER;
    rep=tmp;
    //if(rep<0) {cout<<"+="<<endl;exit(1);};
    return *this;
  }
 
template<>
Z<MON_PREMIER,short int > 
Z<MON_PREMIER,short int >::operator+(const  Z<MON_PREMIER,short int> &b) const
  {
    Z<MON_PREMIER,short int> res; int tmp=rep;
    tmp+=b.rep; if (tmp>MON_PREMIER) tmp-=MON_PREMIER; 
    res.rep=tmp;
    //if(res.rep<0) {cout<<"+"<<endl;exit(1);};
return res;}

template<>
Z<MON_PREMIER,short int>&  
Z<MON_PREMIER,short int>::operator-=(const Z<MON_PREMIER,short int>  z) 
  { 
    //cout<<",kkln,nk"<<rep<<endl;
    //cout<<"desert"<<z.rep<<endl;
    int tmp=rep;
    tmp -= z.rep; 
    if(tmp<0) tmp+=MON_PREMIER;
    rep=tmp;
    //if(rep<0) {cout<<"+="<<endl;exit(1);};
    return *this;
  }

template<>
Z<MON_PREMIER,short int > 
Z<MON_PREMIER,short int >::operator-(const  Z<MON_PREMIER,short int> &b) const
  {
    Z<MON_PREMIER,short int> res; int tmp=rep;
    tmp=(tmp-b.rep); if (tmp<0) tmp+=MON_PREMIER; 
    res.rep=tmp;
    //if(res.rep<0) {cout<<"+"<<endl;exit(1);};
return res;}

#endif 

template<>
Z<MON_PREMIER,short int>&  Z<MON_PREMIER,short int>::operator*=(const Z<MON_PREMIER,short int> z)
{int tmp=rep; rep=(short int)((tmp * z.rep)% MON_PREMIER); return *this;}


template<>
inline
Z<MON_PREMIER,short int> Z<MON_PREMIER,short int>::operator*( const  Z<MON_PREMIER,short int> &b) const
{
  int tmp=rep;
  Z<MON_PREMIER,short int> res; 
  res.rep=(tmp*b.rep)%MON_PREMIER; 
  return res;
}
 
template<>
Z<MON_PREMIER,short int>& Z<MON_PREMIER,short int>::operator/=(const Z<MON_PREMIER,short int> z) 
{ 
  static int flag=1;
  static int *inverse=NULL;
#ifdef PETITFERMAT
  static int expon[MON_PREMIER];
#endif //
  //  int r=z.rep, a=z.rep, b=p,u1=0,u2=0, q;
  if (flag)
    {
      inverse=(int*)malloc(MON_PREMIER*sizeof(int));
      
      cout<<"je malloc tout ca"<<endl;
      for(int i=0;i<MON_PREMIER;i++)
    inverse[i]=-1;
    //  {
//        int inversez;
//        for(inversez=0;((inversez*i % p)!=1);inversez++);
//        inverse[i]=inversez;
//      }
      flag=0;
    }
  //cout<<"z.rep "<<z.rep<<endl;
  if(z==0)
    {std::cerr<<"division by zero"<<endl;exit(-10);}
  if(inverse[(short int)z.rep]==-1) 
    {
      int inversez;
#ifndef PETITFERMAT
 {
   for(inversez=0;((inversez*z.rep % MON_PREMIER)!=1);inversez++);
   int bez_1_1=1,bez_1_2=0;
   int bez_2_1=0,bez_2_2=1;
   
 }
#else
      int comp=1;
      expon[1]=z.rep;
      for(;comp*2<MON_PREMIER;comp*=2)
    expon[2*comp]=(expon[comp]*expon[comp]) %MON_PREMIER;
      inversez=expon[comp];
      for(int total=comp;total<MON_PREMIER-2;total+=comp)
    {
      for(;comp+total>MON_PREMIER-2;comp/=2);
      inversez=(inversez * expon[comp]) %MON_PREMIER;
    }
#endif //
      inverse[z.rep]=inversez;
    }
  flag=rep *inverse[z.rep];
  rep =flag % MON_PREMIER;
  flag=0;
  return *this;
}

//void convert(Z<MON_PREMIER,int> &toto,char * chaine)
//{
//  toto=Z<MON_PREMIER,int>(chaine);
//}
template<typename T>
void convert(Z<MON_PREMIER,T> &toto,char * chaine)
{
  toto=Z<MON_PREMIER,T>(chaine);
}

template<typename T>
int
dpivotL(
        const int  jcol,     /* in */
        const Z<MON_PREMIER,T> u,      /* in - diagonal pivoting threshold */
        int        *usepr,   /* re-use the pivot sequence given by perm_r/iperm_r */
        int        *perm_r,  /* may be modified */
        int        *iperm_r, /* in - inverse of perm_r */
        int        *iperm_c, /* in - used to find diagonal of Pc*A*Pc' */
        int        *pivrow,  /* out */
        GlobalLU_t<Z<MON_PREMIER,T> > *Glu      /* modified - global LU data structures */
       )
{
/*
 * Purpose
 * =======
 *   Performs the numerical pivoting on the current column of L,
 *   and the CDIV operation.
 *
 *   Pivot policy:
 *   (1) Compute thresh = u * max_(i>=j) abs(A_ij);
 *   (2) IF user specifies pivot row k and abs(A_kj) >= thresh THEN
 *           pivot row = k;
 *       ELSE IF abs(A_jj) >= thresh THEN
 *           pivot row = j;
 *       ELSE
 *           pivot row = m;
 * 
 *   Note: If you absolutely want to use a given pivot order, then set u=0.0.
 *
 *   Return value: 0      success;
 *                 i > 0  U(i,i) is exactly zero.
 *
 */
    int          fsupc;     /* first column in the supernode */
    int          nsupc;     /* no of columns in the supernode */
    int          nsupr;     /* no of rows in the supernode */
    int          lptr;      /* points to the starting subscript of the supernode */
    int          pivptr, old_pivptr, diag, diagind;
    Z<MON_PREMIER,T>       pivmax, rtemp, thresh;
    Z<MON_PREMIER,T>       temp;
    Z<MON_PREMIER,T>       *lu_sup_ptr; 
    Z<MON_PREMIER,T>       *lu_col_ptr;
    int          *lsub_ptr;
    int          isub, icol, k, itemp;
    int          *lsub, *xlsub;
    Z<MON_PREMIER,T>       *lusup;
    int          *xlusup;
    extern SuperLUStat_t SuperLUStat;
    flops_t  *ops = SuperLUStat.ops;

    /* Initialize pointers */
    lsub       = Glu->lsub;
    xlsub      = Glu->xlsub;
    lusup      = Glu->lusup;
    xlusup     = Glu->xlusup;
    fsupc      = (Glu->xsup)[(Glu->supno)[jcol]];
    nsupc      = jcol - fsupc;          /* excluding jcol; nsupc >= 0 */
    lptr       = xlsub[fsupc];
    nsupr      = xlsub[fsupc+1] - lptr;
    lu_sup_ptr = &lusup[xlusup[fsupc]]; /* start of the current supernode */
    lu_col_ptr = &lusup[xlusup[jcol]];  /* start of jcol in the supernode */
    lsub_ptr   = &lsub[lptr];   /* start of row indices of the supernode */

#ifdef DEBUG
if ( jcol == MIN_COL ) {
    printf("Before cdiv: col %d\n", jcol);
    for (k = nsupc; k < nsupr; k++) 
    printf("  lu[%d] %f\n", lsub_ptr[k], lu_col_ptr[k]);
}
#endif
    
    /* Determine the largest abs numerical value for partial pivoting;
       Also search for user-specified pivot, and diagonal element. */
    if ( *usepr ) *pivrow = iperm_r[jcol];
    diagind = iperm_c[jcol];
    pivmax =(Z<MON_PREMIER,T>)0;
    pivptr = nsupc;
    diag = EMPTY;
    old_pivptr = nsupc;
    for (isub = nsupc; isub < nsupr; ++isub) {
    rtemp = lu_col_ptr[isub];
    
    if ( rtemp != 0 ) {
      pivmax = rtemp;
      pivptr = isub;
      break;
    }
    if ( *usepr && lsub_ptr[isub] == *pivrow ) old_pivptr = isub;
    if ( lsub_ptr[isub] == diagind ) diag = isub;
    }
    //cout<<"comparaison dans dpivvotL "<<pivmax<<endl;
    /* Test for singularity */
    if ( pivmax ==0) {
    *pivrow = lsub_ptr[nsupc];
    //cout<<"degenere et pivrow "<<*pivrow<<endl;
    perm_r[*pivrow] = jcol;
    *usepr = 0;
    return (jcol+1);
    }

    thresh = u * pivmax;
    
    /* Choose appropriate pivotal element by our policy. */
    if ( *usepr ) {
      //fabs!!!!  
      rtemp = lu_col_ptr[old_pivptr];
      //fabs!!!!!!!!!!!!!!!!! faire gaffe en generic 
      if ( rtemp != (Z<MON_PREMIER,T>)0 )
        pivptr = old_pivptr;
    else
        *usepr = 0;
    }
    if ( *usepr == 0 ) {
    /* Use diagonal pivot? */
    if ( diag >= 0 ) { /* diagonal exists */
        rtemp = lu_col_ptr[diag];
        if ( rtemp != (Z<MON_PREMIER,T>)0  ) pivptr = diag;
        }
    *pivrow = lsub_ptr[pivptr];
    }
    
    /* Record pivot row */
    //    cout<<"piv row "<<*pivrow<<" jcol "<<jcol<<endl;
    perm_r[*pivrow] = jcol;
    
    /* Interchange row subscripts */
    if ( pivptr != nsupc ) {
      itemp = lsub_ptr[pivptr];
      lsub_ptr[pivptr] = lsub_ptr[nsupc];
      lsub_ptr[nsupc] = itemp;

    /* Interchange numerical values as well, for the whole snode, such 
     * that L is indexed the same way as A.
     */
    for (icol = 0; icol <= nsupc; icol++) {
        itemp = pivptr + icol * nsupr;
        temp = lu_sup_ptr[itemp];
        lu_sup_ptr[itemp] = lu_sup_ptr[nsupc + icol*nsupr];
        lu_sup_ptr[nsupc + icol*nsupr] = temp;
    }
    } /* if */

    /* cdiv operation */
    ops[FACT] += nsupr - nsupc;
    temp = (Z<MON_PREMIER,T>)1 / lu_col_ptr[nsupc];
    for (k = nsupc+1; k < nsupr; k++) 
    lu_col_ptr[k] *= temp;

    return 0;
}
int Iszero(const Z<MON_PREMIER,short int> & c)
{return c==(short int)0;}

template<typename T>
inline int Iszero(const Z<MON_PREMIER,T> & c)
{return c==(T)0;}  


#endif
#endif //ALREADY_zpspecial