void dLCP::transfer_i_to_C (int i)
{
  {
    if (m_nC > 0) {
      // ell,Dell were computed by solve1(). note, ell = D \ L1solve (L,A(i,C))
      {
        const int nC = m_nC;
        dReal *const Ltgt = m_L + nC*m_nskip, *ell = m_ell;
        for (int j=0; j<nC; ++j) Ltgt[j] = ell[j];
      }
      const int nC = m_nC;
      m_d[nC] = dRecip (AROW(i)[i] - dDot(m_ell,m_Dell,nC));
    }
    else {
      m_d[0] = dRecip (AROW(i)[i]);
    }

    swapProblem (m_A,m_x,m_b,m_w,m_lo,m_hi,m_p,m_state,m_findex,m_n,m_nC,i,m_nskip,1);

    const int nC = m_nC;
    m_C[nC] = nC;
    m_nC = nC + 1; // nC value is outdated after this line
  }

# ifdef DEBUG_LCP
  checkFactorization (m_A,m_L,m_d,m_nC,m_C,m_nskip);
  if (i < (m_n-1)) checkPermutations (i+1,m_n,m_nC,m_nN,m_p,m_C);
# endif
}

void dLCP::transfer_i_from_N_to_C (int i)
{
  int j;
  if (nC > 0) {
    dReal *aptr = AROW(i);
#   ifdef NUB_OPTIMIZATIONS
    // if nub>0, initial part of aptr unpermuted
    for (j=0; j<nub; j++) Dell[j] = aptr[j];
    for (j=nub; j<nC; j++) Dell[j] = aptr[C[j]];
#   else
    for (j=0; j<nC; j++) Dell[j] = aptr[C[j]];
#   endif
    dSolveL1 (L,Dell,nC,nskip);
    for (j=0; j<nC; j++) ell[j] = Dell[j] * d[j];
    for (j=0; j<nC; j++) L[nC*nskip+j] = ell[j];
    d[nC] = dRecip (AROW(i)[i] - dDot(ell,Dell,nC));
  }
  else {
    d[0] = dRecip (AROW(i)[i]);
  }
  swapProblem (A,x,b,w,lo,hi,p,state,findex,n,nC,i,nskip,1);
  C[nC] = nC;
  nN--;
  nC++;

  // @@@ TO DO LATER
  // if we just finish here then we'll go back and re-solve for
  // delta_x. but actually we can be more efficient and incrementally
  // update delta_x here. but if we do this, we wont have ell and Dell
  // to use in updating the factorization later.

# ifdef DEBUG_LCP
  checkFactorization (A,L,d,nC,C,nskip);
# endif
}

void dLCP::transfer_i_from_N_to_C (int i)
{
  {
    if (m_nC > 0) {
      {
        dReal *const aptr = AROW(i);
        dReal *Dell = m_Dell;
        const int *C = m_C;
#   ifdef NUB_OPTIMIZATIONS
        // if nub>0, initial part of aptr unpermuted
        const int nub = m_nub;
        int j=0;
        for ( ; j<nub; ++j) Dell[j] = aptr[j];
        const int nC = m_nC;
        for ( ; j<nC; ++j) Dell[j] = aptr[C[j]];
#   else
        const int nC = m_nC;
        for (int j=0; j<nC; ++j) Dell[j] = aptr[C[j]];
#   endif
      }
      dSolveL1 (m_L,m_Dell,m_nC,m_nskip);
      {
        const int nC = m_nC;
        dReal *const Ltgt = m_L + nC*m_nskip;
        dReal *ell = m_ell, *Dell = m_Dell, *d = m_d;
        for (int j=0; j<nC; ++j) Ltgt[j] = ell[j] = Dell[j] * d[j];
      }
      const int nC = m_nC;
      dReal Aii_dDot = AROW(i)[i] - dDot(m_ell, m_Dell, nC);
      if(dFabs(Aii_dDot) < 1e-16) {
          Aii_dDot += 1e-6;
      }
      m_d[nC] = dRecip (Aii_dDot);
    }
    else {
        if(dFabs(AROW(i)[i]) < 1e-16) {
            AROW(i)[i] += 1e-6;
        }
        m_d[0] = dRecip (AROW(i)[i]);
    }

    swapProblem (m_A,m_x,m_b,m_w,m_lo,m_hi,m_p,m_state,m_findex,m_n,m_nC,i,m_nskip,1);

    const int nC = m_nC;
    m_C[nC] = nC;
    m_nN--;
    m_nC = nC + 1; // nC value is outdated after this line
  }

  // @@@ TO DO LATER
  // if we just finish here then we'll go back and re-solve for
  // delta_x. but actually we can be more efficient and incrementally
  // update delta_x here. but if we do this, we wont have ell and Dell
  // to use in updating the factorization later.

# ifdef DEBUG_LCP
  checkFactorization (m_A,m_L,m_d,m_nC,m_C,m_nskip);
# endif
}

void dLCP::pN_plusequals_ANi (dReal *p, int i, int sign)
{
  int k;
  for (k=0; k<n; k++) if (N[k] && k >= i) dDebug (0,"N assumption violated");
  if (sign > 0) {
    for (k=0; k<n; k++) if (N[k]) p[k] += AROW(i)[k];
  }
  else {
    for (k=0; k<n; k++) if (N[k]) p[k] -= AROW(i)[k];
  }
}

dLCP::dLCP (int _n, int _nub, dReal *_Adata, dReal *_x, dReal *_b, dReal *_w,
	    dReal *_lo, dReal *_hi, dReal *_L, dReal *_d,
	    dReal *_Dell, dReal *_ell, dReal *_tmp,
	    int *_state, int *_findex, int *_p, int *_C, dReal **Arows)
{
  dUASSERT (_findex==0,"slow dLCP object does not support findex array");

  n = _n;
  nub = _nub;
  Adata = _Adata;
  A = 0;
  x = _x;
  b = _b;
  w = _w;
  lo = _lo;
  hi = _hi;
  nskip = dPAD(n);
  dSetZero (x,n);
  last_i_for_solve1 = -1;

  int i,j;
  C.setSize (n);
  N.setSize (n);
  for (int i=0; i<n; i++) {
    C[i] = 0;
    N[i] = 0;
  }

# ifdef ROWPTRS
  // make matrix row pointers
  A = Arows;
  for (i=0; i<n; i++) A[i] = Adata + i*nskip;
# else
  A = Adata;
# endif

  // lets make A symmetric
  for (i=0; i<n; i++) {
    for (j=i+1; j<n; j++) AROW(i)[j] = AROW(j)[i];
  }

  // if nub>0, put all indexes 0..nub-1 into C and solve for x
  if (nub > 0) {
    for (i=0; i<nub; i++) memcpy (_L+i*nskip,AROW(i),(i+1)*sizeof(dReal));
    dFactorLDLT (_L,_d,nub,nskip);
    memcpy (x,b,nub*sizeof(dReal));
    dSolveLDLT (_L,_d,x,nub,nskip);
    dSetZero (_w,nub);
    for (i=0; i<nub; i++) C[i] = 1;
  }
}

static void checkFactorization (ATYPE A, dReal *_L, dReal *_d,
                                int nC, int *C, int nskip)
{
  int i,j;
  if (nC==0) return;

  // get A1=A, copy the lower triangle to the upper triangle, get A2=A[C,C]
  dMatrix A1 (nC,nC);
  for (i=0; i<nC; i++) {
    for (j=0; j<=i; j++) A1(i,j) = A1(j,i) = AROW(i)[j];
  }
  dMatrix A2 = A1.select (nC,C,nC,C);

  // printf ("A1=\n"); A1.print(); printf ("\n");
  // printf ("A2=\n"); A2.print(); printf ("\n");

  // compute A3 = L*D*L'
  dMatrix L (nC,nC,_L,nskip,1);
  dMatrix D (nC,nC);
  for (i=0; i<nC; i++) D(i,i) = 1/_d[i];
  L.clearUpperTriangle();
  for (i=0; i<nC; i++) L(i,i) = 1;
  dMatrix A3 = L * D * L.transpose();

  // printf ("L=\n"); L.print(); printf ("\n");
  // printf ("D=\n"); D.print(); printf ("\n");
  // printf ("A3=\n"); A2.print(); printf ("\n");

  // compare A2 and A3
  dReal diff = A2.maxDifference (A3);
  if (diff > 1e-8)
    dDebug (0,"L*D*L' check, maximum difference = %.6e\n",diff);
}

void dLCP::pN_equals_ANC_times_qC (dReal *p, dReal *q)
{
  dReal sum;
  for (int ii=0; ii<n; ii++) if (N[ii]) {
    sum = 0;
    for (int jj=0; jj<n; jj++) if (C[jj]) sum += AROW(ii)[jj] * q[jj];
    p[ii] = sum;
  }
}

void dLCP::transfer_i_to_C (int i)
{
  int j;
  if (nC > 0) {
    // ell,Dell were computed by solve1(). note, ell = D \ L1solve (L,A(i,C))
    for (j=0; j<nC; j++) L[nC*nskip+j] = ell[j];
    d[nC] = dRecip (AROW(i)[i] - dDot(ell,Dell,nC));
  }
  else {
    d[0] = dRecip (AROW(i)[i]);
  }
  swapProblem (A,x,b,w,lo,hi,p,state,findex,n,nC,i,nskip,1);
  C[nC] = nC;
  nC++;

# ifdef DEBUG_LCP
  checkFactorization (A,L,d,nC,C,nskip);
  if (i < (n-1)) checkPermutations (i+1,n,nC,nN,p,C);
# endif
}

void dLCP::solve1 (dReal *a, int i, int dir, int only_transfer)
{
  // the `Dell' and `ell' that are computed here are saved. if index i is
  // later added to the factorization then they can be reused.
  //
  // @@@ question: do we need to solve for entire delta_x??? yes, but
  //     only if an x goes below 0 during the step.

  if (m_nC > 0) {
    {
      dReal *Dell = m_Dell;
      int *C = m_C;
      dReal *aptr = AROW(i);
#   ifdef NUB_OPTIMIZATIONS
      // if nub>0, initial part of aptr[] is guaranteed unpermuted
      const int nub = m_nub;
      int j=0;
      for ( ; j<nub; ++j) Dell[j] = aptr[j];
      const int nC = m_nC;
      for ( ; j<nC; ++j) Dell[j] = aptr[C[j]];
#   else
      const int nC = m_nC;
      for (int j=0; j<nC; ++j) Dell[j] = aptr[C[j]];
#   endif
    }
    dSolveL1 (m_L,m_Dell,m_nC,m_nskip);
    {
      dReal *ell = m_ell, *Dell = m_Dell, *d = m_d;
      const int nC = m_nC;
      for (int j=0; j<nC; ++j) ell[j] = Dell[j] * d[j];
    }

    if (!only_transfer) {
      dReal *tmp = m_tmp, *ell = m_ell;
      {
        const int nC = m_nC;
        for (int j=0; j<nC; ++j) tmp[j] = ell[j];
      }
      dSolveL1T (m_L,tmp,m_nC,m_nskip);
      if (dir > 0) {
        int *C = m_C;
        dReal *tmp = m_tmp;
        const int nC = m_nC;
        for (int j=0; j<nC; ++j) a[C[j]] = -tmp[j];
      } else {
        int *C = m_C;
        dReal *tmp = m_tmp;
        const int nC = m_nC;
        for (int j=0; j<nC; ++j) a[C[j]] = tmp[j];
      }
    }
  }
}

 dReal AiN_times_qN (int i, dReal *q) { return dDot (AROW(i)+nC,q+nC,nN); }

void _fastcall AMemLeaks(ParamBlk *parm)
{
	char *pArrayName;
	Locator lArrayLoc;
	int nErrorNo;
	INTEGER(vMem);
	STRING(vMemInfo);
	char *pMemInfo;
	LPDBGALLOCINFO pDbg = gpDbgInfo;

	if (!pDbg)
	{
		RET_INTEGER(0);
		return;
	}

	if (!NULLTERMINATE(p1))
		RAISEERROR(E_INSUFMEMORY);

	LOCKHAND(p1);	
	pArrayName = HANDTOPTR(p1);

	if (!ALLOCHAND(vMemInfo,VFP2C_ERROR_MESSAGE_LEN))
	{
		nErrorNo = E_INSUFMEMORY;
		goto ErrorOut;
	}
	LOCKHAND(vMemInfo);
	pMemInfo = HANDTOPTR(vMemInfo);

	if (nErrorNo = DimensionEx(pArrayName,&lArrayLoc,1,4))
		goto ErrorOut;

	while (pDbg)
	{
		if (nErrorNo = Dimension(pArrayName,++AROW(lArrayLoc),4))
			goto ErrorOut;

		vMem.ev_long = (int)pDbg->pPointer;
		ADIM(lArrayLoc) = 1;
		if (nErrorNo = STORE(lArrayLoc,vMem))
			goto ErrorOut;

		vMem.ev_long = pDbg->nSize;
        ADIM(lArrayLoc) = 2;
		if (nErrorNo = STORE(lArrayLoc,vMem))
			goto ErrorOut;

		if (pDbg->pProgInfo)
			vMemInfo.ev_length = strncpyex(pMemInfo,pDbg->pProgInfo,VFP2C_ERROR_MESSAGE_LEN);
		else
			vMemInfo.ev_length = 0;
		ADIM(lArrayLoc) = 3;
		if (nErrorNo = STORE(lArrayLoc,vMemInfo))
			goto ErrorOut;

		vMemInfo.ev_length = min(pDbg->nSize,VFP2C_ERROR_MESSAGE_LEN);
		memcpy(pMemInfo,pDbg->pPointer,vMemInfo.ev_length);
		ADIM(lArrayLoc) = 4;
		if (nErrorNo = STORE(lArrayLoc,vMemInfo))
			goto ErrorOut;

		pDbg = pDbg->next;
	}

	UNLOCKHAND(p1);
	UNLOCKHAND(vMemInfo);
	FREEHAND(vMemInfo);
	RET_AROWS(lArrayLoc);
	return;
	
	ErrorOut:
		UNLOCKHAND(p1);
		if (VALIDHAND(vMemInfo))
		{
			UNLOCKHAND(vMemInfo);
			FREEHAND(vMemInfo);
		}
		RAISEERROR(nErrorNo);
}

 dReal AiN_times_qN (int i, dReal *q) const { return dDot (AROW(i)+m_nC, q+m_nC, m_nN); }

dLCP::dLCP (int _n, int _nskip, int _nub, dReal *_Adata, dReal *_x, dReal *_b, dReal *_w,
            dReal *_lo, dReal *_hi, dReal *_L, dReal *_d,
            dReal *_Dell, dReal *_ell, dReal *_tmp,
            bool *_state, int *_findex, int *_p, int *_C, dReal **Arows):
  m_n(_n), m_nskip(_nskip), m_nub(_nub), m_nC(0), m_nN(0),
# ifdef ROWPTRS
  m_A(Arows),
#else
  m_A(_Adata),
#endif
  m_x(_x), m_b(_b), m_w(_w), m_lo(_lo), m_hi(_hi),
  m_L(_L), m_d(_d), m_Dell(_Dell), m_ell(_ell), m_tmp(_tmp),
  m_state(_state), m_findex(_findex), m_p(_p), m_C(_C)
{
  {
    dSetZero (m_x,m_n);
  }

  {
# ifdef ROWPTRS
    // make matrix row pointers
    dReal *aptr = _Adata;
    ATYPE A = m_A;
    const int n = m_n, nskip = m_nskip;
    for (int k=0; k<n; aptr+=nskip, ++k) A[k] = aptr;
# endif
  }

  {
    int *p = m_p;
    const int n = m_n;
    for (int k=0; k<n; ++k) p[k]=k;		// initially unpermuted
  }

  /*
  // for testing, we can do some random swaps in the area i > nub
  {
    const int n = m_n;
    const int nub = m_nub;
    if (nub < n) {
    for (int k=0; k<100; k++) {
      int i1,i2;
      do {
        i1 = dRandInt(n-nub)+nub;
        i2 = dRandInt(n-nub)+nub;
      }
      while (i1 > i2); 
      //printf ("--> %d %d\n",i1,i2);
      swapProblem (m_A,m_x,m_b,m_w,m_lo,m_hi,m_p,m_state,m_findex,n,i1,i2,m_nskip,0);
    }
  }
  */

  // permute the problem so that *all* the unbounded variables are at the
  // start, i.e. look for unbounded variables not included in `nub'. we can
  // potentially push up `nub' this way and get a bigger initial factorization.
  // note that when we swap rows/cols here we must not just swap row pointers,
  // as the initial factorization relies on the data being all in one chunk.
  // variables that have findex >= 0 are *not* considered to be unbounded even
  // if lo=-inf and hi=inf - this is because these limits may change during the
  // solution process.

  {
    int *findex = m_findex;
    dReal *lo = m_lo, *hi = m_hi;
    const int n = m_n;
    for (int k = m_nub; k<n; ++k) {
      if (findex && findex[k] >= 0) continue;
      if (lo[k]==-dInfinity && hi[k]==dInfinity) {
        swapProblem (m_A,m_x,m_b,m_w,lo,hi,m_p,m_state,findex,n,m_nub,k,m_nskip,0);
        m_nub++;
      }
    }
  }

  // if there are unbounded variables at the start, factorize A up to that
  // point and solve for x. this puts all indexes 0..nub-1 into C.
  if (m_nub > 0) {
    const int nub = m_nub;
    {
      dReal *Lrow = m_L;
      const int nskip = m_nskip;
      for (int j=0; j<nub; Lrow+=nskip, ++j) memcpy(Lrow,AROW(j),(j+1)*sizeof(dReal));
    }
    dFactorLDLT (m_L,m_d,nub,m_nskip);
    memcpy (m_x,m_b,nub*sizeof(dReal));
    dSolveLDLT (m_L,m_d,m_x,nub,m_nskip);
    dSetZero (m_w,nub);
    {
      int *C = m_C;
      for (int k=0; k<nub; ++k) C[k] = k;
    }
    m_nC = nub;
  }

  // permute the indexes > nub such that all findex variables are at the end
  if (m_findex) {
    const int nub = m_nub;
    int *findex = m_findex;
    int num_at_end = 0;
//.........这里部分代码省略.........

 dReal Aii (int i) const  { return AROW(i)[i]; }

 dReal AiC_times_qC (int i, dReal *q) const { return dDot (AROW(i), q, m_nC); }

void _fastcall AMemBlocks(ParamBlk *parm)
{
	char *pArrayName;
	Locator lArrayLoc;
	int nErrorNo;
	PROCESS_HEAP_ENTRY pEntry;
	INTEGER(vAddress);
	INTEGER(vSize);
	INTEGER(vOverhead);
	DWORD nLastError;

	if (!fpHeapWalk)
		RAISEERROR(E_NOENTRYPOINT);

	if (!NULLTERMINATE(p1))
		RAISEERROR(E_INSUFMEMORY);

	LOCKHAND(p1);	
	pArrayName = HANDTOPTR(p1);

	if (nErrorNo = DimensionEx(pArrayName,&lArrayLoc,1,3))
		goto ErrorOut;

	pEntry.lpData = NULL;

	if (!fpHeapWalk(ghHeap,&pEntry))
	{
		nLastError = GetLastError();
		UNLOCKHAND(p1);
		if (nLastError == ERROR_NO_MORE_ITEMS)
		{
			RET_INTEGER(0);
			return;
		}
		else
		{
			SAVEWIN32ERROR(HeapWalk,nLastError);
			RET_INTEGER(-1);
		}
	}

	do 
	{
		AROW(lArrayLoc)++;

		if ((nErrorNo = Dimension(pArrayName,AROW(lArrayLoc),3)))
			break;

		ADIM(lArrayLoc) = 1;
		vAddress.ev_long = (long)pEntry.lpData;
		if (nErrorNo = STORE(lArrayLoc,vAddress))
			break;

		ADIM(lArrayLoc) = 2;
		vSize.ev_long = pEntry.cbData;
		if (nErrorNo = STORE(lArrayLoc,vSize))
			break;

		ADIM(lArrayLoc) = 3;
		vOverhead.ev_long = pEntry.cbOverhead;
		if (nErrorNo = STORE(lArrayLoc,vOverhead))
			break;

	} while (fpHeapWalk(ghHeap,&pEntry));
	
	nLastError = GetLastError();

	if (nErrorNo)
		goto ErrorOut;

	UNLOCKHAND(p1);

    if (nLastError == ERROR_NO_MORE_ITEMS)
	{
		RET_AROWS(lArrayLoc);
		return;
	}
	else
	{
		SAVEWIN32ERROR(HeapWalk,nLastError);
		RET_INTEGER(-1);
		return;
	}

	ErrorOut:
		UNLOCKHAND(p1);
		RAISEERROR(nErrorNo);
}

void dLCP::solve1 (dReal *a, int i, int dir, int only_transfer)
{

  ALLOCA (dReal,AA,n*nskip*sizeof(dReal));
#ifdef dUSE_MALLOC_FOR_ALLOCA
    if (AA == NULL) {
      dMemoryFlag = d_MEMORY_OUT_OF_MEMORY;
      return;
    }
#endif
  ALLOCA (dReal,dd,n*sizeof(dReal));
#ifdef dUSE_MALLOC_FOR_ALLOCA
    if (dd == NULL) {
      UNALLOCA(AA);
      dMemoryFlag = d_MEMORY_OUT_OF_MEMORY;
      return;
    }
#endif
  ALLOCA (dReal,bb,n*sizeof(dReal));
#ifdef dUSE_MALLOC_FOR_ALLOCA
    if (bb == NULL) {
      UNALLOCA(AA);
      UNALLOCA(dd);
      dMemoryFlag = d_MEMORY_OUT_OF_MEMORY;
      return;
    }
#endif

  int ii,jj,AAi,AAj;

  last_i_for_solve1 = i;
  AAi = 0;
  for (ii=0; ii<n; ii++) if (C[ii]) {
    AAj = 0;
    for (jj=0; jj<n; jj++) if (C[jj]) {
      AA[AAi*nskip+AAj] = AROW(ii)[jj];
      AAj++;
    }
    bb[AAi] = AROW(i)[ii];
    AAi++;
  }
  if (AAi==0) {
      UNALLOCA (AA);
      UNALLOCA (dd);
      UNALLOCA (bb);
      return;
  }

  dFactorLDLT (AA,dd,AAi,nskip);
  dSolveLDLT (AA,dd,bb,AAi,nskip);

  AAi=0;
  if (dir > 0) {
    for (ii=0; ii<n; ii++) if (C[ii]) a[ii] = -bb[AAi++];
  }
  else {
    for (ii=0; ii<n; ii++) if (C[ii]) a[ii] = bb[AAi++];
  }

  UNALLOCA (AA);
  UNALLOCA (dd);
  UNALLOCA (bb);
}

 dReal Aii (int i) { return AROW(i)[i]; }

dLCP::dLCP (int _n, int _nub, dReal *_Adata, dReal *_x, dReal *_b, dReal *_w,
	    dReal *_lo, dReal *_hi, dReal *_L, dReal *_d,
	    dReal *_Dell, dReal *_ell, dReal *_tmp,
	    int *_state, int *_findex, int *_p, int *_C, dReal **Arows)
{
  n = _n;
  nub = _nub;
  Adata = _Adata;
  A = 0;
  x = _x;
  b = _b;
  w = _w;
  lo = _lo;
  hi = _hi;
  L = _L;
  d = _d;
  Dell = _Dell;
  ell = _ell;
  tmp = _tmp;
  state = _state;
  findex = _findex;
  p = _p;
  C = _C;
  nskip = dPAD(n);
  dSetZero (x,n);

  int k;

# ifdef ROWPTRS
  // make matrix row pointers
  A = Arows;
  for (k=0; k<n; k++) A[k] = Adata + k*nskip;
# else
  A = Adata;
# endif

  nC = 0;
  nN = 0;
  for (k=0; k<n; k++) p[k]=k;		// initially unpermuted

  /*
  // for testing, we can do some random swaps in the area i > nub
  if (nub < n) {
    for (k=0; k<100; k++) {
      int i1,i2;
      do {
	i1 = dRandInt(n-nub)+nub;
	i2 = dRandInt(n-nub)+nub;
      }
      while (i1 > i2); 
      //printf ("--> %d %d\n",i1,i2);
      swapProblem (A,x,b,w,lo,hi,p,state,findex,n,i1,i2,nskip,0);
    }
  }
  */

  // permute the problem so that *all* the unbounded variables are at the
  // start, i.e. look for unbounded variables not included in `nub'. we can
  // potentially push up `nub' this way and get a bigger initial factorization.
  // note that when we swap rows/cols here we must not just swap row pointers,
  // as the initial factorization relies on the data being all in one chunk.
  // variables that have findex >= 0 are *not* considered to be unbounded even
  // if lo=-inf and hi=inf - this is because these limits may change during the
  // solution process.

  for (k=nub; k<n; k++) {
    if (findex && findex[k] >= 0) continue;
    if (lo[k]==-dInfinity && hi[k]==dInfinity) {
      swapProblem (A,x,b,w,lo,hi,p,state,findex,n,nub,k,nskip,0);
      nub++;
    }
  }

  // if there are unbounded variables at the start, factorize A up to that
  // point and solve for x. this puts all indexes 0..nub-1 into C.
  if (nub > 0) {
    for (k=0; k<nub; k++) memcpy (L+k*nskip,AROW(k),(k+1)*sizeof(dReal));
    dFactorLDLT (L,d,nub,nskip);
    memcpy (x,b,nub*sizeof(dReal));
    dSolveLDLT (L,d,x,nub,nskip);
    dSetZero (w,nub);
    for (k=0; k<nub; k++) C[k] = k;
    nC = nub;
  }

  // permute the indexes > nub such that all findex variables are at the end
  if (findex) {
    int num_at_end = 0;
    for (k=n-1; k >= nub; k--) {
      if (findex[k] >= 0) {
	swapProblem (A,x,b,w,lo,hi,p,state,findex,n,k,n-1-num_at_end,nskip,1);
	num_at_end++;
      }
    }
  }

  // print info about indexes
  /*
  for (k=0; k<n; k++) {
    if (k<nub) printf ("C");
//.........这里部分代码省略.........

 dReal AiC_times_qC (int i, dReal *q) { return dDot (AROW(i),q,nC); }