include/borderbasix/dpivmacrev_z.cc

Go to the documentation of this file.

#ifndef ALREADY_dpivmacrev
#define ALREADY_dpivmacrev

template<>
int
dpivotL(
        const int  jcol,     /* in */
        const macrev_Z 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<macrev_Z> *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;
    macrev_Z       pivmax, rtemp, thresh;
    macrev_Z       temp;
    macrev_Z       *lu_sup_ptr; 
    macrev_Z       *lu_col_ptr;
    int          *lsub_ptr;
    int          isub, icol, k, itemp;
    int          *lsub, *xlsub;
    macrev_Z       *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 =(macrev_Z)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 != (macrev_Z)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 != (macrev_Z)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 = (macrev_Z)1 / lu_col_ptr[nsupc];
    for (k = nsupc+1; k < nsupr; k++) 
        lu_col_ptr[k] *= temp;

    return 0;
}

void
dSetRWork(int m, int panel_size, macrev_Z *dworkptr,
         macrev_Z **dense, macrev_Z **tempv)
{
    macrev_Z zero(0);
    int i;
    int maxsuper = sp_ienv(3),
        rowblk   = sp_ienv(4);
    *dense = dworkptr;
    *tempv = *dense + panel_size*m;
    for(i=0;i<m*panel_size;i++)
      mpz_init((*dense)[i].rep);
    for(i=0;i<NUM_TEMPV(m,panel_size,maxsuper,rowblk);i++)
      mpz_init((*tempv)[i].rep);
    //dfill (*dense, m * panel_size, zero);
    //dfill (*tempv, NUM_TEMPV(m,panel_size,maxsuper,rowblk), zero);     
}

void dLUWorkFree(int m,int panel_size, int *iwork, macrev_Z *dwork
                 , GlobalLU_t<macrev_Z > *Glu)
{
  int i;
  int maxsuper = sp_ienv(3),
    rowblk   = sp_ienv(4);
  if ( Glu->MemModel == SYSTEM ) {
        SUPERLU_FREE (iwork);
        for(i=0;i<m*panel_size+NUM_TEMPV(m,panel_size,maxsuper,rowblk);i++)
          mpz_clear((dwork[i].rep));
        SUPERLU_FREE (dwork);
    } else {
        stack.used -= (stack.size - stack.top2);
        stack.top2 = stack.size;
/*      dStackCompress(Glu);  */
    }
    
    SUPERLU_FREE (expanders);
}

void* dexpand (
         int *prev_len,   /* length used from previous call */
         MemType type,    /* which part of the memory to expand */
         int len_to_copy, /* size of the memory to be copied to new store */
         int keep_prev,   /* = 1: use prev_len;
                             = 0: compute new_len to expand */
         GlobalLU_t<macrev_Z > *Glu  /* modified - global LU data structures */
        )
{
    double    EXPAND = 1.5;
    double    alpha;
    void     *new_mem, *old_mem;
    int      new_len, tries, lword, extra, bytes_to_copy;

    alpha = EXPAND;

    if ( no_expand == 0 || keep_prev ) /* First time allocate requested */
        new_len = *prev_len;
    else {
        new_len = (int)(alpha * *prev_len);
    }
    
    if ( type == LSUB || type == USUB ) lword = sizeof(int);
    else lword = sizeof(macrev_Z);

    if ( Glu->MemModel == SYSTEM ) {
        new_mem = (void *) SUPERLU_MALLOC(new_len * lword);
        if(new_mem && lword==sizeof(macrev_Z))
          for(extra=0;extra<new_len;extra++)
            mpz_init(((macrev_Z*)new_mem)[extra].rep);      
/*      new_mem = (void *) calloc(new_len, lword); */
        if ( no_expand != 0 ) {
            tries = 0;
            if ( keep_prev ) {
                if ( !new_mem ) return (NULL);
            } else {
                while ( !new_mem ) 
                  {
                    if ( ++tries > 10 ) return (NULL);
                    alpha = Reduce(alpha);
                    new_len = (int)(alpha * *prev_len);
                    new_mem = (void *) SUPERLU_MALLOC(new_len * lword);
/*                  new_mem = (void *) calloc(new_len, lword); */
                  }
                if(tries  && lword==sizeof(macrev_Z))
                  for(int extra=0;extra<new_len;extra++)
                    mpz_init(((macrev_Z*)new_mem)[extra].rep);
            }
            if ( type == LSUB || type == USUB ) {
                copy_mem_int(len_to_copy, expanders[type].mem, new_mem);
            } else {
                copy_mem_double<macrev_Z >(len_to_copy, 
                                           expanders[type].mem, new_mem);
                //cout<<"'j'en detruit ici"<<endl;
                for(unsigned int i=0;i<len_to_copy;i++)
                  mpz_clear(((macrev_Z*)expanders[type].mem)[i].rep);
            }
            
            SUPERLU_FREE (expanders[type].mem);
        }
        expanders[type].mem = (void *) new_mem;
        
    } else { /* MemModel == USER */
        if ( no_expand == 0 ) {
            new_mem = duser_malloc(new_len * lword, HEAD);
            if ( NotDoubleAlign(new_mem) &&
                (type == LUSUP || type == UCOL) ) {
                old_mem = new_mem;
                new_mem = (void *)DoubleAlign(new_mem);
                extra = (char*)new_mem - (char*)old_mem;
#ifdef DEBUG            
                printf("expand(): not aligned, extra %d\n", extra);
#endif          
                stack.top1 += extra;
                stack.used += extra;
            }
            expanders[type].mem = (void *) new_mem;
        }
        else {
            tries = 0;
            extra = (new_len - *prev_len) * lword;
            if ( keep_prev ) {
                if ( StackFull(extra) ) return (NULL);
            } else {
                while ( StackFull(extra) ) {
                    if ( ++tries > 10 ) return (NULL);
                    alpha = Reduce(alpha);
                    new_len = (int)(alpha * *prev_len);
                    extra = (new_len - *prev_len) * lword;          
                }
            }

            if ( type != USUB ) {
                new_mem = (void*)((char*)expanders[type + 1].mem + extra);
                bytes_to_copy = (char*)stack.array + stack.top1
                    - (char*)expanders[type + 1].mem;
                user_bcopy((char*)expanders[type+1].mem,
                           (char*)new_mem, bytes_to_copy);

                if ( type < USUB ) {
                    Glu->usub = (int*)
                      (expanders[USUB].mem =
                        (void*)((char*)expanders[USUB].mem + extra));
                }
                if ( type < LSUB ) {
                    Glu->lsub =(int*)(expanders[LSUB].mem =
                        (void*)((char*)expanders[LSUB].mem + extra));
                }
                if ( type < UCOL ) {
                    Glu->ucol =(macrev_Z*)( expanders[UCOL].mem =
                        (void*)((char*)expanders[UCOL].mem + extra));
                }
                stack.top1 += extra;
                stack.used += extra;
                if ( type == UCOL ) {
                    stack.top1 += extra;   /* Add same amount for USUB */
                    stack.used += extra;
                }
                
            } /* if ... */

        } /* else ... */
    }

    expanders[type].size = new_len;
    *prev_len = new_len;
    if ( no_expand ) ++no_expand;
    
    return (void *) expanders[type].mem;
    
} /* dexpand */

template<>
harewell<macrev_Z >::harewell(const int nrow,const int ncol)
{
  //Pour eviter le undefined behavior on met tous lespointeurs non initialises a NULL
  this->nrow=nrow;
  this->ncol=ncol;
  this->rowind=(int*)malloc(sizeof(int)*(this->ncol));
  this->colptr=(int*)malloc(sizeof(int)*(1+this->ncol));
  this->nzval=MAC_REV_MALLOC<macrev_Z >(sizeof(macrev_Z)*(this->ncol));
  this->nnz=this->ncol;
  for(int i=0;i<this->ncol;i++)
    {
      this->rowind[i]=i;
      this->colptr[i]=i;
      //    ((T*)this->nzval)[i]=0;
    }
  this->colptr[this->ncol]=this->nnz;
}

template<>
inline void harewell<macrev_Z >::destroystore()
{
  free(rowind);
  rowind=NULL;
  free(colptr);
  colptr=NULL;
  //  cout<<"nnz "<<nnz<<endl;
  MAC_REV_FREE<macrev_Z>(nzval,nnz*sizeof(macrev_Z));
  nzval=NULL;
  nnz=0;
}
/*
void destroystore(SuperMatrix<macrev_Z >& m)
{
  //cout<<"appel de destroystrore"<<endl;
  switch (m.Stype)
    {
    case NC:
      {
        for(int i=0;i<((NCformat*)m.Store)->nnz;i++)
          mpz_clear(((((macrev_Z*)((NCformat*)m.Store)->nzval))[i].rep));
        Destroy_CompCol_Matrix(&m);
        break;
      }
    case SC:
      {
        for(int i=0;i<((SCformat*)m.Store)->nnz;i++)
          mpz_clear(((((macrev_Z*)((SCformat*)m.Store)->nzval))[i]).rep);
        Destroy_SuperNode_Matrix(&m);
        break;
      }
    case NCP:
      {
        //cout<<"compcol permuted"<<endl;
        Destroy_CompCol_Permuted(&m);
        break;
      }
    default:
      {
        ABORT("Used Type Storage not yet supported \n");
      }
    };
  
    }*/

void remiseformatNC(SuperMatrix<macrev_Z > &M,const vector<macrev_Z > &nzval
                    ,const vector<int> &rowind,const vector<int> &colptr
                    ,int coeff)
{
  ((NCformat *)M.Store)->nnz=coeff;
  destroystore(M);
  M.Stype=NC;
  M.Store=(void *)malloc(sizeof(NCformat));
  NCformat *Mstore=(NCformat *)M.Store;
  Mstore->nnz=nzval.size();
  Mstore->nzval=(void *)malloc(sizeof(macrev_Z)*(nzval.size()));
  Mstore->rowind=(int*)malloc(sizeof(int)*(nzval.size()));
  Mstore->colptr=(int*)malloc(sizeof(int)*(colptr.size()));
  for(unsigned int i=0;i<nzval.size();++i)
    {
      Mstore->rowind[i]=rowind[i];
      mpz_init(((macrev_Z*)(Mstore->nzval))[i].rep);
      ((macrev_Z*)(Mstore->nzval))[i]=nzval[i];
    }
  for(unsigned int i=0;i<colptr.size();++i)
    Mstore->colptr[i]=colptr[i];
  return;
}


#endif //ALREADY_dpivmacrev

---