1/* $Id: CoinFactorization.hpp 1448 2011-06-19 15:34:41Z stefan $ */
2// Copyright (C) 2002, 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/*
7 Authors
8
9 John Forrest
10
11 */
12#ifndef CoinFactorization_H
13#define CoinFactorization_H
14//#define COIN_ONE_ETA_COPY 100
15
16#include <iostream>
17#include <string>
18#include <cassert>
19#include <cstdio>
20#include <cmath>
21#include "CoinTypes.hpp"
22#include "CoinIndexedVector.hpp"
23
24class CoinPackedMatrix;
25/** This deals with Factorization and Updates
26
27 This class started with a parallel simplex code I was writing in the
28 mid 90's. The need for parallelism led to many complications and
29 I have simplified as much as I could to get back to this.
30
31 I was aiming at problems where I might get speed-up so I was looking at dense
32 problems or ones with structure. This led to permuting input and output
33 vectors and to increasing the number of rows each rank-one update. This is
34 still in as a minor overhead.
35
36 I have also put in handling for hyper-sparsity. I have taken out
37 all outer loop unrolling, dense matrix handling and most of the
38 book-keeping for slacks. Also I always use FTRAN approach to updating
39 even if factorization fairly dense. All these could improve performance.
40
41 I blame some of the coding peculiarities on the history of the code
42 but mostly it is just because I can't do elegant code (or useful
43 comments).
44
45 I am assuming that 32 bits is enough for number of rows or columns, but CoinBigIndex
46 may be redefined to get 64 bits.
47 */
48
49
50class CoinFactorization {
51 friend void CoinFactorizationUnitTest( const std::string & mpsDir );
52
53public:
54
55 /**@name Constructors and destructor and copy */
56 //@{
57 /// Default constructor
58 CoinFactorization ( );
59 /// Copy constructor
60 CoinFactorization ( const CoinFactorization &other);
61
62 /// Destructor
63 ~CoinFactorization ( );
64 /// Delete all stuff (leaves as after CoinFactorization())
65 void almostDestructor();
66 /// Debug show object (shows one representation)
67 void show_self ( ) const;
68 /// Debug - save on file - 0 if no error
69 int saveFactorization (const char * file ) const;
70 /** Debug - restore from file - 0 if no error on file.
71 If factor true then factorizes as if called from ClpFactorization
72 */
73 int restoreFactorization (const char * file , bool factor=false) ;
74 /// Debug - sort so can compare
75 void sort ( ) const;
76 /// = copy
77 CoinFactorization & operator = ( const CoinFactorization & other );
78 //@}
79
80 /**@name Do factorization */
81 //@{
82 /** When part of LP - given by basic variables.
83 Actually does factorization.
84 Arrays passed in have non negative value to say basic.
85 If status is okay, basic variables have pivot row - this is only needed
86 If status is singular, then basic variables have pivot row
87 and ones thrown out have -1
88 returns 0 -okay, -1 singular, -2 too many in basis, -99 memory */
89 int factorize ( const CoinPackedMatrix & matrix,
90 int rowIsBasic[], int columnIsBasic[] ,
91 double areaFactor = 0.0 );
92 /** When given as triplets.
93 Actually does factorization. maximumL is guessed maximum size of L part of
94 final factorization, maximumU of U part. These are multiplied by
95 areaFactor which can be computed by user or internally.
96 Arrays are copied in. I could add flag to delete arrays to save a
97 bit of memory.
98 If status okay, permutation has pivot rows - this is only needed
99 If status is singular, then basic variables have pivot row
100 and ones thrown out have -1
101 returns 0 -okay, -1 singular, -99 memory */
102 int factorize ( int numberRows,
103 int numberColumns,
104 CoinBigIndex numberElements,
105 CoinBigIndex maximumL,
106 CoinBigIndex maximumU,
107 const int indicesRow[],
108 const int indicesColumn[], const double elements[] ,
109 int permutation[],
110 double areaFactor = 0.0);
111 /** Two part version for maximum flexibility
112 This part creates arrays for user to fill.
113 estimateNumberElements is safe estimate of number
114 returns 0 -okay, -99 memory */
115 int factorizePart1 ( int numberRows,
116 int numberColumns,
117 CoinBigIndex estimateNumberElements,
118 int * indicesRow[],
119 int * indicesColumn[],
120 CoinFactorizationDouble * elements[],
121 double areaFactor = 0.0);
122 /** This is part two of factorization
123 Arrays belong to factorization and were returned by part 1
124 If status okay, permutation has pivot rows - this is only needed
125 If status is singular, then basic variables have pivot row
126 and ones thrown out have -1
127 returns 0 -okay, -1 singular, -99 memory */
128 int factorizePart2 (int permutation[],int exactNumberElements);
129 /// Condition number - product of pivots after factorization
130 double conditionNumber() const;
131
132 //@}
133
134 /**@name general stuff such as permutation or status */
135 //@{
136 /// Returns status
137 inline int status ( ) const {
138 return status_;
139 }
140 /// Sets status
141 inline void setStatus ( int value)
142 { status_=value; }
143 /// Returns number of pivots since factorization
144 inline int pivots ( ) const {
145 return numberPivots_;
146 }
147 /// Sets number of pivots since factorization
148 inline void setPivots ( int value )
149 { numberPivots_=value; }
150 /// Returns address of permute region
151 inline int *permute ( ) const {
152 return permute_.array();
153 }
154 /// Returns address of pivotColumn region (also used for permuting)
155 inline int *pivotColumn ( ) const {
156 return pivotColumn_.array();
157 }
158 /// Returns address of pivot region
159 inline CoinFactorizationDouble *pivotRegion ( ) const {
160 return pivotRegion_.array();
161 }
162 /// Returns address of permuteBack region
163 inline int *permuteBack ( ) const {
164 return permuteBack_.array();
165 }
166 /** Returns address of pivotColumnBack region (also used for permuting)
167 Now uses firstCount to save memory allocation */
168 inline int *pivotColumnBack ( ) const {
169 //return firstCount_.array();
170 return pivotColumnBack_.array();
171 }
172 /// Start of each row in L
173 inline CoinBigIndex * startRowL() const
174 { return startRowL_.array();}
175
176 /// Start of each column in L
177 inline CoinBigIndex * startColumnL() const
178 { return startColumnL_.array();}
179
180 /// Index of column in row for L
181 inline int * indexColumnL() const
182 { return indexColumnL_.array();}
183
184 /// Row indices of L
185 inline int * indexRowL() const
186 { return indexRowL_.array();}
187
188 /// Elements in L (row copy)
189 inline CoinFactorizationDouble * elementByRowL() const
190 { return elementByRowL_.array();}
191
192 /// Number of Rows after iterating
193 inline int numberRowsExtra ( ) const {
194 return numberRowsExtra_;
195 }
196 /// Set number of Rows after factorization
197 inline void setNumberRows(int value)
198 { numberRows_ = value; }
199 /// Number of Rows after factorization
200 inline int numberRows ( ) const {
201 return numberRows_;
202 }
203 /// Number in L
204 inline CoinBigIndex numberL() const
205 { return numberL_;}
206
207 /// Base of L
208 inline CoinBigIndex baseL() const
209 { return baseL_;}
210 /// Maximum of Rows after iterating
211 inline int maximumRowsExtra ( ) const {
212 return maximumRowsExtra_;
213 }
214 /// Total number of columns in factorization
215 inline int numberColumns ( ) const {
216 return numberColumns_;
217 }
218 /// Total number of elements in factorization
219 inline int numberElements ( ) const {
220 return totalElements_;
221 }
222 /// Length of FT vector
223 inline int numberForrestTomlin ( ) const {
224 return numberInColumn_.array()[numberColumnsExtra_];
225 }
226 /// Number of good columns in factorization
227 inline int numberGoodColumns ( ) const {
228 return numberGoodU_;
229 }
230 /// Whether larger areas needed
231 inline double areaFactor ( ) const {
232 return areaFactor_;
233 }
234 inline void areaFactor ( double value ) {
235 areaFactor_=value;
236 }
237 /// Returns areaFactor but adjusted for dense
238 double adjustedAreaFactor() const;
239 /// Allows change of pivot accuracy check 1.0 == none >1.0 relaxed
240 inline void relaxAccuracyCheck(double value)
241 { relaxCheck_ = value;}
242 inline double getAccuracyCheck() const
243 { return relaxCheck_;}
244 /// Level of detail of messages
245 inline int messageLevel ( ) const {
246 return messageLevel_ ;
247 }
248 void messageLevel ( int value );
249 /// Maximum number of pivots between factorizations
250 inline int maximumPivots ( ) const {
251 return maximumPivots_ ;
252 }
253 void maximumPivots ( int value );
254
255 /// Gets dense threshold
256 inline int denseThreshold() const
257 { return denseThreshold_;}
258 /// Sets dense threshold
259 inline void setDenseThreshold(int value)
260 { denseThreshold_ = value;}
261 /// Pivot tolerance
262 inline double pivotTolerance ( ) const {
263 return pivotTolerance_ ;
264 }
265 void pivotTolerance ( double value );
266 /// Zero tolerance
267 inline double zeroTolerance ( ) const {
268 return zeroTolerance_ ;
269 }
270 void zeroTolerance ( double value );
271#ifndef COIN_FAST_CODE
272 /// Whether slack value is +1 or -1
273 inline double slackValue ( ) const {
274 return slackValue_ ;
275 }
276 void slackValue ( double value );
277#endif
278 /// Returns maximum absolute value in factorization
279 double maximumCoefficient() const;
280 /// true if Forrest Tomlin update, false if PFI
281 inline bool forrestTomlin() const
282 { return doForrestTomlin_;}
283 inline void setForrestTomlin(bool value)
284 { doForrestTomlin_=value;}
285 /// True if FT update and space
286 inline bool spaceForForrestTomlin() const
287 {
288 CoinBigIndex start = startColumnU_.array()[maximumColumnsExtra_];
289 CoinBigIndex space = lengthAreaU_ - ( start + numberRowsExtra_ );
290 return (space>=0)&&doForrestTomlin_;
291 }
292 //@}
293
294 /**@name some simple stuff */
295 //@{
296
297 /// Returns number of dense rows
298 inline int numberDense() const
299 { return numberDense_;}
300
301 /// Returns number in U area
302 inline CoinBigIndex numberElementsU ( ) const {
303 return lengthU_;
304 }
305 /// Setss number in U area
306 inline void setNumberElementsU(CoinBigIndex value)
307 { lengthU_ = value; }
308 /// Returns length of U area
309 inline CoinBigIndex lengthAreaU ( ) const {
310 return lengthAreaU_;
311 }
312 /// Returns number in L area
313 inline CoinBigIndex numberElementsL ( ) const {
314 return lengthL_;
315 }
316 /// Returns length of L area
317 inline CoinBigIndex lengthAreaL ( ) const {
318 return lengthAreaL_;
319 }
320 /// Returns number in R area
321 inline CoinBigIndex numberElementsR ( ) const {
322 return lengthR_;
323 }
324 /// Number of compressions done
325 inline CoinBigIndex numberCompressions() const
326 { return numberCompressions_;}
327 /// Number of entries in each row
328 inline int * numberInRow() const
329 { return numberInRow_.array();}
330 /// Number of entries in each column
331 inline int * numberInColumn() const
332 { return numberInColumn_.array();}
333 /// Elements of U
334 inline CoinFactorizationDouble * elementU() const
335 { return elementU_.array();}
336 /// Row indices of U
337 inline int * indexRowU() const
338 { return indexRowU_.array();}
339 /// Start of each column in U
340 inline CoinBigIndex * startColumnU() const
341 { return startColumnU_.array();}
342 /// Maximum number of Columns after iterating
343 inline int maximumColumnsExtra()
344 { return maximumColumnsExtra_;}
345 /** L to U bias
346 0 - U bias, 1 - some U bias, 2 some L bias, 3 L bias
347 */
348 inline int biasLU() const
349 { return biasLU_;}
350 inline void setBiasLU(int value)
351 { biasLU_=value;}
352 /** Array persistence flag
353 If 0 then as now (delete/new)
354 1 then only do arrays if bigger needed
355 2 as 1 but give a bit extra if bigger needed
356 */
357 inline int persistenceFlag() const
358 { return persistenceFlag_;}
359 void setPersistenceFlag(int value);
360 //@}
361
362 /**@name rank one updates which do exist */
363 //@{
364
365 /** Replaces one Column to basis,
366 returns 0=OK, 1=Probably OK, 2=singular, 3=no room
367 If checkBeforeModifying is true will do all accuracy checks
368 before modifying factorization. Whether to set this depends on
369 speed considerations. You could just do this on first iteration
370 after factorization and thereafter re-factorize
371 partial update already in U */
372 int replaceColumn ( CoinIndexedVector * regionSparse,
373 int pivotRow,
374 double pivotCheck ,
375 bool checkBeforeModifying=false,
376 double acceptablePivot=1.0e-8);
377 /** Combines BtranU and delete elements
378 If deleted is NULL then delete elements
379 otherwise store where elements are
380 */
381 void replaceColumnU ( CoinIndexedVector * regionSparse,
382 CoinBigIndex * deleted,
383 int internalPivotRow);
384 //@}
385
386 /**@name various uses of factorization (return code number elements)
387 which user may want to know about */
388 //@{
389 /** Updates one column (FTRAN) from regionSparse2
390 Tries to do FT update
391 number returned is negative if no room
392 regionSparse starts as zero and is zero at end.
393 Note - if regionSparse2 packed on input - will be packed on output
394 */
395 int updateColumnFT ( CoinIndexedVector * regionSparse,
396 CoinIndexedVector * regionSparse2);
397 /** This version has same effect as above with FTUpdate==false
398 so number returned is always >=0 */
399 int updateColumn ( CoinIndexedVector * regionSparse,
400 CoinIndexedVector * regionSparse2,
401 bool noPermute=false) const;
402 /** Updates one column (FTRAN) from region2
403 Tries to do FT update
404 number returned is negative if no room.
405 Also updates region3
406 region1 starts as zero and is zero at end */
407 int updateTwoColumnsFT ( CoinIndexedVector * regionSparse1,
408 CoinIndexedVector * regionSparse2,
409 CoinIndexedVector * regionSparse3,
410 bool noPermuteRegion3=false) ;
411 /** Updates one column (BTRAN) from regionSparse2
412 regionSparse starts as zero and is zero at end
413 Note - if regionSparse2 packed on input - will be packed on output
414 */
415 int updateColumnTranspose ( CoinIndexedVector * regionSparse,
416 CoinIndexedVector * regionSparse2) const;
417 /** makes a row copy of L for speed and to allow very sparse problems */
418 void goSparse();
419 /** get sparse threshold */
420 inline int sparseThreshold ( ) const
421 { return sparseThreshold_;}
422 /** set sparse threshold */
423 void sparseThreshold ( int value );
424 //@}
425 /// *** Below this user may not want to know about
426
427 /**@name various uses of factorization (return code number elements)
428 which user may not want to know about (left over from my LP code) */
429 //@{
430 /// Get rid of all memory
431 inline void clearArrays()
432 { gutsOfDestructor();}
433 //@}
434
435 /**@name various updates - none of which have been written! */
436 //@{
437
438 /** Adds given elements to Basis and updates factorization,
439 can increase size of basis. Returns rank */
440 int add ( CoinBigIndex numberElements,
441 int indicesRow[],
442 int indicesColumn[], double elements[] );
443
444 /** Adds one Column to basis,
445 can increase size of basis. Returns rank */
446 int addColumn ( CoinBigIndex numberElements,
447 int indicesRow[], double elements[] );
448
449 /** Adds one Row to basis,
450 can increase size of basis. Returns rank */
451 int addRow ( CoinBigIndex numberElements,
452 int indicesColumn[], double elements[] );
453
454 /// Deletes one Column from basis, returns rank
455 int deleteColumn ( int Row );
456 /// Deletes one Row from basis, returns rank
457 int deleteRow ( int Row );
458
459 /** Replaces one Row in basis,
460 At present assumes just a singleton on row is in basis
461 returns 0=OK, 1=Probably OK, 2=singular, 3 no space */
462 int replaceRow ( int whichRow, int numberElements,
463 const int indicesColumn[], const double elements[] );
464 /// Takes out all entries for given rows
465 void emptyRows(int numberToEmpty, const int which[]);
466 //@}
467 /**@name used by ClpFactorization */
468 /// See if worth going sparse
469 void checkSparse();
470 /// For statistics
471 inline bool collectStatistics() const
472 { return collectStatistics_;}
473 /// For statistics
474 inline void setCollectStatistics(bool onOff) const
475 { collectStatistics_ = onOff;}
476 /// The real work of constructors etc 0 just scalars, 1 bit normal
477 void gutsOfDestructor(int type=1);
478 /// 1 bit - tolerances etc, 2 more, 4 dummy arrays
479 void gutsOfInitialize(int type);
480 void gutsOfCopy(const CoinFactorization &other);
481
482 /// Reset all sparsity etc statistics
483 void resetStatistics();
484
485
486 //@}
487
488 /**@name used by factorization */
489 /// Gets space for a factorization, called by constructors
490 void getAreas ( int numberRows,
491 int numberColumns,
492 CoinBigIndex maximumL,
493 CoinBigIndex maximumU );
494
495 /** PreProcesses raw triplet data.
496 state is 0 - triplets, 1 - some counts etc , 2 - .. */
497 void preProcess ( int state,
498 int possibleDuplicates = -1 );
499 /// Does most of factorization
500 int factor ( );
501protected:
502 /** Does sparse phase of factorization
503 return code is <0 error, 0= finished */
504 int factorSparse ( );
505 /** Does sparse phase of factorization (for smaller problems)
506 return code is <0 error, 0= finished */
507 int factorSparseSmall ( );
508 /** Does sparse phase of factorization (for larger problems)
509 return code is <0 error, 0= finished */
510 int factorSparseLarge ( );
511 /** Does dense phase of factorization
512 return code is <0 error, 0= finished */
513 int factorDense ( );
514
515 /// Pivots when just one other row so faster?
516 bool pivotOneOtherRow ( int pivotRow,
517 int pivotColumn );
518 /// Does one pivot on Row Singleton in factorization
519 bool pivotRowSingleton ( int pivotRow,
520 int pivotColumn );
521 /// Does one pivot on Column Singleton in factorization
522 bool pivotColumnSingleton ( int pivotRow,
523 int pivotColumn );
524
525 /** Gets space for one Column with given length,
526 may have to do compression (returns True if successful),
527 also moves existing vector,
528 extraNeeded is over and above present */
529 bool getColumnSpace ( int iColumn,
530 int extraNeeded );
531
532 /** Reorders U so contiguous and in order (if there is space)
533 Returns true if it could */
534 bool reorderU();
535 /** getColumnSpaceIterateR. Gets space for one extra R element in Column
536 may have to do compression (returns true)
537 also moves existing vector */
538 bool getColumnSpaceIterateR ( int iColumn, double value,
539 int iRow);
540 /** getColumnSpaceIterate. Gets space for one extra U element in Column
541 may have to do compression (returns true)
542 also moves existing vector.
543 Returns -1 if no memory or where element was put
544 Used by replaceRow (turns off R version) */
545 CoinBigIndex getColumnSpaceIterate ( int iColumn, double value,
546 int iRow);
547 /** Gets space for one Row with given length,
548 may have to do compression (returns True if successful),
549 also moves existing vector */
550 bool getRowSpace ( int iRow, int extraNeeded );
551
552 /** Gets space for one Row with given length while iterating,
553 may have to do compression (returns True if successful),
554 also moves existing vector */
555 bool getRowSpaceIterate ( int iRow,
556 int extraNeeded );
557 /// Checks that row and column copies look OK
558 void checkConsistency ( );
559 /// Adds a link in chain of equal counts
560 inline void addLink ( int index, int count ) {
561 int *nextCount = nextCount_.array();
562 int *firstCount = firstCount_.array();
563 int *lastCount = lastCount_.array();
564 int next = firstCount[count];
565 lastCount[index] = -2 - count;
566 if ( next < 0 ) {
567 //first with that count
568 firstCount[count] = index;
569 nextCount[index] = -1;
570 } else {
571 firstCount[count] = index;
572 nextCount[index] = next;
573 lastCount[next] = index;
574 }}
575 /// Deletes a link in chain of equal counts
576 inline void deleteLink ( int index ) {
577 int *nextCount = nextCount_.array();
578 int *firstCount = firstCount_.array();
579 int *lastCount = lastCount_.array();
580 int next = nextCount[index];
581 int last = lastCount[index];
582 if ( last >= 0 ) {
583 nextCount[last] = next;
584 } else {
585 int count = -last - 2;
586
587 firstCount[count] = next;
588 }
589 if ( next >= 0 ) {
590 lastCount[next] = last;
591 }
592 nextCount[index] = -2;
593 lastCount[index] = -2;
594 return;
595 }
596 /// Separate out links with same row/column count
597 void separateLinks(int count,bool rowsFirst);
598 /// Cleans up at end of factorization
599 void cleanup ( );
600
601 /// Updates part of column (FTRANL)
602 void updateColumnL ( CoinIndexedVector * region, int * indexIn ) const;
603 /// Updates part of column (FTRANL) when densish
604 void updateColumnLDensish ( CoinIndexedVector * region, int * indexIn ) const;
605 /// Updates part of column (FTRANL) when sparse
606 void updateColumnLSparse ( CoinIndexedVector * region, int * indexIn ) const;
607 /// Updates part of column (FTRANL) when sparsish
608 void updateColumnLSparsish ( CoinIndexedVector * region, int * indexIn ) const;
609
610 /// Updates part of column (FTRANR) without FT update
611 void updateColumnR ( CoinIndexedVector * region ) const;
612 /** Updates part of column (FTRANR) with FT update.
613 Also stores update after L and R */
614 void updateColumnRFT ( CoinIndexedVector * region, int * indexIn );
615
616 /// Updates part of column (FTRANU)
617 void updateColumnU ( CoinIndexedVector * region, int * indexIn) const;
618
619 /// Updates part of column (FTRANU) when sparse
620 void updateColumnUSparse ( CoinIndexedVector * regionSparse,
621 int * indexIn) const;
622 /// Updates part of column (FTRANU) when sparsish
623 void updateColumnUSparsish ( CoinIndexedVector * regionSparse,
624 int * indexIn) const;
625 /// Updates part of column (FTRANU)
626 int updateColumnUDensish ( double * COIN_RESTRICT region,
627 int * COIN_RESTRICT regionIndex) const;
628 /// Updates part of 2 columns (FTRANU) real work
629 void updateTwoColumnsUDensish (
630 int & numberNonZero1,
631 double * COIN_RESTRICT region1,
632 int * COIN_RESTRICT index1,
633 int & numberNonZero2,
634 double * COIN_RESTRICT region2,
635 int * COIN_RESTRICT index2) const;
636 /// Updates part of column PFI (FTRAN) (after rest)
637 void updateColumnPFI ( CoinIndexedVector * regionSparse) const;
638 /// Permutes back at end of updateColumn
639 void permuteBack ( CoinIndexedVector * regionSparse,
640 CoinIndexedVector * outVector) const;
641
642 /// Updates part of column transpose PFI (BTRAN) (before rest)
643 void updateColumnTransposePFI ( CoinIndexedVector * region) const;
644 /** Updates part of column transpose (BTRANU),
645 assumes index is sorted i.e. region is correct */
646 void updateColumnTransposeU ( CoinIndexedVector * region,
647 int smallestIndex) const;
648 /** Updates part of column transpose (BTRANU) when sparsish,
649 assumes index is sorted i.e. region is correct */
650 void updateColumnTransposeUSparsish ( CoinIndexedVector * region,
651 int smallestIndex) const;
652 /** Updates part of column transpose (BTRANU) when densish,
653 assumes index is sorted i.e. region is correct */
654 void updateColumnTransposeUDensish ( CoinIndexedVector * region,
655 int smallestIndex) const;
656 /** Updates part of column transpose (BTRANU) when sparse,
657 assumes index is sorted i.e. region is correct */
658 void updateColumnTransposeUSparse ( CoinIndexedVector * region) const;
659 /** Updates part of column transpose (BTRANU) by column
660 assumes index is sorted i.e. region is correct */
661 void updateColumnTransposeUByColumn ( CoinIndexedVector * region,
662 int smallestIndex) const;
663
664 /// Updates part of column transpose (BTRANR)
665 void updateColumnTransposeR ( CoinIndexedVector * region ) const;
666 /// Updates part of column transpose (BTRANR) when dense
667 void updateColumnTransposeRDensish ( CoinIndexedVector * region ) const;
668 /// Updates part of column transpose (BTRANR) when sparse
669 void updateColumnTransposeRSparse ( CoinIndexedVector * region ) const;
670
671 /// Updates part of column transpose (BTRANL)
672 void updateColumnTransposeL ( CoinIndexedVector * region ) const;
673 /// Updates part of column transpose (BTRANL) when densish by column
674 void updateColumnTransposeLDensish ( CoinIndexedVector * region ) const;
675 /// Updates part of column transpose (BTRANL) when densish by row
676 void updateColumnTransposeLByRow ( CoinIndexedVector * region ) const;
677 /// Updates part of column transpose (BTRANL) when sparsish by row
678 void updateColumnTransposeLSparsish ( CoinIndexedVector * region ) const;
679 /// Updates part of column transpose (BTRANL) when sparse (by Row)
680 void updateColumnTransposeLSparse ( CoinIndexedVector * region ) const;
681public:
682 /** Replaces one Column to basis for PFI
683 returns 0=OK, 1=Probably OK, 2=singular, 3=no room.
684 In this case region is not empty - it is incoming variable (updated)
685 */
686 int replaceColumnPFI ( CoinIndexedVector * regionSparse,
687 int pivotRow, double alpha);
688protected:
689 /** Returns accuracy status of replaceColumn
690 returns 0=OK, 1=Probably OK, 2=singular */
691 int checkPivot(double saveFromU, double oldPivot) const;
692 /********************************* START LARGE TEMPLATE ********/
693#ifdef INT_IS_8
694#define COINFACTORIZATION_BITS_PER_INT 64
695#define COINFACTORIZATION_SHIFT_PER_INT 6
696#define COINFACTORIZATION_MASK_PER_INT 0x3f
697#else
698#define COINFACTORIZATION_BITS_PER_INT 32
699#define COINFACTORIZATION_SHIFT_PER_INT 5
700#define COINFACTORIZATION_MASK_PER_INT 0x1f
701#endif
702 template <class T> inline bool
703 pivot ( int pivotRow,
704 int pivotColumn,
705 CoinBigIndex pivotRowPosition,
706 CoinBigIndex pivotColumnPosition,
707 CoinFactorizationDouble work[],
708 unsigned int workArea2[],
709 int increment2,
710 T markRow[] ,
711 int largeInteger)
712{
713 int *indexColumnU = indexColumnU_.array();
714 CoinBigIndex *startColumnU = startColumnU_.array();
715 int *numberInColumn = numberInColumn_.array();
716 CoinFactorizationDouble *elementU = elementU_.array();
717 int *indexRowU = indexRowU_.array();
718 CoinBigIndex *startRowU = startRowU_.array();
719 int *numberInRow = numberInRow_.array();
720 CoinFactorizationDouble *elementL = elementL_.array();
721 int *indexRowL = indexRowL_.array();
722 int *saveColumn = saveColumn_.array();
723 int *nextRow = nextRow_.array();
724 int *lastRow = lastRow_.array() ;
725
726 //store pivot columns (so can easily compress)
727 int numberInPivotRow = numberInRow[pivotRow] - 1;
728 CoinBigIndex startColumn = startColumnU[pivotColumn];
729 int numberInPivotColumn = numberInColumn[pivotColumn] - 1;
730 CoinBigIndex endColumn = startColumn + numberInPivotColumn + 1;
731 int put = 0;
732 CoinBigIndex startRow = startRowU[pivotRow];
733 CoinBigIndex endRow = startRow + numberInPivotRow + 1;
734
735 if ( pivotColumnPosition < 0 ) {
736 for ( pivotColumnPosition = startRow; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
737 int iColumn = indexColumnU[pivotColumnPosition];
738 if ( iColumn != pivotColumn ) {
739 saveColumn[put++] = iColumn;
740 } else {
741 break;
742 }
743 }
744 } else {
745 for (CoinBigIndex i = startRow ; i < pivotColumnPosition ; i++ ) {
746 saveColumn[put++] = indexColumnU[i];
747 }
748 }
749 assert (pivotColumnPosition<endRow);
750 assert (indexColumnU[pivotColumnPosition]==pivotColumn);
751 pivotColumnPosition++;
752 for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
753 saveColumn[put++] = indexColumnU[pivotColumnPosition];
754 }
755 //take out this bit of indexColumnU
756 int next = nextRow[pivotRow];
757 int last = lastRow[pivotRow];
758
759 nextRow[last] = next;
760 lastRow[next] = last;
761 nextRow[pivotRow] = numberGoodU_; //use for permute
762 lastRow[pivotRow] = -2;
763 numberInRow[pivotRow] = 0;
764 //store column in L, compress in U and take column out
765 CoinBigIndex l = lengthL_;
766
767 if ( l + numberInPivotColumn > lengthAreaL_ ) {
768 //need more memory
769 if ((messageLevel_&4)!=0)
770 printf("more memory needed in middle of invert\n");
771 return false;
772 }
773 //l+=currentAreaL_->elementByColumn-elementL;
774 CoinBigIndex lSave = l;
775
776 CoinBigIndex * startColumnL = startColumnL_.array();
777 startColumnL[numberGoodL_] = l; //for luck and first time
778 numberGoodL_++;
779 startColumnL[numberGoodL_] = l + numberInPivotColumn;
780 lengthL_ += numberInPivotColumn;
781 if ( pivotRowPosition < 0 ) {
782 for ( pivotRowPosition = startColumn; pivotRowPosition < endColumn; pivotRowPosition++ ) {
783 int iRow = indexRowU[pivotRowPosition];
784 if ( iRow != pivotRow ) {
785 indexRowL[l] = iRow;
786 elementL[l] = elementU[pivotRowPosition];
787 markRow[iRow] = static_cast<T>(l - lSave);
788 l++;
789 //take out of row list
790 CoinBigIndex start = startRowU[iRow];
791 CoinBigIndex end = start + numberInRow[iRow];
792 CoinBigIndex where = start;
793
794 while ( indexColumnU[where] != pivotColumn ) {
795 where++;
796 } /* endwhile */
797#if DEBUG_COIN
798 if ( where >= end ) {
799 abort ( );
800 }
801#endif
802 indexColumnU[where] = indexColumnU[end - 1];
803 numberInRow[iRow]--;
804 } else {
805 break;
806 }
807 }
808 } else {
809 CoinBigIndex i;
810
811 for ( i = startColumn; i < pivotRowPosition; i++ ) {
812 int iRow = indexRowU[i];
813
814 markRow[iRow] = static_cast<T>(l - lSave);
815 indexRowL[l] = iRow;
816 elementL[l] = elementU[i];
817 l++;
818 //take out of row list
819 CoinBigIndex start = startRowU[iRow];
820 CoinBigIndex end = start + numberInRow[iRow];
821 CoinBigIndex where = start;
822
823 while ( indexColumnU[where] != pivotColumn ) {
824 where++;
825 } /* endwhile */
826#if DEBUG_COIN
827 if ( where >= end ) {
828 abort ( );
829 }
830#endif
831 indexColumnU[where] = indexColumnU[end - 1];
832 numberInRow[iRow]--;
833 assert (numberInRow[iRow]>=0);
834 }
835 }
836 assert (pivotRowPosition<endColumn);
837 assert (indexRowU[pivotRowPosition]==pivotRow);
838 CoinFactorizationDouble pivotElement = elementU[pivotRowPosition];
839 CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;
840
841 pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
842 pivotRowPosition++;
843 for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
844 int iRow = indexRowU[pivotRowPosition];
845
846 markRow[iRow] = static_cast<T>(l - lSave);
847 indexRowL[l] = iRow;
848 elementL[l] = elementU[pivotRowPosition];
849 l++;
850 //take out of row list
851 CoinBigIndex start = startRowU[iRow];
852 CoinBigIndex end = start + numberInRow[iRow];
853 CoinBigIndex where = start;
854
855 while ( indexColumnU[where] != pivotColumn ) {
856 where++;
857 } /* endwhile */
858#if DEBUG_COIN
859 if ( where >= end ) {
860 abort ( );
861 }
862#endif
863 indexColumnU[where] = indexColumnU[end - 1];
864 numberInRow[iRow]--;
865 assert (numberInRow[iRow]>=0);
866 }
867 markRow[pivotRow] = static_cast<T>(largeInteger);
868 //compress pivot column (move pivot to front including saved)
869 numberInColumn[pivotColumn] = 0;
870 //use end of L for temporary space
871 int *indexL = &indexRowL[lSave];
872 CoinFactorizationDouble *multipliersL = &elementL[lSave];
873
874 //adjust
875 int j;
876
877 for ( j = 0; j < numberInPivotColumn; j++ ) {
878 multipliersL[j] *= pivotMultiplier;
879 }
880 //zero out fill
881 CoinBigIndex iErase;
882 for ( iErase = 0; iErase < increment2 * numberInPivotRow;
883 iErase++ ) {
884 workArea2[iErase] = 0;
885 }
886 CoinBigIndex added = numberInPivotRow * numberInPivotColumn;
887 unsigned int *temp2 = workArea2;
888 int * nextColumn = nextColumn_.array();
889
890 //pack down and move to work
891 int jColumn;
892 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
893 int iColumn = saveColumn[jColumn];
894 CoinBigIndex startColumn = startColumnU[iColumn];
895 CoinBigIndex endColumn = startColumn + numberInColumn[iColumn];
896 int iRow = indexRowU[startColumn];
897 CoinFactorizationDouble value = elementU[startColumn];
898 double largest;
899 CoinBigIndex put = startColumn;
900 CoinBigIndex positionLargest = -1;
901 CoinFactorizationDouble thisPivotValue = 0.0;
902
903 //compress column and find largest not updated
904 bool checkLargest;
905 int mark = markRow[iRow];
906
907 if ( mark == largeInteger+1 ) {
908 largest = fabs ( value );
909 positionLargest = put;
910 put++;
911 checkLargest = false;
912 } else {
913 //need to find largest
914 largest = 0.0;
915 checkLargest = true;
916 if ( mark != largeInteger ) {
917 //will be updated
918 work[mark] = value;
919 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
920 int bit = mark & COINFACTORIZATION_MASK_PER_INT;
921
922 temp2[word] = temp2[word] | ( 1 << bit ); //say already in counts
923 added--;
924 } else {
925 thisPivotValue = value;
926 }
927 }
928 CoinBigIndex i;
929 for ( i = startColumn + 1; i < endColumn; i++ ) {
930 iRow = indexRowU[i];
931 value = elementU[i];
932 int mark = markRow[iRow];
933
934 if ( mark == largeInteger+1 ) {
935 //keep
936 indexRowU[put] = iRow;
937 elementU[put] = value;
938 if ( checkLargest ) {
939 double absValue = fabs ( value );
940
941 if ( absValue > largest ) {
942 largest = absValue;
943 positionLargest = put;
944 }
945 }
946 put++;
947 } else if ( mark != largeInteger ) {
948 //will be updated
949 work[mark] = value;
950 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
951 int bit = mark & COINFACTORIZATION_MASK_PER_INT;
952
953 temp2[word] = temp2[word] | ( 1 << bit ); //say already in counts
954 added--;
955 } else {
956 thisPivotValue = value;
957 }
958 }
959 //slot in pivot
960 elementU[put] = elementU[startColumn];
961 indexRowU[put] = indexRowU[startColumn];
962 if ( positionLargest == startColumn ) {
963 positionLargest = put; //follow if was largest
964 }
965 put++;
966 elementU[startColumn] = thisPivotValue;
967 indexRowU[startColumn] = pivotRow;
968 //clean up counts
969 startColumn++;
970 numberInColumn[iColumn] = put - startColumn;
971 int * numberInColumnPlus = numberInColumnPlus_.array();
972 numberInColumnPlus[iColumn]++;
973 startColumnU[iColumn]++;
974 //how much space have we got
975 int next = nextColumn[iColumn];
976 CoinBigIndex space;
977
978 space = startColumnU[next] - put - numberInColumnPlus[next];
979 //assume no zero elements
980 if ( numberInPivotColumn > space ) {
981 //getColumnSpace also moves fixed part
982 if ( !getColumnSpace ( iColumn, numberInPivotColumn ) ) {
983 return false;
984 }
985 //redo starts
986 positionLargest = positionLargest + startColumnU[iColumn] - startColumn;
987 startColumn = startColumnU[iColumn];
988 put = startColumn + numberInColumn[iColumn];
989 }
990 double tolerance = zeroTolerance_;
991
992 int *nextCount = nextCount_.array();
993 for ( j = 0; j < numberInPivotColumn; j++ ) {
994 value = work[j] - thisPivotValue * multipliersL[j];
995 double absValue = fabs ( value );
996
997 if ( absValue > tolerance ) {
998 work[j] = 0.0;
999 assert (put<lengthAreaU_);
1000 elementU[put] = value;
1001 indexRowU[put] = indexL[j];
1002 if ( absValue > largest ) {
1003 largest = absValue;
1004 positionLargest = put;
1005 }
1006 put++;
1007 } else {
1008 work[j] = 0.0;
1009 added--;
1010 int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
1011 int bit = j & COINFACTORIZATION_MASK_PER_INT;
1012
1013 if ( temp2[word] & ( 1 << bit ) ) {
1014 //take out of row list
1015 iRow = indexL[j];
1016 CoinBigIndex start = startRowU[iRow];
1017 CoinBigIndex end = start + numberInRow[iRow];
1018 CoinBigIndex where = start;
1019
1020 while ( indexColumnU[where] != iColumn ) {
1021 where++;
1022 } /* endwhile */
1023#if DEBUG_COIN
1024 if ( where >= end ) {
1025 abort ( );
1026 }
1027#endif
1028 indexColumnU[where] = indexColumnU[end - 1];
1029 numberInRow[iRow]--;
1030 } else {
1031 //make sure won't be added
1032 int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
1033 int bit = j & COINFACTORIZATION_MASK_PER_INT;
1034
1035 temp2[word] = temp2[word] | ( 1 << bit ); //say already in counts
1036 }
1037 }
1038 }
1039 numberInColumn[iColumn] = put - startColumn;
1040 //move largest
1041 if ( positionLargest >= 0 ) {
1042 value = elementU[positionLargest];
1043 iRow = indexRowU[positionLargest];
1044 elementU[positionLargest] = elementU[startColumn];
1045 indexRowU[positionLargest] = indexRowU[startColumn];
1046 elementU[startColumn] = value;
1047 indexRowU[startColumn] = iRow;
1048 }
1049 //linked list for column
1050 if ( nextCount[iColumn + numberRows_] != -2 ) {
1051 //modify linked list
1052 deleteLink ( iColumn + numberRows_ );
1053 addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
1054 }
1055 temp2 += increment2;
1056 }
1057 //get space for row list
1058 unsigned int *putBase = workArea2;
1059 int bigLoops = numberInPivotColumn >> COINFACTORIZATION_SHIFT_PER_INT;
1060 int i = 0;
1061
1062 // do linked lists and update counts
1063 while ( bigLoops ) {
1064 bigLoops--;
1065 int bit;
1066 for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
1067 unsigned int *putThis = putBase;
1068 int iRow = indexL[i];
1069
1070 //get space
1071 int number = 0;
1072 int jColumn;
1073
1074 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
1075 unsigned int test = *putThis;
1076
1077 putThis += increment2;
1078 test = 1 - ( ( test >> bit ) & 1 );
1079 number += test;
1080 }
1081 int next = nextRow[iRow];
1082 CoinBigIndex space;
1083
1084 space = startRowU[next] - startRowU[iRow];
1085 number += numberInRow[iRow];
1086 if ( space < number ) {
1087 if ( !getRowSpace ( iRow, number ) ) {
1088 return false;
1089 }
1090 }
1091 // now do
1092 putThis = putBase;
1093 next = nextRow[iRow];
1094 number = numberInRow[iRow];
1095 CoinBigIndex end = startRowU[iRow] + number;
1096 int saveIndex = indexColumnU[startRowU[next]];
1097
1098 //add in
1099 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
1100 unsigned int test = *putThis;
1101
1102 putThis += increment2;
1103 test = 1 - ( ( test >> bit ) & 1 );
1104 indexColumnU[end] = saveColumn[jColumn];
1105 end += test;
1106 }
1107 //put back next one in case zapped
1108 indexColumnU[startRowU[next]] = saveIndex;
1109 markRow[iRow] = static_cast<T>(largeInteger+1);
1110 number = end - startRowU[iRow];
1111 numberInRow[iRow] = number;
1112 deleteLink ( iRow );
1113 addLink ( iRow, number );
1114 }
1115 putBase++;
1116 } /* endwhile */
1117 int bit;
1118
1119 for ( bit = 0; i < numberInPivotColumn; i++, bit++ ) {
1120 unsigned int *putThis = putBase;
1121 int iRow = indexL[i];
1122
1123 //get space
1124 int number = 0;
1125 int jColumn;
1126
1127 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
1128 unsigned int test = *putThis;
1129
1130 putThis += increment2;
1131 test = 1 - ( ( test >> bit ) & 1 );
1132 number += test;
1133 }
1134 int next = nextRow[iRow];
1135 CoinBigIndex space;
1136
1137 space = startRowU[next] - startRowU[iRow];
1138 number += numberInRow[iRow];
1139 if ( space < number ) {
1140 if ( !getRowSpace ( iRow, number ) ) {
1141 return false;
1142 }
1143 }
1144 // now do
1145 putThis = putBase;
1146 next = nextRow[iRow];
1147 number = numberInRow[iRow];
1148 CoinBigIndex end = startRowU[iRow] + number;
1149 int saveIndex;
1150
1151 saveIndex = indexColumnU[startRowU[next]];
1152
1153 //add in
1154 for ( jColumn = 0; jColumn < numberInPivotRow; jColumn++ ) {
1155 unsigned int test = *putThis;
1156
1157 putThis += increment2;
1158 test = 1 - ( ( test >> bit ) & 1 );
1159
1160 indexColumnU[end] = saveColumn[jColumn];
1161 end += test;
1162 }
1163 indexColumnU[startRowU[next]] = saveIndex;
1164 markRow[iRow] = static_cast<T>(largeInteger+1);
1165 number = end - startRowU[iRow];
1166 numberInRow[iRow] = number;
1167 deleteLink ( iRow );
1168 addLink ( iRow, number );
1169 }
1170 markRow[pivotRow] = static_cast<T>(largeInteger+1);
1171 //modify linked list for pivots
1172 deleteLink ( pivotRow );
1173 deleteLink ( pivotColumn + numberRows_ );
1174 totalElements_ += added;
1175 return true;
1176}
1177
1178 /********************************* END LARGE TEMPLATE ********/
1179 //@}
1180////////////////// data //////////////////
1181protected:
1182
1183 /**@name data */
1184 //@{
1185 /// Pivot tolerance
1186 double pivotTolerance_;
1187 /// Zero tolerance
1188 double zeroTolerance_;
1189#ifndef COIN_FAST_CODE
1190 /// Whether slack value is +1 or -1
1191 double slackValue_;
1192#else
1193#ifndef slackValue_
1194#define slackValue_ -1.0
1195#endif
1196#endif
1197 /// How much to multiply areas by
1198 double areaFactor_;
1199 /// Relax check on accuracy in replaceColumn
1200 double relaxCheck_;
1201 /// Number of Rows in factorization
1202 int numberRows_;
1203 /// Number of Rows after iterating
1204 int numberRowsExtra_;
1205 /// Maximum number of Rows after iterating
1206 int maximumRowsExtra_;
1207 /// Number of Columns in factorization
1208 int numberColumns_;
1209 /// Number of Columns after iterating
1210 int numberColumnsExtra_;
1211 /// Maximum number of Columns after iterating
1212 int maximumColumnsExtra_;
1213 /// Number factorized in U (not row singletons)
1214 int numberGoodU_;
1215 /// Number factorized in L
1216 int numberGoodL_;
1217 /// Maximum number of pivots before factorization
1218 int maximumPivots_;
1219 /// Number pivots since last factorization
1220 int numberPivots_;
1221 /// Number of elements in U (to go)
1222 /// or while iterating total overall
1223 CoinBigIndex totalElements_;
1224 /// Number of elements after factorization
1225 CoinBigIndex factorElements_;
1226 /// Pivot order for each Column
1227 CoinIntArrayWithLength pivotColumn_;
1228 /// Permutation vector for pivot row order
1229 CoinIntArrayWithLength permute_;
1230 /// DePermutation vector for pivot row order
1231 CoinIntArrayWithLength permuteBack_;
1232 /// Inverse Pivot order for each Column
1233 CoinIntArrayWithLength pivotColumnBack_;
1234 /// Status of factorization
1235 int status_;
1236
1237 /** 0 - no increasing rows - no permutations,
1238 1 - no increasing rows but permutations
1239 2 - increasing rows
1240 - taken out as always 2 */
1241 //int increasingRows_;
1242
1243 /// Number of trials before rejection
1244 int numberTrials_;
1245 /// Start of each Row as pointer
1246 CoinBigIndexArrayWithLength startRowU_;
1247
1248 /// Number in each Row
1249 CoinIntArrayWithLength numberInRow_;
1250
1251 /// Number in each Column
1252 CoinIntArrayWithLength numberInColumn_;
1253
1254 /// Number in each Column including pivoted
1255 CoinIntArrayWithLength numberInColumnPlus_;
1256
1257 /** First Row/Column with count of k,
1258 can tell which by offset - Rows then Columns */
1259 CoinIntArrayWithLength firstCount_;
1260
1261 /// Next Row/Column with count
1262 CoinIntArrayWithLength nextCount_;
1263
1264 /// Previous Row/Column with count
1265 CoinIntArrayWithLength lastCount_;
1266
1267 /// Next Column in memory order
1268 CoinIntArrayWithLength nextColumn_;
1269
1270 /// Previous Column in memory order
1271 CoinIntArrayWithLength lastColumn_;
1272
1273 /// Next Row in memory order
1274 CoinIntArrayWithLength nextRow_;
1275
1276 /// Previous Row in memory order
1277 CoinIntArrayWithLength lastRow_;
1278
1279 /// Columns left to do in a single pivot
1280 CoinIntArrayWithLength saveColumn_;
1281
1282 /// Marks rows to be updated
1283 CoinIntArrayWithLength markRow_;
1284
1285 /// Detail in messages
1286 int messageLevel_;
1287
1288 /// Larger of row and column size
1289 int biggerDimension_;
1290
1291 /// Base address for U (may change)
1292 CoinIntArrayWithLength indexColumnU_;
1293
1294 /// Pivots for L
1295 CoinIntArrayWithLength pivotRowL_;
1296
1297 /// Inverses of pivot values
1298 CoinFactorizationDoubleArrayWithLength pivotRegion_;
1299
1300 /// Number of slacks at beginning of U
1301 int numberSlacks_;
1302
1303 /// Number in U
1304 int numberU_;
1305
1306 /// Maximum space used in U
1307 CoinBigIndex maximumU_;
1308
1309 /// Base of U is always 0
1310 //int baseU_;
1311
1312 /// Length of U
1313 CoinBigIndex lengthU_;
1314
1315 /// Length of area reserved for U
1316 CoinBigIndex lengthAreaU_;
1317
1318/// Elements of U
1319 CoinFactorizationDoubleArrayWithLength elementU_;
1320
1321/// Row indices of U
1322 CoinIntArrayWithLength indexRowU_;
1323
1324/// Start of each column in U
1325 CoinBigIndexArrayWithLength startColumnU_;
1326
1327/// Converts rows to columns in U
1328 CoinBigIndexArrayWithLength convertRowToColumnU_;
1329
1330 /// Number in L
1331 CoinBigIndex numberL_;
1332
1333/// Base of L
1334 CoinBigIndex baseL_;
1335
1336 /// Length of L
1337 CoinBigIndex lengthL_;
1338
1339 /// Length of area reserved for L
1340 CoinBigIndex lengthAreaL_;
1341
1342 /// Elements of L
1343 CoinFactorizationDoubleArrayWithLength elementL_;
1344
1345 /// Row indices of L
1346 CoinIntArrayWithLength indexRowL_;
1347
1348 /// Start of each column in L
1349 CoinBigIndexArrayWithLength startColumnL_;
1350
1351 /// true if Forrest Tomlin update, false if PFI
1352 bool doForrestTomlin_;
1353
1354 /// Number in R
1355 int numberR_;
1356
1357 /// Length of R stuff
1358 CoinBigIndex lengthR_;
1359
1360 /// length of area reserved for R
1361 CoinBigIndex lengthAreaR_;
1362
1363 /// Elements of R
1364 CoinFactorizationDouble *elementR_;
1365
1366 /// Row indices for R
1367 int *indexRowR_;
1368
1369 /// Start of columns for R
1370 CoinBigIndexArrayWithLength startColumnR_;
1371
1372 /// Dense area
1373 double * denseArea_;
1374
1375 /// Dense permutation
1376 int * densePermute_;
1377
1378 /// Number of dense rows
1379 int numberDense_;
1380
1381 /// Dense threshold
1382 int denseThreshold_;
1383
1384 /// First work area
1385 CoinFactorizationDoubleArrayWithLength workArea_;
1386
1387 /// Second work area
1388 CoinUnsignedIntArrayWithLength workArea2_;
1389
1390 /// Number of compressions done
1391 CoinBigIndex numberCompressions_;
1392
1393 /// Below are all to collect
1394 mutable double ftranCountInput_;
1395 mutable double ftranCountAfterL_;
1396 mutable double ftranCountAfterR_;
1397 mutable double ftranCountAfterU_;
1398 mutable double btranCountInput_;
1399 mutable double btranCountAfterU_;
1400 mutable double btranCountAfterR_;
1401 mutable double btranCountAfterL_;
1402
1403 /// We can roll over factorizations
1404 mutable int numberFtranCounts_;
1405 mutable int numberBtranCounts_;
1406
1407 /// While these are average ratios collected over last period
1408 double ftranAverageAfterL_;
1409 double ftranAverageAfterR_;
1410 double ftranAverageAfterU_;
1411 double btranAverageAfterU_;
1412 double btranAverageAfterR_;
1413 double btranAverageAfterL_;
1414
1415 /// For statistics
1416 mutable bool collectStatistics_;
1417
1418 /// Below this use sparse technology - if 0 then no L row copy
1419 int sparseThreshold_;
1420
1421 /// And one for "sparsish"
1422 int sparseThreshold2_;
1423
1424 /// Start of each row in L
1425 CoinBigIndexArrayWithLength startRowL_;
1426
1427 /// Index of column in row for L
1428 CoinIntArrayWithLength indexColumnL_;
1429
1430 /// Elements in L (row copy)
1431 CoinFactorizationDoubleArrayWithLength elementByRowL_;
1432
1433 /// Sparse regions
1434 mutable CoinIntArrayWithLength sparse_;
1435 /** L to U bias
1436 0 - U bias, 1 - some U bias, 2 some L bias, 3 L bias
1437 */
1438 int biasLU_;
1439 /** Array persistence flag
1440 If 0 then as now (delete/new)
1441 1 then only do arrays if bigger needed
1442 2 as 1 but give a bit extra if bigger needed
1443 */
1444 int persistenceFlag_;
1445 //@}
1446};
1447// Dense coding
1448#ifdef COIN_HAS_LAPACK
1449#define DENSE_CODE 1
1450/* Type of Fortran integer translated into C */
1451#ifndef ipfint
1452//typedef ipfint FORTRAN_INTEGER_TYPE ;
1453typedef int ipfint;
1454typedef const int cipfint;
1455#endif
1456#endif
1457#endif
1458// Extra for ugly include
1459#ifdef UGLY_COIN_FACTOR_CODING
1460#define FAC_UNSET (FAC_SET+1)
1461{
1462 goodPivot=false;
1463 //store pivot columns (so can easily compress)
1464 CoinBigIndex startColumnThis = startColumn[iPivotColumn];
1465 CoinBigIndex endColumn = startColumnThis + numberDoColumn + 1;
1466 int put = 0;
1467 CoinBigIndex startRowThis = startRow[iPivotRow];
1468 CoinBigIndex endRow = startRowThis + numberDoRow + 1;
1469 if ( pivotColumnPosition < 0 ) {
1470 for ( pivotColumnPosition = startRowThis; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
1471 int iColumn = indexColumn[pivotColumnPosition];
1472 if ( iColumn != iPivotColumn ) {
1473 saveColumn[put++] = iColumn;
1474 } else {
1475 break;
1476 }
1477 }
1478 } else {
1479 for (CoinBigIndex i = startRowThis ; i < pivotColumnPosition ; i++ ) {
1480 saveColumn[put++] = indexColumn[i];
1481 }
1482 }
1483 assert (pivotColumnPosition<endRow);
1484 assert (indexColumn[pivotColumnPosition]==iPivotColumn);
1485 pivotColumnPosition++;
1486 for ( ; pivotColumnPosition < endRow; pivotColumnPosition++ ) {
1487 saveColumn[put++] = indexColumn[pivotColumnPosition];
1488 }
1489 //take out this bit of indexColumn
1490 int next = nextRow[iPivotRow];
1491 int last = lastRow[iPivotRow];
1492
1493 nextRow[last] = next;
1494 lastRow[next] = last;
1495 nextRow[iPivotRow] = numberGoodU_; //use for permute
1496 lastRow[iPivotRow] = -2;
1497 numberInRow[iPivotRow] = 0;
1498 //store column in L, compress in U and take column out
1499 CoinBigIndex l = lengthL_;
1500 // **** HORRID coding coming up but a goto seems best!
1501 {
1502 if ( l + numberDoColumn > lengthAreaL_ ) {
1503 //need more memory
1504 if ((messageLevel_&4)!=0)
1505 printf("more memory needed in middle of invert\n");
1506 goto BAD_PIVOT;
1507 }
1508 //l+=currentAreaL_->elementByColumn-elementL;
1509 CoinBigIndex lSave = l;
1510
1511 CoinBigIndex * startColumnL = startColumnL_.array();
1512 startColumnL[numberGoodL_] = l; //for luck and first time
1513 numberGoodL_++;
1514 startColumnL[numberGoodL_] = l + numberDoColumn;
1515 lengthL_ += numberDoColumn;
1516 if ( pivotRowPosition < 0 ) {
1517 for ( pivotRowPosition = startColumnThis; pivotRowPosition < endColumn; pivotRowPosition++ ) {
1518 int iRow = indexRow[pivotRowPosition];
1519 if ( iRow != iPivotRow ) {
1520 indexRowL[l] = iRow;
1521 elementL[l] = element[pivotRowPosition];
1522 markRow[iRow] = l - lSave;
1523 l++;
1524 //take out of row list
1525 CoinBigIndex start = startRow[iRow];
1526 CoinBigIndex end = start + numberInRow[iRow];
1527 CoinBigIndex where = start;
1528
1529 while ( indexColumn[where] != iPivotColumn ) {
1530 where++;
1531 } /* endwhile */
1532#if DEBUG_COIN
1533 if ( where >= end ) {
1534 abort ( );
1535 }
1536#endif
1537 indexColumn[where] = indexColumn[end - 1];
1538 numberInRow[iRow]--;
1539 } else {
1540 break;
1541 }
1542 }
1543 } else {
1544 CoinBigIndex i;
1545
1546 for ( i = startColumnThis; i < pivotRowPosition; i++ ) {
1547 int iRow = indexRow[i];
1548
1549 markRow[iRow] = l - lSave;
1550 indexRowL[l] = iRow;
1551 elementL[l] = element[i];
1552 l++;
1553 //take out of row list
1554 CoinBigIndex start = startRow[iRow];
1555 CoinBigIndex end = start + numberInRow[iRow];
1556 CoinBigIndex where = start;
1557
1558 while ( indexColumn[where] != iPivotColumn ) {
1559 where++;
1560 } /* endwhile */
1561#if DEBUG_COIN
1562 if ( where >= end ) {
1563 abort ( );
1564 }
1565#endif
1566 indexColumn[where] = indexColumn[end - 1];
1567 numberInRow[iRow]--;
1568 assert (numberInRow[iRow]>=0);
1569 }
1570 }
1571 assert (pivotRowPosition<endColumn);
1572 assert (indexRow[pivotRowPosition]==iPivotRow);
1573 CoinFactorizationDouble pivotElement = element[pivotRowPosition];
1574 CoinFactorizationDouble pivotMultiplier = 1.0 / pivotElement;
1575
1576 pivotRegion_.array()[numberGoodU_] = pivotMultiplier;
1577 pivotRowPosition++;
1578 for ( ; pivotRowPosition < endColumn; pivotRowPosition++ ) {
1579 int iRow = indexRow[pivotRowPosition];
1580
1581 markRow[iRow] = l - lSave;
1582 indexRowL[l] = iRow;
1583 elementL[l] = element[pivotRowPosition];
1584 l++;
1585 //take out of row list
1586 CoinBigIndex start = startRow[iRow];
1587 CoinBigIndex end = start + numberInRow[iRow];
1588 CoinBigIndex where = start;
1589
1590 while ( indexColumn[where] != iPivotColumn ) {
1591 where++;
1592 } /* endwhile */
1593#if DEBUG_COIN
1594 if ( where >= end ) {
1595 abort ( );
1596 }
1597#endif
1598 indexColumn[where] = indexColumn[end - 1];
1599 numberInRow[iRow]--;
1600 assert (numberInRow[iRow]>=0);
1601 }
1602 markRow[iPivotRow] = FAC_SET;
1603 //compress pivot column (move pivot to front including saved)
1604 numberInColumn[iPivotColumn] = 0;
1605 //use end of L for temporary space
1606 int *indexL = &indexRowL[lSave];
1607 CoinFactorizationDouble *multipliersL = &elementL[lSave];
1608
1609 //adjust
1610 int j;
1611
1612 for ( j = 0; j < numberDoColumn; j++ ) {
1613 multipliersL[j] *= pivotMultiplier;
1614 }
1615 //zero out fill
1616 CoinBigIndex iErase;
1617 for ( iErase = 0; iErase < increment2 * numberDoRow;
1618 iErase++ ) {
1619 workArea2[iErase] = 0;
1620 }
1621 CoinBigIndex added = numberDoRow * numberDoColumn;
1622 unsigned int *temp2 = workArea2;
1623 int * nextColumn = nextColumn_.array();
1624
1625 //pack down and move to work
1626 int jColumn;
1627 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
1628 int iColumn = saveColumn[jColumn];
1629 CoinBigIndex startColumnThis = startColumn[iColumn];
1630 CoinBigIndex endColumn = startColumnThis + numberInColumn[iColumn];
1631 int iRow = indexRow[startColumnThis];
1632 CoinFactorizationDouble value = element[startColumnThis];
1633 double largest;
1634 CoinBigIndex put = startColumnThis;
1635 CoinBigIndex positionLargest = -1;
1636 CoinFactorizationDouble thisPivotValue = 0.0;
1637
1638 //compress column and find largest not updated
1639 bool checkLargest;
1640 int mark = markRow[iRow];
1641
1642 if ( mark == FAC_UNSET ) {
1643 largest = fabs ( value );
1644 positionLargest = put;
1645 put++;
1646 checkLargest = false;
1647 } else {
1648 //need to find largest
1649 largest = 0.0;
1650 checkLargest = true;
1651 if ( mark != FAC_SET ) {
1652 //will be updated
1653 workArea[mark] = value;
1654 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
1655 int bit = mark & COINFACTORIZATION_MASK_PER_INT;
1656
1657 temp2[word] = temp2[word] | ( 1 << bit ); //say already in counts
1658 added--;
1659 } else {
1660 thisPivotValue = value;
1661 }
1662 }
1663 CoinBigIndex i;
1664 for ( i = startColumnThis + 1; i < endColumn; i++ ) {
1665 iRow = indexRow[i];
1666 value = element[i];
1667 int mark = markRow[iRow];
1668
1669 if ( mark == FAC_UNSET ) {
1670 //keep
1671 indexRow[put] = iRow;
1672 element[put] = value;
1673 if ( checkLargest ) {
1674 double absValue = fabs ( value );
1675
1676 if ( absValue > largest ) {
1677 largest = absValue;
1678 positionLargest = put;
1679 }
1680 }
1681 put++;
1682 } else if ( mark != FAC_SET ) {
1683 //will be updated
1684 workArea[mark] = value;
1685 int word = mark >> COINFACTORIZATION_SHIFT_PER_INT;
1686 int bit = mark & COINFACTORIZATION_MASK_PER_INT;
1687
1688 temp2[word] = temp2[word] | ( 1 << bit ); //say already in counts
1689 added--;
1690 } else {
1691 thisPivotValue = value;
1692 }
1693 }
1694 //slot in pivot
1695 element[put] = element[startColumnThis];
1696 indexRow[put] = indexRow[startColumnThis];
1697 if ( positionLargest == startColumnThis ) {
1698 positionLargest = put; //follow if was largest
1699 }
1700 put++;
1701 element[startColumnThis] = thisPivotValue;
1702 indexRow[startColumnThis] = iPivotRow;
1703 //clean up counts
1704 startColumnThis++;
1705 numberInColumn[iColumn] = put - startColumnThis;
1706 int * numberInColumnPlus = numberInColumnPlus_.array();
1707 numberInColumnPlus[iColumn]++;
1708 startColumn[iColumn]++;
1709 //how much space have we got
1710 int next = nextColumn[iColumn];
1711 CoinBigIndex space;
1712
1713 space = startColumn[next] - put - numberInColumnPlus[next];
1714 //assume no zero elements
1715 if ( numberDoColumn > space ) {
1716 //getColumnSpace also moves fixed part
1717 if ( !getColumnSpace ( iColumn, numberDoColumn ) ) {
1718 goto BAD_PIVOT;
1719 }
1720 //redo starts
1721 positionLargest = positionLargest + startColumn[iColumn] - startColumnThis;
1722 startColumnThis = startColumn[iColumn];
1723 put = startColumnThis + numberInColumn[iColumn];
1724 }
1725 double tolerance = zeroTolerance_;
1726
1727 int *nextCount = nextCount_.array();
1728 for ( j = 0; j < numberDoColumn; j++ ) {
1729 value = workArea[j] - thisPivotValue * multipliersL[j];
1730 double absValue = fabs ( value );
1731
1732 if ( absValue > tolerance ) {
1733 workArea[j] = 0.0;
1734 element[put] = value;
1735 indexRow[put] = indexL[j];
1736 if ( absValue > largest ) {
1737 largest = absValue;
1738 positionLargest = put;
1739 }
1740 put++;
1741 } else {
1742 workArea[j] = 0.0;
1743 added--;
1744 int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
1745 int bit = j & COINFACTORIZATION_MASK_PER_INT;
1746
1747 if ( temp2[word] & ( 1 << bit ) ) {
1748 //take out of row list
1749 iRow = indexL[j];
1750 CoinBigIndex start = startRow[iRow];
1751 CoinBigIndex end = start + numberInRow[iRow];
1752 CoinBigIndex where = start;
1753
1754 while ( indexColumn[where] != iColumn ) {
1755 where++;
1756 } /* endwhile */
1757#if DEBUG_COIN
1758 if ( where >= end ) {
1759 abort ( );
1760 }
1761#endif
1762 indexColumn[where] = indexColumn[end - 1];
1763 numberInRow[iRow]--;
1764 } else {
1765 //make sure won't be added
1766 int word = j >> COINFACTORIZATION_SHIFT_PER_INT;
1767 int bit = j & COINFACTORIZATION_MASK_PER_INT;
1768
1769 temp2[word] = temp2[word] | ( 1 << bit ); //say already in counts
1770 }
1771 }
1772 }
1773 numberInColumn[iColumn] = put - startColumnThis;
1774 //move largest
1775 if ( positionLargest >= 0 ) {
1776 value = element[positionLargest];
1777 iRow = indexRow[positionLargest];
1778 element[positionLargest] = element[startColumnThis];
1779 indexRow[positionLargest] = indexRow[startColumnThis];
1780 element[startColumnThis] = value;
1781 indexRow[startColumnThis] = iRow;
1782 }
1783 //linked list for column
1784 if ( nextCount[iColumn + numberRows_] != -2 ) {
1785 //modify linked list
1786 deleteLink ( iColumn + numberRows_ );
1787 addLink ( iColumn + numberRows_, numberInColumn[iColumn] );
1788 }
1789 temp2 += increment2;
1790 }
1791 //get space for row list
1792 unsigned int *putBase = workArea2;
1793 int bigLoops = numberDoColumn >> COINFACTORIZATION_SHIFT_PER_INT;
1794 int i = 0;
1795
1796 // do linked lists and update counts
1797 while ( bigLoops ) {
1798 bigLoops--;
1799 int bit;
1800 for ( bit = 0; bit < COINFACTORIZATION_BITS_PER_INT; i++, bit++ ) {
1801 unsigned int *putThis = putBase;
1802 int iRow = indexL[i];
1803
1804 //get space
1805 int number = 0;
1806 int jColumn;
1807
1808 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
1809 unsigned int test = *putThis;
1810
1811 putThis += increment2;
1812 test = 1 - ( ( test >> bit ) & 1 );
1813 number += test;
1814 }
1815 int next = nextRow[iRow];
1816 CoinBigIndex space;
1817
1818 space = startRow[next] - startRow[iRow];
1819 number += numberInRow[iRow];
1820 if ( space < number ) {
1821 if ( !getRowSpace ( iRow, number ) ) {
1822 goto BAD_PIVOT;
1823 }
1824 }
1825 // now do
1826 putThis = putBase;
1827 next = nextRow[iRow];
1828 number = numberInRow[iRow];
1829 CoinBigIndex end = startRow[iRow] + number;
1830 int saveIndex = indexColumn[startRow[next]];
1831
1832 //add in
1833 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
1834 unsigned int test = *putThis;
1835
1836 putThis += increment2;
1837 test = 1 - ( ( test >> bit ) & 1 );
1838 indexColumn[end] = saveColumn[jColumn];
1839 end += test;
1840 }
1841 //put back next one in case zapped
1842 indexColumn[startRow[next]] = saveIndex;
1843 markRow[iRow] = FAC_UNSET;
1844 number = end - startRow[iRow];
1845 numberInRow[iRow] = number;
1846 deleteLink ( iRow );
1847 addLink ( iRow, number );
1848 }
1849 putBase++;
1850 } /* endwhile */
1851 int bit;
1852
1853 for ( bit = 0; i < numberDoColumn; i++, bit++ ) {
1854 unsigned int *putThis = putBase;
1855 int iRow = indexL[i];
1856
1857 //get space
1858 int number = 0;
1859 int jColumn;
1860
1861 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
1862 unsigned int test = *putThis;
1863
1864 putThis += increment2;
1865 test = 1 - ( ( test >> bit ) & 1 );
1866 number += test;
1867 }
1868 int next = nextRow[iRow];
1869 CoinBigIndex space;
1870
1871 space = startRow[next] - startRow[iRow];
1872 number += numberInRow[iRow];
1873 if ( space < number ) {
1874 if ( !getRowSpace ( iRow, number ) ) {
1875 goto BAD_PIVOT;
1876 }
1877 }
1878 // now do
1879 putThis = putBase;
1880 next = nextRow[iRow];
1881 number = numberInRow[iRow];
1882 CoinBigIndex end = startRow[iRow] + number;
1883 int saveIndex;
1884
1885 saveIndex = indexColumn[startRow[next]];
1886
1887 //add in
1888 for ( jColumn = 0; jColumn < numberDoRow; jColumn++ ) {
1889 unsigned int test = *putThis;
1890
1891 putThis += increment2;
1892 test = 1 - ( ( test >> bit ) & 1 );
1893
1894 indexColumn[end] = saveColumn[jColumn];
1895 end += test;
1896 }
1897 indexColumn[startRow[next]] = saveIndex;
1898 markRow[iRow] = FAC_UNSET;
1899 number = end - startRow[iRow];
1900 numberInRow[iRow] = number;
1901 deleteLink ( iRow );
1902 addLink ( iRow, number );
1903 }
1904 markRow[iPivotRow] = FAC_UNSET;
1905 //modify linked list for pivots
1906 deleteLink ( iPivotRow );
1907 deleteLink ( iPivotColumn + numberRows_ );
1908 totalElements_ += added;
1909 goodPivot= true;
1910 // **** UGLY UGLY UGLY
1911 }
1912 BAD_PIVOT:
1913
1914 ;
1915}
1916#undef FAC_UNSET
1917#endif
1918