1 /* $Id: CoinOslFactorization3.cpp 2083 2019-01-06 19:38:09Z unxusr $ */
2 /*
3 Copyright (C) 1987, 2009, International Business Machines
4 Corporation and others. All Rights Reserved.
5
6 This code is licensed under the terms of the Eclipse Public License (EPL).
7 */
8 #if COIN_BIG_INDEX == 0
9 #include "CoinOslFactorization.hpp"
10 #include "CoinOslC.h"
11 #include "CoinFinite.hpp"
12 #define GO_DENSE 70
13 #define GO_DENSE_RATIO 1.8
14 int c_ekkclco(const EKKfactinfo *fact, int *hcoli,
15 int *mrstrt, int *hinrow, int xnewro);
16
17 void c_ekkclcp(const int *hcol, const double *dels, const int *mrstrt,
18 int *hrow, double *dels2, int *mcstrt,
19 int *hincol, int itype, int nnrow, int nncol,
20 int ninbas);
21
22 int c_ekkcmfc(EKKfactinfo *fact,
23 EKKHlink *rlink, EKKHlink *clink,
24 EKKHlink *mwork, void *maction_void,
25 int nnetas,
26 int *nsingp, int *xrejctp,
27 int *xnewrop, int xnewco,
28 int *ncompactionsp);
29
30 int c_ekkcmfy(EKKfactinfo *fact,
31 EKKHlink *rlink, EKKHlink *clink,
32 EKKHlink *mwork, void *maction_void,
33 int nnetas,
34 int *nsingp, int *xrejctp,
35 int *xnewrop, int xnewco,
36 int *ncompactionsp);
37
38 int c_ekkcmfd(EKKfactinfo *fact,
39 int *mcol,
40 EKKHlink *rlink, EKKHlink *clink,
41 int *maction,
42 int nnetas,
43 int *nnentlp, int *nnentup,
44 int *nsingp);
45 int c_ekkford(const EKKfactinfo *fact, const int *hinrow, const int *hincol,
46 int *hpivro, int *hpivco,
47 EKKHlink *rlink, EKKHlink *clink);
48 void c_ekkrowq(int *hrow, int *hcol, double *dels,
49 int *mrstrt,
50 const int *hinrow, int nnrow, int ninbas);
51 int c_ekkrwco(const EKKfactinfo *fact, double *dluval, int *hcoli, int *mrstrt, int *hinrow, int xnewro);
52
53 int c_ekkrwcs(const EKKfactinfo *fact, double *dluval, int *hcoli, int *mrstrt,
54 const int *hinrow, const EKKHlink *mwork,
55 int nfirst);
56
57 void c_ekkrwct(const EKKfactinfo *fact, double *dluval, int *hcoli, int *mrstrt,
58 const int *hinrow, const EKKHlink *mwork,
59 const EKKHlink *rlink,
60 const short *msort, double *dsort,
61 int nlast, int xnewro);
62
63 int c_ekkshff(EKKfactinfo *fact,
64 EKKHlink *clink, EKKHlink *rlink,
65 int xnewro);
66
67 void c_ekkshfv(EKKfactinfo *fact, EKKHlink *rlink, EKKHlink *clink,
68 int xnewro);
69 int c_ekktria(EKKfactinfo *fact,
70 EKKHlink *rlink,
71 EKKHlink *clink,
72 int *nsingp,
73 int *xnewcop, int *xnewrop,
74 int *nlrowtp,
75 const int ninbas);
76 #if 0
77 static void c_ekkafpv(int *hentry, int *hcoli,
78 double *dluval, int *mrstrt,
79 int *hinrow, int nentry)
80 {
81 int j;
82 int nel, krs;
83 int koff;
84 int irow;
85 int ientry;
86 int * index;
87
88 for (ientry = 0; ientry < nentry; ++ientry) {
89 #ifdef INTEL
90 int * els_long,maxaij_long;
91 #endif
92 double * els;
93 irow = UNSHIFT_INDEX(hentry[ientry]);
94 nel = hinrow[irow];
95 krs = mrstrt[irow];
96 index=&hcoli[krs];
97 els=&dluval[krs];
98 #ifdef INTEL
99 els_long=reinterpret_cast<int *> (els);
100 maxaij_long=0;
101 #else
102 double maxaij = 0.f;
103 #endif
104 koff = 0;
105 j=0;
106 if ((nel&1)!=0) {
107 #ifdef INTEL
108 maxaij_long = els_long[1] & 0x7fffffff;
109 #else
110 maxaij=fabs(els[0]);
111 #endif
112 j=1;
113 }
114
115 while (j<nel) {
116 #ifdef INTEL
117 UNROLL_LOOP_BODY2({
118 int d_long = els_long[1+(j<<1)] & 0x7fffffff;
119 if (maxaij_long < d_long) {
120 maxaij_long = d_long;
121 koff=j;
122 }
123 j++;
124 });
125 #else
126 UNROLL_LOOP_BODY2({
127 double d = fabs(els[j]);
128 if (maxaij < d) {
129 maxaij = d;
130 koff=j;
131 }
132 j++;
133 });
134 #endif
135 }
136
137 SWAP(int, index[koff], index[0]);
138 SWAP(double, els[koff], els[0]);
139 }
140 } /* c_ekkafpv */
141 #endif
142
143 /* Uwe H. Suhl, March 1987 */
144 /* This routine processes col singletons during the LU-factorization. */
145 /* Return codes (checked version 1.11): */
146 /* 0: ok */
147 /* 6: pivot element too small */
148 /* -43: system error at label 420 (ipivot not found) */
149 /*
150 * This routine processes singleton columns during factorization of the
151 * nucleus. It is very similar to the first part of c_ekktria,
152 * but is more complex, because we now have to maintain the length
153 * lists.
154 * The differences are:
155 * (1) here we use the length list for length 1 rather than a queue.
156 * This routine is only called if it is known that there is a singleton
157 * column.
158 *
159 * (2) here we maintain hrowi by moving the last entry into the pivot
160 * column entry; that means we don't have to search for the pivot row
161 * entry like we do in c_ekktria.
162 *
163 * (3) here the hlink data structure is in use for the length lists,
164 * so we maintain it as we shorten rows and removing columns altogether.
165 *
166 */
c_ekkcsin(EKKfactinfo * fact,EKKHlink * rlink,EKKHlink * clink,int * nsingp)167 int c_ekkcsin(EKKfactinfo *fact,
168 EKKHlink *rlink, EKKHlink *clink,
169
170 int *nsingp)
171 {
172 #if 1
173 int *hcoli = fact->xecadr;
174 double *dluval = fact->xeeadr;
175 //double *dvalpv = fact->kw3adr;
176 int *mrstrt = fact->xrsadr;
177 int *hrowi = fact->xeradr;
178 int *mcstrt = fact->xcsadr;
179 int *hinrow = fact->xrnadr;
180 int *hincol = fact->xcnadr;
181 int *hpivro = fact->krpadr;
182 int *hpivco = fact->kcpadr;
183 #endif
184 const int nrow = fact->nrow;
185 const double drtpiv = fact->drtpiv;
186
187 int j, k, kc, kce, kcs, nzj;
188 double pivot;
189 int kipis, kipie;
190 int jpivot;
191 #ifndef NDEBUG
192 int kpivot = -1;
193 #else
194 int kpivot = -1;
195 #endif
196
197 bool small_pivot = false;
198
199 /* next singleton column.
200 * Note that when the pivot column itself was removed from the
201 * list, the column in the list after it (if any) moves to the
202 * head of the list.
203 * Also, if any column from the pivot row was reduced to length 1,
204 * then it will have been added to the list and now be in front.
205 */
206 for (jpivot = hpivco[1]; jpivot > 0; jpivot = hpivco[1]) {
207 const int ipivot = hrowi[mcstrt[jpivot]]; /* (2) */
208 assert(ipivot);
209 /* The pivot row is being eliminated (3) */
210 C_EKK_REMOVE_LINK(hpivro, hinrow, rlink, ipivot);
211
212 /* Loop over nonzeros in pivot row: */
213 kipis = mrstrt[ipivot];
214 kipie = kipis + hinrow[ipivot] - 1;
215 for (k = kipis; k <= kipie; ++k) {
216 j = hcoli[k];
217
218 /*
219 * We're eliminating column jpivot,
220 * so we're eliminating the row it occurs in,
221 * so every column in this row is becoming one shorter.
222 *
223 * I don't know why we don't do the same for rejected columns.
224 *
225 * if xrejct is false, then no column has ever been rejected
226 * and this test wouldn't have to be made.
227 * However, that means this whole loop would have to be copied.
228 */
229 if (!(clink[j].pre > nrow)) {
230 C_EKK_REMOVE_LINK(hpivco, hincol, clink, j); /* (3) */
231 }
232 --hincol[j];
233
234 kcs = mcstrt[j];
235 kce = kcs + hincol[j];
236 for (kc = kcs; kc <= kce; ++kc) {
237 if (ipivot == hrowi[kc]) {
238 break;
239 }
240 }
241 /* ASSERT !(kc>kce) */
242
243 /* (2) */
244 hrowi[kc] = hrowi[kce];
245 hrowi[kce] = 0;
246
247 if (j == jpivot) {
248 /* remember the slot corresponding to the pivot column */
249 kpivot = k;
250 } else {
251 /*
252 * We just reduced the length of the column.
253 * If we haven't eliminated all of its elements completely,
254 * then we have to put it back in its new length list.
255 *
256 * If the column was rejected, we only put it back in a length
257 * list when it has been reduced to a singleton column,
258 * because it would just be rejected again.
259 */
260 nzj = hincol[j];
261 if (!(nzj <= 0) && !(clink[j].pre > nrow && nzj != 1)) {
262 C_EKK_ADD_LINK(hpivco, nzj, clink, j); /* (3) */
263 }
264 }
265 }
266 assert(kpivot > 0);
267
268 /* store pivot sequence number */
269 ++fact->npivots;
270 rlink[ipivot].pre = -fact->npivots;
271 clink[jpivot].pre = -fact->npivots;
272
273 /* compute how much room we'll need later */
274 fact->nuspike += hinrow[ipivot];
275
276 /* check the pivot */
277 pivot = dluval[kpivot];
278 if (fabs(pivot) < drtpiv) {
279 /* pivot element too small */
280 small_pivot = true;
281 rlink[ipivot].pre = -nrow - 1;
282 clink[jpivot].pre = -nrow - 1;
283 ++(*nsingp);
284 }
285
286 /* swap the pivoted column entry with the first entry in the row */
287 dluval[kpivot] = dluval[kipis];
288 dluval[kipis] = pivot;
289 hcoli[kpivot] = hcoli[kipis];
290 hcoli[kipis] = jpivot;
291 }
292
293 return (small_pivot);
294 } /* c_ekkcsin */
295
296 /* Uwe H. Suhl, March 1987 */
297 /* This routine processes row singletons during the computation of */
298 /* an LU-decomposition for the nucleus. */
299 /* Return codes (checked version 1.16): */
300 /* -5: not enough space in row file */
301 /* -6: not enough space in column file */
302 /* 7: pivot element too small */
303 /* -52: system error at label 220 (ipivot not found) */
304 /* -53: system error at label 400 (jpivot not found) */
c_ekkrsin(EKKfactinfo * fact,EKKHlink * rlink,EKKHlink * clink,EKKHlink * mwork,int nfirst,int * nsingp,int * xnewcop,int * xnewrop,int * nnentup,int * kmxetap,int * ncompactionsp,int * nnentlp)305 int c_ekkrsin(EKKfactinfo *fact,
306 EKKHlink *rlink, EKKHlink *clink,
307 EKKHlink *mwork, int nfirst,
308 int *nsingp,
309 int *xnewcop, int *xnewrop,
310 int *nnentup,
311 int *kmxetap, int *ncompactionsp,
312 int *nnentlp)
313
314 {
315 #if 1
316 int *hcoli = fact->xecadr;
317 double *dluval = fact->xeeadr;
318 //double *dvalpv = fact->kw3adr;
319 int *mrstrt = fact->xrsadr;
320 int *hrowi = fact->xeradr;
321 int *mcstrt = fact->xcsadr;
322 int *hinrow = fact->xrnadr;
323 int *hincol = fact->xcnadr;
324 int *hpivro = fact->krpadr;
325 int *hpivco = fact->kcpadr;
326 #endif
327 const int nrow = fact->nrow;
328 const double drtpiv = fact->drtpiv;
329
330 int xnewro = *xnewrop;
331 int xnewco = *xnewcop;
332 int kmxeta = *kmxetap;
333 int nnentu = *nnentup;
334 int ncompactions = *ncompactionsp;
335 int nnentl = *nnentlp;
336
337 int i, j, k, kc, kr, npr, nzi;
338 double pivot;
339 int kjpis, kjpie, knprs, knpre;
340 double elemnt, maxaij;
341 int ipivot, epivco, lstart;
342 #ifndef NDEBUG
343 int kpivot = -1;
344 #else
345 int kpivot = -1;
346 #endif
347 int irtcod = 0;
348 const int nnetas = fact->nnetas;
349
350 lstart = nnetas - nnentl + 1;
351
352 for (ipivot = hpivro[1]; ipivot > 0; ipivot = hpivro[1]) {
353 const int jpivot = hcoli[mrstrt[ipivot]];
354
355 kjpis = mcstrt[jpivot];
356 kjpie = kjpis + hincol[jpivot];
357 for (k = kjpis; k < kjpie; ++k) {
358 i = hrowi[k];
359
360 /*
361 * We're eliminating row ipivot,
362 * so we're eliminating the column it occurs in,
363 * so every row in this column is becoming one shorter.
364 *
365 * No exception is made for rejected rows.
366 */
367 C_EKK_REMOVE_LINK(hpivro, hinrow, rlink, i);
368 }
369
370 /* The pivot column is being eliminated */
371 /* I don't know why there is an exception for rejected columns */
372 if (!(clink[jpivot].pre > nrow)) {
373 C_EKK_REMOVE_LINK(hpivco, hincol, clink, jpivot);
374 }
375
376 epivco = hincol[jpivot] - 1;
377 kjpie = kjpis + epivco;
378 for (kc = kjpis; kc <= kjpie; ++kc) {
379 if (ipivot == hrowi[kc]) {
380 break;
381 }
382 }
383 /* ASSERT !(kc>kjpie) */
384
385 /* move the last column entry into this deleted one to keep */
386 /* the entries compact */
387 hrowi[kc] = hrowi[kjpie];
388
389 hrowi[kjpie] = 0;
390
391 /* store pivot sequence number */
392 ++fact->npivots;
393 rlink[ipivot].pre = -fact->npivots;
394 clink[jpivot].pre = -fact->npivots;
395
396 /* Check if row or column files have to be compressed */
397 if (!(xnewro + epivco < lstart)) {
398 if (!(nnentu + epivco < lstart)) {
399 return (-5);
400 }
401 {
402 int iput = c_ekkrwcs(fact, dluval, hcoli, mrstrt, hinrow, mwork, nfirst);
403 kmxeta += xnewro - iput;
404 xnewro = iput - 1;
405 ++ncompactions;
406 }
407 }
408 if (!(xnewco + epivco < lstart)) {
409 if (!(nnentu + epivco < lstart)) {
410 return (-5);
411 }
412 xnewco = c_ekkclco(fact, hrowi, mcstrt, hincol, xnewco);
413 ++ncompactions;
414 }
415
416 /* This column has no more entries in it */
417 hincol[jpivot] = 0;
418
419 /* Perform numerical part of elimination. */
420 pivot = dluval[mrstrt[ipivot]];
421 if (fabs(pivot) < drtpiv) {
422 irtcod = 7;
423 rlink[ipivot].pre = -nrow - 1;
424 clink[jpivot].pre = -nrow - 1;
425 ++(*nsingp);
426 }
427
428 /* If epivco is 0, then we can treat this like a singleton column (?)*/
429 if (!(epivco <= 0)) {
430 ++fact->xnetal;
431 mcstrt[fact->xnetal] = lstart - 1;
432 hpivco[fact->xnetal] = ipivot;
433
434 /* Loop over nonzeros in pivot column. */
435 kjpis = mcstrt[jpivot];
436 kjpie = kjpis + epivco;
437 nnentl += epivco;
438 nnentu -= epivco;
439 for (kc = kjpis; kc < kjpie; ++kc) {
440 npr = hrowi[kc];
441 /* zero out the row entries as we go along */
442 hrowi[kc] = 0;
443
444 /* each row in the column is getting shorter */
445 --hinrow[npr];
446
447 /* find the entry in this row for the pivot column */
448 knprs = mrstrt[npr];
449 knpre = knprs + hinrow[npr];
450 for (kr = knprs; kr <= knpre; ++kr) {
451 if (jpivot == hcoli[kr])
452 break;
453 }
454 /* ASSERT !(kr>knpre) */
455
456 elemnt = dluval[kr];
457 /* move the last pivot column entry into this one */
458 /* to keep entries compact */
459 dluval[kr] = dluval[knpre];
460 hcoli[kr] = hcoli[knpre];
461
462 /*
463 * c_ekkmltf put the largest entries in front, and
464 * we want to maintain that property.
465 * There is only a problem if we just pivoted out the first
466 * entry, and there is more than one entry in the list.
467 */
468 if (!(kr != knprs || hinrow[npr] <= 1)) {
469 maxaij = 0.f;
470 for (k = knprs; k <= knpre; ++k) {
471 if (!(fabs(dluval[k]) <= maxaij)) {
472 maxaij = fabs(dluval[k]);
473 kpivot = k;
474 }
475 }
476 assert(kpivot > 0);
477 maxaij = dluval[kpivot];
478 dluval[kpivot] = dluval[knprs];
479 dluval[knprs] = maxaij;
480
481 j = hcoli[kpivot];
482 hcoli[kpivot] = hcoli[knprs];
483 hcoli[knprs] = j;
484 }
485
486 /* store elementary row transformation */
487 --lstart;
488 dluval[lstart] = -elemnt / pivot;
489 hrowi[lstart] = SHIFT_INDEX(npr);
490
491 /* Only add the row back in a length list if it isn't empty */
492 nzi = hinrow[npr];
493 if (!(nzi <= 0)) {
494 C_EKK_ADD_LINK(hpivro, nzi, rlink, npr);
495 }
496 }
497 ++fact->nuspike;
498 }
499 }
500
501 *xnewrop = xnewro;
502 *xnewcop = xnewco;
503 *kmxetap = kmxeta;
504 *nnentup = nnentu;
505 *ncompactionsp = ncompactions;
506 *nnentlp = nnentl;
507
508 return (irtcod);
509 } /* c_ekkrsin */
510
c_ekkfpvt(const EKKfactinfo * fact,EKKHlink * rlink,EKKHlink * clink,int * nsingp,int * xrejctp,int * xipivtp,int * xjpivtp)511 int c_ekkfpvt(const EKKfactinfo *fact,
512 EKKHlink *rlink, EKKHlink *clink,
513 int *nsingp, int *xrejctp,
514 int *xipivtp, int *xjpivtp)
515 {
516 double zpivlu = fact->zpivlu;
517 #if 1
518 int *hcoli = fact->xecadr;
519 double *dluval = fact->xeeadr;
520 //double *dvalpv = fact->kw3adr;
521 int *mrstrt = fact->xrsadr;
522 int *hrowi = fact->xeradr;
523 int *mcstrt = fact->xcsadr;
524 int *hinrow = fact->xrnadr;
525 int *hincol = fact->xcnadr;
526 int *hpivro = fact->krpadr;
527 int *hpivco = fact->kcpadr;
528 #endif
529 int i, j, k, ke, kk, ks, nz, nz1, kce, kcs, kre, krs;
530 double minsze;
531 int marcst, mincst, mincnt, trials, nentri;
532 int jpivot = -1;
533 bool rjectd;
534 int ipivot;
535 const int nrow = fact->nrow;
536 int irtcod = 0;
537
538 /* this used to be initialized in c_ekklfct */
539 const int xtrial = 1;
540
541 trials = 0;
542 ipivot = 0;
543 mincst = COIN_INT_MAX;
544 mincnt = COIN_INT_MAX;
545 for (nz = 2; nz <= nrow; ++nz) {
546 nz1 = nz - 1;
547 if (mincnt <= nz) {
548 goto L900;
549 }
550
551 /* Search rows for a pivot */
552 for (i = hpivro[nz]; !(i <= 0); i = rlink[i].suc) {
553
554 ks = mrstrt[i];
555 ke = ks + nz - 1;
556 /* Determine magnitude of minimal acceptable element */
557 minsze = fabs(dluval[ks]) * zpivlu;
558 for (k = ks; k <= ke; ++k) {
559 /* Consider a column only if it passes the stability test */
560 if (!(fabs(dluval[k]) < minsze)) {
561 j = hcoli[k];
562 marcst = nz1 * hincol[j];
563 if (!(marcst >= mincst)) {
564 mincst = marcst;
565 mincnt = hincol[j];
566 ipivot = i;
567 jpivot = j;
568 if (mincnt <= nz + 1) {
569 goto L900;
570 }
571 }
572 }
573 }
574 ++trials;
575
576 if (trials >= xtrial) {
577 goto L900;
578 }
579 }
580
581 /* Search columns for a pivot */
582 j = hpivco[nz];
583 while (!(j <= 0)) {
584 /* XSEARD = XSEARD + 1 */
585 rjectd = false;
586 kcs = mcstrt[j];
587 kce = kcs + nz - 1;
588 for (k = kcs; k <= kce; ++k) {
589 i = hrowi[k];
590 nentri = hinrow[i];
591 marcst = nz1 * nentri;
592 if (!(marcst >= mincst)) {
593 /* Determine magnitude of minimal acceptable element */
594 minsze = fabs(dluval[mrstrt[i]]) * zpivlu;
595 krs = mrstrt[i];
596 kre = krs + nentri - 1;
597 for (kk = krs; kk <= kre; ++kk) {
598 if (hcoli[kk] == j)
599 break;
600 }
601 /* ASSERT (kk <= kre) */
602
603 /* perform stability test */
604 if (!(fabs(dluval[kk]) < minsze)) {
605 mincst = marcst;
606 mincnt = nentri;
607 ipivot = i;
608 jpivot = j;
609 rjectd = false;
610 if (mincnt <= nz) {
611 goto L900;
612 }
613 } else {
614 if (ipivot == 0) {
615 rjectd = true;
616 }
617 }
618 }
619 }
620 ++trials;
621 if (trials >= xtrial && ipivot > 0) {
622 goto L900;
623 }
624 if (rjectd) {
625 int jsuc = clink[j].suc;
626 ++(*xrejctp);
627 C_EKK_REMOVE_LINK(hpivco, hincol, clink, j);
628 clink[j].pre = nrow + 1;
629 j = jsuc;
630 } else {
631 j = clink[j].suc;
632 }
633 }
634 }
635
636 /* FLAG REJECTED ROWS (should this be columns ?) */
637 for (j = 1; j <= nrow; ++j) {
638 if (hinrow[j] == 0) {
639 rlink[j].pre = -nrow - 1;
640 ++(*nsingp);
641 }
642 }
643 irtcod = 10;
644
645 L900:
646 *xipivtp = ipivot;
647 *xjpivtp = jpivot;
648 return (irtcod);
649 } /* c_ekkfpvt */
c_ekkprpv(EKKfactinfo * fact,EKKHlink * rlink,EKKHlink * clink,int xrejct,int ipivot,int jpivot)650 void c_ekkprpv(EKKfactinfo *fact,
651 EKKHlink *rlink, EKKHlink *clink,
652
653 int xrejct,
654 int ipivot, int jpivot)
655 {
656 #if 1
657 int *hcoli = fact->xecadr;
658 double *dluval = fact->xeeadr;
659 //double *dvalpv = fact->kw3adr;
660 int *mrstrt = fact->xrsadr;
661 int *hrowi = fact->xeradr;
662 int *mcstrt = fact->xcsadr;
663 int *hinrow = fact->xrnadr;
664 int *hincol = fact->xcnadr;
665 int *hpivro = fact->krpadr;
666 int *hpivco = fact->kcpadr;
667 #endif
668 int i, k;
669 int kc;
670 double pivot;
671 int kipis = mrstrt[ipivot];
672 int kipie = kipis + hinrow[ipivot] - 1;
673
674 #ifndef NDEBUG
675 int kpivot = -1;
676 #else
677 int kpivot = -1;
678 #endif
679 const int nrow = fact->nrow;
680
681 /* Update data structures */
682 {
683 int kjpis = mcstrt[jpivot];
684 int kjpie = kjpis + hincol[jpivot];
685 for (k = kjpis; k < kjpie; ++k) {
686 i = hrowi[k];
687 C_EKK_REMOVE_LINK(hpivro, hinrow, rlink, i);
688 }
689 }
690
691 for (k = kipis; k <= kipie; ++k) {
692 int j = hcoli[k];
693
694 if ((xrejct == 0) || !(clink[j].pre > nrow)) {
695 C_EKK_REMOVE_LINK(hpivco, hincol, clink, j);
696 }
697
698 --hincol[j];
699 int kcs = mcstrt[j];
700 int kce = kcs + hincol[j];
701
702 for (kc = kcs; kc < kce; kc++) {
703 if (hrowi[kc] == ipivot)
704 break;
705 }
706 assert(kc < kce || hrowi[kce] == ipivot);
707 hrowi[kc] = hrowi[kce];
708 hrowi[kce] = 0;
709 if (j == jpivot) {
710 kpivot = k;
711 }
712 }
713 assert(kpivot > 0);
714
715 /* Store the pivot sequence number */
716 ++fact->npivots;
717 rlink[ipivot].pre = -fact->npivots;
718 clink[jpivot].pre = -fact->npivots;
719
720 pivot = dluval[kpivot];
721 dluval[kpivot] = dluval[kipis];
722 dluval[kipis] = pivot;
723 hcoli[kpivot] = hcoli[kipis];
724 hcoli[kipis] = jpivot;
725
726 } /* c_ekkprpv */
727
728 /*
729 * c_ekkclco is almost exactly like c_ekkrwco.
730 */
c_ekkclco(const EKKfactinfo * fact,int * hcoli,int * mrstrt,int * hinrow,int xnewro)731 int c_ekkclco(const EKKfactinfo *fact, int *hcoli, int *mrstrt, int *hinrow, int xnewro)
732 {
733 #if 0
734 int *hcoli = fact->xecadr;
735 double *dluval = fact->xeeadr;
736 double *dvalpv = fact->kw3adr;
737 int *mrstrt = fact->xrsadr;
738 int *hrowi = fact->xeradr;
739 int *mcstrt = fact->xcsadr;
740 int *hinrow = fact->xrnadr;
741 int *hincol = fact->xcnadr;
742 int *hpivro = fact->krpadr;
743 int *hpivco = fact->kcpadr;
744 #endif
745 int i, k, nz, kold;
746 int kstart;
747 const int nrow = fact->nrow;
748
749 for (i = 1; i <= nrow; ++i) {
750 nz = hinrow[i];
751 if (0 < nz) {
752 /* save the last column entry of row i in hinrow */
753 /* and replace that entry with -i */
754 k = mrstrt[i] + nz - 1;
755 hinrow[i] = hcoli[k];
756 hcoli[k] = -i;
757 }
758 }
759
760 kstart = 0;
761 kold = 0;
762 for (k = 1; k <= xnewro; ++k) {
763 if (hcoli[k] != 0) {
764 ++kstart;
765
766 /* if this is the last entry for the row... */
767 if (hcoli[k] < 0) {
768 /* restore the entry */
769 i = -hcoli[k];
770 hcoli[k] = hinrow[i];
771
772 /* update mrstart and hinrow */
773 mrstrt[i] = kold + 1;
774 hinrow[i] = kstart - kold;
775 kold = kstart;
776 }
777
778 hcoli[kstart] = hcoli[k];
779 }
780 }
781
782 /* INSERTED INCASE CALLED FROM YTRIAN JJHF */
783 mrstrt[nrow + 1] = kstart + 1;
784
785 return (kstart);
786 } /* c_ekkclco */
787
788 #undef MACTION_T
789 #define COIN_OSL_CMFC
790 #define MACTION_T short int
c_ekkcmfc(EKKfactinfo * fact,EKKHlink * rlink,EKKHlink * clink,EKKHlink * mwork,void * maction_void,int nnetas,int * nsingp,int * xrejctp,int * xnewrop,int xnewco,int * ncompactionsp)791 int c_ekkcmfc(EKKfactinfo *fact,
792 EKKHlink *rlink, EKKHlink *clink,
793 EKKHlink *mwork, void *maction_void,
794 int nnetas,
795 int *nsingp, int *xrejctp,
796 int *xnewrop, int xnewco,
797 int *ncompactionsp)
798
799 #include "CoinOslC.h"
800 #undef COIN_OSL_CMFC
801 #undef MACTION_T
802 static int c_ekkidmx(int n, const double *dx)
803 {
804 int ret_val;
805 int i;
806 double dmax;
807 --dx;
808
809 /* Function Body */
810
811 if (n < 1) {
812 return (0);
813 }
814
815 if (n == 1) {
816 return (1);
817 }
818
819 ret_val = 1;
820 dmax = fabs(dx[1]);
821 for (i = 2; i <= n; ++i) {
822 if (fabs(dx[i]) > dmax) {
823 ret_val = i;
824 dmax = fabs(dx[i]);
825 }
826 }
827
828 return ret_val;
829 } /* c_ekkidmx */
830 /* Return codes in IRTCOD/IRTCOD are */
831 /* 4: numerical problems */
832 /* 5: not enough space in row file */
833 /* 6: not enough space in column file */
c_ekkcmfd(EKKfactinfo * fact,int * mcol,EKKHlink * rlink,EKKHlink * clink,int * maction,int nnetas,int * nnentlp,int * nnentup,int * nsingp)834 int c_ekkcmfd(EKKfactinfo *fact,
835 int *mcol,
836 EKKHlink *rlink, EKKHlink *clink,
837 int *maction,
838 int nnetas,
839 int *nnentlp, int *nnentup,
840 int *nsingp)
841 {
842 int *hcoli = fact->xecadr;
843 double *dluval = fact->xeeadr;
844 int *mrstrt = fact->xrsadr;
845 int *hrowi = fact->xeradr;
846 int *mcstrt = fact->xcsadr;
847 int *hinrow = fact->xrnadr;
848 int *hincol = fact->xcnadr;
849 int *hpivro = fact->krpadr;
850 int *hpivco = fact->kcpadr;
851 int nnentl = *nnentlp;
852 int nnentu = *nnentup;
853 int storeZero = fact->ndenuc;
854
855 int mkrs[8];
856 double dpivyy[8];
857
858 /* Local variables */
859 int i, j;
860 double d0, dx;
861 int nz, ndo, krs;
862 int kend, jcol;
863 int irow, iput, jrow, krxs;
864 int mjcol[8];
865 double pivot;
866 int count;
867 int ilast, isort;
868 double dpivx, dsave;
869 double dpivxx[8];
870 double multip;
871 int lstart, ndense, krlast, kcount, idense, ipivot,
872 jdense, kchunk, jpivot;
873
874 const int nrow = fact->nrow;
875
876 int irtcod = 0;
877
878 lstart = nnetas - nnentl + 1;
879
880 /* put list of columns in last HROWI */
881 /* fix row order once for all */
882 ndense = nrow - fact->npivots;
883 iput = ndense + 1;
884 for (i = 1; i <= nrow; ++i) {
885 if (hpivro[i] > 0) {
886 irow = hpivro[i];
887 for (j = 1; j <= nrow; ++j) {
888 --iput;
889 maction[iput] = irow;
890 irow = rlink[irow].suc;
891 if (irow == 0) {
892 break;
893 }
894 }
895 }
896 }
897 if (iput != 1) {
898 ++(*nsingp);
899 } else {
900 /* Use HCOLI just for last row */
901 ilast = maction[1];
902 krlast = mrstrt[ilast];
903 /* put list of columns in last HCOLI */
904 iput = 0;
905 for (i = 1; i <= nrow; ++i) {
906 if (clink[i].pre >= 0) {
907 hcoli[krlast + iput] = i;
908 ++iput;
909 }
910 }
911 if (iput != ndense) {
912 ++(*nsingp);
913 } else {
914 ndo = ndense / 8;
915 /* do most */
916 for (kcount = 1; kcount <= ndo; ++kcount) {
917 idense = ndense;
918 isort = 8;
919 for (count = ndense; count >= ndense - 7; --count) {
920 ipivot = maction[count];
921 krs = mrstrt[ipivot];
922 --isort;
923 mkrs[isort] = krs;
924 }
925 isort = 8;
926 for (count = ndense; count >= ndense - 7; --count) {
927 /* Find a pivot element */
928 --isort;
929 ipivot = maction[count];
930 krs = mkrs[isort];
931 jcol = c_ekkidmx(idense, &dluval[krs]) - 1;
932 pivot = dluval[krs + jcol];
933 --idense;
934 mcol[count] = jcol;
935 mjcol[isort] = mcol[count];
936 dluval[krs + jcol] = dluval[krs + idense];
937 if (fabs(pivot) < fact->zeroTolerance) {
938 pivot = 0.;
939 dpivx = 0.;
940 } else {
941 dpivx = 1. / pivot;
942 }
943 dluval[krs + idense] = pivot;
944 dpivxx[isort] = dpivx;
945 for (j = isort - 1; j >= 0; --j) {
946 krxs = mkrs[j];
947 multip = -dluval[krxs + jcol] * dpivx;
948 dluval[krxs + jcol] = dluval[krxs + idense];
949 /* for moment skip if zero */
950 if (fabs(multip) > fact->zeroTolerance) {
951 for (i = 0; i < idense; ++i) {
952 dluval[krxs + i] += multip * dluval[krs + i];
953 }
954 } else {
955 multip = 0.;
956 }
957 dluval[krxs + idense] = multip;
958 }
959 }
960 /* sort all U in rows already done */
961 for (i = 7; i >= 0; --i) {
962 /* **** this is important bit */
963 krs = mkrs[i];
964 for (j = i - 1; j >= 0; --j) {
965 jcol = mjcol[j];
966 dsave = dluval[krs + jcol];
967 dluval[krs + jcol] = dluval[krs + idense + j];
968 dluval[krs + idense + j] = dsave;
969 }
970 }
971 /* leave IDENSE as it is */
972 if (ndense <= 400) {
973 for (jrow = ndense - 8; jrow >= 1; --jrow) {
974 irow = maction[jrow];
975 krxs = mrstrt[irow];
976 for (j = 7; j >= 0; --j) {
977 jcol = mjcol[j];
978 dsave = dluval[krxs + jcol];
979 dluval[krxs + jcol] = dluval[krxs + idense + j];
980 dluval[krxs + idense + j] = dsave;
981 }
982 for (j = 7; j >= 0; --j) {
983 krs = mkrs[j];
984 jdense = idense + j;
985 dpivx = dpivxx[j];
986 multip = -dluval[krxs + jdense] * dpivx;
987 if (fabs(multip) <= fact->zeroTolerance) {
988 multip = 0.;
989 }
990 dpivyy[j] = multip;
991 dluval[krxs + jdense] = multip;
992 for (i = idense; i < jdense; ++i) {
993 dluval[krxs + i] += multip * dluval[krs + i];
994 }
995 }
996 for (i = 0; i < idense; ++i) {
997 dx = dluval[krxs + i];
998 d0 = dpivyy[0] * dluval[mkrs[0] + i];
999 dx += dpivyy[1] * dluval[mkrs[1] + i];
1000 d0 += dpivyy[2] * dluval[mkrs[2] + i];
1001 dx += dpivyy[3] * dluval[mkrs[3] + i];
1002 d0 += dpivyy[4] * dluval[mkrs[4] + i];
1003 dx += dpivyy[5] * dluval[mkrs[5] + i];
1004 d0 += dpivyy[6] * dluval[mkrs[6] + i];
1005 dx += dpivyy[7] * dluval[mkrs[7] + i];
1006 dluval[krxs + i] = d0 + dx;
1007 }
1008 }
1009 } else {
1010 for (jrow = ndense - 8; jrow >= 1; --jrow) {
1011 irow = maction[jrow];
1012 krxs = mrstrt[irow];
1013 for (j = 7; j >= 0; --j) {
1014 jcol = mjcol[j];
1015 dsave = dluval[krxs + jcol];
1016 dluval[krxs + jcol] = dluval[krxs + idense + j];
1017 dluval[krxs + idense + j] = dsave;
1018 }
1019 for (j = 7; j >= 0; --j) {
1020 krs = mkrs[j];
1021 jdense = idense + j;
1022 dpivx = dpivxx[j];
1023 multip = -dluval[krxs + jdense] * dpivx;
1024 if (fabs(multip) <= fact->zeroTolerance) {
1025 multip = 0.;
1026 }
1027 dluval[krxs + jdense] = multip;
1028 for (i = idense; i < jdense; ++i) {
1029 dluval[krxs + i] += multip * dluval[krs + i];
1030 }
1031 }
1032 }
1033 for (kchunk = 0; kchunk < idense; kchunk += 400) {
1034 kend = CoinMin(idense - 1, kchunk + 399);
1035 for (jrow = ndense - 8; jrow >= 1; --jrow) {
1036 irow = maction[jrow];
1037 krxs = mrstrt[irow];
1038 for (j = 7; j >= 0; --j) {
1039 dpivyy[j] = dluval[krxs + idense + j];
1040 }
1041 for (i = kchunk; i <= kend; ++i) {
1042 dx = dluval[krxs + i];
1043 d0 = dpivyy[0] * dluval[mkrs[0] + i];
1044 dx += dpivyy[1] * dluval[mkrs[1] + i];
1045 d0 += dpivyy[2] * dluval[mkrs[2] + i];
1046 dx += dpivyy[3] * dluval[mkrs[3] + i];
1047 d0 += dpivyy[4] * dluval[mkrs[4] + i];
1048 dx += dpivyy[5] * dluval[mkrs[5] + i];
1049 d0 += dpivyy[6] * dluval[mkrs[6] + i];
1050 dx += dpivyy[7] * dluval[mkrs[7] + i];
1051 dluval[krxs + i] = d0 + dx;
1052 }
1053 }
1054 }
1055 }
1056 /* resort all U in rows already done */
1057 for (i = 7; i >= 0; --i) {
1058 krs = mkrs[i];
1059 for (j = 0; j < i; ++j) {
1060 jcol = mjcol[j];
1061 dsave = dluval[krs + jcol];
1062 dluval[krs + jcol] = dluval[krs + idense + j];
1063 dluval[krs + idense + j] = dsave;
1064 }
1065 }
1066 ndense += -8;
1067 }
1068 idense = ndense;
1069 /* do remainder */
1070 for (count = ndense; count >= 1; --count) {
1071 /* Find a pivot element */
1072 ipivot = maction[count];
1073 krs = mrstrt[ipivot];
1074 jcol = c_ekkidmx(idense, &dluval[krs]) - 1;
1075 pivot = dluval[krs + jcol];
1076 --idense;
1077 mcol[count] = jcol;
1078 dluval[krs + jcol] = dluval[krs + idense];
1079 if (fabs(pivot) < fact->zeroTolerance) {
1080 dluval[krs + idense] = 0.;
1081 } else {
1082 dpivx = 1. / pivot;
1083 dluval[krs + idense] = pivot;
1084 for (jrow = idense; jrow >= 1; --jrow) {
1085 irow = maction[jrow];
1086 krxs = mrstrt[irow];
1087 multip = -dluval[krxs + jcol] * dpivx;
1088 dluval[krxs + jcol] = dluval[krxs + idense];
1089 /* for moment skip if zero */
1090 if (fabs(multip) > fact->zeroTolerance) {
1091 dluval[krxs + idense] = multip;
1092 for (i = 0; i < idense; ++i) {
1093 dluval[krxs + i] += multip * dluval[krs + i];
1094 }
1095 } else {
1096 dluval[krxs + idense] = 0.;
1097 }
1098 }
1099 }
1100 }
1101 /* now create in form for OSL */
1102 ndense = nrow - fact->npivots;
1103 idense = ndense;
1104 for (count = ndense; count >= 1; --count) {
1105 /* Find a pivot element */
1106 ipivot = maction[count];
1107 krs = mrstrt[ipivot];
1108 --idense;
1109 jcol = mcol[count];
1110 jpivot = hcoli[krlast + jcol];
1111 ++fact->npivots;
1112 pivot = dluval[krs + idense];
1113 if (pivot == 0.) {
1114 hinrow[ipivot] = 0;
1115 rlink[ipivot].pre = -nrow - 1;
1116 ++(*nsingp);
1117 irtcod = 10;
1118 } else {
1119 rlink[ipivot].pre = -fact->npivots;
1120 clink[jpivot].pre = -fact->npivots;
1121 hincol[jpivot] = 0;
1122 ++fact->xnetal;
1123 mcstrt[fact->xnetal] = lstart - 1;
1124 hpivco[fact->xnetal] = ipivot;
1125 for (jrow = idense; jrow >= 1; --jrow) {
1126 irow = maction[jrow];
1127 krxs = mrstrt[irow];
1128 multip = dluval[krxs + idense];
1129 /* for moment skip if zero */
1130 if (multip != 0. || storeZero) {
1131 /* Store elementary row transformation */
1132 ++nnentl;
1133 --nnentu;
1134 --lstart;
1135 dluval[lstart] = multip;
1136 hrowi[lstart] = SHIFT_INDEX(irow);
1137 }
1138 }
1139 hcoli[krlast + jcol] = hcoli[krlast + idense];
1140 /* update pivot row and last row HCOLI */
1141 dluval[krs + idense] = dluval[krs];
1142 hcoli[krlast + idense] = hcoli[krlast];
1143 nz = 1;
1144 dluval[krs] = pivot;
1145 hcoli[krs] = jpivot;
1146 if (!storeZero) {
1147 for (i = 1; i <= idense; ++i) {
1148 if (fabs(dluval[krs + i]) > fact->zeroTolerance) {
1149 ++nz;
1150 hcoli[krs + nz - 1] = hcoli[krlast + i];
1151 dluval[krs + nz - 1] = dluval[krs + i];
1152 }
1153 }
1154 hinrow[ipivot] = nz;
1155 } else {
1156 for (i = 1; i <= idense; ++i) {
1157 ++nz;
1158 hcoli[krs + nz - 1] = hcoli[krlast + i];
1159 dluval[krs + nz - 1] = dluval[krs + i];
1160 }
1161 hinrow[ipivot] = nz;
1162 }
1163 }
1164 }
1165 }
1166 }
1167
1168 *nnentlp = nnentl;
1169 *nnentup = nnentu;
1170
1171 return (irtcod);
1172 } /* c_ekkcmfd */
1173 /* ***C_EKKCMFC */
1174
1175 /*
1176 * Generate a variant of c_ekkcmfc that uses an maction array of type
1177 * int rather than short.
1178 */
1179 #undef MACTION_T
1180 #define C_EKKCMFY
1181 #define COIN_OSL_CMFC
1182 #define MACTION_T int
c_ekkcmfy(EKKfactinfo * fact,EKKHlink * rlink,EKKHlink * clink,EKKHlink * mwork,void * maction_void,int nnetas,int * nsingp,int * xrejctp,int * xnewrop,int xnewco,int * ncompactionsp)1183 int c_ekkcmfy(EKKfactinfo *fact,
1184 EKKHlink *rlink, EKKHlink *clink,
1185 EKKHlink *mwork, void *maction_void,
1186 int nnetas,
1187 int *nsingp, int *xrejctp,
1188 int *xnewrop, int xnewco,
1189 int *ncompactionsp)
1190
1191 #include "CoinOslC.h"
1192 #undef COIN_OSL_CMFC
1193 #undef C_EKKCMFY
1194 #undef MACTION_T
1195 int c_ekkford(const EKKfactinfo *fact, const int *hinrow, const int *hincol,
1196 int *hpivro, int *hpivco,
1197 EKKHlink *rlink, EKKHlink *clink)
1198 {
1199 int i, iri, nzi;
1200 const int nrow = fact->nrow;
1201 int nsing = 0;
1202
1203 /* Uwe H. Suhl, August 1986 */
1204 /* Builds linked lists of rows and cols of nucleus for efficient */
1205 /* pivot searching. */
1206
1207 memset(hpivro + 1, 0, nrow * sizeof(int));
1208 memset(hpivco + 1, 0, nrow * sizeof(int));
1209 for (i = 1; i <= nrow; ++i) {
1210 //hpivro[i] = 0;
1211 //hpivco[i] = 0;
1212 assert(rlink[i].suc == 0);
1213 assert(clink[i].suc == 0);
1214 }
1215
1216 /* Generate double linked list of rows having equal numbers of */
1217 /* nonzeros in each row. Skip pivotal rows. */
1218 for (i = 1; i <= nrow; ++i) {
1219 if (!(rlink[i].pre < 0)) {
1220 nzi = hinrow[i];
1221 if (nzi <= 0) {
1222 ++nsing;
1223 rlink[i].pre = -nrow - 1;
1224 } else {
1225 iri = hpivro[nzi];
1226 hpivro[nzi] = i;
1227 rlink[i].suc = iri;
1228 rlink[i].pre = 0;
1229 if (iri != 0) {
1230 rlink[iri].pre = i;
1231 }
1232 }
1233 }
1234 }
1235
1236 /* Generate double linked list of cols having equal numbers of */
1237 /* nonzeros in each col. Skip pivotal cols. */
1238 for (i = 1; i <= nrow; ++i) {
1239 if (!(clink[i].pre < 0)) {
1240 nzi = hincol[i];
1241 if (nzi <= 0) {
1242 ++nsing;
1243 clink[i].pre = -nrow - 1;
1244 } else {
1245 iri = hpivco[nzi];
1246 hpivco[nzi] = i;
1247 clink[i].suc = iri;
1248 clink[i].pre = 0;
1249 if (iri != 0) {
1250 clink[iri].pre = i;
1251 }
1252 }
1253 }
1254 }
1255
1256 return (nsing);
1257 } /* c_ekkford */
1258
1259 /* c version of OSL from 36100 */
1260
1261 /* Assumes that a basis exists in correct form */
1262 /* Calls Uwe's routines (approximately) */
1263 /* Then if OK shuffles U into column order */
1264 /* Return codes: */
1265
1266 /* 0: everything ok */
1267 /* 1: everything ok but performance would be better if more space */
1268 /* would be make available */
1269 /* 4: growth rate of element in U too big */
1270 /* 5: not enough space in row file */
1271 /* 6: not enough space in column file */
1272 /* 7: pivot too small - col sing */
1273 /* 8: pivot too small - row sing */
1274 /* 10: matrix is singular */
1275
1276 /* I suspect c_ekklfct never returns 1 */
1277 /*
1278 * layout of data
1279 *
1280 * dluval/hcoli: (L^-1)B - hole - L factors
1281 *
1282 * The L factors are written from high to low, starting from nnetas.
1283 * There are nnentl factors in L. lstart the next entry to use for the
1284 * L factors. Eventually, (L^-1)B turns into U.
1285
1286 * The ninbas coefficients of matrix B are originally in the start of
1287 * dluval/hcoli. As L transforms it, rows may have to be expanded.
1288 * If there is room, they are copied to the start of the hole,
1289 * otherwise the first part of this area is compacted, and hopefully
1290 * there is then room.
1291 * There are nnentu coefficients in (L^-1)B.
1292 * nnentu + nnentl >= ninbas.
1293 * nnentu + nnentl == ninbas if there has been no fill-in.
1294 * nnentu is decreased when the pivot eliminates elements
1295 * (in which case there is a corresponding increase in nnentl),
1296 * and if pivoting happens to cancel out factors (in which case
1297 * there is no corresponding increase in L).
1298 * nnentu is increased if there is fill-in (no decrease in L).
1299 * If nnentu + nnentl >= nnetas, then we've run out of room.
1300 * It is not the case that the elements of (L^-1)B are all in the
1301 * first nnentu positions of dluval/hcoli, but that is of course
1302 * the lower bound on the number of positions needed to store it.
1303 * nuspik is roughly the sum of the row lengths of the rows that were pivoted
1304 * out. singleton rows in c_ekktria do not change nuspik, but
1305 * c_ekkrsin does increment it for each singleton row.
1306 * That is, there are nuspik elements that in the upper part of (L^-1)B,
1307 * and (nnentu - nuspik) elements left in B.
1308 */
1309 /*
1310 * As part of factorization, we test candidate pivots for numerical
1311 * stability; if the largest element in a row/col is much larger than
1312 * the smallest, this generally causes problems. To easily determine
1313 * what the largest element is, we ensure that it is always in front.
1314 * This establishes this property; later on we take steps to preserve it.
1315 */
c_ekkmltf(const EKKfactinfo * fact,double * dluval,int * hcoli,const int * mrstrt,const int * hinrow,const EKKHlink * rlink)1316 static void c_ekkmltf(const EKKfactinfo *fact, double *dluval, int *hcoli,
1317 const int *mrstrt, const int *hinrow,
1318 const EKKHlink *rlink)
1319 {
1320 #if 0
1321 int *hcoli = fact->xecadr;
1322 double *dluval = fact->xeeadr;
1323 double *dvalpv = fact->kw3adr;
1324 int *mrstrt = fact->xrsadr;
1325 int *hrowi = fact->xeradr;
1326 int *mcstrt = fact->xcsadr;
1327 int *hinrow = fact->xrnadr;
1328 int *hincol = fact->xcnadr;
1329 int *hpivro = fact->krpadr;
1330 int *hpivco = fact->kcpadr;
1331 #endif
1332 int i, j, k;
1333 int koff = -1;
1334 const int nrow = fact->nrow;
1335
1336 for (i = 1; i <= nrow; ++i) {
1337 /* ignore rows that have already been pivoted */
1338 /* if it is a singleton row, the property trivially holds */
1339 if (!(rlink[i].pre < 0 || hinrow[i] <= 1)) {
1340 const int krs = mrstrt[i];
1341 const int kre = krs + hinrow[i] - 1;
1342
1343 double maxaij = 0.f;
1344
1345 /* this assumes that at least one of the dluvals is non-zero. */
1346 for (k = krs; k <= kre; ++k) {
1347 if (!(fabs(dluval[k]) <= maxaij)) {
1348 maxaij = fabs(dluval[k]);
1349 koff = k;
1350 }
1351 }
1352 assert(koff > 0);
1353 maxaij = dluval[koff];
1354 j = hcoli[koff];
1355
1356 dluval[koff] = dluval[krs];
1357 hcoli[koff] = hcoli[krs];
1358
1359 dluval[krs] = maxaij;
1360 hcoli[krs] = j;
1361 }
1362 }
1363 } /* c_ekkmltf */
c_ekklfct(register EKKfactinfo * fact)1364 int c_ekklfct(register EKKfactinfo *fact)
1365 {
1366 const int nrow = fact->nrow;
1367 int ninbas = fact->xcsadr[nrow + 1] - 1;
1368 int ifvsol = fact->ifvsol;
1369 int *hcoli = fact->xecadr;
1370 double *dluval = fact->xeeadr;
1371 int *mrstrt = fact->xrsadr;
1372 int *hrowi = fact->xeradr;
1373 int *mcstrt = fact->xcsadr;
1374 int *hinrow = fact->xrnadr;
1375 int *hincol = fact->xcnadr;
1376 int *hpivro = fact->krpadr;
1377 int *hpivco = fact->kcpadr;
1378
1379 EKKHlink *rlink = fact->kp1adr;
1380 EKKHlink *clink = fact->kp2adr;
1381 EKKHlink *mwork = (reinterpret_cast< EKKHlink * >(fact->kw1adr)) - 1;
1382
1383 int nsing, kdnspt, xnewro, xnewco;
1384 int i;
1385 int xrejct;
1386 int irtcod;
1387 const int nnetas = fact->nnetas;
1388
1389 int ncompactions;
1390 double save_drtpiv = fact->drtpiv;
1391 double save_zpivlu = fact->zpivlu;
1392 if (ifvsol > 0 && fact->invok < 0) {
1393 fact->zpivlu = CoinMin(0.9, fact->zpivlu * 10.);
1394 fact->drtpiv = 1.0e-8;
1395 }
1396
1397 rlink--;
1398 clink--;
1399
1400 /* Function Body */
1401 hcoli[nnetas] = 1;
1402 hrowi[nnetas] = 1;
1403 dluval[nnetas] = 0.0;
1404 /* set amount of work */
1405 xrejct = 0;
1406 nsing = 0;
1407 kdnspt = nnetas + 1;
1408 fact->ndenuc = 0;
1409 /* Triangularize */
1410 irtcod = c_ekktria(fact, rlink, clink,
1411 &nsing,
1412 &xnewco, &xnewro,
1413 &ncompactions, ninbas);
1414 fact->nnentl = ninbas - fact->nnentu;
1415
1416 if (irtcod < 0) {
1417 /* no space or system error */
1418 goto L8000;
1419 }
1420
1421 if (irtcod != 0 && fact->invok >= 0) {
1422 goto L8500; /* 7 or 8 - pivot too small */
1423 }
1424 #if 0
1425 /* is this necessary ? */
1426 lstart = nnetas - fact->nnentl + 1;
1427 for (i = lstart; i <= nnetas; ++i) {
1428 hrowi[i] = (hcoli[i] << 3);
1429 }
1430 #endif
1431
1432 /* See if finished */
1433 if (!(fact->npivots >= nrow)) {
1434 int nsing1;
1435
1436 /* No - do nucleus */
1437
1438 nsing1 = c_ekkford(fact, hinrow, hincol, hpivro, hpivco, rlink, clink);
1439 nsing += nsing1;
1440 if (nsing1 != 0 && fact->invok >= 0) {
1441 irtcod = 7;
1442 goto L8500;
1443 }
1444 c_ekkmltf(fact, dluval, hcoli, mrstrt, hinrow, rlink);
1445
1446 {
1447 bool callcmfy = false;
1448
1449 if (nrow > 32767) {
1450 int count = 0;
1451 for (i = 1; i <= nrow; ++i) {
1452 count = CoinMax(count, hinrow[i]);
1453 }
1454 if (count + nrow - fact->npivots > 32767) {
1455 /* will have to use I*4 version of CMFC */
1456 /* no changes to pointer params */
1457 callcmfy = true;
1458 }
1459 }
1460
1461 irtcod = (callcmfy ? c_ekkcmfy : c_ekkcmfc)(fact,
1462 rlink, clink,
1463 mwork, &mwork[nrow + 1],
1464 nnetas,
1465 &nsing, &xrejct,
1466 &xnewro, xnewco,
1467 &ncompactions);
1468
1469 /* irtcod one of 0,-5,7,10 */
1470 }
1471
1472 if (irtcod < 0) {
1473 goto L8000;
1474 }
1475 kdnspt = nnetas - fact->nnentl;
1476 }
1477
1478 /* return if error */
1479 if (nsing > 0 || irtcod == 10) {
1480 irtcod = 99;
1481 }
1482 /* irtcod one of 0,7,99 */
1483
1484 if (irtcod != 0) {
1485 goto L8500;
1486 }
1487 ++fact->xnetal;
1488 mcstrt[fact->xnetal] = nnetas - fact->nnentl;
1489
1490 /* give message if tight on memory */
1491 if (ncompactions > 2) {
1492 if (1) {
1493 int etasize = CoinMax(4 * fact->nnentu + (nnetas - fact->nnentl) + 1000, fact->eta_size);
1494 fact->eta_size = CoinMin(static_cast< int >(1.2 * fact->eta_size), etasize);
1495 if (fact->maxNNetas > 0 && fact->eta_size > fact->maxNNetas) {
1496 fact->eta_size = fact->maxNNetas;
1497 }
1498 } /* endif */
1499 }
1500 /* Shuffle U and multiply L by 8 (if assembler) */
1501 {
1502 int jrtcod = c_ekkshff(fact, clink, rlink,
1503 xnewro);
1504
1505 /* nR_etas is the number of R transforms;
1506 * it is incremented only in c_ekketsj.
1507 */
1508 fact->nR_etas = 0;
1509 /*fact->R_etas_start = mcstrt+nrow+fact->nnentl+3;*/
1510 fact->R_etas_start[1] = /*kdnspt - 1*/ 0; /* magic */
1511 fact->R_etas_index = &fact->xeradr[kdnspt - 1];
1512 fact->R_etas_element = &fact->xeeadr[kdnspt - 1];
1513
1514 if (jrtcod != 0) {
1515 irtcod = jrtcod;
1516 /* irtcod == 2 */
1517 }
1518 }
1519 goto L8500;
1520
1521 /* Fatal error */
1522 L8000:
1523
1524 if (1) {
1525 if (fact->maxNNetas != fact->eta_size && nnetas) {
1526 /* return and get more space */
1527
1528 /* double eta_size, unless that exceeds max (if there is one) */
1529 fact->eta_size = fact->eta_size << 1;
1530 if (fact->maxNNetas > 0 && fact->eta_size > fact->maxNNetas) {
1531 fact->eta_size = fact->maxNNetas;
1532 }
1533 return (5);
1534 }
1535 }
1536 /*c_ekkmesg_no_i1(121, -irtcod);*/
1537 irtcod = 3;
1538
1539 L8500:
1540 /* restore pivot tolerance */
1541 fact->drtpiv = save_drtpiv;
1542 fact->zpivlu = save_zpivlu;
1543 #ifndef NDEBUG
1544 if (fact->rows_ok) {
1545 int *hinrow = fact->xrnadr;
1546 if (!fact->xe2adr) {
1547 for (int i = 1; i <= fact->nrow; i++) {
1548 assert(hinrow[i] >= 0 && hinrow[i] <= fact->nrow);
1549 }
1550 }
1551 }
1552 #endif
1553 return (irtcod);
1554 } /* c_ekklfct */
1555
1556 /*
1557 summary of return codes
1558
1559 c_ekktria:
1560 7 small pivot
1561 -5 no memory
1562
1563 c_ekkcsin:
1564 returns true if small pivot
1565
1566 c_ekkrsin:
1567 -5 no memory
1568 7 small pivot
1569
1570 c_ekkfpvt:
1571 10: no pivots found (singular)
1572
1573 c_ekkcmfd:
1574 10: zero pivot (not just small)
1575
1576 c_ekkcmfc:
1577 -5: no memory
1578 any non-zero code from c_ekkcsin, c_ekkrsin, c_ekkfpvt, c_ekkprpv, c_ekkcmfd
1579
1580 c_ekkshff:
1581 2: singular
1582
1583 c_ekklfct:
1584 any positive code from c_ekktria, c_ekkcmfc, c_ekkshff (2,7,10)
1585 *except* 10, which is changed to 99.
1586 all negative return codes are changed to 5 or 3
1587 (5 == ran out of memory but could get more,
1588 3 == ran out of memory, no luck)
1589 so: 2,3,5,7,99
1590
1591 c_ekklfct1:
1592 1: c_ekksmem_invert failed
1593 2: c_ekkslcf/c_ekkslct ran out of room
1594 any return code from c_ekklfct, except 2 and 5
1595 */
1596
c_ekkrowq(int * hrow,int * hcol,double * dels,int * mrstrt,const int * hinrow,int nnrow,int ninbas)1597 void c_ekkrowq(int *hrow, int *hcol, double *dels,
1598 int *mrstrt,
1599 const int *hinrow, int nnrow, int ninbas)
1600 {
1601 int i, k, iak, jak;
1602 double daik;
1603 int iloc;
1604 double dsave;
1605 int isave, jsave;
1606
1607 /* Order matrix rowwise using MRSTRT, DELS, HCOL */
1608
1609 k = 1;
1610 /* POSITION AFTER END OF ROW */
1611 for (i = 1; i <= nnrow; ++i) {
1612 k += hinrow[i];
1613 mrstrt[i] = k;
1614 }
1615
1616 for (k = ninbas; k >= 1; --k) {
1617 iak = hrow[k];
1618 if (iak != 0) {
1619 daik = dels[k];
1620 jak = hcol[k];
1621 hrow[k] = 0;
1622 while (1) {
1623 --mrstrt[iak];
1624
1625 iloc = mrstrt[iak];
1626
1627 dsave = dels[iloc];
1628 isave = hrow[iloc];
1629 jsave = hcol[iloc];
1630 dels[iloc] = daik;
1631 hrow[iloc] = 0;
1632 hcol[iloc] = jak;
1633
1634 if (isave == 0)
1635 break;
1636 daik = dsave;
1637 iak = isave;
1638 jak = jsave;
1639 }
1640 }
1641 }
1642 } /* c_ekkrowq */
1643
c_ekkrwco(const EKKfactinfo * fact,double * dluval,int * hcoli,int * mrstrt,int * hinrow,int xnewro)1644 int c_ekkrwco(const EKKfactinfo *fact, double *dluval,
1645 int *hcoli, int *mrstrt, int *hinrow, int xnewro)
1646 {
1647 int i, k, nz, kold;
1648 int kstart;
1649 const int nrow = fact->nrow;
1650
1651 for (i = 1; i <= nrow; ++i) {
1652 nz = hinrow[i];
1653 if (0 < nz) {
1654 /* save the last column entry of row i in hinrow */
1655 /* and replace that entry with -i */
1656 k = mrstrt[i] + nz - 1;
1657 hinrow[i] = hcoli[k];
1658 hcoli[k] = -i;
1659 }
1660 }
1661
1662 kstart = 0;
1663 kold = 0;
1664 for (k = 1; k <= xnewro; ++k) {
1665 if (hcoli[k] != 0) {
1666 ++kstart;
1667
1668 /* if this is the last entry for the row... */
1669 if (hcoli[k] < 0) {
1670 /* restore the entry */
1671 i = -hcoli[k];
1672 hcoli[k] = hinrow[i];
1673
1674 /* update mrstart and hinrow */
1675 /* ACTUALLY, hinrow should already be accurate */
1676 mrstrt[i] = kold + 1;
1677 hinrow[i] = kstart - kold;
1678 kold = kstart;
1679 }
1680
1681 /* move the entry */
1682 dluval[kstart] = dluval[k];
1683 hcoli[kstart] = hcoli[k];
1684 }
1685 }
1686
1687 return (kstart);
1688 } /* c_ekkrwco */
1689
c_ekkrwcs(const EKKfactinfo * fact,double * dluval,int * hcoli,int * mrstrt,const int * hinrow,const EKKHlink * mwork,int nfirst)1690 int c_ekkrwcs(const EKKfactinfo *fact, double *dluval, int *hcoli, int *mrstrt,
1691 const int *hinrow, const EKKHlink *mwork,
1692 int nfirst)
1693 {
1694 #if 0
1695 int *hcoli = fact->xecadr;
1696 double *dluval = fact->xeeadr;
1697 double *dvalpv = fact->kw3adr;
1698 int *mrstrt = fact->xrsadr;
1699 int *hrowi = fact->xeradr;
1700 int *mcstrt = fact->xcsadr;
1701 int *hinrow = fact->xrnadr;
1702 int *hincol = fact->xcnadr;
1703 int *hpivro = fact->krpadr;
1704 int *hpivco = fact->kcpadr;
1705 #endif
1706 int i, k, k1, k2, nz;
1707 int irow, iput;
1708 const int nrow = fact->nrow;
1709
1710 /* Compress row file */
1711
1712 iput = 1;
1713 irow = nfirst;
1714 for (i = 1; i <= nrow; ++i) {
1715 nz = hinrow[irow];
1716 k1 = mrstrt[irow];
1717 if (k1 != iput) {
1718 mrstrt[irow] = iput;
1719 k2 = k1 + nz - 1;
1720 for (k = k1; k <= k2; ++k) {
1721 dluval[iput] = dluval[k];
1722 hcoli[iput] = hcoli[k];
1723 ++iput;
1724 }
1725 } else {
1726 iput += nz;
1727 }
1728 irow = mwork[irow].suc;
1729 }
1730
1731 return (iput);
1732 } /* c_ekkrwcs */
c_ekkrwct(const EKKfactinfo * fact,double * dluval,int * hcoli,int * mrstrt,const int * hinrow,const EKKHlink * mwork,const EKKHlink * rlink,const short * msort,double * dsort,int nlast,int xnewro)1733 void c_ekkrwct(const EKKfactinfo *fact, double *dluval, int *hcoli, int *mrstrt,
1734 const int *hinrow, const EKKHlink *mwork,
1735 const EKKHlink *rlink,
1736 const short *msort, double *dsort,
1737 int nlast, int xnewro)
1738 {
1739 #if 0
1740 int *hcoli = fact->xecadr;
1741 double *dluval = fact->xeeadr;
1742 double *dvalpv = fact->kw3adr;
1743 int *mrstrt = fact->xrsadr;
1744 int *hrowi = fact->xeradr;
1745 int *mcstrt = fact->xcsadr;
1746 int *hinrow = fact->xrnadr;
1747 int *hincol = fact->xcnadr;
1748 int *hpivro = fact->krpadr;
1749 int *hpivco = fact->kcpadr;
1750 #endif
1751 int i, k, k1, nz, icol;
1752 int kmax;
1753 int irow, iput;
1754 int ilook;
1755 const int nrow = fact->nrow;
1756
1757 iput = xnewro;
1758 irow = nlast;
1759 kmax = nrow - fact->npivots;
1760 for (i = 1; i <= nrow; ++i) {
1761 nz = hinrow[irow];
1762 k1 = mrstrt[irow] - 1;
1763 if (rlink[irow].pre < 0) {
1764 /* pivoted on already */
1765 iput -= nz;
1766 if (k1 != iput) {
1767 mrstrt[irow] = iput + 1;
1768 for (k = nz; k >= 1; --k) {
1769 dluval[iput + k] = dluval[k1 + k];
1770 hcoli[iput + k] = hcoli[k1 + k];
1771 }
1772 }
1773 } else {
1774 /* not pivoted - going dense */
1775 iput -= kmax;
1776 mrstrt[irow] = iput + 1;
1777 c_ekkdzero(kmax, &dsort[1]);
1778 for (k = 1; k <= nz; ++k) {
1779 icol = hcoli[k1 + k];
1780 ilook = msort[icol];
1781 dsort[ilook] = dluval[k1 + k];
1782 }
1783 c_ekkdcpy(kmax,
1784 (dsort + 1), (dluval + iput + 1));
1785 }
1786 irow = mwork[irow].pre;
1787 }
1788 } /* c_ekkrwct */
1789 /* takes Uwe's modern structures and puts them back 20 years */
c_ekkshff(EKKfactinfo * fact,EKKHlink * clink,EKKHlink * rlink,int xnewro)1790 int c_ekkshff(EKKfactinfo *fact,
1791 EKKHlink *clink, EKKHlink *rlink,
1792 int xnewro)
1793 {
1794 int *hpivro = fact->krpadr;
1795
1796 int i, j;
1797 int nbas, icol;
1798 int ipiv;
1799 const int nrow = fact->nrow;
1800 int nsing;
1801
1802 for (i = 1; i <= nrow; ++i) {
1803 j = -rlink[i].pre;
1804 rlink[i].pre = j;
1805 if (j > 0 && j <= nrow) {
1806 hpivro[j] = i;
1807 }
1808 j = -clink[i].pre;
1809 clink[i].pre = j;
1810 }
1811 /* hpivro[j] is now (hopefully) the row that was pivoted on step j */
1812 /* rlink[i].pre is the step in which row i was pivoted */
1813
1814 nbas = 0;
1815 nsing = 0;
1816 /* Decide if permutation wanted */
1817 fact->first_dense = nrow - fact->ndenuc + 1 + 1;
1818 fact->last_dense = nrow;
1819
1820 /* rlink[].suc is dead at this point */
1821
1822 /*
1823 * replace the the basis index
1824 * with the pivot (or permuted) index generated by factorization.
1825 * This eventually goes into mpermu.
1826 */
1827 for (icol = 1; icol <= nrow; ++icol) {
1828 int ibasis = icol;
1829 ipiv = clink[ibasis].pre;
1830
1831 if (0 < ipiv && ipiv <= nrow) {
1832 rlink[ibasis].suc = ipiv;
1833 ++nbas;
1834 }
1835 }
1836
1837 nsing = nrow - nbas;
1838 if (nsing > 0) {
1839 abort();
1840 }
1841
1842 /* if we reach here, then rlink[1..nrow].suc == clink[1..nrow].pre */
1843
1844 /* switch off sparse update if any dense section */
1845 {
1846 const int notMuchRoom = (fact->nnentu + xnewro + 10 > fact->nnetas - fact->nnentl);
1847
1848 /* must be same as in c_ekkshfv */
1849 if (fact->ndenuc || notMuchRoom || nrow < C_EKK_GO_SPARSE) {
1850 #if PRINT_DEBUG
1851 if (fact->if_sparse_update) {
1852 printf("**** Switching off sparse update - dense - c_ekkshff\n");
1853 }
1854 #endif
1855 fact->if_sparse_update = 0;
1856 }
1857 }
1858
1859 /* hpivro[1..nrow] is not read by c_ekkshfv */
1860 c_ekkshfv(fact,
1861 rlink, clink,
1862 xnewro);
1863
1864 return (0);
1865 } /* c_ekkshff */
1866 /* sorts on indices dragging elements with */
c_ekk_sort2(int * key,double * array2,int number)1867 static void c_ekk_sort2(int *key, double *array2, int number)
1868 {
1869 int minsize = 10;
1870 int n = number;
1871 int sp;
1872 int *v = key;
1873 int *m, t;
1874 int *ls[32], *rs[32];
1875 int *l, *r, c;
1876 double it;
1877 int j;
1878 /*check already sorted */
1879 #ifndef LONG_MAX
1880 #define LONG_MAX 0x7fffffff;
1881 #endif
1882 int last = -LONG_MAX;
1883 for (j = 0; j < number; j++) {
1884 if (key[j] >= last) {
1885 last = key[j];
1886 } else {
1887 break;
1888 } /* endif */
1889 } /* endfor */
1890 if (j == number) {
1891 return;
1892 } /* endif */
1893 sp = 0;
1894 ls[sp] = v;
1895 rs[sp] = v + (n - 1);
1896 while (sp >= 0) {
1897 if (rs[sp] - ls[sp] > minsize) {
1898 l = ls[sp];
1899 r = rs[sp];
1900 m = l + (r - l) / 2;
1901 if (*l > *m) {
1902 t = *l;
1903 *l = *m;
1904 *m = t;
1905 it = array2[l - v];
1906 array2[l - v] = array2[m - v];
1907 array2[m - v] = it;
1908 }
1909 if (*m > *r) {
1910 t = *m;
1911 *m = *r;
1912 *r = t;
1913 it = array2[m - v];
1914 array2[m - v] = array2[r - v];
1915 array2[r - v] = it;
1916 if (*l > *m) {
1917 t = *l;
1918 *l = *m;
1919 *m = t;
1920 it = array2[l - v];
1921 array2[l - v] = array2[m - v];
1922 array2[m - v] = it;
1923 }
1924 }
1925 c = *m;
1926 while (r - l > 1) {
1927 while (*(++l) < c)
1928 ;
1929 while (*(--r) > c)
1930 ;
1931 t = *l;
1932 *l = *r;
1933 *r = t;
1934 it = array2[l - v];
1935 array2[l - v] = array2[r - v];
1936 array2[r - v] = it;
1937 }
1938 l = r - 1;
1939 if (l < m) {
1940 ls[sp + 1] = ls[sp];
1941 rs[sp + 1] = l;
1942 ls[sp] = r;
1943 } else {
1944 ls[sp + 1] = r;
1945 rs[sp + 1] = rs[sp];
1946 rs[sp] = l;
1947 }
1948 sp++;
1949 } else
1950 sp--;
1951 }
1952 for (l = v, m = v + (n - 1); l < m; l++) {
1953 if (*l > *(l + 1)) {
1954 c = *(l + 1);
1955 it = array2[(l - v) + 1];
1956 for (r = l; r >= v && *r > c; r--) {
1957 *(r + 1) = *r;
1958 array2[(r - v) + 1] = array2[(r - v)];
1959 }
1960 *(r + 1) = c;
1961 array2[(r - v) + 1] = it;
1962 }
1963 }
1964 }
1965 /* For each row compute reciprocal of pivot element and take out of */
1966 /* Also use HLINK(1 to permute column numbers */
1967 /* and HPIVRO to permute row numbers */
1968 /* Sort into column order as was stored by row */
1969 /* If Assembler then shift row numbers in L by 3 */
1970 /* Put column numbers in U for L-U update */
1971 /* and multiply U elements by - reciprocal of pivot element */
1972 /* and set up backward pointers for pivot rows */
c_ekkshfv(EKKfactinfo * fact,EKKHlink * rlink,EKKHlink * clink,int xnewro)1973 void c_ekkshfv(EKKfactinfo *fact,
1974 EKKHlink *rlink, EKKHlink *clink,
1975 int xnewro)
1976 {
1977 int *hcoli = fact->xecadr;
1978 double *dluval = fact->xeeadr;
1979 double *dvalpv = fact->kw3adr;
1980 int *mrstrt = fact->xrsadr;
1981 int *hrowi = fact->xeradr;
1982 int *mcstrt = fact->xcsadr;
1983 int *hinrow = fact->xrnadr;
1984 int *hincol = fact->xcnadr;
1985 int *hpivro = fact->krpadr;
1986 int *hpivco = fact->kcpadr;
1987 double *dpermu = fact->kadrpm;
1988 double *de2val = fact->xe2adr ? fact->xe2adr - 1 : 0;
1989 int nnentu = fact->nnentu;
1990 int xnetal = fact->xnetal;
1991
1992 int numberSlacks; /* numberSlacks not read */
1993
1994 int i, j, k, kk, nel;
1995 int nroom;
1996 bool need_more_space;
1997 int ndenuc = fact->ndenuc;
1998 int if_sparse_update = fact->if_sparse_update;
1999 int nnentl = fact->nnentl;
2000 int nnetas = fact->nnetas;
2001
2002 int *ihlink = (reinterpret_cast< int * >(clink)) + 1; /* can't use rlink for simple loop below */
2003
2004 const int nrow = fact->nrow;
2005 const int maxinv = fact->maxinv;
2006
2007 /* this is not just a temporary - c_ekkbtrn etc use this */
2008 int *mpermu = (reinterpret_cast< int * >(dpermu + nrow)) + 1;
2009
2010 int *temp = ihlink + nrow;
2011 int *temp2 = temp + nrow;
2012 const int notMuchRoom = (nnentu + xnewro + 10 > nnetas - nnentl);
2013
2014 /* compress hlink and make simpler */
2015 for (i = 1; i <= nrow; ++i) {
2016 mpermu[i] = rlink[i].pre;
2017 ihlink[i] = rlink[i].suc;
2018 }
2019 /* mpermu[i] == the step in which row i was pivoted */
2020 /* ihlink[i] == the step in which col i was pivoted */
2021
2022 /* must be same as in c_ekkshff */
2023 if (fact->ndenuc || notMuchRoom || nrow < C_EKK_GO_SPARSE) {
2024 int ninbas;
2025
2026 /* CHANGE COLUMN NUMBERS AND FILL IN RECIPROCALS */
2027 /* ALSO RECOMPUTE NUMBER IN COLUMN */
2028 /* initialize with a fake pivot in each column */
2029 c_ekkscpy_0_1(nrow, 1, &hincol[1]);
2030
2031 if (notMuchRoom) {
2032 fact->eta_size = static_cast< int >(1.05 * fact->eta_size);
2033
2034 /* eta_size can be no larger than maxNNetas */
2035 if (fact->maxNNetas > 0 && fact->eta_size > fact->maxNNetas) {
2036 fact->eta_size = fact->maxNNetas;
2037 }
2038 } /* endif */
2039
2040 /* For each row compute reciprocal of pivot element and take out of U */
2041 /* Also use ihlink to permute column numbers */
2042 /* the rows are not stored compactly or in order,
2043 * so we have to find out where the last one is stored */
2044 ninbas = 0;
2045 for (i = 1; i <= nrow; ++i) {
2046 int jpiv = mpermu[i];
2047 int nin = hinrow[i];
2048 int krs = mrstrt[i];
2049 int kre = krs + nin;
2050
2051 temp[jpiv] = krs;
2052 temp2[jpiv] = nin;
2053
2054 ninbas = CoinMax(kre, ninbas);
2055
2056 /* c_ekktria etc ensure that the first row entry is the pivot */
2057 dvalpv[jpiv] = 1. / dluval[krs];
2058 hcoli[krs] = 0; /* probably needed for c_ekkrowq */
2059 /* room for the pivot has already been allocated, so hincol ok */
2060
2061 for (kk = krs + 1; kk < kre; ++kk) {
2062 int j = ihlink[hcoli[kk]];
2063 hcoli[kk] = j; /* permute the col index */
2064 hrowi[kk] = jpiv; /* permute the row index */
2065 ++hincol[j];
2066 }
2067 }
2068 /* temp [mpermu[i]] == mrstrt[i] */
2069 /* temp2[mpermu[i]] == hinrow[i] */
2070
2071 ninbas--; /* ???? */
2072 c_ekkscpy(nrow, &temp[1], &mrstrt[1]);
2073 c_ekkscpy(nrow, &temp2[1], &hinrow[1]);
2074
2075 /* now mrstrt, hinrow, hcoli and hrowi have been permuted */
2076
2077 /* Sort into column order as was stored by row */
2078 /* There will be an empty entry in front of each each column,
2079 * because we initialized hincol to 1s, and c_ekkrowq fills in
2080 * entries from the back */
2081 c_ekkrowq(hcoli, hrowi, dluval, mcstrt, hincol, nrow, ninbas);
2082
2083 /* The shuffle zeroed out column pointers */
2084 /* Put them back for L-U update */
2085 /* Also multiply U elements by - reciprocal of pivot element */
2086 /* Also decrement mcstrt/hincol to give "real" sizes */
2087 for (i = 1; i <= nrow; ++i) {
2088 int kx = --mcstrt[i];
2089 nel = --hincol[i];
2090 hrowi[kx] = nel;
2091 dluval[kx] = dvalpv[i];
2092 #ifndef NO_SHIFT
2093 for (int j = kx + 1; j <= kx + nel; j++)
2094 hrowi[j] = SHIFT_INDEX(hrowi[j]);
2095 #endif
2096 }
2097
2098 /* sort dense part */
2099 for (i = nrow - ndenuc + 1; i <= nrow; i++) {
2100 int kx = mcstrt[i] + 1; /* "real" entries start after pivot */
2101 int nel = hincol[i];
2102 c_ekk_sort2(&hrowi[kx], &dluval[kx], nel);
2103 }
2104
2105 /* Recompute number in U */
2106 nnentu = mcstrt[nrow] + hincol[nrow];
2107 mcstrt[nrow + 4] = nnentu + 1; /* magic - AND DEAD */
2108
2109 /* as not much room switch off fast etas */
2110 mrstrt[1] = 0; /* magic */
2111 fact->rows_ok = false;
2112 i = nrow + maxinv + 5; /* DEAD */
2113 } else {
2114 /* *************************************** */
2115 /* enough memory to do a bit faster */
2116 /* For each row compute reciprocal of pivot element and */
2117 /* take out of U */
2118 /* Also use HLINK(1 to permute column numbers */
2119 int ninbas = 0;
2120 int ilast; /* last available entry */
2121 int spareSpace;
2122 double *dluval2;
2123 /*int * hlink2 = ihlink+nrow;
2124 int * mrstrt2 = hlink2+nrow;*/
2125 /* mwork has order of row copy */
2126 EKKHlink *mwork = (reinterpret_cast< EKKHlink * >(fact->kw1adr)) - 1;
2127 fact->rows_ok = true;
2128
2129 if (if_sparse_update) {
2130 ilast = nnetas - nnentl;
2131 } else {
2132 /* missing out nnentl stuff */
2133 ilast = nnetas;
2134 }
2135 spareSpace = ilast - nnentu;
2136 need_more_space = false;
2137 /* save clean row copy if enough room */
2138 nroom = (spareSpace) / nrow;
2139 if (nrow < 10000) {
2140 if (nroom < 10) {
2141 need_more_space = true;
2142 }
2143 } else {
2144 if (nroom < 5 && !if_sparse_update) {
2145 need_more_space = true;
2146 }
2147 }
2148 if (nroom > CoinMin(50, maxinv)) {
2149 need_more_space = false;
2150 }
2151 if (need_more_space) {
2152 if (if_sparse_update) {
2153 int i1 = fact->eta_size + 10 * nrow;
2154 fact->eta_size = static_cast< int >(1.2 * fact->eta_size);
2155 if (i1 > fact->eta_size) {
2156 fact->eta_size = i1;
2157 }
2158 } else {
2159 fact->eta_size = static_cast< int >(1.05 * fact->eta_size);
2160 }
2161 } else {
2162 if (nroom < 11) {
2163 if (if_sparse_update) {
2164 int i1 = fact->eta_size + (11 - nroom) * nrow;
2165 fact->eta_size = static_cast< int >(1.2 * fact->eta_size);
2166 if (i1 > fact->eta_size) {
2167 fact->eta_size = i1;
2168 }
2169 }
2170 }
2171 }
2172 if (fact->maxNNetas > 0 && fact->eta_size > fact->maxNNetas) {
2173 fact->eta_size = fact->maxNNetas;
2174 }
2175 {
2176 /* we can swap de2val and dluval to save copying */
2177 int *eta_last = mpermu + nrow * 2 + 3;
2178 int *eta_next = eta_last + nrow + 2;
2179 int last = 0;
2180 eta_last[0] = -1;
2181 if (nnentl) {
2182 /* went into c_ekkcmfc - if not then in order */
2183 int next;
2184 /*next=mwork[((nrow+1)<<1)+1];*/
2185 next = mwork[nrow + 1].pre;
2186 #ifdef DEBUG
2187 j = mrstrt[next];
2188 #endif
2189 for (i = 1; i <= nrow; ++i) {
2190 int iperm = mpermu[next];
2191 eta_next[last] = iperm;
2192 eta_last[iperm] = last;
2193 temp[iperm] = mrstrt[next];
2194 temp2[iperm] = hinrow[next];
2195 #ifdef DEBUG
2196 if (mrstrt[next] != j)
2197 abort();
2198 j = mrstrt[next] + hinrow[next];
2199 #endif
2200 /*next= mwork[(next<<1)+2];*/
2201 next = mwork[next].suc;
2202 last = iperm;
2203 }
2204 } else {
2205 #ifdef DEBUG
2206 j = 0;
2207 #endif
2208 for (i = 1; i <= nrow; ++i) {
2209 int iperm = mpermu[i];
2210 eta_next[last] = iperm;
2211 eta_last[iperm] = last;
2212 temp[iperm] = mrstrt[i];
2213 temp2[iperm] = hinrow[i];
2214 last = iperm;
2215 #ifdef DEBUG
2216 if (mrstrt[i] <= j)
2217 abort();
2218 if (i > 1 && mrstrt[i] != j + hinrow[i - 1])
2219 abort();
2220 j = mrstrt[i];
2221 #endif
2222 }
2223 }
2224 eta_next[last] = nrow + 1;
2225 eta_last[nrow + 1] = last;
2226 eta_next[nrow + 1] = nrow + 2;
2227 c_ekkscpy(nrow, &temp[1], &mrstrt[1]);
2228 c_ekkscpy(nrow, &temp2[1], &hinrow[1]);
2229 i = eta_last[nrow + 1];
2230 ninbas = mrstrt[i] + hinrow[i] - 1;
2231 #ifdef DEBUG
2232 if (spareSpace < ninbas) {
2233 abort();
2234 }
2235 #endif
2236 c_ekkizero(nrow, &hincol[1]);
2237 #ifdef DEBUG
2238 for (i = nrow; i > 0; i--) {
2239 int krs = mrstrt[i];
2240 int jpiv = hcoli[krs];
2241 if (ihlink[jpiv] != i)
2242 abort();
2243 }
2244 #endif
2245 for (i = 1; i <= ninbas; ++i) {
2246 k = hcoli[i];
2247 k = ihlink[k];
2248 #ifdef DEBUG
2249 if (k <= 0 || k > nrow)
2250 abort();
2251 #endif
2252 hcoli[i] = k;
2253 hincol[k]++;
2254 }
2255 #ifdef DEBUG
2256 for (i = nrow; i > 0; i--) {
2257 int krs = mrstrt[i];
2258 int jpiv = hcoli[krs];
2259 if (jpiv != i)
2260 abort();
2261 if (krs > ninbas)
2262 abort();
2263 }
2264 #endif
2265 /* Sort into column order as was stored by row */
2266 k = 1;
2267 /* Position */
2268 for (kk = 1; kk <= nrow; ++kk) {
2269 nel = hincol[kk];
2270 mcstrt[kk] = k;
2271 hrowi[k] = nel - 1;
2272 k += hincol[kk];
2273 hincol[kk] = 0;
2274 }
2275 if (de2val) {
2276 dluval2 = de2val;
2277 } else {
2278 dluval2 = dluval + ninbas;
2279 }
2280 nnentu = k - 1;
2281 mcstrt[nrow + 4] = nnentu + 1;
2282 /* create column copy */
2283 for (i = nrow; i > 0; i--) {
2284 int krs = mrstrt[i];
2285 int kre = krs + hinrow[i];
2286 hinrow[i]--;
2287 mrstrt[i]++;
2288 {
2289 int kx = mcstrt[i];
2290 /*nel = hincol[i];
2291 if (hrowi[kx]!=nel) abort();
2292 hrowi[kx] = nel-1;*/
2293 dluval2[kx] = 1.0 / dluval[krs];
2294 /*hincol[i]=0;*/
2295 for (kk = krs + 1; kk < kre; ++kk) {
2296 int j = hcoli[kk];
2297 int iput = hincol[j] + 1;
2298 hincol[j] = iput;
2299 iput += mcstrt[j];
2300 hrowi[iput] = SHIFT_INDEX(i);
2301 dluval2[iput] = dluval[kk];
2302 }
2303 }
2304 }
2305 if (de2val) {
2306 double *a = dluval;
2307 double *address;
2308 /* move first down */
2309 i = eta_next[0];
2310 {
2311 int krs = mrstrt[i];
2312 nel = hinrow[i];
2313 for (j = 1; j <= nel; j++) {
2314 hcoli[j] = hcoli[j + krs - 1];
2315 dluval[j] = dluval[j + krs - 1];
2316 }
2317 }
2318 mrstrt[i] = 1;
2319 /****** swap dluval and de2val !!!! ******/
2320 /* should work even for dspace */
2321 /* move L part across */
2322 address = fact->xeeadr + 1;
2323 fact->xeeadr = fact->xe2adr - 1;
2324 fact->xe2adr = address;
2325 if (nnentl) {
2326 int n = xnetal - nrow - maxinv - 5;
2327 int j1, j2;
2328 int *mcstrt2 = mcstrt + nrow + maxinv + 4;
2329 j2 = mcstrt2[1];
2330 j1 = mcstrt2[n + 1] + 1;
2331 #if 0
2332 memcpy(de2val+j1,dluval+j1,(j2-j1+1)*sizeof(double));
2333 #else
2334 c_ekkdcpy(j2 - j1 + 1,
2335 (dluval + j1), (de2val + j1));
2336 #endif
2337 }
2338 dluval = de2val;
2339 de2val = a;
2340 } else {
2341 /* copy down dluval */
2342 #if 0
2343 memcpy(&dluval[1],&dluval2[1],ninbas*sizeof(double));
2344 #else
2345 c_ekkdcpy(ninbas,
2346 (dluval2 + 1), (dluval + 1));
2347 #endif
2348 }
2349 /* sort dense part */
2350 for (i = nrow - ndenuc + 1; i <= nrow; i++) {
2351 int kx = mcstrt[i] + 1;
2352 int nel = hincol[i];
2353 c_ekk_sort2(&hrowi[kx], &dluval[kx], nel);
2354 }
2355 }
2356 mrstrt[nrow + 1] = ilast + 1;
2357 }
2358 /* Find first non slack */
2359 for (i = 1; i <= nrow; ++i) {
2360 int kcs = mcstrt[i];
2361 if (hincol[i] != 0 || dluval[kcs] != SLACK_VALUE) {
2362 break;
2363 }
2364 }
2365 numberSlacks = i - 1;
2366 {
2367 /* set slacks to 1 */
2368 int *array = fact->krpadr + (fact->nrowmx + 2);
2369 int nSet = (numberSlacks) >> 5;
2370 int n2 = (fact->nrowmx + 32) >> 5;
2371 int i;
2372 memset(array, 0xff, nSet * sizeof(int));
2373 memset(array + nSet, 0, (n2 - nSet) * sizeof(int));
2374 for (i = nSet << 5; i <= numberSlacks; i++)
2375 c_ekk_Set(array, i);
2376 c_ekk_Unset(array, fact->nrow + 1); /* make sure off end not slack */
2377 #ifndef NDEBUG
2378 for (i = 1; i <= numberSlacks; i++)
2379 assert(c_ekk_IsSet(array, i));
2380 for (; i <= fact->nrow; i++)
2381 assert(!c_ekk_IsSet(array, i));
2382 #endif
2383 }
2384
2385 /* and set up backward pointers */
2386 /* clean up HPIVCO for fancy assembler stuff */
2387 /* xnetal was initialized to nrow + maxinv + 4 in c_ekktria, and grows */
2388 c_ekkscpy_0_1(maxinv + 1, 1, &hpivco[nrow + 4]); /* magic */
2389
2390 hpivco[xnetal] = 1;
2391 /* shuffle down for gaps so can get rid of hpivco for L */
2392 {
2393 const int lstart = nrow + maxinv + 5;
2394 int n = xnetal - lstart; /* number of L entries */
2395 int add, iel;
2396 int *hpivco_L = &hpivco[lstart];
2397 int *mcstrt_L = &mcstrt[lstart];
2398 if (nnentl) {
2399 /* elements of L were stored in descending order in dluval/hcoli */
2400 int kle = mcstrt_L[0];
2401 int kls = mcstrt_L[n] + 1;
2402
2403 if (if_sparse_update) {
2404 int i2, iel;
2405 int *mrstrt2 = &mrstrt[nrow];
2406
2407 /* need row copy of L */
2408 /* hpivro is spare for counts; just used as a temp buffer */
2409 c_ekkizero(nrow, &hpivro[1]);
2410
2411 /* permute L indices; count L row lengths */
2412 for (iel = kls; iel <= kle; ++iel) {
2413 int jrow = mpermu[UNSHIFT_INDEX(hrowi[iel])];
2414 hpivro[jrow]++;
2415 hrowi[iel] = SHIFT_INDEX(jrow);
2416 }
2417 {
2418 int ibase = nnetas - nnentl + 1;
2419 int firstDoRow = 0;
2420 for (i = 1; i <= nrow; i++) {
2421 mrstrt2[i] = ibase;
2422 if (hpivro[i] && !firstDoRow) {
2423 firstDoRow = i;
2424 }
2425 ibase += hpivro[i];
2426 hpivro[i] = mrstrt2[i];
2427 }
2428 if (!firstDoRow) {
2429 firstDoRow = nrow + 1;
2430 }
2431 mrstrt2[i] = ibase;
2432 fact->firstDoRow = firstDoRow;
2433 }
2434 i2 = mcstrt_L[n];
2435 for (i = n - 1; i >= 0; --i) {
2436 int i1 = mcstrt_L[i];
2437 int ipiv = hpivco_L[i];
2438 ipiv = mpermu[ipiv];
2439 hpivco_L[i] = ipiv;
2440 for (iel = i2; iel < i1; iel++) {
2441 int irow = UNSHIFT_INDEX(hrowi[iel + 1]);
2442 int iput = hpivro[irow];
2443 hpivro[irow] = iput + 1;
2444 hcoli[iput] = ipiv;
2445 de2val[iput] = dluval[iel + 1];
2446 }
2447 i2 = i1;
2448 }
2449 } else {
2450 /* just permute row numbers */
2451
2452 for (j = 0; j < n; ++j) {
2453 hpivco_L[j] = mpermu[hpivco_L[j]];
2454 }
2455 for (iel = kls; iel <= kle; ++iel) {
2456 int jrow = mpermu[UNSHIFT_INDEX(hrowi[iel])];
2457 hrowi[iel] = SHIFT_INDEX(jrow);
2458 }
2459 }
2460
2461 add = hpivco_L[n - 1] - hpivco_L[0] - n + 1;
2462 if (add) {
2463 int i;
2464 int last = hpivco_L[n - 1];
2465 int laststart = mcstrt_L[n];
2466 int base = hpivco_L[0] - 1;
2467 /* adjust so numbers match */
2468 mcstrt_L -= base;
2469 hpivco_L -= base;
2470 mcstrt_L[last] = laststart;
2471 for (i = n - 1; i >= 0; i--) {
2472 int ipiv = hpivco_L[i + base];
2473 while (ipiv < last) {
2474 mcstrt_L[last - 1] = laststart;
2475 hpivco_L[last - 1] = last;
2476 last--;
2477 }
2478 laststart = mcstrt_L[i + base];
2479 mcstrt_L[last - 1] = laststart;
2480 hpivco_L[last - 1] = last;
2481 last--;
2482 }
2483 xnetal += add;
2484 }
2485 }
2486 //int lstart=fact->lstart;
2487 //const int * COIN_RESTRICT hpivco = fact->kcpadr;
2488 fact->firstLRow = hpivco[lstart];
2489 }
2490 fact->nnentu = nnentu;
2491 fact->xnetal = xnetal;
2492 /* now we have xnetal * we can set up pointers */
2493 clp_setup_pointers(fact);
2494
2495 /* this is the array used in c_ekkbtrn; it is passed to c_ekkbtju as hpivco.
2496 * this gets modified by F-T as we pivot columns in and out.
2497 */
2498 {
2499 /* do new hpivco */
2500 int *hpivco_new = fact->kcpadr + 1;
2501 int *back = &fact->kcpadr[2 * nrow + maxinv + 4];
2502 /* set zeroth to stop illegal read */
2503 back[0] = 1;
2504
2505 hpivco_new[nrow + 1] = nrow + 1; /* deliberate loop for dense tests */
2506 hpivco_new[0] = 1;
2507
2508 for (i = 1; i <= nrow; i++) {
2509 hpivco_new[i] = i + 1;
2510 back[i + 1] = i;
2511 }
2512 back[1] = 0;
2513
2514 fact->first_dense = CoinMax(fact->first_dense, 4);
2515 fact->numberSlacks = numberSlacks;
2516 fact->lastSlack = numberSlacks;
2517 fact->firstNonSlack = hpivco_new[numberSlacks];
2518 }
2519
2520 /* also zero out permute region and nonzero */
2521 c_ekkdzero(nrow, (dpermu + 1));
2522
2523 if (if_sparse_update) {
2524 char *nonzero = reinterpret_cast< char * >(&mpermu[nrow + 1]); /* used in c_ekkbtrn */
2525 /*c_ekkizero(nrow,(int *)nonzero);*/
2526 c_ekkczero(nrow, nonzero);
2527 /*memset(nonzero,0,nrow*sizeof(int));*/ /* for faster method */
2528 }
2529 for (i = 1; i <= nrow; ++i) {
2530 hpivro[mpermu[i]] = i;
2531 }
2532
2533 } /* c_ekkshfv */
2534
c_ekkclcp1(const int * hcol,const int * mrstrt,int * hrow,int * mcstrt,int * hincol,int nnrow,int nncol,int ninbas)2535 static void c_ekkclcp1(const int *hcol, const int *mrstrt,
2536 int *hrow, int *mcstrt,
2537 int *hincol, int nnrow, int nncol,
2538 int ninbas)
2539 {
2540 int i, j, kc, kr, kre, krs, icol;
2541 int iput;
2542
2543 /* Create columnwise storage of row indices */
2544
2545 kc = 1;
2546 for (j = 1; j <= nncol; ++j) {
2547 mcstrt[j] = kc;
2548 kc += hincol[j];
2549 hincol[j] = 0;
2550 }
2551 mcstrt[nncol + 1] = ninbas + 1;
2552
2553 for (i = 1; i <= nnrow; ++i) {
2554 krs = mrstrt[i];
2555 kre = mrstrt[i + 1] - 1;
2556 for (kr = krs; kr <= kre; ++kr) {
2557 icol = hcol[kr];
2558 iput = hincol[icol];
2559 hincol[icol] = iput + 1;
2560 iput += mcstrt[icol];
2561 hrow[iput] = i;
2562 }
2563 }
2564 } /* c_ekkclcp */
c_ekkclcp2(const int * hcol,const double * dels,const int * mrstrt,int * hrow,double * dels2,int * mcstrt,int * hincol,int nnrow,int nncol,int ninbas)2565 inline void c_ekkclcp2(const int *hcol, const double *dels, const int *mrstrt,
2566 int *hrow, double *dels2, int *mcstrt,
2567 int *hincol, int nnrow, int nncol,
2568 int ninbas)
2569 {
2570 int i, j, kc, kr, kre, krs, icol;
2571 int iput;
2572
2573 /* Create columnwise storage of row indices */
2574
2575 kc = 1;
2576 for (j = 1; j <= nncol; ++j) {
2577 mcstrt[j] = kc;
2578 kc += hincol[j];
2579 hincol[j] = 0;
2580 }
2581 mcstrt[nncol + 1] = ninbas + 1;
2582
2583 for (i = 1; i <= nnrow; ++i) {
2584 krs = mrstrt[i];
2585 kre = mrstrt[i + 1] - 1;
2586 for (kr = krs; kr <= kre; ++kr) {
2587 icol = hcol[kr];
2588 iput = hincol[icol];
2589 hincol[icol] = iput + 1;
2590 iput += mcstrt[icol];
2591 hrow[iput] = i;
2592 dels2[iput] = dels[kr];
2593 }
2594 }
2595 } /* c_ekkclcp */
c_ekkslcf(register const EKKfactinfo * fact)2596 int c_ekkslcf(register const EKKfactinfo *fact)
2597 {
2598 int *hrow = fact->xeradr;
2599 int *hcol = fact->xecadr;
2600 double *dels = fact->xeeadr;
2601 int *hinrow = fact->xrnadr;
2602 int *hincol = fact->xcnadr;
2603 int *mrstrt = fact->xrsadr;
2604 int *mcstrt = fact->xcsadr;
2605 const int nrow = fact->nrow;
2606 int ninbas;
2607 /* space for etas */
2608 const int nnetas = fact->nnetas;
2609 ninbas = mcstrt[nrow + 1] - 1;
2610
2611 /* Now sort */
2612 if (ninbas << 1 > nnetas) {
2613 /* Put it in row order */
2614 int i, k;
2615 c_ekkrowq(hrow, hcol, dels, mrstrt, hinrow, nrow, ninbas);
2616 k = 1;
2617 for (i = 1; i <= nrow; ++i) {
2618 mrstrt[i] = k;
2619 k += hinrow[i];
2620 }
2621 mrstrt[nrow + 1] = k;
2622
2623 /* make a column copy without the extra values */
2624 c_ekkclcp1(hcol, mrstrt, hrow, mcstrt, hincol, nrow, nrow, ninbas);
2625 } else {
2626 /* Move elements up memory */
2627 c_ekkdcpy(ninbas,
2628 (dels + 1), (dels + ninbas + 1));
2629
2630 /* make a row copy with the extra values */
2631 c_ekkclcp2(hrow, &dels[ninbas], mcstrt, hcol, dels, mrstrt, hinrow, nrow, nrow, ninbas);
2632 }
2633 return (ninbas);
2634 } /* c_ekkslcf */
2635 /* Uwe H. Suhl, September 1986 */
2636 /* Removes lower and upper triangular factors from the matrix. */
2637 /* Code for routine: 102 */
2638 /* Return codes: */
2639 /* 0: ok */
2640 /* -5: not enough space in row file */
2641 /* 7: pivot too small - col sing */
2642 /*
2643 * This selects singleton columns and rows for the LU factorization.
2644 * Singleton columns require no
2645 *
2646 * (1) Note that columns are processed using a queue, not a stack;
2647 * this produces better pivots.
2648 *
2649 * (2) At most nrows elements are ever entered into the queue.
2650 *
2651 * (3) When pivoting singleton columns, every column that is part of
2652 * the pivot row is shortened by one, including the singleton column
2653 * itself; the hincol entries are updated appropriately.
2654 * Thus, pivoting on a singleton column may create other singleton columns
2655 * (but not singleton rows).
2656 * The dual property is true for rows.
2657 *
2658 * (4) Row entries (hrowi) are not changed when pivoting singleton columns.
2659 * Singleton columns that are created as a result of pivoting the
2660 * rows of other singleton columns will therefore have row entries
2661 * corresponding to those pivoted rows. Since we need to find the
2662 * row entry for the row being pivoted, we have to
2663 * search its row entries for the one whose hlink entry indicates
2664 * that it has not yet been pivoted.
2665 *
2666 * (9) As a result of pivoting columns, sections in hrowi corresponding to
2667 * pivoted columns are no longer needed, and entries in sections
2668 * for non-pivoted columns may have entries corresponding to pivoted rows.
2669 * This is why hrowi needs to be compacted.
2670 *
2671 * (5) When the row_pre and col_pre fields of the hlink struct contain
2672 * negative values, they indicate that the row has been pivoted, and
2673 * the negative of that value is the pivot order.
2674 * That is the only use for these fields in this routine.
2675 *
2676 * (6) This routine assumes that hlink is initialized to zeroes.
2677 * Under this assumption, the following is an invariant in this routine:
2678 *
2679 * (clink[i].pre < 0) ==> (hincol[i]==0)
2680 *
2681 * The converse is not true; see (15).
2682 *
2683 * The dual is also true, but only while pivoting singletong rows,
2684 * since we don't update hinrow while pivoting columns;
2685 * THESE VALUES ARE USED LATER, BUT I DON'T UNDERSTAND HOW YET.
2686 *
2687 * (7) hpivco is used for two purposes. The low end is used to implement the
2688 * queue when pivoting columns; the high end is used to hold eta-matrix
2689 * entries.
2690 *
2691 * (8) As a result of pivoting columns, for all i:1<=i<=nrow, either
2692 * hinrow[i] has not changed
2693 * or
2694 * hinrow[i] = 0
2695 * This is another way of saying that pivoting singleton columns cannot
2696 * create singleton rows.
2697 * The dual holds for hincol after pivoting rows.
2698 *
2699 * (10) In constrast to (4), while pivoting rows we
2700 * do not let the hcoli get out-of-date. That is because as part of
2701 * the process of numerical pivoting we have to find the row entries
2702 * for all the rows in the pivot column, so we may as well keep the
2703 * entries up to date. This is done by moving the last column entry
2704 * for each row into the entry that was used for the pivot column.
2705 *
2706 * (11) When pivoting a column, we must find the pivot row entry in
2707 * its row table. Sometimes we search for other things at the same time.
2708 * The same is true for pivoting columns. This search should never
2709 * fail.
2710 *
2711 * (12) Information concerning the eta matrices is stored in the high
2712 * ends of arrays that are also used to store information concerning
2713 * the basis; these arrays are: hpivco, mcstrt, dluval and hcoli.
2714 * Information is only stored in these arrays as a part of pivoting
2715 * singleton rows, since the only thing that needs to be saved as
2716 * a part of pivoting singleton columns is which rows and columns were chosen,
2717 * and this is stored in hlink.
2718 * Since they have to share the same array, the eta information grows
2719 * downward instead of upward. Eventually, eta information may grow
2720 * down to the top of the basis information. As pivoting proceeds,
2721 * more and more of this information is no longer needed, so when this
2722 * happens we can try compacting the arrays to see if we can recover
2723 * enough space. lstart points at the bottom entry in the arrays,
2724 * xnewro/xnewco at the top of the basis information, and each time we
2725 * pivot a singleton row we know that we will need exactly as many new
2726 * entries as there are rows in the pivot column, so we can easily
2727 * determine if we need more room. The variable maxinv may be used
2728 * to reserve extra room when inversion starts.
2729 *
2730 * (13) Eta information is stored in a fashion that is similar to how
2731 * matrices are stored. There is one entry in hpivco and mcstrt for
2732 * each eta (other than the initial ones for singleton columns and
2733 * for singleton rows that turn out to be singleton columns),
2734 * in the order they were chosen. hpivco records the pivot row,
2735 * and mcstrt points at the first entry in the other two arrays
2736 * for this row. dluval contains the actual eta values for the column,
2737 * and hcoli the rows these values were in.
2738 * These entries in mcstrt and hpivco grow upward; they start above
2739 * the entries used to store basis information.
2740 * (Actually, I don't see why they need to start maxinv+4 entries past the top).
2741 *
2742 * (14) c_ekkrwco assumes that invalidated hrowi/hcoli entries contain 0.
2743 *
2744 * (15) When pivoting singleton columns, it may possibly happen
2745 * that a row with all singleton column entries is created.
2746 * In this case, all of the columns will be enqueued, and pivoting
2747 * on any of them eliminates the rest, without their being chosen
2748 * as pivots. The dual holds for singleton rows.
2749 * DOES THIS INDICATE A SINGULARITY?
2750 *
2751 * (15) There are some aspects of the implementation that I find odd.
2752 * hrowi is not set to 0 for pivot rows while pivoting singleton columns,
2753 * which would make sense to me. Things don't work if this isn't done,
2754 * so the information is used somehow later on. Also, the information
2755 * for the pivot column is shifted to the front of the pivot row
2756 * when pivoting singleton columns; this is also necessary for reasons
2757 * I don't understand.
2758 */
c_ekktria(EKKfactinfo * fact,EKKHlink * rlink,EKKHlink * clink,int * nsingp,int * xnewcop,int * xnewrop,int * ncompactionsp,const int ninbas)2759 int c_ekktria(EKKfactinfo *fact,
2760 EKKHlink *rlink,
2761 EKKHlink *clink,
2762 int *nsingp,
2763 int *xnewcop, int *xnewrop,
2764 int *ncompactionsp,
2765 const int ninbas)
2766 {
2767 const int nrow = fact->nrow;
2768 const int maxinv = fact->maxinv;
2769 int *hcoli = fact->xecadr;
2770 double *dluval = fact->xeeadr;
2771 int *mrstrt = fact->xrsadr;
2772 int *hrowi = fact->xeradr;
2773 int *mcstrt = fact->xcsadr;
2774 int *hinrow = fact->xrnadr;
2775 int *hincol = fact->xcnadr;
2776 int *stack = fact->krpadr; /* normally hpivro */
2777 int *hpivco = fact->kcpadr;
2778 const double drtpiv = fact->drtpiv;
2779 CoinZeroN(reinterpret_cast< int * >(rlink + 1), static_cast< int >(nrow * (sizeof(EKKHlink) / sizeof(int))));
2780 CoinZeroN(reinterpret_cast< int * >(clink + 1), static_cast< int >(nrow * (sizeof(EKKHlink) / sizeof(int))));
2781
2782 fact->npivots = 0;
2783 /* Use NUSPIK to keep sum of deactivated row counts */
2784 fact->nuspike = 0;
2785 int xnetal = nrow + maxinv + 4;
2786 int xnewro = mrstrt[nrow] + hinrow[nrow] - 1;
2787 int xnewco = xnewro;
2788 int kmxeta = ninbas;
2789 int ncompactions = 0;
2790
2791 int i, j, k, kc, kce, kcs, npr;
2792 double pivot;
2793 int kipis, kipie, kjpis, kjpie, knprs, knpre;
2794 int ipivot, jpivot, stackc, stackr;
2795 #ifndef NDEBUG
2796 int kpivot = -1;
2797 #else
2798 int kpivot = -1;
2799 #endif
2800 int epivco, kstart, maxstk;
2801 int irtcod = 0;
2802 int lastSlack = 0;
2803
2804 int lstart = fact->nnetas + 1;
2805 /*int nnentu = ninbas; */
2806 int lstart_minus_nnentu = lstart - ninbas;
2807 /* do initial column singletons - as can do faster */
2808 for (jpivot = 1; jpivot <= nrow; ++jpivot) {
2809 if (hincol[jpivot] == 1) {
2810 ipivot = hrowi[mcstrt[jpivot]];
2811 if (ipivot > lastSlack) {
2812 lastSlack = ipivot;
2813 } else {
2814 /* so we can't put a structural over a slack */
2815 break;
2816 }
2817 kipis = mrstrt[ipivot];
2818 #if 1
2819 assert(hcoli[kipis] == jpivot);
2820 #else
2821 if (hcoli[kipis] != jpivot) {
2822 kpivot = kipis + 1;
2823 while (hcoli[kpivot] != jpivot)
2824 kpivot++;
2825 #ifdef DEBUG
2826 kipie = kipis + hinrow[ipivot];
2827 if (kpivot >= kipie) {
2828 abort();
2829 }
2830 #endif
2831 pivot = dluval[kpivot];
2832 dluval[kpivot] = dluval[kipis];
2833 dluval[kipis] = pivot;
2834 hcoli[kpivot] = hcoli[kipis];
2835 hcoli[kipis] = jpivot;
2836 }
2837 #endif
2838 if (dluval[kipis] == SLACK_VALUE) {
2839 /* record the new pivot row and column */
2840 ++fact->npivots;
2841 rlink[ipivot].pre = -fact->npivots;
2842 clink[jpivot].pre = -fact->npivots;
2843 hincol[jpivot] = 0;
2844 fact->nuspike += hinrow[ipivot];
2845 } else {
2846 break;
2847 }
2848 } else {
2849 break;
2850 }
2851 }
2852 /* Fill queue with other column singletons and clean up */
2853 maxstk = 0;
2854 for (j = 1; j <= nrow; ++j) {
2855 if (hincol[j]) {
2856 int n = 0;
2857 kcs = mcstrt[j];
2858 kce = mcstrt[j + 1];
2859 for (k = kcs; k < kce; ++k) {
2860 if (!(rlink[hrowi[k]].pre < 0)) {
2861 n++;
2862 }
2863 }
2864 hincol[j] = n;
2865 if (n == 1) {
2866 /* we just created a new singleton column - enqueue it */
2867 ++maxstk;
2868 stack[maxstk] = j;
2869 }
2870 }
2871 }
2872 stackc = 0; /* (1) */
2873
2874 while (!(stackc >= maxstk)) { /* (1) */
2875 /* dequeue the next entry */
2876 ++stackc;
2877 jpivot = stack[stackc];
2878
2879 /* (15) */
2880 if (hincol[jpivot] != 0) {
2881
2882 for (k = mcstrt[jpivot]; rlink[hrowi[k]].pre < 0; k++) {
2883 /* (4) */
2884 }
2885 ipivot = hrowi[k];
2886
2887 /* All the columns in this row are being shortened. */
2888 kipis = mrstrt[ipivot];
2889 kipie = kipis + hinrow[ipivot];
2890 for (k = kipis; k < kipie; ++k) {
2891 j = hcoli[k];
2892 --hincol[j]; /* (3) (6) */
2893
2894 if (j == jpivot) {
2895 kpivot = k; /* (11) */
2896 } else if (hincol[j] == 1) {
2897 /* we just created a new singleton column - enqueue it */
2898 ++maxstk;
2899 stack[maxstk] = j;
2900 }
2901 }
2902
2903 /* record the new pivot row and column */
2904 ++fact->npivots;
2905 rlink[ipivot].pre = -fact->npivots;
2906 clink[jpivot].pre = -fact->npivots;
2907
2908 fact->nuspike += hinrow[ipivot];
2909
2910 /* check the pivot */
2911 assert(kpivot > 0);
2912 pivot = dluval[kpivot];
2913 if (fabs(pivot) < drtpiv) {
2914 irtcod = 7;
2915 ++(*nsingp);
2916 rlink[ipivot].pre = -nrow - 1;
2917 clink[jpivot].pre = -nrow - 1;
2918 }
2919
2920 /* swap the pivot column entry with the first one. */
2921 /* I don't know why. */
2922 dluval[kpivot] = dluval[kipis];
2923 dluval[kipis] = pivot;
2924 hcoli[kpivot] = hcoli[kipis];
2925 hcoli[kipis] = jpivot;
2926 }
2927 }
2928 /* (8) */
2929
2930 /* The entire basis may already be triangular */
2931 if (fact->npivots < nrow) {
2932
2933 /* (9) */
2934 kstart = 0;
2935 for (j = 1; j <= nrow; ++j) {
2936 if (!(clink[j].pre < 0)) {
2937 kcs = mcstrt[j];
2938 kce = mcstrt[j + 1];
2939
2940 mcstrt[j] = kstart + 1;
2941
2942 for (k = kcs; k < kce; ++k) {
2943 if (!(rlink[hrowi[k]].pre < 0)) {
2944 ++kstart;
2945 hrowi[kstart] = hrowi[k];
2946 }
2947 }
2948 hincol[j] = kstart - mcstrt[j] + 1;
2949 }
2950 }
2951 xnewco = kstart;
2952
2953 /* Fill stack with initial row singletons that haven't been pivoted away */
2954 stackr = 0;
2955 for (i = 1; i <= nrow; ++i) {
2956 if (!(rlink[i].pre < 0) && (hinrow[i] == 1)) {
2957 ++stackr;
2958 stack[stackr] = i;
2959 }
2960 }
2961
2962 while (!(stackr <= 0)) {
2963 ipivot = stack[stackr];
2964 assert(ipivot);
2965 --stackr;
2966
2967 #if 1
2968 assert(rlink[ipivot].pre >= 0);
2969 #else
2970 /* This test is probably unnecessary: rlink[i].pre < 0 ==> hinrow[i]==0 */
2971 if (rlink[ipivot].pre < 0) {
2972 continue;
2973 }
2974 #endif
2975
2976 /* (15) */
2977 if (hinrow[ipivot] != 0) {
2978
2979 /* This is a singleton row, which means it has exactly one column */
2980 jpivot = hcoli[mrstrt[ipivot]];
2981
2982 kjpis = mcstrt[jpivot];
2983 epivco = hincol[jpivot] - 1;
2984 hincol[jpivot] = 0; /* this column is being pivoted away */
2985
2986 /* (11) */
2987 kjpie = kjpis + epivco;
2988 for (k = kjpis; k <= kjpie; ++k) {
2989 if (ipivot == hrowi[k])
2990 break;
2991 }
2992 /* ASSERT (k <= kjpie) */
2993
2994 /* move the last row entry for the pivot column into the pivot row's entry */
2995 /* I don't know why */
2996 hrowi[k] = hrowi[kjpie];
2997
2998 /* invalidate the (old) last row entry of the pivot column */
2999 /* I don't know why */
3000 hrowi[kjpie] = 0;
3001
3002 /* (12) */
3003 if (!(xnewro + epivco < lstart)) {
3004 int kstart;
3005
3006 if (!(epivco < lstart_minus_nnentu)) {
3007 irtcod = -5;
3008 break;
3009 }
3010 kstart = c_ekkrwco(fact, dluval, hcoli, mrstrt, hinrow, xnewro);
3011 ++ncompactions;
3012 kmxeta += (xnewro - kstart) << 1;
3013 xnewro = kstart;
3014 }
3015 if (!(xnewco + epivco < lstart)) {
3016 if (!(epivco < lstart_minus_nnentu)) {
3017 irtcod = -5;
3018 break;
3019 }
3020 xnewco = c_ekkclco(fact, hrowi, mcstrt, hincol, xnewco);
3021 ++ncompactions;
3022
3023 /* HINCOL MAY HAVE CHANGED ??? (JJHF) */
3024 epivco = hincol[jpivot];
3025 }
3026
3027 /* record the new pivot row and column */
3028 ++fact->npivots;
3029 rlink[ipivot].pre = -fact->npivots;
3030 clink[jpivot].pre = -fact->npivots;
3031
3032 /* no update for nuspik */
3033
3034 /* check the pivot */
3035 pivot = dluval[mrstrt[ipivot]];
3036 if (fabs(pivot) < drtpiv) {
3037 /* If the pivot is too small, reject it, but keep going */
3038 irtcod = 7;
3039 rlink[ipivot].pre = -nrow - 1;
3040 clink[jpivot].pre = -nrow - 1;
3041 }
3042
3043 /* Perform numerical part of elimination. */
3044 if (!(epivco <= 0)) {
3045 ++xnetal;
3046 mcstrt[xnetal] = lstart - 1;
3047 hpivco[xnetal] = ipivot;
3048 pivot = -1.f / pivot;
3049
3050 kcs = mcstrt[jpivot];
3051 kce = kcs + epivco - 1;
3052 hincol[jpivot] = 0;
3053
3054 for (kc = kcs; kc <= kce; ++kc) {
3055 npr = hrowi[kc];
3056
3057 /* why bother? */
3058 hrowi[kc] = 0;
3059
3060 --hinrow[npr]; /* (3) */
3061 if (hinrow[npr] == 1) {
3062 /* this may create new singleton rows */
3063 ++stackr;
3064 stack[stackr] = npr;
3065 }
3066
3067 /* (11) */
3068 knprs = mrstrt[npr];
3069 knpre = knprs + hinrow[npr];
3070 for (k = knprs; k <= knpre; ++k) {
3071 if (jpivot == hcoli[k]) {
3072 kpivot = k;
3073 break;
3074 }
3075 }
3076 /* ASSERT (kpivot <= knpre) */
3077
3078 {
3079 /* (10) */
3080 double elemnt = dluval[kpivot];
3081 dluval[kpivot] = dluval[knpre];
3082 hcoli[kpivot] = hcoli[knpre];
3083
3084 hcoli[knpre] = 0; /* (14) */
3085
3086 /* store elementary row transformation */
3087 --lstart;
3088 dluval[lstart] = elemnt * pivot;
3089 hcoli[lstart] = npr;
3090 }
3091 }
3092 }
3093 }
3094 }
3095 }
3096 /* (8) */
3097
3098 *xnewcop = xnewco;
3099 *xnewrop = xnewro;
3100 fact->xnetal = xnetal;
3101 fact->nnentu = lstart - lstart_minus_nnentu;
3102 fact->kmxeta = kmxeta;
3103 *ncompactionsp = ncompactions;
3104
3105 return (irtcod);
3106 } /* c_ekktria */
3107 #endif
3108
3109 /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
3110 */
3111