122 #define btLCP_FAST // use fast btLCP object
126 #define BTATYPE btScalar **
127 #define BTAROW(i) (m_A[i])
134 #define BTNUB_OPTIMIZATIONS
151 btScalar Z11,m11,Z21,m21,p1,q1,p2,*ex;
155 for (i=0; i < n; i+=2) {
163 for (j=i-2; j >= 0; j -= 2) {
204 Z21 = ex[1] - Z21 - p1*Z11;
223 btScalar Z11,m11,Z12,m12,Z21,m21,Z22,m22,p1,q1,p2,q2,*ex;
227 for (i=0; i < n; i+=2) {
237 for (j=i-2; j >= 0; j -= 2) {
292 Z12 = ex[lskip1] - Z12;
295 Z21 = ex[1] - Z21 - p1*Z11;
297 Z22 = ex[1+lskip1] - Z22 - p1*Z12;
307 btScalar sum,*ell,*dee,dd,p1,p2,q1,q2,Z11,m11,Z21,m21,Z22,m22;
310 for (i=0; i<=n-2; i += 2) {
320 for (j=i-6; j >= 0; j -= 6) {
423 Z21 = ell[nskip1] - Z21;
424 Z22 = ell[1+nskip1] - Z22;
451 for (j=i-6; j >= 0; j -= 6) {
529 btScalar Z11,Z21,Z31,Z41,p1,q1,p2,p3,p4,*ex;
531 int lskip2,lskip3,i,j;
536 for (i=0; i <= n-4; i+=4) {
546 for (j=i-12; j >= 0; j -= 12) {
706 Z21 = ex[1] - Z21 - p1*Z11;
710 Z31 = ex[2] - Z31 - p1*Z11 - p2*Z21;
715 Z41 = ex[3] - Z41 - p1*Z11 - p2*Z21 - p3*Z31;
727 for (j=i-12; j >= 0; j -= 12) {
822 btScalar Z11,m11,Z21,m21,Z31,m31,Z41,m41,p1,q1,p2,p3,p4,*ex;
824 int lskip2,lskip3,i,j;
826 L = L + (n-1)*(lskip1+1);
833 for (i=0; i <= n-4; i+=4) {
843 for (j=i-4; j >= 0; j -= 4) {
936 Z21 = ex[-1] - Z21 - p1*Z11;
940 Z31 = ex[-2] - Z31 - p1*Z11 - p2*Z21;
945 Z41 = ex[-3] - Z41 - p1*Z11 - p2*Z21 - p3*Z31;
957 for (j=i-4; j >= 0; j -= 4) {
1012 for (
int i=0; i<n; i++) {
1019 btAssert (L && d && b && n > 0 && nskip >= n);
1036 int do_fast_row_swaps)
1038 btAssert (A && n > 0 && i1 >= 0 && i2 >= 0 && i1 < n && i2 < n &&
1039 nskip >= n && i1 < i2);
1044 for (
int i=i1+1; i<i2; ++i) {
1049 A_i1[i2] = A_i1[i1];
1050 A_i1[i1] = A_i2[i1];
1051 A_i2[i1] = A_i2[i2];
1053 if (do_fast_row_swaps) {
1059 for (
int k = 0; k <= i2; ++k) {
1066 for (
int j=i2+1; j<n; ++j) {
1075 for (
int k = 0; k < i1; ++k) {
1081 for (
int i=i1+1; i<i2; A_i+=nskip, ++i) {
1088 A_i1[i1] = A_i2[i2];
1092 for (
int j=i2+1; j<n; A_j+=nskip, ++j) {
1104 btScalar *hi,
int *p,
bool *state,
int *findex,
1105 int n,
int i1,
int i2,
int nskip,
1106 int do_fast_row_swaps)
1111 btAssert (n>0 && i1 >=0 && i2 >= 0 && i1 < n && i2 < n && nskip >= n && i1 <= i2);
1141 state[i1] = state[i2];
1146 findex[i1] = findex[i2];
1204 bool *_state,
int *_findex,
int *_p,
int *_C,
btScalar **Arows);
1221 void solve1 (
btScalar *a,
int i,
int dir=1,
int only_transfer=0);
1229 bool *_state,
int *_findex,
int *_p,
int *_C,
btScalar **Arows):
1230 m_n(_n), m_nskip(_nskip), m_nub(_nub), m_nC(0), m_nN(0),
1236 m_x(_x), m_b(_b), m_w(_w), m_lo(_lo), m_hi(_hi),
1237 m_L(_L), m_d(_d), m_Dell(_Dell), m_ell(_ell), m_tmp(_tmp),
1238 m_state(_state), m_findex(_findex), m_p(_p), m_C(_C)
1250 for (
int k=0; k<n; aptr+=nskip, ++k) A[k] = aptr;
1257 for (
int k=0; k<n; ++k) p[k]=k;
1292 for (
int k =
m_nub; k<n; ++k) {
1293 if (findex && findex[k] >= 0)
continue;
1295 btSwapProblem (
m_A,
m_x,
m_b,
m_w,lo,hi,
m_p,
m_state,findex,n,
m_nub,k,
m_nskip,0);
1304 const int nub =
m_nub;
1308 for (
int j=0; j<nub; Lrow+=nskip, ++j) memcpy(Lrow,
BTAROW(j),(j+1)*
sizeof(
btScalar));
1316 for (
int k=0; k<nub; ++k) C[k] = k;
1323 const int nub =
m_nub;
1326 for (
int k=
m_n-1; k >= nub; k--) {
1327 if (findex[k] >= 0) {
1328 btSwapProblem (
m_A,
m_x,
m_b,
m_w,
m_lo,
m_hi,
m_p,
m_state,findex,
m_n,k,
m_n-1-num_at_end,
m_nskip,1);
1356 const int nC =
m_nC;
1358 for (
int j=0; j<nC; ++j) Ltgt[j] = ell[j];
1360 const int nC =
m_nC;
1367 btSwapProblem (
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);
1369 const int nC =
m_nC;
1385 # ifdef BTNUB_OPTIMIZATIONS
1387 const int nub =
m_nub;
1389 for ( ; j<nub; ++j) Dell[j] = aptr[j];
1390 const int nC =
m_nC;
1391 for ( ; j<nC; ++j) Dell[j] = aptr[C[j]];
1393 const int nC =
m_nC;
1394 for (
int j=0; j<nC; ++j) Dell[j] = aptr[C[j]];
1399 const int nC =
m_nC;
1402 for (
int j=0; j<nC; ++j) Ltgt[j] = ell[j] = Dell[j] * d[j];
1404 const int nC =
m_nC;
1411 btSwapProblem (
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);
1413 const int nC =
m_nC;
1429 btAssert(A && n > 0 && nskip >= n && r >= 0 && r < n);
1430 if (r >= n-1)
return;
1433 const size_t move_size = (n-r-1)*
sizeof(
btScalar);
1435 for (
int i=0; i<r; Adst+=nskip,++i) {
1437 memmove (Adst,Asrc,move_size);
1441 const size_t cpy_size = r*
sizeof(
btScalar);
1443 for (
int i=r; i<(n-1); ++i) {
1445 memcpy (Adst,Asrc,cpy_size);
1451 const size_t cpy_size = (n-r-1)*
sizeof(
btScalar);
1452 btScalar *Adst = A + r * (nskip + 1);
1453 for (
int i=r; i<(n-1); ++i) {
1454 btScalar *Asrc = Adst + (nskip + 1);
1455 memcpy (Adst,Asrc,cpy_size);
1466 btAssert (L && d && a && n > 0 && nskip >= n);
1476 for (
int j=1; j<n; ++j) {
1487 btScalar alphanew = alpha1 + (W11*W11)*dee;
1493 alphanew = alpha2 - (W21*W21)*dee;
1498 btScalar k2 = W21*gamma1*W11 - W21;
1500 for (
int p=1; p<n; ll+=nskip, ++p) {
1503 W1[p] = Wp - W11*ell;
1504 W2[p] = k1*Wp + k2*ell;
1509 for (
int j=1; j<n; ll+=nskip+1, ++j) {
1514 btScalar alphanew = alpha1 + (k1*k1)*dee;
1520 alphanew = alpha2 - (k2*k2)*dee;
1528 for (
int p=j+1; p<n; l+=nskip, ++p) {
1533 Wp = W2[p] - k2 * ell;
1542 #define _BTGETA(i,j) (A[i][j])
1544 #define BTGETA(i,j) ((i > j) ? _BTGETA(i,j) : _BTGETA(j,i))
1548 return nskip * 2 *
sizeof(
btScalar);
1555 btAssert(A && p && L && d && n1 > 0 && n2 > 0 && r >= 0 && r < n2 &&
1556 n1 >= n2 && nskip >= n1);
1558 for (
int i=0; i<n2; ++i)
1568 scratch.
resize(nskip * 2+n2);
1572 const int p_0 = p[0];
1573 for (
int i=0; i<n2; ++i) {
1574 a[i] = -
BTGETA(p[i],p_0);
1583 for (
int i=0; i<r; ++Lcurr, ++i) {
1585 t[i] = *Lcurr / d[i];
1591 const int *pp_r = p + r, p_r = *pp_r;
1592 const int n2_minus_r = n2-r;
1593 for (
int i=0; i<n2_minus_r; Lcurr+=nskip,++i) {
1598 btLDLTAddTL (L + r*nskip+r, d+r, a, n2-r, nskip, scratch);
1604 if (r < (n2-1)) memmove (d+r,d+r+1,(n2-r-1)*
sizeof(
btScalar));
1615 const int nC =
m_nC;
1617 for ( ; j<nC; ++j) {
1624 if (last_idx == -1) {
1625 for (k=j+1 ; k<nC; ++k) {
1636 if (j < (nC-1)) memmove (C+j,C+j+1,(nC-j-1)*
sizeof(
int));
1642 btSwapProblem (
m_A,
m_x,
m_b,
m_w,
m_lo,
m_hi,
m_p,
m_state,
m_findex,
m_n,i,nC-1,
m_nskip,1);
1658 const int nC =
m_nC;
1660 const int nN =
m_nN;
1661 for (
int i=0; i<nN; ++i) {
1669 const int nC =
m_nC;
1673 const int nN =
m_nN;
1674 for (
int j=0; j<nN; ++j) ptgt[j] += aptr[j];
1677 const int nN =
m_nN;
1678 for (
int j=0; j<nN; ++j) ptgt[j] -= aptr[j];
1684 const int nC =
m_nC;
1685 for (
int i=0; i<nC; ++i) {
1692 const int nC =
m_nC;
1693 btScalar *ptgt = p + nC, *qsrc = q + nC;
1694 const int nN =
m_nN;
1695 for (
int i=0; i<nN; ++i) {
1696 ptgt[i] += s*qsrc[i];
1713 # ifdef BTNUB_OPTIMIZATIONS
1715 const int nub =
m_nub;
1717 for ( ; j<nub; ++j) Dell[j] = aptr[j];
1718 const int nC =
m_nC;
1719 for ( ; j<nC; ++j) Dell[j] = aptr[C[j]];
1721 const int nC =
m_nC;
1722 for (
int j=0; j<nC; ++j) Dell[j] = aptr[C[j]];
1728 const int nC =
m_nC;
1729 for (
int j=0; j<nC; ++j) ell[j] = Dell[j] * d[j];
1732 if (!only_transfer) {
1735 const int nC =
m_nC;
1736 for (
int j=0; j<nC; ++j) tmp[j] = ell[j];
1742 const int nC =
m_nC;
1743 for (
int j=0; j<nC; ++j) a[C[j]] = -tmp[j];
1747 const int nC =
m_nC;
1748 for (
int j=0; j<nC; ++j) a[C[j]] = tmp[j];
1763 for (
int j=0; j<n; ++j) x[p[j]] = tmp[j];
1770 for (
int j=0; j<n; ++j) w[p[j]] = tmp[j];
1774 #endif // btLCP_FAST
1786 btAssert (n>0 && A && x && b && lo && hi && nub >= 0 && nub <= n);
1792 for (
int k=0; k<n; ++k)
1793 btAssert (lo[k] <= 0 && hi[k] >= 0);
1812 const int nskip = (n);
1832 btLCP lcp(n,nskip,nub,A,x,b,w,lo,hi,&scratchMem.
L[0],&scratchMem.
d[0],&scratchMem.
Dell[0],&scratchMem.
ell[0],&scratchMem.
delta_w[0],&scratchMem.
state[0],findex,&scratchMem.
p[0],&scratchMem.
C[0],&scratchMem.
Arows[0]);
1833 int adj_nub = lcp.
getNub();
1844 bool hit_first_friction_index =
false;
1845 for (
int i=adj_nub; i<n; ++i)
1861 if (!hit_first_friction_index && findex && findex[i] >= 0) {
1863 for (
int j=0; j<n; ++j) scratchMem.
delta_w[scratchMem.
p[j]] = x[j];
1866 for (
int k=i; k<n; ++k) {
1873 hi[k] =
btFabs (hi[k] * wfk);
1877 hit_first_friction_index =
true;
1895 if (lo[i]==0 && w[i] >= 0) {
1897 scratchMem.
state[i] =
false;
1899 else if (hi[i]==0 && w[i] <= 0) {
1901 scratchMem.
state[i] =
true;
1966 const int numN = lcp.
numN();
1967 for (
int k=0; k < numN; ++k) {
1968 const int indexN_k = lcp.
indexN(k);
1969 if (!scratchMem.
state[indexN_k] ? scratchMem.
delta_w[indexN_k] < 0 : scratchMem.
delta_w[indexN_k] > 0) {
1971 if (lo[indexN_k] == 0 && hi[indexN_k] == 0)
continue;
1983 const int numC = lcp.
numC();
1984 for (
int k=adj_nub; k < numC; ++k) {
1985 const int indexC_k = lcp.
indexC(k);
1987 btScalar s2 = (lo[indexC_k]-x[indexC_k]) / scratchMem.
delta_x[indexC_k];
1995 btScalar s2 = (hi[indexC_k]-x[indexC_k]) / scratchMem.
delta_x[indexC_k];
2029 w[i] += s * scratchMem.
delta_w[i];
2040 scratchMem.
state[i] =
false;
2045 scratchMem.
state[i] =
true;
2054 scratchMem.
state[si] =
false;
2059 scratchMem.
state[si] =
true;
2064 if (cmd <= 3)
break;
static T sum(const btAlignedObjectArray< T > &items)
void btRemoveRowCol(btScalar *A, int n, int nskip, int r)
void btLDLTAddTL(btScalar *L, btScalar *d, const btScalar *a, int n, int nskip, btAlignedObjectArray< btScalar > &scratch)
void btFactorLDLT(btScalar *A, btScalar *d, int n, int nskip1)
btLCP(int _n, int _nskip, int _nub, btScalar *_Adata, btScalar *_x, btScalar *_b, btScalar *_w, btScalar *_lo, btScalar *_hi, btScalar *_L, btScalar *_d, btScalar *_Dell, btScalar *_ell, btScalar *_tmp, bool *_state, int *_findex, int *_p, int *_C, btScalar **Arows)
btAlignedObjectArray< btScalar * > Arows
void pN_plusequals_ANi(btScalar *p, int i, int sign=1)
btAlignedObjectArray< btScalar > delta_w
btScalar AiN_times_qN(int i, btScalar *q) const
void transfer_i_from_N_to_C(int i)
void btLDLTRemove(btScalar **A, const int *p, btScalar *L, btScalar *d, int n1, int n2, int r, int nskip, btAlignedObjectArray< btScalar > &scratch)
btScalar *const *const *const m_w
btScalar Aii(int i) const
btScalar *const *const m_d
static void btSwapProblem(BTATYPE A, btScalar *x, btScalar *b, btScalar *w, btScalar *lo, btScalar *hi, int *p, bool *state, int *findex, int n, int i1, int i2, int nskip, int do_fast_row_swaps)
void solve1(btScalar *a, int i, int dir=1, int only_transfer=0)
btScalar *const *const m_ell
btScalar *const *const *const m_tmp
void btVectorScale(btScalar *a, const btScalar *d, int n)
btScalar *const *const m_b
btScalar *const *const *const *const *const m_hi
size_t btEstimateLDLTAddTLTmpbufSize(int nskip)
void transfer_i_to_N(int i)
void pN_equals_ANC_times_qC(btScalar *p, btScalar *q)
btAlignedObjectArray< btScalar > m_scratch
static void btSwapRowsAndCols(BTATYPE A, int n, int i1, int i2, int nskip, int do_fast_row_swaps)
btAlignedObjectArray< btScalar > Dell
void pC_plusequals_s_times_qC(btScalar *p, btScalar s, btScalar *q)
void btSetZero(T *a, int n)
void btSolveL1T(const btScalar *L, btScalar *B, int n, int lskip1)
void btSolveLDLT(const btScalar *L, const btScalar *d, btScalar *b, int n, int nskip)
static void btSolveL1_2(const btScalar *L, btScalar *B, int n, int lskip1)
btAlignedObjectArray< int > p
btAlignedObjectArray< int > C
static void btSolveL1_1(const btScalar *L, btScalar *B, int n, int lskip1)
bool btSolveDantzigLCP(int n, btScalar *A, btScalar *x, btScalar *b, btScalar *outer_w, int nub, btScalar *lo, btScalar *hi, int *findex, btDantzigScratchMemory &scratchMem)
void resize(int newsize, const T &fillData=T())
void transfer_i_from_C_to_N(int i, btAlignedObjectArray< btScalar > &scratch)
btScalar AiC_times_qC(int i, btScalar *q) const
btAlignedObjectArray< btScalar > ell
btScalar *const *const *const *const m_lo
btAlignedObjectArray< btScalar > delta_x
void btSolveL1(const btScalar *L, btScalar *B, int n, int lskip1)
btAlignedObjectArray< btScalar > L
btAlignedObjectArray< btScalar > d
void pN_plusequals_s_times_qN(btScalar *p, btScalar s, btScalar *q)
btAlignedObjectArray< bool > state
int *const *const *const m_C
void transfer_i_to_C(int i)
float btScalar
The btScalar type abstracts floating point numbers, to easily switch between double and single floati...
btScalar btLargeDot(const btScalar *a, const btScalar *b, int n)
btScalar btFabs(btScalar x)