自学教程：C++ AROW函数代码示例

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); }