include/borderbasix/dpivmacrev_z.cc File Reference

Go to the source code of this file.

Functions

Function Documentation

void* dexpand ( int *  prev_len,
MemType  type,
int  len_to_copy,
int  keep_prev,
GlobalLU_t< macrev_Z > *  Glu 
)

Definition at line 194 of file dpivmacrev_z.cc.

00198                                 : use prev_len;
00199                              = 0: compute new_len to expand */
00200          GlobalLU_t<macrev_Z > *Glu  /* modified - global LU data structures */
00201         )
00202 {
00203     double    EXPAND = 1.5;
00204     double    alpha;
00205     void     *new_mem, *old_mem;
00206     int      new_len, tries, lword, extra, bytes_to_copy;
00207 
00208     alpha = EXPAND;
00209 
00210     if ( no_expand == 0 || keep_prev ) /* First time allocate requested */
00211         new_len = *prev_len;
00212     else {
00213         new_len = (int)(alpha * *prev_len);
00214     }
00215     
00216     if ( type == LSUB || type == USUB ) lword = sizeof(int);
00217     else lword = sizeof(macrev_Z);
00218 
00219     if ( Glu->MemModel == SYSTEM ) {
00220         new_mem = (void *) SUPERLU_MALLOC(new_len * lword);
00221         if(new_mem && lword==sizeof(macrev_Z))
00222           for(extra=0;extra<new_len;extra++)
00223             mpz_init(((macrev_Z*)new_mem)[extra].rep);      
00224 /*      new_mem = (void *) calloc(new_len, lword); */
00225         if ( no_expand != 0 ) {
00226             tries = 0;
00227             if ( keep_prev ) {
00228                 if ( !new_mem ) return (NULL);
00229             } else {
00230                 while ( !new_mem ) 
00231                   {
00232                     if ( ++tries > 10 ) return (NULL);
00233                     alpha = Reduce(alpha);
00234                     new_len = (int)(alpha * *prev_len);
00235                     new_mem = (void *) SUPERLU_MALLOC(new_len * lword);
00236 /*                  new_mem = (void *) calloc(new_len, lword); */
00237                   }
00238                 if(tries  && lword==sizeof(macrev_Z))
00239                   for(int extra=0;extra<new_len;extra++)
00240                     mpz_init(((macrev_Z*)new_mem)[extra].rep);
00241             }
00242             if ( type == LSUB || type == USUB ) {
00243                 copy_mem_int(len_to_copy, expanders[type].mem, new_mem);
00244             } else {
00245                 copy_mem_double<macrev_Z >(len_to_copy, 
00246                                            expanders[type].mem, new_mem);
00247                 //cout<<"'j'en detruit ici"<<endl;
00248                 for(unsigned int i=0;i<len_to_copy;i++)
00249                   mpz_clear(((macrev_Z*)expanders[type].mem)[i].rep);
00250             }
00251             
00252             SUPERLU_FREE (expanders[type].mem);
00253         }
00254         expanders[type].mem = (void *) new_mem;
00255         
00256     } else { /* MemModel == USER */
00257         if ( no_expand == 0 ) {
00258             new_mem = duser_malloc(new_len * lword, HEAD);
00259             if ( NotDoubleAlign(new_mem) &&
00260                 (type == LUSUP || type == UCOL) ) {
00261                 old_mem = new_mem;
00262                 new_mem = (void *)DoubleAlign(new_mem);
00263                 extra = (char*)new_mem - (char*)old_mem;
00264 #ifdef DEBUG            
00265                 printf("expand(): not aligned, extra %d\n", extra);
00266 #endif          
00267                 stack.top1 += extra;
00268                 stack.used += extra;
00269             }
00270             expanders[type].mem = (void *) new_mem;
00271         }
00272         else {
00273             tries = 0;
00274             extra = (new_len - *prev_len) * lword;
00275             if ( keep_prev ) {
00276                 if ( StackFull(extra) ) return (NULL);
00277             } else {
00278                 while ( StackFull(extra) ) {
00279                     if ( ++tries > 10 ) return (NULL);
00280                     alpha = Reduce(alpha);
00281                     new_len = (int)(alpha * *prev_len);
00282                     extra = (new_len - *prev_len) * lword;          
00283                 }
00284             }
00285 
00286             if ( type != USUB ) {
00287                 new_mem = (void*)((char*)expanders[type + 1].mem + extra);
00288                 bytes_to_copy = (char*)stack.array + stack.top1
00289                     - (char*)expanders[type + 1].mem;
00290                 user_bcopy((char*)expanders[type+1].mem,
00291                            (char*)new_mem, bytes_to_copy);
00292 
00293                 if ( type < USUB ) {
00294                     Glu->usub = (int*)
00295                       (expanders[USUB].mem =
00296                         (void*)((char*)expanders[USUB].mem + extra));
00297                 }
00298                 if ( type < LSUB ) {
00299                     Glu->lsub =(int*)(expanders[LSUB].mem =
00300                         (void*)((char*)expanders[LSUB].mem + extra));
00301                 }
00302                 if ( type < UCOL ) {
00303                     Glu->ucol =(macrev_Z*)( expanders[UCOL].mem =
00304                         (void*)((char*)expanders[UCOL].mem + extra));
00305                 }
00306                 stack.top1 += extra;
00307                 stack.used += extra;
00308                 if ( type == UCOL ) {
00309                     stack.top1 += extra;   /* Add same amount for USUB */
00310                     stack.used += extra;
00311                 }
00312                 
00313             } /* if ... */
00314 
00315         } /* else ... */
00316     }
00317 
00318     expanders[type].size = new_len;
00319     *prev_len = new_len;
00320     if ( no_expand ) ++no_expand;
00321     
00322     return (void *) expanders[type].mem;
00323     
00324 } /* dexpand */

void dLUWorkFree ( int  m,
int  panel_size,
int *  iwork,
macrev_Z *  dwork,
GlobalLU_t< macrev_Z > *  Glu 
)

Definition at line 174 of file dpivmacrev_z.cc.

00176 {
00177   int i;
00178   int maxsuper = sp_ienv(3),
00179     rowblk   = sp_ienv(4);
00180   if ( Glu->MemModel == SYSTEM ) {
00181         SUPERLU_FREE (iwork);
00182         for(i=0;i<m*panel_size+NUM_TEMPV(m,panel_size,maxsuper,rowblk);i++)
00183           mpz_clear((dwork[i].rep));
00184         SUPERLU_FREE (dwork);
00185     } else {
00186         stack.used -= (stack.size - stack.top2);
00187         stack.top2 = stack.size;
00188 /*      dStackCompress(Glu);  */
00189     }
00190     
00191     SUPERLU_FREE (expanders);
00192 }

int dpivotL ( const int  jcol,
const macrev_Z  u,
int *  usepr,
int *  perm_r,
int *  iperm_r,
int *  iperm_c,
int *  pivrow,
GlobalLU_t< macrev_Z > *  Glu 
) [inline]

Definition at line 6 of file dpivmacrev_z.cc.

00016 {
00017 /*
00018  * Purpose
00019  * =======
00020  *   Performs the numerical pivoting on the current column of L,
00021  *   and the CDIV operation.
00022  *
00023  *   Pivot policy:
00024  *   (1) Compute thresh = u * max_(i>=j) abs(A_ij);
00025  *   (2) IF user specifies pivot row k and abs(A_kj) >= thresh THEN
00026  *           pivot row = k;
00027  *       ELSE IF abs(A_jj) >= thresh THEN
00028  *           pivot row = j;
00029  *       ELSE
00030  *           pivot row = m;
00031  * 
00032  *   Note: If you absolutely want to use a given pivot order, then set u=0.0.
00033  *
00034  *   Return value: 0      success;
00035  *                 i > 0  U(i,i) is exactly zero.
00036  *
00037  */
00038     int          fsupc;     /* first column in the supernode */
00039     int          nsupc;     /* no of columns in the supernode */
00040     int          nsupr;     /* no of rows in the supernode */
00041     int          lptr;      /* points to the starting subscript of the supernode */
00042     int          pivptr, old_pivptr, diag, diagind;
00043     macrev_Z       pivmax, rtemp, thresh;
00044     macrev_Z       temp;
00045     macrev_Z       *lu_sup_ptr; 
00046     macrev_Z       *lu_col_ptr;
00047     int          *lsub_ptr;
00048     int          isub, icol, k, itemp;
00049     int          *lsub, *xlsub;
00050     macrev_Z       *lusup;
00051     int          *xlusup;
00052     extern SuperLUStat_t SuperLUStat;
00053     flops_t  *ops = SuperLUStat.ops;
00054 
00055     /* Initialize pointers */
00056     lsub       = Glu->lsub;
00057     xlsub      = Glu->xlsub;
00058     lusup      = Glu->lusup;
00059     xlusup     = Glu->xlusup;
00060     fsupc      = (Glu->xsup)[(Glu->supno)[jcol]];
00061     nsupc      = jcol - fsupc;          /* excluding jcol; nsupc >= 0 */
00062     lptr       = xlsub[fsupc];
00063     nsupr      = xlsub[fsupc+1] - lptr;
00064     lu_sup_ptr = &lusup[xlusup[fsupc]]; /* start of the current supernode */
00065     lu_col_ptr = &lusup[xlusup[jcol]];  /* start of jcol in the supernode */
00066     lsub_ptr   = &lsub[lptr];   /* start of row indices of the supernode */
00067 
00068 #ifdef DEBUG
00069 if ( jcol == MIN_COL ) {
00070     printf("Before cdiv: col %d\n", jcol);
00071     for (k = nsupc; k < nsupr; k++) 
00072         printf("  lu[%d] %f\n", lsub_ptr[k], lu_col_ptr[k]);
00073 }
00074 #endif
00075     
00076     /* Determine the largest abs numerical value for partial pivoting;
00077        Also search for user-specified pivot, and diagonal element. */
00078     if ( *usepr ) *pivrow = iperm_r[jcol];
00079     diagind = iperm_c[jcol];
00080     pivmax =(macrev_Z)0;
00081     pivptr = nsupc;
00082     diag = EMPTY;
00083     old_pivptr = nsupc;
00084     for (isub = nsupc; isub < nsupr; ++isub) {
00085         rtemp = lu_col_ptr[isub];
00086         
00087         if ( rtemp != 0 ) {
00088             pivmax = rtemp;
00089             pivptr = isub;
00090             break;
00091         }
00092         if ( *usepr && lsub_ptr[isub] == *pivrow ) old_pivptr = isub;
00093         if ( lsub_ptr[isub] == diagind ) diag = isub;
00094     }
00095     //cout<<"comparaison dans dpivvotL "<<pivmax<<endl;
00096     /* Test for singularity */
00097     if ( pivmax ==0) {
00098         *pivrow = lsub_ptr[nsupc];
00099         //cout<<"degenere et pivrow "<<*pivrow<<endl;
00100         perm_r[*pivrow] = jcol;
00101         *usepr = 0;
00102         return (jcol+1);
00103     }
00104 
00105     thresh = u * pivmax;
00106     
00107     /* Choose appropriate pivotal element by our policy. */
00108     if ( *usepr ) {
00109       //fabs!!!!  
00110       rtemp = lu_col_ptr[old_pivptr];
00111       //fabs!!!!!!!!!!!!!!!!! faire gaffe en generic 
00112       if ( rtemp != (macrev_Z)0 )
00113             pivptr = old_pivptr;
00114         else
00115             *usepr = 0;
00116     }
00117     if ( *usepr == 0 ) {
00118         /* Use diagonal pivot? */
00119         if ( diag >= 0 ) { /* diagonal exists */
00120             rtemp = lu_col_ptr[diag];
00121             if ( rtemp != (macrev_Z)0  ) pivptr = diag;
00122         }
00123         *pivrow = lsub_ptr[pivptr];
00124     }
00125     
00126     /* Record pivot row */
00127     //    cout<<"piv row "<<*pivrow<<" jcol "<<jcol<<endl;
00128     perm_r[*pivrow] = jcol;
00129     
00130     /* Interchange row subscripts */
00131     if ( pivptr != nsupc ) {
00132       itemp = lsub_ptr[pivptr];
00133       lsub_ptr[pivptr] = lsub_ptr[nsupc];
00134       lsub_ptr[nsupc] = itemp;
00135 
00136         /* Interchange numerical values as well, for the whole snode, such 
00137          * that L is indexed the same way as A.
00138          */
00139         for (icol = 0; icol <= nsupc; icol++) {
00140             itemp = pivptr + icol * nsupr;
00141             temp = lu_sup_ptr[itemp];
00142             lu_sup_ptr[itemp] = lu_sup_ptr[nsupc + icol*nsupr];
00143             lu_sup_ptr[nsupc + icol*nsupr] = temp;
00144         }
00145     } /* if */
00146 
00147     /* cdiv operation */
00148     ops[FACT] += nsupr - nsupc;
00149     temp = (macrev_Z)1 / lu_col_ptr[nsupc];
00150     for (k = nsupc+1; k < nsupr; k++) 
00151         lu_col_ptr[k] *= temp;
00152 
00153     return 0;
00154 }

void dSetRWork ( int  m,
int  panel_size,
macrev_Z *  dworkptr,
macrev_Z **  dense,
macrev_Z **  tempv 
)

Definition at line 157 of file dpivmacrev_z.cc.

00159 {
00160     macrev_Z zero(0);
00161     int i;
00162     int maxsuper = sp_ienv(3),
00163         rowblk   = sp_ienv(4);
00164     *dense = dworkptr;
00165     *tempv = *dense + panel_size*m;
00166     for(i=0;i<m*panel_size;i++)
00167       mpz_init((*dense)[i].rep);
00168     for(i=0;i<NUM_TEMPV(m,panel_size,maxsuper,rowblk);i++)
00169       mpz_init((*tempv)[i].rep);
00170     //dfill (*dense, m * panel_size, zero);
00171     //dfill (*tempv, NUM_TEMPV(m,panel_size,maxsuper,rowblk), zero);     
00172 }

void remiseformatNC ( SuperMatrix< macrev_Z > &  M,
const vector< macrev_Z > &  nzval,
const vector< int > &  rowind,
const vector< int > &  colptr,
int  coeff 
)

Definition at line 391 of file dpivmacrev_z.cc.

00394 {
00395   ((NCformat *)M.Store)->nnz=coeff;
00396   destroystore(M);
00397   M.Stype=NC;
00398   M.Store=(void *)malloc(sizeof(NCformat));
00399   NCformat *Mstore=(NCformat *)M.Store;
00400   Mstore->nnz=nzval.size();
00401   Mstore->nzval=(void *)malloc(sizeof(macrev_Z)*(nzval.size()));
00402   Mstore->rowind=(int*)malloc(sizeof(int)*(nzval.size()));
00403   Mstore->colptr=(int*)malloc(sizeof(int)*(colptr.size()));
00404   for(unsigned int i=0;i<nzval.size();++i)
00405     {
00406       Mstore->rowind[i]=rowind[i];
00407       mpz_init(((macrev_Z*)(Mstore->nzval))[i].rep);
00408       ((macrev_Z*)(Mstore->nzval))[i]=nzval[i];
00409     }
00410   for(unsigned int i=0;i<colptr.size();++i)
00411     Mstore->colptr[i]=colptr[i];
00412   return;
00413 }

Generated on 6 Dec 2012 for borderbasix by  doxygen 1.6.1