2 #ifndef COIN_OSL_C_INCLUDE
5 #define COIN_OSL_C_INCLUDE
10 #define C_EKK_GO_SPARSE 200
14 #if __BYTE_ORDER == __LITTLE_ENDIAN
35 int * mpt,
int first_nonzero);
38 int * mpt,
int ipivrw,
int * spare);
42 int *mpt2,
double dalpha,
int orig_nincol,
43 int npivot,
int *nuspikp,
44 const int ipivrw,
int * spare);
47 double * dpermu,
int * mpt,
int numberNonZero);
50 double *dwork1,
int *mpt,
int *nincolp);
52 double * dpermu1,
int * mpt1,
int *nincolp,
53 double *dwork1_ft,
int *mpt_ft,
int *nincolp_ft);
57 inline void c_ekkscpy(
int n,
const int *marr1,
int *marr2)
59 inline void c_ekkdcpy(
int n,
const double *marr1,
double *marr2)
65 void c_ekkzero(
int length,
int n,
void * array);
76 #define c_ekkscpy_0_1(s,ival,array) CoinFillN(array,s,ival)
77 #define c_ekks1cpy( n,marr1,marr2) CoinMemcpyN(marr1,n, marr2)
81 int *
clp_int(
int number_entries);
85 #define SLACK_VALUE -1.0
86 #define C_EKK_REMOVE_LINK(hpiv,hin,link,ipivot) \
88 int ipre = link[ipivot].pre; \
89 int isuc = link[ipivot].suc; \
91 link[ipre].suc = isuc; \
94 hpiv[hin[ipivot]] = isuc; \
97 link[isuc].pre = ipre; \
101 #define C_EKK_ADD_LINK(hpiv,nzi,link, npr) \
103 int ifiri = hpiv[nzi]; \
105 link[npr].suc = ifiri; \
108 link[ifiri].pre = npr; \
114 #define SHIFT_INDEX(limit) (limit)
115 #define UNSHIFT_INDEX(limit) (limit)
116 #define SHIFT_REF(arr,ind) (arr)[ind]
120 #define SHIFT_INDEX(limit) ((limit)<<3)
121 #define UNSHIFT_INDEX(limit) ((unsigned int)(limit)>>3)
122 #define SHIFT_REF(arr,ind) (*(double*)((char*)(arr) + (ind)))
127 #define NOT_ZERO(x) (((*((reinterpret_cast<unsigned char *>(&x))+7)) & 0x7F) != 0)
129 #define NOT_ZERO(x) ((x) != 0.0)
132 #define SWAP(type,_x,_y) { type _tmp = (_x); (_x) = (_y); (_y) = _tmp;}
134 #define UNROLL_LOOP_BODY1(code) \
136 #define UNROLL_LOOP_BODY2(code) \
138 #define UNROLL_LOOP_BODY4(code) \
139 {{code} {code} {code} {code}}
149 int *hcoli = fact->xecadr;
150 double *dluval = fact->xeeadr;
151 double *dvalpv = fact->kw3adr;
152 int *mrstrt = fact->xrsadr;
153 int *hrowi = fact->xeradr;
154 int *mcstrt = fact->xcsadr;
155 int *hinrow = fact->xrnadr;
156 int *hincol = fact->xcnadr;
157 int *hpivro = fact->krpadr;
158 int *hpivco = fact->kcpadr;
160 int nnentl = fact->nnentl;
161 int nnentu = fact->nnentu;
162 int kmxeta = fact->kmxeta;
163 int xnewro = *xnewrop;
164 int ncompactions = *ncompactionsp;
166 MACTION_T *maction =
reinterpret_cast<MACTION_T*
>(maction_void);
171 int jj, kk, kr, nz, jj1, jj2, kce, kcs, kqq, npr;
175 int knpr, irow, iadd32, ibase;
181 int kipie, kcpiv, knprs, knpre;
183 double multip, elemnt;
184 int ipivot, jpivot, epivro, epivco, lstart, ifdens, nfirst;
185 int nzpivj, kfill, kstart;
189 int noRoomForDense=0;
190 int if_sparse_update=fact->if_sparse_update;
193 const int nrow = fact->nrow;
199 lstart = nnetas - nnentl + 1;
200 for (i = lstart; i <= nnetas; ++i) {
205 for (i = 1; i <= nrow; ++i) {
207 mwork[i].pre = i - 1;
208 mwork[i].suc = i + 1;
217 for (count = 1; count <= nrow; ++count) {
220 if (! (hpivco[1] <= 0)) {
221 int small_pivot = c_ekkcsin(fact,
227 if (fact->invok >= 0) {
231 if (fact->npivots >= nrow) {
237 if (! (hpivro[1] <= 0)) {
238 irtcod = c_ekkrsin(fact,
245 &kmxeta, &ncompactions,
248 if (irtcod < 0 || fact->invok >= 0) {
255 if (fact->npivots >= nrow) {
258 lstart = nnetas - nnentl + 1;
262 irtcod = c_ekkfpvt(fact,
264 nsingp, xrejctp, &ipivot, &jpivot);
270 c_ekkprpv(fact, rlink, clink,
271 *xrejctp, ipivot, jpivot);
273 epivco = hincol[jpivot];
275 mcstrt[fact->xnetal] = lstart - 1;
276 hpivco[fact->xnetal] = ipivot;
277 epivro = hinrow[ipivot];
279 kipis = mrstrt[ipivot];
280 pivot = dluval[kipis];
282 kipie = kipis + epivr1;
286 double size = nrow - fact->npivots;
287 if (size > GO_DENSE && (nnentu - fact->nuspike) * GO_DENSE_RATIO > size * size) {
298 for (k = kipis; k <= kipie; ++k) {
300 dvalpv[k - kipis + 1] = dluval[k];
301 maction[irow] =
static_cast<MACTION_T
>(k - kipis + 1);
305 kcpiv = mcstrt[jpivot] - 1;
306 for (nzpivj = 1; nzpivj <= epivco; ++nzpivj) {
319 knpre = knprs + enpr;
325 if (jpivot == hcoli[knpr]) {
332 multip = -dluval[knpr] * dpivx;
335 dluval[knpr] = dluval[knpre];
336 hcoli[knpr] = hcoli[knpre];
341 kfill = epivr1 - (knpre - knprs + 1);
342 nres = ((knpre - knprs + 1) & 1) + knprs;
349 if (maction[j] == 0) {
353 maction[j] =
static_cast<MACTION_T
>(-maction[j]);
354 dluval[knprs] += multip * dvalpv[jj];
355 d1 = fabs(dluval[knprs]);
358 j2 = hcoli[nres + 1];
360 for (kr = nres; kr < knpre; kr += 2) {
365 maction[j1] =
static_cast<MACTION_T
>(-maction[j1]);
366 dluval[kr] += multip * dvalpv[jj1];
367 cancel = cancel || ! (fact->zeroTolerance < d1);
368 d1 = fabs(dluval[kr]);
374 maction[j2] =
static_cast<MACTION_T
>(-maction[j2]);
375 dluval[kr + 1] += multip * dvalpv[jj2];
376 cancel = cancel || ! (fact->zeroTolerance < d1);
377 d1 = fabs(dluval[kr + 1]);
382 cancel = cancel || ! (fact->zeroTolerance < d1);
399 kfill = epivr1 - (knpre - knprs + 1);
402 for (kr = knprs; kr <= knpre; kr++) {
410 maction[j1] = -maction[j1];
411 dluval[kr] += multip * dvalpv[jj1];
412 d1 = fabs(dluval[kr]);
413 cancel = cancel || ! (fact->zeroTolerance < d1);
423 for (kr = knprs; kr <= knpre; ++kr) {
425 if (fabs(dluval[kr]) > fact->zeroTolerance) {
427 dluval[kstart] = dluval[kr];
435 kce = kcs + hincol[j];
436 for (kk = kcs; kk <= kce; ++kk) {
437 if (hrowi[kk] == npr) {
438 hrowi[kk] = hrowi[kce];
450 for (k = kipis; k <= kipie; ++k) {
451 maction[hcoli[k]] =
static_cast<MACTION_T
>(-maction[hcoli[k]]);
455 naft = mwork[npr].suc;
456 kqq = mrstrt[naft] - knpre - 1;
462 if (! (xnewro + nznpr + 1 < lstart)) {
463 if (! (nnentu + nznpr + 1 < lstart)) {
473 int iput = c_ekkrwcs(fact,dluval, hcoli, mrstrt, hinrow, mwork, nfirst);
474 kmxeta += xnewro - iput ;
479 kipis = mrstrt[ipivot] + 1;
480 kipie = kipis + epivr1 - 1;
487 knpre = knprs + enpr - 1;
495 npre = mwork[npr].pre;
500 mwork[naft].pre = npre;
501 mwork[npre].suc = naft;
503 mwork[nfirst].pre = npr;
504 mwork[nlast].suc = npr;
505 mwork[npr].pre = nlast;
506 mwork[npr].suc = nfirst;
509 mrstrt[npr] = kstart + 1;
510 nmove = knpre - knprs + 1;
511 ibase = kstart + 1 - knprs;
512 for (kr = knprs; kr <= knpre; ++kr) {
513 dluval[ibase + kr] = dluval[kr];
514 hcoli[ibase + kr] = hcoli[kr];
535 if (kstart + (nrow << 1) + 100 < lstart) {
536 ileft = ((nrow - fact->npivots + 32) & -32);
537 if (kstart + ileft * ileft + 32 < lstart) {
539 xnewro =
CoinMax(kstart,xnewro);
540 xnewro = (xnewro & -32) + ileft;
542 xnewro = ((kstart + 31) & -32);
550 }
else if (! (nnentu + kqq + 2 < lstart)) {
555 for (kr = kipis; kr <= kipie; ++kr) {
559 elemnt = multip * dvalpv[jj];
560 if (fabs(elemnt) > fact->zeroTolerance) {
562 dluval[kstart] = elemnt;
571 if (xnewco + 1 >= lstart) {
572 if (xnewco + nz + 1 >= lstart) {
574 if (nnentu + nz + 1 < lstart) {
575 xnewco = c_ekkclco(fact,hrowi, mcstrt, hincol, xnewco);
578 kcpiv = mcstrt[jpivot] - 1;
589 mcstrt[j] = xnewco + 1;
590 ibase = mcstrt[j] - kcs;
591 for (kk = kcs; kk <= kce; ++kk) {
592 hrowi[ibase + kk] = hrowi[kk];
595 kce = xnewco + kce - kcs + 1;
600 }
else if (hrowi[kce + 1] != 0) {
602 if (xnewco + nz + 1 >= lstart) {
604 if (nnentu + nz + 1 < lstart) {
605 xnewco = c_ekkclco(fact,hrowi, mcstrt, hincol, xnewco);
608 kcpiv = mcstrt[jpivot] - 1;
619 mcstrt[j] = xnewco + 1;
620 ibase = mcstrt[j] - kcs;
621 for (kk = kcs; kk <= kce; ++kk) {
622 hrowi[ibase + kk] = hrowi[kk];
625 kce = xnewco + kce - kcs + 1;
629 hrowi[kce + 1] = npr;
633 maction[j] =
static_cast<MACTION_T
>(-maction[j]);
640 hinrow[npr] = kstart - mrstrt[npr] + 1;
642 if (! (xnewco + 1 < lstart)) {
643 xnewco = c_ekkclco(fact,hrowi, mcstrt, hincol, xnewco);
646 kcpiv = mcstrt[jpivot] - 1;
648 if (! (xnewro + 1 < lstart)) {
649 int iput = c_ekkrwcs(fact,dluval, hcoli, mrstrt, hinrow, mwork, nfirst);
650 kmxeta += xnewro - iput ;
654 kipis = mrstrt[ipivot] + 1;
655 kipie = kipis + epivr1 - 1;
661 dluval[lstart] = multip;
664 #define INLINE_AFPV 3
685 double d = fabs(els[j]);
695 double maxaij=fabs(els[0]);
696 for (j=1;j<nel;j++) {
697 double d = fabs(els[j]);
715 double d = fabs(els[j]);
724 SWAP(
int, index[koff], index[0]);
725 SWAP(
double, els[koff], els[0]);
730 int nzi = hinrow[npr];
739 int nn = mcstrt[fact->xnetal] - lstart + 1;
740 c_ekkafpv(hrowi+lstart, hcoli, dluval, mrstrt, hinrow, nn);
744 for (k = kipis; k <= kipie; ++k) {
745 maction[hcoli[k]] = 0;
749 for (k = kipis; k <= kipie; ++k) {
753 ! ((clink[j].pre > nrow && nzj != 1))) {
758 for (k = kipis; k <= kipie; ++k) {
766 fact->nuspike += hinrow[ipivot];
771 int ndense = nrow - fact->npivots;
772 if (! (xnewro + ndense * ndense >= lstart)) {
775 c_ekkizero( nrow, reinterpret_cast<int *> (maction+1));
777 for (i = 1; i <= nrow; ++i) {
778 if (clink[i].pre >= 0) {
780 maction[i] =
static_cast<short int>(iput);
785 for (i = 1; i <= nrow; ++i) {
786 if (rlink[i].pre >= 0) {
787 nspare = nspare + ndense - hinrow[i];
790 if (iput != nrow - fact->npivots) {
792 c_ekkizero( nrow, reinterpret_cast<int *> (maction+1));
795 int iput = c_ekkrwcs(fact,dluval, hcoli, mrstrt, hinrow, mwork, nfirst);
796 kmxeta += xnewro - iput ;
801 if (xnewro + nspare + ndense * ndense >= lstart) {
802 c_ekkizero( nrow, reinterpret_cast<int *> (maction+1));
806 c_ekkrwct(fact,dluval, hcoli, mrstrt, hinrow, mwork,
807 rlink, maction, dvalpv,
810 if (nnentu + nnentl > nrow * 5 &&
811 (ndense*ndense)>(nnentu+nnentl)>>2 &&
813 fact->ndenuc = ndense;
815 irtcod = c_ekkcmfd(fact,
816 (reinterpret_cast<int*>(dvalpv)+1),
818 (reinterpret_cast<int*>(maction+1))+1,
832 if (!noRoomForDense) {
833 int etasize =
CoinMax(4*fact->nnentu+(nnetas-fact->nnentl)+1000,fact->eta_size);
834 noRoomForDense=ndense;
835 fact->eta_size=
CoinMin(static_cast<int>(1.2*fact->eta_size),etasize);
836 if (fact->maxNNetas>0&&fact->eta_size>
838 fact->eta_size=fact->maxNNetas;
849 int iput = c_ekkrwcs(fact,dluval, hcoli, mrstrt, hinrow, mwork, nfirst);
850 kmxeta += xnewro - iput;
857 mwork[nrow+1].pre = nfirst;
858 mwork[nrow+1].suc = nlast;
860 fact->nnentl = nnentl;
861 fact->nnentu = nnentu;
862 fact->kmxeta = kmxeta;
864 *ncompactionsp = ncompactions;
int c_ekkbtrn_ipivrw(register const EKKfactinfo *fact, double *dwork1, int *mpt, int ipivrw, int *spare)
void clp_free(void *oldArray)
#define SHIFT_INDEX(limit)
void c_ekkczero(int n, char *marray)
double * clp_double(int number_entries)
T CoinMax(register const T x1, register const T x2)
Return the larger (according to operator<() of the arguments.
#define C_EKK_ADD_LINK(hpiv, nzi, link, npr)
void c_ekkdcpy(int n, const double *marr1, double *marr2)
void c_ekkizero(int n, int *marray)
void CoinZeroN(register T *to, const int size)
This helper function fills an array with zero.
int c_ekklfct(register EKKfactinfo *fact)
int * clp_int(int number_entries)
void clp_memory(int type)
int c_ekk_IsSet(const int *array, int bit)
int c_ekketsj(registerEKKfactinfo *fact, double *dwork1, int *mpt2, double dalpha, int orig_nincol, int npivot, int *nuspikp, const int ipivrw, int *spare)
void c_ekk_Unset(int *array, int bit)
void c_ekkftrn2(register EKKfactinfo *fact, double *dwork1, double *dpermu1, int *mpt1, int *nincolp, double *dwork1_ft, int *mpt_ft, int *nincolp_ft)
#define SWAP(type, _x, _y)
void CoinMemcpyN(register const T *from, const int size, register T *to)
This helper function copies an array to another location.
void * clp_malloc(int number_entries)
void c_ekkzero(int length, int n, void *array)
#define UNROLL_LOOP_BODY2(code)
void c_ekkscpy(int n, const int *marr1, int *marr2)
int c_ekkftrn(register const EKKfactinfo *fact, double *dwork1, double *dpermu, int *mpt, int numberNonZero)
int c_ekkslcf(register const EKKfactinfo *fact)
void c_ekk_Set(int *array, int bit)
void clp_setup_pointers(EKKfactinfo *fact)
T CoinMin(register const T x1, register const T x2)
Return the smaller (according to operator<() of the arguments.
int c_ekkbtrn(register const EKKfactinfo *fact, double *dwork1, int *mpt, int first_nonzero)
#define UNROLL_LOOP_BODY4(code)
void c_ekkdzero(int n, double *marray)
int c_ekkftrn_ft(register EKKfactinfo *fact, double *dwork1, int *mpt, int *nincolp)