psymbfact_util.c

SuperLU 2.2版本。对大型、稀疏、非对称的线性系统的直接求解

  }
  else {
    nelts = xsub[vtxXp_lid+1] - xsub[vtxXp_lid];
    if (nelts == 0) nelts = 1;
    len_texp = xsub[fstVtx_nextLvl_lid] - xsub[vtxXp_lid+1] +
      4 * nelts;
    prev_len = xsub[vtxXp_lid];
  }
  
  if (prev_len + len_texp >= prev_xsub_nextLvl) {
    /* not enough memory */
    min_new_len = prev_len + len_texp + (sz_prev_mem - prev_xsub_nextLvl);
    if (mem_error = 
	psymbfact_LUXpandMem (iam, n, vtxXp, next, min_new_len, 
			      mem_type, rout_type, 0, Pslu_freeable, Llu_symbfact,
			      VInfo, PS))
      return (mem_error);
    if ( mem_type == LSUB ) 
      new_mem = Llu_symbfact->lsub;
    else if ( mem_type == USUB ) 
      new_mem = Llu_symbfact->usub;
  }
  else 
    new_mem = prev_mem;

  if (mem_type == LSUB && PS->estimLSz < (prev_len + len_texp))
    PS->estimLSz = prev_len + len_texp;
  if (mem_type == USUB && PS->estimUSz < (prev_len + len_texp))
    PS->estimUSz = prev_len;

  /* expand the space */
  if (rout_type == LL_SYMB) {
    i = xsub[vtxXp_lid] + len_texp;
    vtx_lid = fstVtx_nextLvl_lid - 1;
    for (; vtx_lid > vtxXp_lid; vtx_lid --) {
      j = xsub[vtx_lid];  
      nel = 0;
      while (j < xsub[vtx_lid+1] && prev_mem[j] != EMPTY) {
	nel ++; j ++;
      }
      j = xsub[vtx_lid] + nel - 1;  
      k = i - (xsub[vtx_lid+1] - xsub[vtx_lid]) + nel - 1;
      if (k+1 < i)  new_mem[k+1] = EMPTY; 
      while (j >= xsub[vtx_lid]) {
	new_mem[k] = prev_mem[j]; k--; j--;
      }
      k = i;
      i -= (xsub[vtx_lid+1] - xsub[vtx_lid]);
      xsub[vtx_lid+1] = k;
    }
    xsub[vtx_lid+1] = i;
    k = *p_next;
    if (k < xsub[vtx_lid+1])
      new_mem[k] = EMPTY;
  }

  if (rout_type == RL_SYMB) {
    *p_next -= xsub[vtxXp_lid];
    i = xsub[fstVtxLvl_loc_lid] + len_texp;
    vtx_lid = fstVtx_nextLvl_lid - 1;
    for (; vtx_lid >= fstVtxLvl_loc_lid; vtx_lid --) {
      nelts = 2 * (xsub[vtx_lid+1] - xsub[vtx_lid]);
      if (nelts == 0) nelts = 2;
      if (nelts > Llu_symbfact->cntelt_vtcs[vtx_lid])
	nelts = Llu_symbfact->cntelt_vtcs[vtx_lid];
      j = xsub[vtx_lid];  
      nel = 0;
      while (j < xsub[vtx_lid+1] && prev_mem[j] != EMPTY) {
	nel ++; j ++;
      }
      j = xsub[vtx_lid] + nel - 1;  
      k = i - nelts + nel - 1;
      if (k+1 < i) new_mem[k+1] = EMPTY; 
      while (j >= xsub[vtx_lid]) {
	new_mem[k] = prev_mem[j]; k--; j--;
      }
      k = i;
      i -= nelts;
      xsub[vtx_lid+1] = k;
    }
    *p_next += xsub[vtxXp_lid];
  }  

  if (free_prev_mem && new_mem != prev_mem)
    SUPERLU_FREE (prev_mem);
  Llu_symbfact->no_expcp ++;
  
  return SUCCES_RET;
}

/*
 * Expand the data structures for L and U during the factorization.
 * Return value:   0 - successful return
 *               > 0 - number of bytes allocated when run out of space
 */
/************************************************************************/
int_t psymbfact_LUXpand_RL
/************************************************************************/
(
 int_t iam, 
 int_t n,           /* total number of columns */
 int_t vtxXp,       /* current vertex */
 int_t next,        /* number of elements currently in the factors */
 int_t len_texp,    /* length to expand */
 int_t mem_type,    /* which type of memory to expand  */
 Pslu_freeable_t *Pslu_freeable, 
 Llu_symbfact_t *Llu_symbfact,  /* modified - global LU data structures */
 vtcsInfo_symbfact_t *VInfo,
 psymbfact_stat_t *PS
 )
{
  int_t  *new_mem, *prev_mem, *xsub, mem_error, sz_prev_mem;
  /* size of the memory to be copied to new store starting from the 
     beginning/end of the memory */
  int_t exp, prev_xsub_nextLvl, vtxXp_lid, xsub_nextLvl;
  int_t *globToLoc, nvtcs_loc, maxNvtcsPProc;
  int_t fstVtx_nextLvl, fstVtx_nextLvl_lid;
  int_t i, j, k, vtx_lid, nel;
  int_t fstVtxLvl_loc_lid, prev_len, min_new_len;

#ifdef TEST_SYMB
  printf ("Pe[%d] Expand LU_RL mem_t %d vtxXp %d\n", 
	  iam, mem_type, vtxXp); 
#endif
  globToLoc = Pslu_freeable->globToLoc;
  nvtcs_loc = VInfo->nvtcs_loc;
  maxNvtcsPProc  = Pslu_freeable->maxNvtcsPProc;
  fstVtx_nextLvl = VInfo->fstVtx_nextLvl;
  
  vtxXp_lid = LOCAL_IND( globToLoc[vtxXp] );
  if (fstVtx_nextLvl == n)
    fstVtx_nextLvl_lid = VInfo->nvtcs_loc;
  else
    fstVtx_nextLvl_lid = LOCAL_IND( globToLoc[fstVtx_nextLvl] );  

  if ( mem_type == LSUB ) {
    xsub = Llu_symbfact->xlsub;
    prev_mem = Llu_symbfact->lsub;
    prev_xsub_nextLvl = VInfo->xlsub_nextLvl;
    sz_prev_mem = Llu_symbfact->szLsub;
  } else if ( mem_type == USUB ) {
    xsub = Llu_symbfact->xusub;
    prev_mem = Llu_symbfact->usub;
    prev_xsub_nextLvl = VInfo->xusub_nextLvl;
    sz_prev_mem = Llu_symbfact->szUsub;
  }
  else ABORT("Tries to expand nonexisting memory type.\n");
  
  /* Try to expand the size of xsub in the existing memory */
  prev_len = xsub[vtxXp_lid];
  
  if (prev_len + len_texp >= prev_xsub_nextLvl) {
    /* not enough memory */
    min_new_len = prev_len + len_texp + (sz_prev_mem - prev_xsub_nextLvl);
    if (mem_error = 
	psymbfact_LUXpandMem (iam, n, vtxXp, next, min_new_len, 
			      mem_type, RL_SYMB, 0, Pslu_freeable, Llu_symbfact,
			      VInfo, PS))
      return (mem_error);
    if ( mem_type == LSUB ) 
      new_mem = Llu_symbfact->lsub;
    else if ( mem_type == USUB ) 
      new_mem = Llu_symbfact->usub;
  }
  else 
    new_mem = prev_mem;

  /* expand the space */
  if (mem_type == LSUB && PS->estimLSz < (prev_len + len_texp))
    PS->estimLSz = prev_len + len_texp;
  if (mem_type == USUB && PS->estimUSz < (prev_len + len_texp))
    PS->estimUSz = prev_len;

  i = xsub[vtxXp_lid] + len_texp;
  vtx_lid = fstVtx_nextLvl_lid - 1;
  for (; vtx_lid > vtxXp_lid; vtx_lid --) {
    j = xsub[vtx_lid];  
    nel = 0;
    while (j < xsub[vtx_lid+1] && prev_mem[j] != EMPTY) {
      nel ++; j++;
    }
    j = xsub[vtx_lid] + nel - 1;  
    k = i - Llu_symbfact->cntelt_vtcs[vtx_lid] + nel - 1;
    if (k+1 < i) 
      new_mem[k+1] = EMPTY; 
    while (j >= xsub[vtx_lid]) {
      new_mem[k] = prev_mem[j];
      k--; j--;
    }
    k = i;
    i -= Llu_symbfact->cntelt_vtcs[vtx_lid];
    xsub[vtx_lid+1] = k;
  }
  xsub[vtx_lid+1] = i;
  k = next;
  if (k < xsub[vtx_lid+1])
    new_mem[k] = EMPTY;
  
  if (new_mem != prev_mem)
    SUPERLU_FREE (prev_mem);
  Llu_symbfact->no_expcp ++;
  
  return SUCCES_RET;
}

/*
 * Expand the data structures for L and U pruned during the factorization.
 * Return value: SUCCES_RET - successful return
 *               ERROR_RET - error when run out of space
 */
/************************************************************************/
int_t psymbfact_prLUXpand
/************************************************************************/
(
 int_t iam, 
 int_t min_new_len, /* minimum new length to allocate */ 
 MemType mem_type,  /* which type of memory to expand  */
 Llu_symbfact_t *Llu_symbfact, /* modified L/U pruned structures */
 psymbfact_stat_t *PS
 )
{
  int_t *prev_mem, *new_mem;
  int_t prev_len, new_len, len_tcopy_fbeg;
  
  if ( mem_type == LSUB_PR ) {
    prev_len = Llu_symbfact->szLsubPr;
    prev_mem = Llu_symbfact->lsubPr;
    len_tcopy_fbeg = Llu_symbfact->indLsubPr;
  } else if ( mem_type == USUB_PR ) {
    prev_len = Llu_symbfact->szUsubPr;
    prev_mem = Llu_symbfact->usubPr;
    len_tcopy_fbeg = Llu_symbfact->indUsubPr;
  } else ABORT("Tries to expand nonexisting memory type.\n");
  
#ifdef TEST_SYMB
  printf ("Pe[%d] Expand prmem prev_len %d min_new_l %d len_tfbeg %d\n", 
	  iam, prev_len, min_new_len, len_tcopy_fbeg);
#endif
  
  new_mem = expand (prev_len, min_new_len, prev_mem, 
		    &new_len, len_tcopy_fbeg, 0, PS);
  
  if ( !new_mem ) {
    fprintf(stderr, "Can't expand MemType %d: \n", mem_type);
    return (ERROR_RET);
  }
  
  Llu_symbfact->no_expand_pr ++;
  if ( mem_type == LSUB_PR ) {
    Llu_symbfact->lsubPr  = new_mem;
    Llu_symbfact->szLsubPr = new_len;
  } else if ( mem_type == USUB_PR ) {
    Llu_symbfact->usubPr  = new_mem;
    Llu_symbfact->szUsubPr = new_len;
  } else ABORT("Tries to expand nonexisting memory type.\n");
  
  SUPERLU_FREE (prev_mem);

  return SUCCES_RET;
}