1/* $Id: CoinOslFactorization.cpp 1448 2011-06-19 15:34:41Z stefan $ */
2// Copyright (C) 1987, 2009, International Business Machines
3// Corporation and others. All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#include "CoinUtilsConfig.h"
7
8#include <cassert>
9#include "CoinPragma.hpp"
10#include "CoinOslFactorization.hpp"
11#include "CoinOslC.h"
12#include "CoinIndexedVector.hpp"
13#include "CoinHelperFunctions.hpp"
14#include "CoinPackedMatrix.hpp"
15#include "CoinTypes.hpp"
16#include "CoinFinite.hpp"
17#include <stdio.h>
18static void c_ekksmem(EKKfactinfo *fact,int numberRows,int maximumPivots);
19static void c_ekksmem_copy(EKKfactinfo *fact,const EKKfactinfo * rhsFact);
20static void c_ekksmem_delete(EKKfactinfo *fact);
21//:class CoinOslFactorization. Deals with Factorization and Updates
22// CoinOslFactorization. Constructor
23CoinOslFactorization::CoinOslFactorization ( )
24 : CoinOtherFactorization()
25{
26 gutsOfInitialize();
27}
28
29/// Copy constructor
30CoinOslFactorization::CoinOslFactorization ( const CoinOslFactorization &other)
31 : CoinOtherFactorization(other)
32{
33 gutsOfInitialize();
34 gutsOfCopy(other);
35}
36// Clone
37CoinOtherFactorization *
38CoinOslFactorization::clone() const
39{
40 return new CoinOslFactorization(*this);
41}
42/// The real work of constructors etc
43void CoinOslFactorization::gutsOfDestructor(bool clearFact)
44{
45 delete [] elements_;
46 delete [] pivotRow_;
47 delete [] workArea_;
48 elements_ = NULL;
49 pivotRow_ = NULL;
50 workArea_ = NULL;
51 numberRows_ = 0;
52 numberColumns_ = 0;
53 numberGoodU_ = 0;
54 status_ = -1;
55 maximumRows_=0;
56 maximumSpace_=0;
57 solveMode_=0;
58 if (clearFact)
59 c_ekksmem_delete(&factInfo_);
60}
61void CoinOslFactorization::gutsOfInitialize(bool zapFact)
62{
63 pivotTolerance_ = 1.0e-1;
64 zeroTolerance_ = 1.0e-13;
65#ifndef COIN_FAST_CODE
66 slackValue_ = -1.0;
67#endif
68 maximumPivots_=200;
69 relaxCheck_=1.0;
70 numberRows_ = 0;
71 numberColumns_ = 0;
72 numberGoodU_ = 0;
73 status_ = -1;
74 numberPivots_ = 0;
75 maximumRows_=0;
76 maximumSpace_=0;
77 elements_ = NULL;
78 pivotRow_ = NULL;
79 workArea_ = NULL;
80 solveMode_=0;
81 if (zapFact) {
82 memset(&factInfo_,0,sizeof(factInfo_));
83 factInfo_.maxinv=100;
84 factInfo_.drtpiv=1.0e-10;
85 factInfo_.zeroTolerance=1.0e-12;
86 factInfo_.zpivlu=0.1;
87 factInfo_.areaFactor=1.0;
88 factInfo_.nbfinv=100;
89 }
90}
91// ~CoinOslFactorization. Destructor
92CoinOslFactorization::~CoinOslFactorization ( )
93{
94 gutsOfDestructor();
95}
96// =
97CoinOslFactorization & CoinOslFactorization::operator = ( const CoinOslFactorization & other ) {
98 if (this != &other) {
99 bool noGood = factInfo_.nrowmx!=other.factInfo_.nrowmx&&
100 factInfo_.eta_size!=other.factInfo_.eta_size;
101 gutsOfDestructor(noGood);
102 gutsOfInitialize(noGood);
103 gutsOfCopy(other);
104 }
105 return *this;
106}
107#define WORK_MULT 2
108void CoinOslFactorization::gutsOfCopy(const CoinOslFactorization &other)
109{
110 pivotTolerance_ = other.pivotTolerance_;
111 zeroTolerance_ = other.zeroTolerance_;
112#ifndef COIN_FAST_CODE
113 slackValue_ = other.slackValue_;
114#endif
115 relaxCheck_ = other.relaxCheck_;
116 numberRows_ = other.numberRows_;
117 numberColumns_ = other.numberColumns_;
118 maximumRows_ = other.maximumRows_;
119 maximumSpace_ = other.maximumSpace_;
120 solveMode_ = other.solveMode_;
121 numberGoodU_ = other.numberGoodU_;
122 maximumPivots_ = other.maximumPivots_;
123 numberPivots_ = other.numberPivots_;
124 factorElements_ = other.factorElements_;
125 status_ = other.status_;
126 elements_ = NULL;
127 pivotRow_ = NULL;
128 workArea_ = NULL;
129 c_ekksmem_copy(&factInfo_,&other.factInfo_);
130}
131
132// getAreas. Gets space for a factorization
133//called by constructors
134void
135CoinOslFactorization::getAreas ( int numberOfRows,
136 int numberOfColumns,
137 CoinBigIndex maximumL,
138 CoinBigIndex maximumU )
139{
140
141 numberRows_ = numberOfRows;
142 numberColumns_ = numberOfColumns;
143 CoinBigIndex size = static_cast<CoinBigIndex>(factInfo_.areaFactor*
144 (maximumL+maximumU));
145 factInfo_.zeroTolerance=zeroTolerance_;
146 // If wildly out redo
147 if (maximumRows_>numberRows_+1000) {
148 maximumRows_=0;
149 maximumSpace_=0;
150 factInfo_.last_eta_size=0;
151 }
152 if (size>maximumSpace_) {
153 //delete [] elements_;
154 //elements_ = new CoinFactorizationDouble [size];
155 maximumSpace_ = size;
156 }
157 factInfo_.lastEtaCount = factInfo_.nnentu+factInfo_.nnentl;
158 int oldnnetas=factInfo_.last_eta_size;
159 // If we are going to increase then be on safe side
160 if (size>oldnnetas)
161 size = static_cast<int>(1.1*size);
162 factInfo_.eta_size=CoinMax(size,oldnnetas);
163 //printf("clp size %d, old %d now %d - iteration %d - last count %d - rows %d,%d,%d\n",
164 // size,oldnnetas,factInfo_.eta_size,factInfo_.iterno,factInfo_.lastEtaCount,
165 //numberRows_,factInfo_.nrowmx,factInfo_.nrow);
166 //if (!factInfo_.iterno) {
167 //printf("here\n");
168 //}
169 /** Get solve mode e.g. 0 C++ code, 1 Lapack, 2 choose
170 If 4 set then values pass
171 if 8 set then has iterated
172 */
173 solveMode_ &= 4+8; // clear bottom bits
174 factInfo_.ifvsol= ((solveMode_&4)!=0) ? 1 : 0;
175 if ((solveMode_&8)!=0) {
176 factInfo_.ifvsol=0;
177 factInfo_.invok=1;
178 } else {
179 factInfo_.iter0=factInfo_.iterno;
180 factInfo_.invok=-1;
181 factInfo_.if_sparse_update=0;
182 }
183#if 0
184 if (!factInfo_.if_sparse_update &&
185 factInfo_.iterno>factInfo_.iter0 &&
186 numberRows_>=C_EKK_GO_SPARSE) {
187 printf("count %d rows %d etasize %d\n",
188 factInfo_.lastEtaCount,factInfo_.nrow,factInfo_.eta_size);
189
190 }
191#endif
192 if (!factInfo_.if_sparse_update &&
193 factInfo_.iterno>factInfo_.iter0 &&
194 numberRows_>=C_EKK_GO_SPARSE &&
195 (factInfo_.lastEtaCount>>2)<factInfo_.nrow&&
196 !factInfo_.switch_off_sparse_update) {
197#if PRINT_DEBUG
198 printf("**** Switching on sparse update - etacount\n");
199#endif
200 /* I suspect this can go into c_ekksslvf;
201 * if c_ekkshff decides to switch sparse_update off,
202 * then it problably always switches it off (?)
203 */
204 factInfo_.if_sparse_update=2;
205 }
206 c_ekksmem(&factInfo_,numberRows_,maximumPivots_);
207 if (numberRows_>maximumRows_) {
208 maximumRows_ = numberRows_;
209 //delete [] pivotRow_;
210 //delete [] workArea_;
211 //pivotRow_ = new int [2*maximumRows_+maximumPivots_];
212 //workArea_ = new CoinFactorizationDouble [maximumRows_*WORK_MULT];
213 }
214}
215
216// preProcess.
217void
218CoinOslFactorization::preProcess ()
219{
220 factInfo_.zpivlu=pivotTolerance_;
221 // Go to Fortran
222 int * hcoli=factInfo_.xecadr+1;
223 int * indexRowU = factInfo_.xeradr+1;
224 CoinBigIndex * startColumnU=factInfo_.xcsadr+1;
225 for (int i=0;i<numberRows_;i++) {
226 int start = startColumnU[i];
227 startColumnU[i]++; // to Fortran
228 for (int j=start;j<startColumnU[i+1];j++) {
229 indexRowU[j]++; // to Fortran
230 hcoli[j]=i+1; // to Fortran
231 }
232 }
233 startColumnU[numberRows_]++; // to Fortran
234
235 /* can do in column order - no zeros or duplicates */
236#ifndef NDEBUG
237 int ninbas =
238#endif
239 c_ekkslcf(&factInfo_);
240 assert (ninbas>0);
241}
242
243//Does factorization
244int
245CoinOslFactorization::factor ( )
246{
247 /* Uwe's factorization (sort of) */
248 int irtcod = c_ekklfct(&factInfo_);
249
250 /* Check return code */
251 /* 0 - Fine , 1 - Backtrack, 2 - Singularities on initial, 3-Fatal */
252 /* now 5-Need more memory */
253
254 status_= 0;
255 if (factInfo_.eta_size>factInfo_.last_eta_size) {
256 factInfo_.areaFactor *= factInfo_.eta_size;
257 factInfo_.areaFactor /= factInfo_.last_eta_size;
258#ifdef CLP_INVESTIGATE
259 printf("areaFactor increased to %g\n",factInfo_.areaFactor);
260#endif
261 }
262 if (irtcod==5) {
263 status_=-99;
264 assert (factInfo_.eta_size>factInfo_.last_eta_size) ;
265#ifdef CLP_INVESTIGATE
266 printf("need more memory\n");
267#endif
268 } else if (irtcod) {
269 status_=-1;
270 //printf("singular %d\n",irtcod);
271 }
272 return status_;
273}
274// Makes a non-singular basis by replacing variables
275void
276CoinOslFactorization::makeNonSingular(int * sequence, int numberColumns)
277{
278 const EKKHlink *rlink = factInfo_.kp1adr;
279 const EKKHlink *clink = factInfo_.kp2adr;
280 int nextRow=0;
281 //int * mark = reinterpret_cast<int *>(factInfo_.kw1adr);
282 //int nr=0;
283 //int nc=0;
284#if 0
285 for (int i=0;i<numberRows_;i++) {
286 //mark[i]=-1;
287 if (rlink[i].pre>=0||rlink[i].pre==-(numberRows_+1)) {
288 nr++;
289 printf("%d rl %d cl %d\n",i,rlink[i].pre,clink[i].pre);
290 }
291 if (clink[i].pre>=0||clink[i].pre==-(numberRows_+1)) {
292 nc++;
293 printf("%d rl %d cl %d\n",i,rlink[i].pre,clink[i].pre);
294 }
295 }
296#endif
297 //printf("nr %d nc %d\n",nr,nc);
298 bool goodPass=true;
299 int numberDone=0;
300 for (int i=0;i<numberRows_;i++) {
301 int cRow =(-clink[i].pre)-1;
302 if (cRow==numberRows_||cRow<0) {
303 // throw out
304 for (;nextRow<numberRows_;nextRow++) {
305 int rRow =(-rlink[nextRow].pre)-1;
306 if (rRow==numberRows_||rRow<0)
307 break;
308 }
309 if (nextRow<numberRows_) {
310 sequence[i]=nextRow+numberColumns;
311 nextRow++;
312 numberDone++;
313 } else {
314 goodPass=false;
315 assert(numberDone);
316 //printf("BAD singular at row %d\n",i);
317 break;
318 }
319 }
320 }
321#ifndef NDEBUG
322 if (goodPass) {
323 for (;nextRow<numberRows_;nextRow++) {
324 int rRow =(-rlink[nextRow].pre)-1;
325 assert (!(rRow==numberRows_||rRow<0));
326 }
327 }
328#endif
329}
330// Does post processing on valid factorization - putting variables on correct rows
331void
332CoinOslFactorization::postProcess(const int * sequence, int * pivotVariable)
333{
334 factInfo_.iterin=factInfo_.iterno;
335 factInfo_.npivots=0;
336 numberPivots_=0;
337 const int * permute3 = factInfo_.mpermu+1;
338 assert (permute3==reinterpret_cast<const int *>
339 (factInfo_.kadrpm+numberRows_+1));
340 // this is ridiculous - must be better way
341 int * permute2 = reinterpret_cast<int *>(factInfo_.kw1adr);
342 const int * permute = reinterpret_cast<const int *>(factInfo_.kp2adr);
343 for (int i=0;i<numberRows_;i++) {
344 permute2[permute[i]-1]=i;
345 }
346 for (int i=0;i<numberRows_;i++) {
347 // the row is i
348 // the column is whatever matches k3[i] in k1
349 int look=permute3[i]-1;
350 int j=permute2[look];
351 int k = sequence[j];
352 pivotVariable[i]=k;
353 }
354#ifdef CLP_REUSE_ETAS
355 int * start = factInfo_.xcsadr+1;
356 int * putSeq = factInfo_.xrsadr+2*factInfo_.nrowmx+2;
357 int * position = putSeq+factInfo_.maxinv;
358 int * putStart = position+factInfo_.maxinv;
359 memcpy(putStart,start,numberRows_*sizeof(int));
360 int iLast=start[numberRows_-1];
361 putStart[numberRows_]=iLast+factInfo_.xeradr[iLast]+1;
362 factInfo_.save_nnentu=factInfo_.nnentu;
363#endif
364#ifndef NDEBUG
365 {
366 int lstart=numberRows_+factInfo_.maxinv+5;
367 int ndo = factInfo_.xnetal-lstart;
368 double * dluval=factInfo_.xeeadr;
369 int * mcstrt = factInfo_.xcsadr+lstart;
370 if (ndo)
371 assert (dluval[mcstrt[ndo]+1]<1.0e50);
372 }
373#endif
374}
375/* Replaces one Column to basis,
376 returns 0=OK, 1=Probably OK, 2=singular, 3=no room
377 If checkBeforeModifying is true will do all accuracy checks
378 before modifying factorization. Whether to set this depends on
379 speed considerations. You could just do this on first iteration
380 after factorization and thereafter re-factorize
381 partial update already in U */
382int
383CoinOslFactorization::replaceColumn ( CoinIndexedVector * regionSparse,
384 int pivotRow,
385 double pivotCheck ,
386 bool /*checkBeforeModifying*/,
387 double acceptablePivot)
388{
389 if (numberPivots_+1==maximumPivots_)
390 return 3;
391 int *regionIndex = regionSparse->getIndices ( );
392 double *region = regionSparse->denseVector ( );
393 int orig_nincol=0;
394 double saveTolerance = factInfo_.drtpiv;
395 factInfo_.drtpiv=acceptablePivot;
396 int returnCode=c_ekketsj(&factInfo_,region-1,
397 regionIndex,
398 pivotCheck,orig_nincol,
399 numberPivots_,&factInfo_.nuspike,
400 pivotRow+1,
401 reinterpret_cast<int *>(factInfo_.kw1adr));
402 factInfo_.drtpiv=saveTolerance;
403 if (returnCode!=2)
404 numberPivots_++;
405#ifndef NDEBUG
406 {
407 int lstart=numberRows_+factInfo_.maxinv+5;
408 int ndo = factInfo_.xnetal-lstart;
409 double * dluval=factInfo_.xeeadr;
410 int * mcstrt = factInfo_.xcsadr+lstart;
411 if (ndo)
412 assert (dluval[mcstrt[ndo]+1]<1.0e50);
413 }
414#endif
415 return returnCode;
416}
417/* This version has same effect as above with FTUpdate==false
418 so number returned is always >=0 */
419int
420CoinOslFactorization::updateColumn ( CoinIndexedVector * regionSparse,
421 CoinIndexedVector * regionSparse2,
422 bool /*noPermute*/) const
423{
424#ifndef NDEBUG
425 {
426 int lstart=numberRows_+factInfo_.maxinv+5;
427 int ndo = factInfo_.xnetal-lstart;
428 double * dluval=factInfo_.xeeadr;
429 int * mcstrt = factInfo_.xcsadr+lstart;
430 if (ndo)
431 assert (dluval[mcstrt[ndo]+1]<1.0e50);
432 }
433#endif
434 assert (numberRows_==numberColumns_);
435 double *region2 = regionSparse2->denseVector ( );
436 int *regionIndex2 = regionSparse2->getIndices ( );
437 int numberNonZero = regionSparse2->getNumElements ( );
438 double *region = regionSparse->denseVector ( );
439 //const int * permuteIn = factInfo_.mpermu+1;
440 // Stuff is put one up so won't get illegal read
441 assert (!region[numberRows_]);
442 assert (!regionSparse2->packedMode());
443#if 0
444 int first=numberRows_;
445 for (int j=0;j<numberNonZero;j++) {
446 int jRow = regionIndex2[j];
447 int iRow = permuteIn[jRow];
448 region[iRow]=region2[jRow];
449 first=CoinMin(first,iRow);
450 region2[jRow]=0.0;
451 }
452#endif
453 numberNonZero=c_ekkftrn(&factInfo_,
454 region2-1,region,regionIndex2,numberNonZero);
455 regionSparse2->setNumElements(numberNonZero);
456 return 0;
457}
458/* Updates one column (FTRAN) from regionSparse2
459 Tries to do FT update
460 number returned is negative if no room
461 regionSparse starts as zero and is zero at end.
462 Note - if regionSparse2 packed on input - will be packed on output
463*/
464int
465CoinOslFactorization::updateColumnFT ( CoinIndexedVector * regionSparse,
466 CoinIndexedVector * regionSparse2,
467 bool /*noPermute*/)
468{
469 assert (numberRows_==numberColumns_);
470 double *region2 = regionSparse2->denseVector ( );
471 int *regionIndex2 = regionSparse2->getIndices ( );
472 int numberNonZero = regionSparse2->getNumElements ( );
473 assert (regionSparse2->packedMode());
474 // packed mode
475 //int numberNonInOriginal=numberNonZero;
476 //double *dpermu = factInfo_.kadrpm;
477 // Use region instead of dpermu
478 double * save =factInfo_.kadrpm;
479 factInfo_.kadrpm=regionSparse->denseVector()-1;
480 int nuspike=c_ekkftrn_ft(&factInfo_, region2,regionIndex2,
481 &numberNonZero);
482 factInfo_.kadrpm=save;
483 regionSparse2->setNumElements(numberNonZero);
484 //regionSparse2->print();
485 factInfo_.nuspike=nuspike;
486 return nuspike;
487}
488
489
490int
491CoinOslFactorization::updateTwoColumnsFT(CoinIndexedVector * regionSparse1,
492 CoinIndexedVector * regionSparse2,
493 CoinIndexedVector * regionSparse3,
494 bool /*noPermute*/)
495{
496#if 1
497 // probably best to merge on a LU part by part
498 // but can try full merge
499 double *region2 = regionSparse2->denseVector ( );
500 int *regionIndex2 = regionSparse2->getIndices ( );
501 int numberNonZero2 = regionSparse2->getNumElements ( );
502 assert (regionSparse2->packedMode());
503
504 assert (numberRows_==numberColumns_);
505 double *region3 = regionSparse3->denseVector ( );
506 int *regionIndex3 = regionSparse3->getIndices ( );
507 int numberNonZero3 = regionSparse3->getNumElements ( );
508 double *region = regionSparse1->denseVector ( );
509 // Stuff is put one up so won't get illegal read
510 assert (!region[numberRows_]);
511 assert (!regionSparse3->packedMode());
512 // packed mode
513 //double *dpermu = factInfo_.kadrpm;
514#if 0
515 factInfo_.nuspike=c_ekkftrn_ft(&factInfo_, region2,regionIndex2,
516 &numberNonZero2);
517 numberNonZero3=c_ekkftrn(&factInfo_,
518 region3-1,region,regionIndex3,numberNonZero3);
519#else
520 c_ekkftrn2(&factInfo_,region3-1,region,regionIndex3,&numberNonZero3,
521 region2,regionIndex2,&numberNonZero2);
522#endif
523 regionSparse2->setNumElements(numberNonZero2);
524 regionSparse3->setNumElements(numberNonZero3);
525 return factInfo_.nuspike;
526#else
527 // probably best to merge on a LU part by part
528 // but can try full merge
529 int returnCode= updateColumnFT(regionSparse1,
530 regionSparse2);
531 updateColumn(regionSparse1,
532 regionSparse3,
533 noPermute);
534 return returnCode;
535#endif
536}
537
538/* Updates one column (BTRAN) from regionSparse2
539 regionSparse starts as zero and is zero at end
540 Note - if regionSparse2 packed on input - will be packed on output
541*/
542int
543CoinOslFactorization::updateColumnTranspose ( CoinIndexedVector * regionSparse,
544 CoinIndexedVector * regionSparse2) const
545{
546 assert (numberRows_==numberColumns_);
547 double *region2 = regionSparse2->denseVector ( );
548 int *regionIndex2 = regionSparse2->getIndices ( );
549 int numberNonZero = regionSparse2->getNumElements ( );
550 //double *region = regionSparse->denseVector ( );
551 /*int *regionIndex = regionSparse->getIndices ( );*/
552 const int * permuteIn = factInfo_.mpermu+1;
553 factInfo_.packedMode = regionSparse2->packedMode() ? 1 : 0;
554 // Use region instead of dpermu
555 double * save =factInfo_.kadrpm;
556 factInfo_.kadrpm=regionSparse->denseVector()-1;
557 // use internal one for now (address is one off)
558 double * region = factInfo_.kadrpm;
559 if (numberNonZero<2) {
560 if (numberNonZero) {
561 int ipivrw=regionIndex2[0];
562 if (factInfo_.packedMode) {
563 double value=region2[0];
564 region2[0]=0.0;
565 region2[ipivrw]=value;
566 }
567 numberNonZero=c_ekkbtrn_ipivrw(&factInfo_, region2-1,
568 regionIndex2-1,ipivrw+1,
569 reinterpret_cast<int *>(factInfo_.kp1adr));
570 }
571 } else {
572#ifndef NDEBUG
573 {
574 int *mcstrt = factInfo_.xcsadr;
575 int * hpivco_new=factInfo_.kcpadr+1;
576 int nrow=factInfo_.nrow;
577 int i;
578 int ipiv = hpivco_new[0];
579 int last = mcstrt[ipiv];
580 for (i=0;i<nrow-1;i++) {
581 ipiv=hpivco_new[ipiv];
582 assert (mcstrt[ipiv]>last);
583 last=mcstrt[ipiv];
584 }
585 }
586#endif
587 int iSmallest = COIN_INT_MAX;
588 int iPiv=0;
589 const int *mcstrt = factInfo_.xcsadr;
590 // permute and save where nonzeros are
591 if (!factInfo_.packedMode) {
592 if ((numberRows_<200||(numberNonZero<<4)>numberRows_)) {
593 for (int j=0;j<numberNonZero;j++) {
594 int jRow = regionIndex2[j];
595 int iRow = permuteIn[jRow];
596 regionIndex2[j]=iRow;
597 region[iRow]=region2[jRow];
598 region2[jRow]=0.0;
599 }
600 } else {
601 for (int j=0;j<numberNonZero;j++) {
602 int jRow = regionIndex2[j];
603 int iRow = permuteIn[jRow];
604 regionIndex2[j]=iRow;
605 region[iRow]=region2[jRow];
606 if (mcstrt[iRow]<iSmallest) {
607 iPiv=iRow;
608 iSmallest=mcstrt[iRow];
609 }
610 region2[jRow]=0.0;
611 }
612 }
613 } else {
614 for (int j=0;j<numberNonZero;j++) {
615 int jRow = regionIndex2[j];
616 int iRow = permuteIn[jRow];
617 regionIndex2[j]=iRow;
618 region[iRow]=region2[j];
619 region2[j]=0.0;
620 }
621 }
622 assert (iPiv>=0);
623 numberNonZero=c_ekkbtrn(&factInfo_, region2-1,regionIndex2-1,iPiv);
624 }
625 factInfo_.kadrpm=save;
626 factInfo_.packedMode=0;
627 regionSparse2->setNumElements(numberNonZero);
628 return 0;
629}
630// Number of entries in each row
631int *
632CoinOslFactorization::numberInRow() const
633{ return reinterpret_cast<int *> (factInfo_.xrnadr+1);}
634// Number of entries in each column
635int *
636CoinOslFactorization::numberInColumn() const
637{ return reinterpret_cast<int *> (factInfo_.xcnadr+1);}
638// Returns array to put basis starts in
639CoinBigIndex *
640CoinOslFactorization::starts() const
641{ return reinterpret_cast<CoinBigIndex *> (factInfo_.xcsadr+1);}
642// Returns array to put basis elements in
643CoinFactorizationDouble *
644CoinOslFactorization::elements() const
645{ return factInfo_.xeeadr+1;}
646// Returns pivot row
647int *
648CoinOslFactorization::pivotRow() const
649{ return factInfo_.krpadr+1;}
650// Returns work area
651CoinFactorizationDouble *
652CoinOslFactorization::workArea() const
653{ return factInfo_.kw1adr;}
654// Returns int work area
655int *
656CoinOslFactorization::intWorkArea() const
657{ return reinterpret_cast<int *> (factInfo_.kw1adr);}
658// Returns permute back
659int *
660CoinOslFactorization::permuteBack() const
661{ return factInfo_.kcpadr+1;}
662// Returns array to put basis indices in
663int *
664CoinOslFactorization::indices() const
665{ return factInfo_.xeradr+1;}
666// Returns true if wants tableauColumn in replaceColumn
667bool
668CoinOslFactorization::wantsTableauColumn() const
669{ return false;}
670/* Useful information for factorization
671 0 - iteration number
672 whereFrom is 0 for factorize and 1 for replaceColumn
673*/
674#ifdef CLP_REUSE_ETAS
675void
676CoinOslFactorization::setUsefulInformation(const int * info,int whereFrom)
677{
678 factInfo_.iterno=info[0];
679 if (whereFrom) {
680 factInfo_.reintro=-1;
681 if( factInfo_.first_dense>=factInfo_.last_dense) {
682 int * putSeq = factInfo_.xrsadr+2*factInfo_.nrowmx+2;
683 int * position = putSeq+factInfo_.maxinv;
684 //int * putStart = position+factInfo_.maxinv;
685 int iSequence=info[1];
686 if (whereFrom==1) {
687 putSeq[factInfo_.npivots]=iSequence;
688 } else {
689 int i;
690 for (i=factInfo_.npivots-1;i>=0;i--) {
691 if (putSeq[i]==iSequence)
692 break;
693 }
694 if (i>=0) {
695 factInfo_.reintro=position[i];
696 } else {
697 factInfo_.reintro=-1;
698 }
699 factInfo_.nnentu=factInfo_.save_nnentu;
700 }
701 }
702 }
703}
704#else
705void
706CoinOslFactorization::setUsefulInformation(const int * info,int /*whereFrom*/)
707{ factInfo_.iterno=info[0]; }
708#endif
709
710// Get rid of all memory
711void
712CoinOslFactorization::clearArrays()
713{
714 factInfo_.nR_etas=0;
715 factInfo_.nnentu=0;
716 factInfo_.nnentl=0;
717 maximumRows_=0;
718 maximumSpace_=0;
719 factInfo_.last_eta_size=0;
720 gutsOfDestructor(false);
721}
722void
723CoinOslFactorization::maximumPivots ( int value )
724{
725 maximumPivots_ = value;
726}
727#define CLP_FILL 15
728/*#undef NDEBUG*/
729//#define CLP_DEBUG_MALLOC 1000000
730#if CLP_DEBUG_MALLOC
731static int malloc_number=0;
732static int malloc_check=-1;
733static int malloc_counts_on=0;
734struct malloc_struct {
735 void * previous;
736 void * next;
737 int size;
738 int when;
739 int type;
740};
741static double malloc_times=0.0;
742static double malloc_total=0.0;
743static double malloc_current=0.0;
744static double malloc_max=0.0;
745static int malloc_amount[]={0,32,128,256,1024,4096,16384,65536,262144,
746 2000000000};
747static int malloc_n=10;
748double malloc_counts[10]={0,0,0,0,0,0,0,0,0,0};
749typedef struct malloc_struct malloc_struct;
750static malloc_struct startM = {0,0,0,0};
751static malloc_struct endM = {0,0,0,0};
752static int extra=4;
753void clp_memory(int type)
754{
755 if (type==0) {
756 /* switch on */
757 malloc_counts_on=1;
758 startM.next=&endM;
759 endM.previous=&startM;
760 } else {
761 /* summary */
762 double average = malloc_total/malloc_times;
763 int i;
764 malloc_struct * previous = (malloc_struct *) endM.previous;
765 printf("count %g bytes %g - average %g\n",malloc_times,malloc_total,average);
766 printf("current bytes %g - maximum %g\n",malloc_current,malloc_max);
767
768 for ( i=0;i<malloc_n;i++)
769 printf("%g ",malloc_counts[i]);
770 printf("\n");
771 malloc_counts_on=0;
772 if (previous->previous!=&startM) {
773 int n=0;
774 printf("Allocated blocks\n");
775 while (previous->previous!=&startM) {
776 printf("(%d at %d) ",previous->size,previous->when);
777 n++;
778 if ((n%5)==0)
779 printf("\n");
780 previous = (malloc_struct *) previous->previous;
781 }
782 printf("\n - total %d\n",n);
783 }
784 }
785 malloc_number=0;
786 malloc_times=0.0;
787 malloc_total=0.0;
788 malloc_current=0.0;
789 malloc_max=0.0;
790 memset(malloc_counts,0,sizeof(malloc_counts));
791}
792static void clp_adjust(void * temp,int size,int type)
793{
794 malloc_struct * itemp = (malloc_struct *) temp;
795 malloc_struct * first = (malloc_struct *) startM.next;
796 int i;
797 startM.next=temp;
798 first->previous=temp;
799 itemp->previous=&startM;
800 itemp->next=first;
801 malloc_number++;
802 if (malloc_number==malloc_check) {
803 printf("Allocation of %d bytes at %d (type %d) not freed\n",
804 size,malloc_number,type);
805 }
806 itemp->when=malloc_number;
807 itemp->size=size;
808 itemp->type=type;
809 malloc_times ++;
810 malloc_total += size;
811 malloc_current += size;
812 malloc_max=CoinMax(malloc_max,malloc_current);
813 for (i=0;i<malloc_n;i++) {
814 if ((int) size<=malloc_amount[i]) {
815 malloc_counts[i]++;
816 break;
817 }
818 }
819}
820#endif
821/* covers */
822double * clp_double(int number_entries)
823{
824#if CLP_DEBUG_MALLOC==0
825 return reinterpret_cast<double *>( malloc(number_entries*sizeof(double)));
826#else
827 double * temp = reinterpret_cast<double *>( malloc((number_entries+extra)*sizeof(double)));
828 clp_adjust(temp,number_entries*sizeof(double),1);
829#if CLP_DEBUG_MALLOC>1
830 if (number_entries*sizeof(double)>=CLP_DEBUG_MALLOC)
831 printf("WWW %x malloced by double %d - size %d\n",
832 temp+extra,malloc_number,number_entries);
833#endif
834 return temp+extra;
835#endif
836}
837int * clp_int(int number_entries)
838{
839#if CLP_DEBUG_MALLOC==0
840 return reinterpret_cast<int *>( malloc(number_entries*sizeof(int)));
841#else
842 double * temp = reinterpret_cast<double *>( malloc(((number_entries+1)/2+extra)*sizeof(double)));
843 clp_adjust(temp,number_entries*sizeof(int),2);
844#if CLP_DEBUG_MALLOC>1
845 if (number_entries*sizeof(int)>=CLP_DEBUG_MALLOC)
846 printf("WWW %x malloced by int %d - size %d\n",
847 temp+extra,malloc_number,number_entries);
848#endif
849 return reinterpret_cast<int *>( (temp+extra));
850#endif
851}
852void * clp_malloc(int number_entries)
853{
854#if CLP_DEBUG_MALLOC==0
855 return malloc(number_entries);
856#else
857 double * temp = reinterpret_cast<double *>( malloc(number_entries+extra*sizeof(double)));
858 clp_adjust(temp,number_entries,0);
859#if CLP_DEBUG_MALLOC>1
860 if (number_entries>=CLP_DEBUG_MALLOC)
861 printf("WWW %x malloced by void %d - size %d\n",
862 temp+extra,malloc_number,number_entries);
863#endif
864 return (void *) (temp+extra);
865#endif
866}
867void clp_free(void * oldArray)
868{
869#if CLP_DEBUG_MALLOC==0
870 free(oldArray);
871#else
872 if (oldArray) {
873 double * temp = (reinterpret_cast<double *>( oldArray)-extra);
874 malloc_struct * itemp = (malloc_struct *) temp;
875 malloc_struct * next = (malloc_struct *) itemp->next;
876 malloc_struct * previous = (malloc_struct *) itemp->previous;
877 previous->next=next;
878 next->previous=previous;
879 malloc_current -= itemp->size;
880#if CLP_DEBUG_MALLOC>1
881 if (itemp->size>=CLP_DEBUG_MALLOC)
882 printf("WWW %x freed by free %d - old length %d - type %d\n",
883 oldArray,itemp->when,itemp->size,itemp->type);
884#endif
885 free(temp);
886 }
887#endif
888}
889/*#define FIX_ADD 4*nrowmx+5
890 #define FIX_ADD2 4*nrowmx+5*/
891#define FIX_ADD 1*nrowmx+5
892#define FIX_ADD2 1*nrowmx+5
893#define ALIGNMENT 32
894static void * clp_align (void * memory)
895{
896 if (sizeof(int)==sizeof(void *)&&ALIGNMENT) {
897 CoinInt64 k = reinterpret_cast<CoinInt64> (memory);
898 if ((k&(ALIGNMENT-1))!=0) {
899 k &= ~(ALIGNMENT-1);
900 k += ALIGNMENT;
901 memory = reinterpret_cast<void *> (k);
902 }
903 return memory;
904 } else {
905 return memory;
906 }
907}
908void clp_setup_pointers(EKKfactinfo * fact)
909{
910 /* do extra stuff */
911 int nrow=fact->nrow;
912 int nrowmx=fact->nrowmx;
913 int maxinv=fact->maxinv;
914 fact->lstart = nrow + maxinv + 5;
915 /* this is the number of L transforms */
916 fact->xnetalval = fact->xnetal - fact->lstart;
917 fact->mpermu = (reinterpret_cast<int*> (fact->kadrpm+nrow))+1;
918 fact->bitArray = fact->krpadr + ( nrowmx+2);
919 fact->back = fact->kcpadr+2*nrow + maxinv + 4;
920 fact->hpivcoR = fact->kcpadr+nrow+3;
921 fact->nonzero = (reinterpret_cast<char *>( &fact->mpermu[nrow+1]))-1;
922}
923#ifndef NDEBUG
924int ets_count=0;
925int ets_check=-1;
926//static int adjust_count=0;
927//static int adjust_check=-1;
928#endif
929static void clp_adjust_pointers(EKKfactinfo * fact, int adjust)
930{
931#if 0 //ndef NDEBUG
932 adjust_count++;
933 if (adjust_check>=0&&adjust_count>=adjust_check) {
934 printf("trouble\n");
935 }
936#endif
937 if (fact->trueStart) {
938 fact->kadrpm += adjust;
939 fact->krpadr += adjust;
940 fact->kcpadr += adjust;
941 fact->xrsadr += adjust;
942 fact->xcsadr += adjust;
943 fact->xrnadr += adjust;
944 fact->xcnadr += adjust;
945 }
946 if (fact->xeradr) {
947 fact->xeradr += adjust;
948 fact->xecadr += adjust;
949 fact->xeeadr += adjust;
950 }
951}
952/* deals with memory for complicated array
953 0 just do addresses
954 1 just get memory */
955static double *
956clp_alloc_memory(EKKfactinfo * fact,int type, int * length)
957{
958 int nDouble=0;
959 int nInt=0;
960 int nrowmxp;
961 int ntot1;
962 int ntot2;
963 int ntot3;
964 int nrowmx;
965 int * tempI;
966 double * tempD;
967 nrowmx=fact->nrowmx;
968 nrowmxp = nrowmx + 2;
969 ntot1 = nrowmxp;
970 ntot2 = 3*nrowmx+5; /* space for three lists */
971 ntot3 = 2*nrowmx;
972 if ((ntot1<<1)<ntot2) {
973 ntot1=ntot2>>1;
974 }
975 ntot3=CoinMax(ntot3,ntot1);
976 /* Row work regions */
977 /* must be contiguous so allocate as one chunk */
978 /* may only need 2.5 */
979 /* now doing all at once - far too much - reduce later */
980 tempD=fact->kw1adr;
981 tempD+=nrowmxp;
982 tempD = reinterpret_cast<double *>( clp_align(tempD));
983 fact->kw2adr=tempD;
984 tempD+=nrowmxp;
985 tempD = reinterpret_cast<double *>( clp_align(tempD));
986 fact->kw3adr=tempD-1;
987 tempD+=nrowmxp;
988 tempD = reinterpret_cast<double *>( clp_align(tempD));
989 fact->kp1adr=reinterpret_cast<EKKHlink *>(tempD);
990 tempD+=nrowmxp;
991 tempD = reinterpret_cast<double *>( clp_align(tempD));
992 fact->kp2adr=reinterpret_cast<EKKHlink *>(tempD);
993 //tempD+=ntot3;
994 tempD+=nrowmxp;
995 tempD = reinterpret_cast<double *>( clp_align(tempD));
996 /*printf("zz %x %x\n",tempD,fact->kadrpm);*/
997 fact->kadrpm = tempD;
998 /* seems a lot */
999 tempD += ((6*nrowmx +8)*(sizeof(int))/sizeof(double));
1000 /* integer arrays */
1001 tempI = reinterpret_cast<int *>( tempD);
1002 tempI = reinterpret_cast<int *>( clp_align(tempI));
1003 fact->xrsadr = tempI;
1004#ifdef CLP_REUSE_ETAS
1005 tempI +=( 3*(nrowmx+fact->maxinv+1));
1006#else
1007 tempI +=( (nrowmx<<1)+fact->maxinv+1);
1008#endif
1009 tempI = reinterpret_cast<int *>( clp_align(tempI));
1010 fact->xcsadr = tempI;
1011#if 1 //def CLP_REUSE_ETAS
1012 tempI += ( 2*nrowmx+8+2*fact->maxinv);
1013#else
1014 tempI += ( 2*nrowmx+8+fact->maxinv);
1015#endif
1016 tempI += FIX_ADD+FIX_ADD2;
1017 tempI = reinterpret_cast<int *>( clp_align(tempI));
1018 fact->xrnadr = tempI;
1019 tempI += nrowmx;
1020 tempI = reinterpret_cast<int *>( clp_align(tempI));
1021 fact->xcnadr = tempI;
1022 tempI += nrowmx;
1023 tempI = reinterpret_cast<int *>( clp_align(tempI));
1024 fact->krpadr = tempI;
1025 tempI += ( nrowmx+1) +((nrowmx+33)>>5);
1026 /*printf("zzz %x %x\n",tempI,fact->kcpadr);*/
1027 tempI = reinterpret_cast<int *>( clp_align(tempI));
1028 fact->kcpadr = tempI;
1029 tempI += 3*nrowmx+8+fact->maxinv;
1030 fact->R_etas_start = fact->xcsadr+nrowmx+fact->maxinv+4;
1031 fact->R_etas_start += FIX_ADD;
1032 nInt = static_cast<int>(tempI-(reinterpret_cast<int *>( fact->trueStart)));
1033 nDouble = static_cast<int>(sizeof(int)*(nInt+1)/sizeof(double));
1034 *length = nDouble;
1035 /*printf("nDouble %d - type %d\n",nDouble,type);*/
1036 nDouble += static_cast<int>((2*ALIGNMENT)/sizeof(double));
1037 if (type) {
1038 /*printf("%d allocated\n",nDouble);*/
1039 tempD = reinterpret_cast<double *>( clp_double(nDouble));
1040#ifndef NDEBUG
1041 memset(tempD,CLP_FILL,nDouble*sizeof(double));
1042#endif
1043 }
1044 return tempD;
1045}
1046static void c_ekksmem(EKKfactinfo *fact,int nrow,int maximumPivots)
1047{
1048 /* space for invert */
1049 int nnetas=fact->eta_size;
1050 fact->nrow=nrow;
1051 if (!(nnetas>fact->last_eta_size||(!fact->xe2adr&&fact->if_sparse_update)||
1052 nrow>fact->nrowmx||maximumPivots>fact->maxinv))
1053 return;
1054 clp_adjust_pointers(fact, +1);
1055 if (nrow>fact->nrowmx||maximumPivots>fact->maxinv) {
1056 int length;
1057 fact->nrowmx=CoinMax(nrow,fact->nrowmx);
1058 fact->maxinv=CoinMax(maximumPivots,fact->maxinv);
1059 clp_free(fact->trueStart);
1060 fact->trueStart=0;
1061 fact->kw1adr=0;
1062 fact->trueStart=clp_alloc_memory(fact,1,&length);
1063 fact->kw1adr=reinterpret_cast<double *>( clp_align(fact->trueStart));
1064 clp_alloc_memory(fact,0,&length);
1065 }
1066 /*if (!fact->iterno) fact->eta_size+=1000000;*//* TEMP*/
1067 if (nnetas>fact->last_eta_size||(!fact->xe2adr&&fact->if_sparse_update)) {
1068 fact->last_eta_size = nnetas;
1069 clp_free(reinterpret_cast<char *>(fact->xe2adr));
1070 /* if malloc fails - we have lost memory - start again */
1071 if (!fact->ndenuc &&fact->if_sparse_update) {
1072 /* allow second copy of elements */
1073 fact->xe2adr = clp_double(nnetas);
1074#ifndef NDEBUG
1075 memset(fact->xe2adr,CLP_FILL,nnetas*sizeof(double));
1076#endif
1077 if (!fact->xe2adr) {
1078 fact->maxNNetas=fact->last_eta_size; /* dont allow any increase */
1079 nnetas=fact->last_eta_size;
1080 fact->eta_size=nnetas;
1081#ifdef PRINT_DEBUG
1082 if (fact->if_sparse_update) {
1083 printf("*** Sparse update off due to memory\n");
1084 }
1085#endif
1086 fact->if_sparse_update=0;
1087 fact->switch_off_sparse_update=1;
1088 }
1089 } else {
1090 fact->xe2adr = 0;
1091 fact->if_sparse_update=0;
1092 }
1093 clp_free(fact->xeradr);
1094 fact->xeradr= clp_int( nnetas);
1095#ifndef NDEBUG
1096 memset(fact->xeradr,CLP_FILL,nnetas*sizeof(int));
1097#endif
1098 if (!fact->xeradr) {
1099 nnetas=0;
1100 }
1101 if (nnetas) {
1102 clp_free(fact->xecadr);
1103 fact->xecadr= clp_int( nnetas);
1104#ifndef NDEBUG
1105 memset(fact->xecadr,CLP_FILL,nnetas*sizeof(int));
1106#endif
1107 if (!fact->xecadr) {
1108 nnetas=0;
1109 }
1110 }
1111 if (nnetas) {
1112 clp_free(fact->xeeadr);
1113 fact->xeeadr= clp_double(nnetas);
1114#ifndef NDEBUG
1115 memset(fact->xeeadr,CLP_FILL,nnetas*sizeof(double));
1116#endif
1117 if (!fact->xeeadr) {
1118 nnetas=0;
1119 }
1120 }
1121 }
1122 if (!nnetas) {
1123 char msg[100];
1124 sprintf(msg,"Unable to allocate factorization memory for %d elements",
1125 nnetas);
1126 throw(msg);
1127 }
1128 /*c_ekklplp->nnetas=nnetas;*/
1129 fact->nnetas=nnetas;
1130 clp_adjust_pointers(fact, -1);
1131}
1132static void c_ekksmem_copy(EKKfactinfo *fact,const EKKfactinfo * rhsFact)
1133{
1134 /* space for invert */
1135 int nrowmx=rhsFact->nrowmx,nnetas=rhsFact->nnetas;
1136 int canReuseEtas= (fact->eta_size==rhsFact->eta_size) ? 1 : 0;
1137 int canReuseArrays = (fact->nrowmx==rhsFact->nrowmx) ? 1 : 0;
1138 clp_adjust_pointers(fact, +1);
1139 clp_adjust_pointers(const_cast<EKKfactinfo *>(rhsFact), +1);
1140 /*memset(fact,0,sizeof(EKKfactinfo));*/
1141 /* copy scalars */
1142 memcpy(&fact->drtpiv,&rhsFact->drtpiv,5*sizeof(double));
1143 memcpy(&fact->nrow,&rhsFact->nrow,((&fact->maxNNetas-&fact->nrow)+1)*
1144 sizeof(int));
1145 if (nrowmx) {
1146 int length;
1147 int kCopyEnd,nCopyEnd,nCopyStart;
1148 if (!canReuseEtas) {
1149 clp_free(fact->xeradr);
1150 clp_free(fact->xecadr);
1151 clp_free(fact->xeeadr);
1152 clp_free(fact->xe2adr);
1153 fact->xeradr = 0;
1154 fact->xecadr = 0;
1155 fact->xeeadr = 0;
1156 fact->xe2adr = 0;
1157 }
1158 if (!canReuseArrays) {
1159 clp_free(fact->trueStart);
1160 fact->trueStart=0;
1161 fact->kw1adr=0;
1162 fact->trueStart=clp_alloc_memory(fact,1,&length);
1163 fact->kw1adr=reinterpret_cast<double *>( clp_align(fact->trueStart));
1164 }
1165 clp_alloc_memory(fact,0,&length);
1166 nnetas=fact->eta_size;
1167 assert (nnetas);
1168 {
1169 int n2 = rhsFact->nR_etas;
1170 int n3 = n2 ? rhsFact->R_etas_start[1+n2]: 0;
1171 int * startR = rhsFact->R_etas_index+n3;
1172 nCopyEnd=static_cast<int>((rhsFact->xeradr+nnetas)-startR);
1173 nCopyStart=rhsFact->nnentu;
1174 nCopyEnd = CoinMin(nCopyEnd+20,nnetas);
1175 kCopyEnd = nnetas-nCopyEnd;
1176 nCopyStart = CoinMin(nCopyStart+20,nnetas);
1177 if (!n2&&!rhsFact->nnentu&&!rhsFact->nnentl) {
1178 nCopyStart=nCopyEnd=0;
1179 }
1180 }
1181 /* copy */
1182 if(nCopyStart||nCopyEnd||true) {
1183#if 1
1184 memcpy(fact->kw1adr,rhsFact->kw1adr,length*sizeof(double));
1185#else
1186 c_ekkscpy((length*sizeof(double))/sizeof(int),
1187 reinterpret_cast<int *>( rhsFact->kw1adr,reinterpret_cast<int *>( fact->kw1adr));
1188#endif
1189 }
1190 /* if malloc fails - we have lost memory - start again */
1191 if (!fact->ndenuc &&fact->if_sparse_update) {
1192 /* allow second copy of elements */
1193 if (!canReuseEtas)
1194 fact->xe2adr = clp_double(nnetas);
1195 if (!fact->xe2adr) {
1196 fact->maxNNetas=nnetas; /* dont allow any increase */
1197#ifdef PRINT_DEBUG
1198 if (fact->if_sparse_update) {
1199 printf("*** Sparse update off due to memory\n");
1200 }
1201#endif
1202 fact->if_sparse_update=0;
1203 } else {
1204#ifndef NDEBUG
1205 memset(fact->xe2adr,CLP_FILL,nnetas*sizeof(double));
1206#endif
1207 }
1208 } else {
1209 clp_free(fact->xe2adr);
1210 fact->xe2adr = 0;
1211 fact->if_sparse_update=0;
1212 }
1213
1214 if (!canReuseEtas)
1215 fact->xeradr= clp_int(nnetas);
1216 if (!fact->xeradr) {
1217 nnetas=0;
1218 } else {
1219#ifndef NDEBUG
1220 memset(fact->xeradr,CLP_FILL,nnetas*sizeof(int));
1221#endif
1222
1223 /* copy */
1224 if(nCopyStart||nCopyEnd) {
1225#if 0
1226 memcpy(fact->xeradr,rhsFact->xeradr,nCopyStart*sizeof(int));
1227 memcpy(fact->xeradr+kCopyEnd,rhsFact->xeradr+kCopyEnd,nCopyEnd*sizeof(int));
1228#else
1229 c_ekkscpy(nCopyStart,rhsFact->xeradr,fact->xeradr);
1230 c_ekkscpy(nCopyEnd,rhsFact->xeradr+kCopyEnd,fact->xeradr+kCopyEnd);
1231#endif
1232 }
1233 }
1234 if (nnetas) {
1235 if (!canReuseEtas)
1236 fact->xecadr= clp_int(nnetas);
1237 if (!fact->xecadr) {
1238 nnetas=0;
1239 } else {
1240#ifndef NDEBUG
1241 memset(fact->xecadr,CLP_FILL,nnetas*sizeof(int));
1242#endif
1243 /* copy */
1244 if (fact->rows_ok&&(nCopyStart||nCopyEnd)) {
1245 int i;
1246 int * hcoliR=rhsFact->xecadr-1;
1247 int * hcoli=fact->xecadr-1;
1248 int * mrstrt=fact->xrsadr;
1249 int * hinrow=fact->xrnadr;
1250#if 0
1251 memcpy(fact->xecadr+kCopyEnd,rhsFact->xecadr+kCopyEnd,
1252 nCopyEnd*sizeof(int));
1253#else
1254 c_ekkscpy(nCopyEnd,rhsFact->xecadr+kCopyEnd,fact->xecadr+kCopyEnd);
1255#endif
1256 if (!fact->xe2adr) {
1257 for (i=0;i<fact->nrow;i++) {
1258 int istart = mrstrt[i];
1259 assert (istart>0&&istart<=nnetas);
1260 assert (hinrow[i]>=0&&hinrow[i]<=fact->nrow);
1261 memcpy(hcoli+istart,hcoliR+istart,hinrow[i]*sizeof(int));
1262 }
1263 } else {
1264 double * de2valR=rhsFact->xe2adr-1;
1265 double * de2val=fact->xe2adr-1;
1266#if 0
1267 memcpy(fact->xe2adr+kCopyEnd,rhsFact->xe2adr+kCopyEnd,
1268 nCopyEnd*sizeof(double));
1269#else
1270 c_ekkdcpy(nCopyEnd,rhsFact->xe2adr+kCopyEnd
1271 ,fact->xe2adr+kCopyEnd);
1272#endif
1273 for (i=0;i<fact->nrow;i++) {
1274 int istart = mrstrt[i];
1275 assert (istart>0&&istart<=nnetas);
1276 assert (hinrow[i]>=0&&hinrow[i]<=fact->nrow);
1277 memcpy(hcoli+istart,hcoliR+istart,hinrow[i]*sizeof(int));
1278 memcpy(de2val+istart,de2valR+istart,hinrow[i]*sizeof(double));
1279#ifndef NDEBUG
1280 {
1281 int j;
1282 for (j=istart;j<istart+hinrow[i];j++)
1283 assert (fabs(de2val[j])<1.0e50);
1284 }
1285#endif
1286 }
1287 }
1288 }
1289 }
1290 }
1291 if (nnetas) {
1292 if (!canReuseEtas)
1293 fact->xeeadr= clp_double(nnetas);
1294 if (!fact->xeeadr) {
1295 nnetas=0;
1296 } else {
1297#ifndef NDEBUG
1298 memset(fact->xeeadr,CLP_FILL,nnetas*sizeof(double));
1299#endif
1300 /* copy */
1301 if(nCopyStart||nCopyEnd) {
1302#if 0
1303 memcpy(fact->xeeadr,rhsFact->xeeadr,nCopyStart*sizeof(double));
1304 memcpy(fact->xeeadr+kCopyEnd,rhsFact->xeeadr+kCopyEnd,nCopyEnd*sizeof(double));
1305#else
1306 c_ekkdcpy(nCopyStart,
1307 rhsFact->xeeadr,fact->xeeadr);
1308 c_ekkdcpy(nCopyEnd,
1309 rhsFact->xeeadr+kCopyEnd,
1310 fact->xeeadr+kCopyEnd);
1311#endif
1312 }
1313 /*fact->R_etas_index = &XERADR1()[kdnspt - 1];
1314 fact->R_etas_element = &XEEADR1()[kdnspt - 1];*/
1315 fact->R_etas_start = fact->xcsadr+
1316 (rhsFact->R_etas_start-rhsFact->xcsadr);
1317 fact->R_etas_index = fact->xeradr+
1318 (rhsFact->R_etas_index-rhsFact->xeradr);
1319 fact->R_etas_element = fact->xeeadr+
1320 (rhsFact->R_etas_element-rhsFact->xeeadr);
1321 }
1322 }
1323 }
1324 assert (nnetas||!nrowmx);
1325 fact->nnetas=nnetas;
1326 clp_adjust_pointers(fact, -1);
1327 clp_setup_pointers(fact);
1328 clp_adjust_pointers(const_cast<EKKfactinfo *>(rhsFact), -1);
1329}
1330static void c_ekksmem_delete(EKKfactinfo *fact)
1331{
1332 clp_adjust_pointers(fact, +1);
1333 clp_free(fact->trueStart);
1334 clp_free(fact->xe2adr);
1335 clp_free(fact->xecadr);
1336 clp_free(fact->xeradr);
1337 clp_free(fact->xeeadr);
1338 fact->eta_size=0;
1339 fact->xrsadr = 0;
1340 fact->xcsadr = 0;
1341 fact->xrnadr = 0;
1342 fact->xcnadr = 0;
1343 fact->krpadr = 0;
1344 fact->kcpadr = 0;
1345 fact->xeradr = 0;
1346 fact->xecadr = 0;
1347 fact->xeeadr = 0;
1348 fact->xe2adr = 0;
1349 fact->trueStart = 0;
1350 fact->kw2adr = 0;
1351 fact->kw3adr = 0;
1352 fact->kp1adr = 0;
1353 fact->kp2adr = 0;
1354 fact->kadrpm = 0;
1355 fact->kw1adr = 0;
1356}
1357int c_ekk_IsSet(const int * array,int bit);
1358void c_ekk_Set(int * array,int bit);
1359void c_ekk_Unset(int * array,int bit);
1360int c_ekk_IsSet(const int * array,int bit)
1361{
1362 int iWord = bit>>5;
1363 int iBit = bit&31;
1364 int word = array[iWord];
1365 return (word&(1<<iBit))!=0;
1366}
1367void c_ekk_Set(int * array,int bit)
1368{
1369 int iWord = bit>>5;
1370 int iBit = bit&31;
1371 array[iWord] |= (1<<iBit);
1372}
1373void c_ekk_Unset(int * array,int bit)
1374{
1375 int iWord = bit>>5;
1376 int iBit = bit&31;
1377 array[iWord] &= ~(1<<iBit);
1378}
1379int CoinOslFactorization::factorize (
1380 const CoinPackedMatrix & matrix,
1381 int rowIsBasic[],
1382 int columnIsBasic[],
1383 double areaFactor )
1384{
1385 setSolveMode(10);
1386 if (areaFactor)
1387 factInfo_.areaFactor = areaFactor;
1388 const int * row = matrix.getIndices();
1389 const CoinBigIndex * columnStart = matrix.getVectorStarts();
1390 const int * columnLength = matrix.getVectorLengths();
1391 const double * element = matrix.getElements();
1392 int numberRows=matrix.getNumRows();
1393 int numberColumns=matrix.getNumCols();
1394 int numberBasic = 0;
1395 CoinBigIndex numberElements=0;
1396 int numberRowBasic=0;
1397
1398 // compute how much in basis
1399
1400 int i;
1401 // Move pivot variables across if they look good
1402 int * pivotTemp = new int [numberRows];
1403
1404 for (i=0;i<numberRows;i++) {
1405 if (rowIsBasic[i]>=0)
1406 pivotTemp[numberRowBasic++]=i;
1407 }
1408
1409 numberBasic = numberRowBasic;
1410
1411 for (i=0;i<numberColumns;i++) {
1412 if (columnIsBasic[i]>=0) {
1413 pivotTemp[numberBasic++]=i;
1414 numberElements += columnLength[i];
1415 }
1416 }
1417 if ( numberBasic > numberRows ) {
1418 return -2; // say too many in basis
1419 }
1420 numberElements = 3 * numberRows + 3 * numberElements + 20000;
1421 setUsefulInformation(&numberRows,0);
1422 getAreas ( numberRows, numberRows, numberElements,
1423 2 * numberElements );
1424 //fill
1425 numberBasic=0;
1426 numberElements=0;
1427 // Fill in counts so we can skip part of preProcess
1428 double * elementU=elements();
1429 int * indexRowU=indices();
1430 int * startColumnU=starts();
1431 int * numberInRow=this->numberInRow();
1432 int * numberInColumn=this->numberInColumn();
1433 CoinZeroN ( numberInRow, numberRows );
1434 CoinZeroN ( numberInColumn, numberRows );
1435 for (i=0;i<numberRowBasic;i++) {
1436 int iRow = pivotTemp[i];
1437 // Change pivotTemp to correct sequence
1438 pivotTemp[i]=iRow+numberColumns;
1439 indexRowU[i]=iRow;
1440 startColumnU[i]=i;
1441 elementU[i]=-1.0;
1442 numberInRow[iRow]=1;
1443 numberInColumn[i]=1;
1444 }
1445 startColumnU[numberRowBasic]=numberRowBasic;
1446 numberElements=numberRowBasic;
1447 numberBasic=numberRowBasic;
1448 for (i=0;i<numberColumns;i++) {
1449 if (columnIsBasic[i]>=0) {
1450 CoinBigIndex j;
1451 for (j=columnStart[i];j<columnStart[i]+columnLength[i];j++) {
1452 int iRow=row[j];
1453 numberInRow[iRow]++;
1454 indexRowU[numberElements]=iRow;
1455 elementU[numberElements++]=element[j];
1456 }
1457 numberInColumn[numberBasic]=columnLength[i];
1458 numberBasic++;
1459 startColumnU[numberBasic]=numberElements;
1460 }
1461 }
1462 preProcess ( );
1463 factor ( );
1464 if (status() == 0) {
1465 int * pivotVariable = new int [numberRows];
1466 postProcess(pivotTemp,pivotVariable);
1467 for (i=0;i<numberRows;i++) {
1468 int iPivot=pivotVariable[i];
1469 if (iPivot<numberColumns) {
1470 assert (columnIsBasic[iPivot]>=0);
1471 columnIsBasic[iPivot]=i;
1472 } else {
1473 iPivot-=numberColumns;
1474 assert (rowIsBasic[iPivot]>=0);
1475 rowIsBasic[iPivot]=i;
1476 }
1477 }
1478 delete [] pivotVariable;
1479 }
1480 delete [] pivotTemp;
1481 return status_;
1482}
1483// Condition number - product of pivots after factorization
1484double
1485CoinOslFactorization::conditionNumber() const
1486{
1487 double condition = 1.0;
1488 const double *dluval = factInfo_.xeeadr+1-1; // stored before
1489 const int *mcstrt = factInfo_.xcsadr+1;
1490 for (int i=0;i<numberRows_;i++) {
1491 const int kx = mcstrt[i];
1492 const double dpiv = dluval[kx];
1493 condition *= dpiv;
1494 }
1495 condition = CoinMax(fabs(condition),1.0e-50);
1496 return 1.0/condition;
1497}
1498