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