source: trunk/Clp/src/AbcMatrix.cpp @ 2271

Last change on this file since 2271 was 2074, checked in by forrest, 5 years ago

changes to try and make more stable

  • Property svn:keywords set to Id
File size: 119.2 KB
Line 
1/* $Id: AbcMatrix.cpp 2074 2014-12-10 09:43:54Z forrest $ */
2// Copyright (C) 2002, International Business Machines
3// Corporation and others, Copyright (C) 2012, FasterCoin.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6
7
8#include <cstdio>
9#include "CoinPragma.hpp"
10#include "CoinIndexedVector.hpp"
11#include "CoinHelperFunctions.hpp"
12#include "CoinAbcHelperFunctions.hpp"
13#include "AbcSimplexFactorization.hpp"
14#include "AbcPrimalColumnDantzig.hpp"
15#include "AbcPrimalColumnSteepest.hpp"
16#include "CoinTime.hpp"
17
18#include "AbcSimplex.hpp"
19#include "AbcSimplexDual.hpp"
20// at end to get min/max!
21#include "AbcMatrix.hpp"
22#include "ClpMessage.hpp"
23#ifdef INTEL_MKL
24#include "mkl_spblas.h"
25#endif
26#if ABC_INSTRUMENT>1
27extern int abcPricing[20];
28extern int abcPricingDense[20];
29#endif
30//=============================================================================
31
32//#############################################################################
33// Constructors / Destructor / Assignment
34//#############################################################################
35
36//-------------------------------------------------------------------
37// Default Constructor
38//-------------------------------------------------------------------
39AbcMatrix::AbcMatrix ()
40  : matrix_(NULL),
41    model_(NULL),
42    rowStart_(NULL),
43    element_(NULL),
44    column_(NULL),
45    numberColumnBlocks_(0),
46    numberRowBlocks_(0),
47#ifdef COUNT_COPY
48    countRealColumn_(NULL),
49    countStartLarge_(NULL),
50    countRow_(NULL),
51    countElement_(NULL),
52    smallestCount_(0),
53    largestCount_(0),
54#endif
55    startFraction_(0.0),
56    endFraction_(1.0),
57    savedBestDj_(0.0),
58    originalWanted_(0),
59    currentWanted_(0),
60    savedBestSequence_(-1),
61    minimumObjectsScan_(-1),
62    minimumGoodReducedCosts_(-1)
63{
64}
65
66//-------------------------------------------------------------------
67// Copy constructor
68//-------------------------------------------------------------------
69AbcMatrix::AbcMatrix (const AbcMatrix & rhs)
70{
71#ifndef COIN_SPARSE_MATRIX
72  matrix_ = new CoinPackedMatrix(*(rhs.matrix_), -1, -1);
73#else
74  matrix_ = new CoinPackedMatrix(*(rhs.matrix_), -0, -0);
75#endif
76  model_=rhs.model_;
77  rowStart_ = NULL;
78  element_ = NULL;
79  column_ = NULL;
80#ifdef COUNT_COPY
81  countRealColumn_ = NULL;
82  countStartLarge_ = NULL;
83  countRow_ = NULL;
84  countElement_ = NULL;
85#endif
86  numberColumnBlocks_ = rhs.numberColumnBlocks_;
87  CoinAbcMemcpy(startColumnBlock_,rhs.startColumnBlock_,numberColumnBlocks_+1);
88  numberRowBlocks_ = rhs.numberRowBlocks_;
89  if (numberRowBlocks_) {
90    assert (model_);
91    int numberRows=model_->numberRows();
92    int numberElements = matrix_->getNumElements();
93    memcpy(blockStart_,rhs.blockStart_,sizeof(blockStart_));
94    rowStart_=CoinCopyOfArray(rhs.rowStart_,numberRows*(numberRowBlocks_+2));
95    element_=CoinCopyOfArray(rhs.element_,numberElements);
96    column_=CoinCopyOfArray(rhs.column_,numberElements);
97#ifdef COUNT_COPY
98    smallestCount_ = rhs.smallestCount_;
99    largestCount_ = rhs.largestCount_;
100    int numberColumns=model_->numberColumns();
101    countRealColumn_=CoinCopyOfArray(rhs.countRealColumn_,numberColumns);
102    memcpy(countStart_,rhs.countStart_,reinterpret_cast<char *>(&countRealColumn_)-
103           reinterpret_cast<char *>(countStart_));
104    int numberLarge = numberColumns-countStart_[MAX_COUNT];
105    countStartLarge_=CoinCopyOfArray(rhs.countStartLarge_,numberLarge+1);
106    numberElements=countStartLarge_[numberLarge];
107    countElement_=CoinCopyOfArray(rhs.countElement_,numberElements);
108    countRow_=CoinCopyOfArray(rhs.countRow_,numberElements);
109#endif
110  }
111}
112
113AbcMatrix::AbcMatrix (const CoinPackedMatrix & rhs)
114{
115#ifndef COIN_SPARSE_MATRIX
116  matrix_ = new CoinPackedMatrix(rhs, -1, -1);
117#else
118  matrix_ = new CoinPackedMatrix(rhs, -0, -0);
119#endif
120  matrix_->cleanMatrix();
121  model_=NULL;
122  rowStart_ = NULL;
123  element_ = NULL;
124  column_ = NULL;
125#ifdef COUNT_COPY
126  countRealColumn_ = NULL;
127  countStartLarge_ = NULL;
128  countRow_ = NULL;
129  countElement_ = NULL;
130  smallestCount_ = 0;
131  largestCount_ = 0;
132#endif
133  numberColumnBlocks_=1;
134  startColumnBlock_[0]=0;
135  startColumnBlock_[1]=0;
136  numberRowBlocks_ = 0;
137  startFraction_ = 0;
138  endFraction_ = 1.0;
139  savedBestDj_ = 0;
140  originalWanted_ = 0;
141  currentWanted_ = 0;
142  savedBestSequence_ = -1;
143  minimumObjectsScan_ = -1;
144  minimumGoodReducedCosts_ = -1;
145}
146#ifdef ABC_SPRINT
147/* Subset constructor (without gaps). */
148AbcMatrix::AbcMatrix (const AbcMatrix & wholeMatrix,
149                      int numberRows, const int * whichRows,
150                      int numberColumns, const int * whichColumns)
151{
152#ifndef COIN_SPARSE_MATRIX
153  matrix_ = new CoinPackedMatrix(*wholeMatrix.matrix_,numberRows,whichRows,
154                                 numberColumns,whichColumns);
155#else
156  matrix_ = new CoinPackedMatrix(rhs, -0, -0);
157  abort();
158#endif
159  matrix_->cleanMatrix();
160  model_=NULL;
161  rowStart_ = NULL;
162  element_ = NULL;
163  column_ = NULL;
164#ifdef COUNT_COPY
165  countRealColumn_ = NULL;
166  countStartLarge_ = NULL;
167  countRow_ = NULL;
168  countElement_ = NULL;
169  smallestCount_ = 0;
170  largestCount_ = 0;
171#endif
172  numberColumnBlocks_=1;
173  startColumnBlock_[0]=0;
174  startColumnBlock_[1]=0;
175  numberRowBlocks_ = 0;
176  startFraction_ = 0;
177  endFraction_ = 1.0;
178  savedBestDj_ = 0;
179  originalWanted_ = 0;
180  currentWanted_ = 0;
181  savedBestSequence_ = -1;
182  minimumObjectsScan_ = -1;
183  minimumGoodReducedCosts_ = -1;
184}
185#endif
186//-------------------------------------------------------------------
187// Destructor
188//-------------------------------------------------------------------
189AbcMatrix::~AbcMatrix ()
190{
191  delete matrix_;
192  delete [] rowStart_;
193  delete [] element_;
194  delete [] column_;
195#ifdef COUNT_COPY
196  delete [] countRealColumn_;
197  delete [] countStartLarge_;
198  delete [] countRow_;
199  delete [] countElement_;
200#endif
201}
202
203//----------------------------------------------------------------
204// Assignment operator
205//-------------------------------------------------------------------
206AbcMatrix &
207AbcMatrix::operator=(const AbcMatrix& rhs)
208{
209  if (this != &rhs) {
210    delete matrix_;
211#ifndef COIN_SPARSE_MATRIX
212    matrix_ = new CoinPackedMatrix(*(rhs.matrix_));
213#else
214    matrix_ = new CoinPackedMatrix(*(rhs.matrix_), -0, -0);
215#endif
216    model_=rhs.model_;
217    delete [] rowStart_;
218    delete [] element_;
219    delete [] column_;
220#ifdef COUNT_COPY
221    delete [] countRealColumn_;
222    delete [] countStartLarge_;
223    delete [] countRow_;
224    delete [] countElement_;
225#endif
226    rowStart_ = NULL;
227    element_ = NULL;
228    column_ = NULL;
229#ifdef COUNT_COPY
230    countRealColumn_ = NULL;
231    countStartLarge_ = NULL;
232    countRow_ = NULL;
233    countElement_ = NULL;
234#endif
235    numberColumnBlocks_ = rhs.numberColumnBlocks_;
236    CoinAbcMemcpy(startColumnBlock_,rhs.startColumnBlock_,numberColumnBlocks_+1);
237    numberRowBlocks_ = rhs.numberRowBlocks_;
238    if (numberRowBlocks_) {
239      assert (model_);
240      int numberRows=model_->numberRows();
241      int numberElements = matrix_->getNumElements();
242      memcpy(blockStart_,rhs.blockStart_,sizeof(blockStart_));
243      rowStart_=CoinCopyOfArray(rhs.rowStart_,numberRows*(numberRowBlocks_+2));
244      element_=CoinCopyOfArray(rhs.element_,numberElements);
245      column_=CoinCopyOfArray(rhs.column_,numberElements);
246#ifdef COUNT_COPY
247      smallestCount_ = rhs.smallestCount_;
248      largestCount_ = rhs.largestCount_;
249      int numberColumns=model_->numberColumns();
250      countRealColumn_=CoinCopyOfArray(rhs.countRealColumn_,numberColumns);
251      memcpy(countStart_,rhs.countStart_,reinterpret_cast<char *>(&countRealColumn_)-
252             reinterpret_cast<char *>(countStart_));
253      int numberLarge = numberColumns-countStart_[MAX_COUNT];
254      countStartLarge_=CoinCopyOfArray(rhs.countStartLarge_,numberLarge+1);
255      numberElements=countStartLarge_[numberLarge];
256      countElement_=CoinCopyOfArray(rhs.countElement_,numberElements);
257      countRow_=CoinCopyOfArray(rhs.countRow_,numberElements);
258#endif
259    }
260    startFraction_ = rhs.startFraction_;
261    endFraction_ = rhs.endFraction_;
262    savedBestDj_ = rhs.savedBestDj_;
263    originalWanted_ = rhs.originalWanted_;
264    currentWanted_ = rhs.currentWanted_;
265    savedBestSequence_ = rhs.savedBestSequence_;
266    minimumObjectsScan_ = rhs.minimumObjectsScan_;
267    minimumGoodReducedCosts_ = rhs.minimumGoodReducedCosts_;
268  }
269  return *this;
270}
271// Sets model
272void 
273AbcMatrix::setModel(AbcSimplex * model)
274{
275  model_=model;
276  int numberColumns=model_->numberColumns();
277  bool needExtension=numberColumns>matrix_->getNumCols();
278  if (needExtension) {
279    CoinBigIndex lastElement = matrix_->getNumElements();
280    matrix_->reserve(numberColumns,lastElement,true);
281    CoinBigIndex * columnStart = matrix_->getMutableVectorStarts();
282    for (int i=numberColumns;i>=0;i--) {
283      if (columnStart[i]==0)
284        columnStart[i]=lastElement;
285      else
286        break;
287    }
288    assert (lastElement==columnStart[numberColumns]);
289  }
290}
291/* Returns a new matrix in reverse order without gaps */
292CoinPackedMatrix *
293AbcMatrix::reverseOrderedCopy() const
294{
295  CoinPackedMatrix * matrix = new CoinPackedMatrix();
296  matrix->setExtraGap(0.0);
297  matrix->setExtraMajor(0.0);
298  matrix->reverseOrderedCopyOf(*matrix_);
299  return matrix;
300}
301/// returns number of elements in column part of basis,
302CoinBigIndex
303AbcMatrix::countBasis( const int * whichColumn,
304                       int & numberColumnBasic)
305{
306  const int *  COIN_RESTRICT columnLength = matrix_->getVectorLengths();
307  CoinBigIndex numberElements = 0;
308  int numberRows=model_->numberRows();
309  // just count - can be over so ignore zero problem
310  for (int i = 0; i < numberColumnBasic; i++) {
311    int iColumn = whichColumn[i]-numberRows;
312    numberElements += columnLength[iColumn];
313  }
314  return numberElements;
315}
316void
317AbcMatrix::fillBasis(const int * COIN_RESTRICT whichColumn,
318                     int & numberColumnBasic,
319                     int * COIN_RESTRICT indexRowU,
320                     int * COIN_RESTRICT start,
321                     int * COIN_RESTRICT rowCount,
322                     int * COIN_RESTRICT columnCount,
323                     CoinFactorizationDouble * COIN_RESTRICT elementU)
324{
325  const int * COIN_RESTRICT columnLength = matrix_->getVectorLengths();
326  CoinBigIndex numberElements = start[0];
327  // fill
328  const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
329  const int * COIN_RESTRICT row = matrix_->getIndices();
330  const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
331  int numberRows=model_->numberRows();
332  for (int i = 0; i < numberColumnBasic; i++) {
333    int iColumn = whichColumn[i]-numberRows;
334    int length = columnLength[iColumn];
335    CoinBigIndex startThis = columnStart[iColumn];
336    columnCount[i] = length;
337    CoinBigIndex endThis = startThis + length;
338    for (CoinBigIndex j = startThis; j < endThis; j++) {
339      int iRow = row[j];
340      indexRowU[numberElements] = iRow;
341      rowCount[iRow]++;
342      assert (elementByColumn[j]);
343      elementU[numberElements++] = elementByColumn[j];
344    }
345    start[i+1] = numberElements;
346  }
347}
348#ifdef ABC_LONG_FACTORIZATION
349void
350AbcMatrix::fillBasis(const int * COIN_RESTRICT whichColumn,
351                     int & numberColumnBasic,
352                     int * COIN_RESTRICT indexRowU,
353                     int * COIN_RESTRICT start,
354                     int * COIN_RESTRICT rowCount,
355                     int * COIN_RESTRICT columnCount,
356                     long double * COIN_RESTRICT elementU)
357{
358  const int * COIN_RESTRICT columnLength = matrix_->getVectorLengths();
359  CoinBigIndex numberElements = start[0];
360  // fill
361  const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
362  const int * COIN_RESTRICT row = matrix_->getIndices();
363  const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
364  int numberRows=model_->numberRows();
365  for (int i = 0; i < numberColumnBasic; i++) {
366    int iColumn = whichColumn[i]-numberRows;
367    int length = columnLength[iColumn];
368    CoinBigIndex startThis = columnStart[iColumn];
369    columnCount[i] = length;
370    CoinBigIndex endThis = startThis + length;
371    for (CoinBigIndex j = startThis; j < endThis; j++) {
372      int iRow = row[j];
373      indexRowU[numberElements] = iRow;
374      rowCount[iRow]++;
375      assert (elementByColumn[j]);
376      elementU[numberElements++] = elementByColumn[j];
377    }
378    start[i+1] = numberElements;
379  }
380}
381#endif 
382#if 0 
383/// Move largest in column to beginning
384void
385AbcMatrix::moveLargestToStart()
386{
387  // get matrix data pointers
388  int *  COIN_RESTRICT row = matrix_->getMutableIndices();
389  const CoinBigIndex *  COIN_RESTRICT columnStart = matrix_->getVectorStarts();
390  double *  COIN_RESTRICT elementByColumn = matrix_->getMutableElements();
391  int numberColumns=model_->numberColumns();
392  CoinBigIndex start = 0;
393  for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
394    CoinBigIndex end = columnStart[iColumn+1];
395    double largest=0.0;
396    int position=-1;
397    for (CoinBigIndex j = start; j < end; j++) {
398      double value = fabs(elementByColumn[j]);
399      if (value>largest) {
400        largest=value;
401        position=j;
402      }
403    }
404    assert (position>=0); // ? empty column
405    if (position>start) {
406      double value=elementByColumn[start];
407      elementByColumn[start]=elementByColumn[position];
408      elementByColumn[position]=value;
409      int iRow=row[start];
410      row[start]=row[position];
411      row[position]=iRow;
412    }
413    start=end;
414  }
415}
416#endif
417// Creates row copy
418void 
419AbcMatrix::createRowCopy()
420{
421#if ABC_PARALLEL
422  if (model_->parallelMode()==0)
423#endif
424    numberRowBlocks_=1;
425#if ABC_PARALLEL
426  else
427    numberRowBlocks_=CoinMin(NUMBER_ROW_BLOCKS,model_->numberCpus());
428#endif
429  int maximumRows=model_->maximumAbcNumberRows();
430  int numberRows=model_->numberRows();
431  int numberColumns=model_->numberColumns();
432  int numberElements = matrix_->getNumElements();
433  assert (!rowStart_);
434  char * whichBlock_=new char [numberColumns];
435  rowStart_=new CoinBigIndex[numberRows*(numberRowBlocks_+2)];
436  element_ = new double [numberElements];
437  column_ = new int [numberElements];
438  const CoinBigIndex *  COIN_RESTRICT columnStart = matrix_->getVectorStarts();
439  memset(blockStart_,0,sizeof(blockStart_));
440  int ecount[10];
441  assert (numberRowBlocks_<16);
442  CoinAbcMemset0(ecount,10);
443  // allocate to blocks (put a bit less in first as will be dealing with slacks) LATER
444  CoinBigIndex start=0;
445  int block=0;
446  CoinBigIndex work=(2*numberColumns+matrix_->getNumElements()+numberRowBlocks_-1)/numberRowBlocks_;
447  CoinBigIndex thisWork=work;
448  for (int iColumn=0;iColumn<numberColumns;iColumn++) {
449    CoinBigIndex end=columnStart[iColumn+1];
450    assert (block>=0&&block<numberRowBlocks_);
451    whichBlock_[iColumn]=static_cast<char>(block);
452    thisWork -= 2+end-start;
453    ecount[block]+=end-start;
454    start=end;
455    blockStart_[block]++;
456    if (thisWork<=0) {
457      block++;
458      thisWork=work;
459    }
460  }
461#if 0
462  printf("Blocks ");
463  for (int i=0;i<numberRowBlocks_;i++)
464    printf("(%d %d) ",blockStart_[i],ecount[i]);
465  printf("\n");
466#endif
467  start=0;
468  for (int i=0;i<numberRowBlocks_;i++) {
469    int next=start+blockStart_[i];
470    blockStart_[i]=start;
471    start=next;
472  }
473  blockStart_[numberRowBlocks_]=start;
474  assert (start==numberColumns);
475  const int *  COIN_RESTRICT row = matrix_->getIndices();
476  const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
477  // counts
478  CoinAbcMemset0(rowStart_,numberRows*(numberRowBlocks_+2));
479  int * COIN_RESTRICT last = rowStart_+numberRows*(numberRowBlocks_+1);
480  for (int iColumn=0;iColumn<numberColumns;iColumn++) {
481    int block=whichBlock_[iColumn];
482    blockStart_[block]++;
483    int base=(block+1)*numberRows;
484    for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn+1];j++) {
485      int iRow=row[j];
486      rowStart_[base+iRow]++;
487      last[iRow]++;
488    }
489  }
490  // back to real starts
491  for (int iBlock=numberRowBlocks_-1;iBlock>=0;iBlock--) 
492    blockStart_[iBlock+1]=blockStart_[iBlock];
493  blockStart_[0]=0;
494  CoinBigIndex put=0;
495  for (int iRow=1;iRow<numberRows;iRow++) {
496    int n=last[iRow-1];
497    last[iRow-1]=put;
498    put += n;
499    rowStart_[iRow]=put;
500  }
501  int n=last[numberRows-1];
502  last[numberRows-1]=put;
503  put += n;
504  assert (put==numberElements);
505  //last[numberRows-1]=put;
506  // starts
507  int base=0;
508  for (int iBlock=0;iBlock<numberRowBlocks_;iBlock++) {
509    for (int iRow=0;iRow<numberRows;iRow++) {
510      rowStart_[base+numberRows+iRow]+=rowStart_[base+iRow];
511    }
512    base+= numberRows;
513  }
514  for (int iColumn=0;iColumn<numberColumns;iColumn++) {
515    int block=whichBlock_[iColumn];
516    int where=blockStart_[block];
517    blockStart_[block]++;
518    int base=block*numberRows;
519    for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn+1];j++) {
520      int iRow=row[j];
521      CoinBigIndex put=rowStart_[base+iRow];
522      rowStart_[base+iRow]++;
523      column_[put]=where+maximumRows;
524      element_[put]=elementByColumn[j];
525    }
526  }
527  // back to real starts etc
528  base=numberRows*numberRowBlocks_;
529  for (int iBlock=numberRowBlocks_-1;iBlock>=0;iBlock--)  {
530    blockStart_[iBlock+1]=blockStart_[iBlock]+maximumRows;
531    CoinAbcMemcpy(rowStart_+base,rowStart_+base-numberRows,numberRows);
532    base -= numberRows;
533  }
534  blockStart_[0]=0;//maximumRows;
535  delete [] whichBlock_;
536  // and move
537  CoinAbcMemcpy(rowStart_,last,numberRows);
538  // All in useful
539  CoinAbcMemcpy(rowStart_+(numberRowBlocks_+1)*numberRows,
540                rowStart_+(numberRowBlocks_)*numberRows,numberRows);
541#ifdef COUNT_COPY
542  // now blocked by element count
543  countRealColumn_=new int [numberColumns];
544  int counts [2*MAX_COUNT];
545  memset(counts,0,sizeof(counts));
546  //memset(countFirst_,0,sizeof(countFirst_));
547  int numberLarge=0;
548  for (int iColumn=0;iColumn<numberColumns;iColumn++) {
549    int n=columnStart[iColumn+1]-columnStart[iColumn];
550    if (n<MAX_COUNT)
551      counts[n]++;
552    else
553      numberLarge++;
554  }
555  CoinBigIndex numberExtra=3; // for alignment
556#define SMALL_COUNT 1
557  for (int i=0;i<MAX_COUNT;i++) {
558    int n=counts[i];
559    if (n>=SMALL_COUNT) {
560      n &= 3;
561      int extra=(4-n)&3;
562      numberExtra+= i*extra;
563    } else {
564      // treat as large
565      numberLarge+=n;
566    }
567  }
568  countElement_= new double [numberElements+numberExtra];
569  countRow_=new int [numberElements+numberExtra];
570  countStartLarge_ = new CoinBigIndex [numberLarge+1];
571  countStartLarge_[numberLarge]=numberElements+numberExtra;
572  //return;
573  CoinInt64 xx = reinterpret_cast<CoinInt64>(countElement_);
574  int iBottom = static_cast<int>(xx & 31);
575  int offset = iBottom>>3;
576  CoinBigIndex firstElementLarge=0;
577  if (offset)
578    firstElementLarge=4-offset;
579  //countStart_[0]=firstElementLarge;
580  int positionLarge=0;
581  smallestCount_=0;
582  largestCount_=0;
583  for (int i=0;i<MAX_COUNT;i++) {
584    int n=counts[i];
585    countFirst_[i]=positionLarge;
586    countStart_[i]=firstElementLarge;
587    if (n>=SMALL_COUNT) {
588      counts[i+MAX_COUNT]=1;
589      if (smallestCount_==0)
590        smallestCount_=i;
591      largestCount_=i;
592      positionLarge+=n;
593      firstElementLarge+=n*i;
594      n &= 3;
595      int extra=(4-n)&3;
596      firstElementLarge+= i*extra;
597    }
598    counts[i]=0;
599  }
600  largestCount_++;
601  countFirst_[MAX_COUNT]=positionLarge;
602  countStart_[MAX_COUNT]=firstElementLarge;
603  numberLarge=0;
604  for (int iColumn=0;iColumn<numberColumns;iColumn++) {
605    CoinBigIndex start=columnStart[iColumn];
606    int n=columnStart[iColumn+1]-start;
607    CoinBigIndex put;
608    int position;
609    if (n<MAX_COUNT&&counts[MAX_COUNT+n]!=0) {
610      int iWhich=counts[n];
611      position = countFirst_[n]+iWhich;
612      int iBlock4=iWhich&(~3);
613      iWhich &= 3;
614      put = countStart_[n]+iBlock4*n;
615      put += iWhich;
616      counts[n]++;
617      for (int i=0;i<n;i++) {
618        countRow_[put]=row[start+i];
619        countElement_[put]=elementByColumn[start+i];
620        put+=4;
621      }
622    } else {
623      position = positionLarge+numberLarge;
624      put=firstElementLarge;
625      countStartLarge_[numberLarge]=put;
626      firstElementLarge+=n;
627      numberLarge++;
628      CoinAbcMemcpy(countRow_+put,row+start,n);
629      CoinAbcMemcpy(countElement_+put,elementByColumn+start,n);
630    }
631    countRealColumn_[position]=iColumn+maximumRows;
632  }
633  countStartLarge_[numberLarge]=firstElementLarge;
634  assert (firstElementLarge<=numberElements+numberExtra);
635#endif
636}
637// Make all useful
638void 
639AbcMatrix::makeAllUseful(CoinIndexedVector & /*spare*/)
640{
641  int numberRows=model_->numberRows();
642  CoinBigIndex * COIN_RESTRICT rowEnd = rowStart_+numberRows;
643  const CoinBigIndex * COIN_RESTRICT rowReallyEnd = rowStart_+2*numberRows;
644  for (int iRow=0;iRow<numberRows;iRow++) {
645    rowEnd[iRow]=rowReallyEnd[iRow];
646  }
647}
648#define SQRT_ARRAY
649// Creates scales for column copy (rowCopy in model may be modified)
650void
651AbcMatrix::scale(int numberAlreadyScaled)
652{
653  int numberRows=model_->numberRows();
654  bool inBranchAndBound=(model_->specialOptions(),0x1000000)!=0;
655  bool doScaling=numberAlreadyScaled>=0;
656  if (!doScaling)
657    numberAlreadyScaled=0;
658  if (numberAlreadyScaled==numberRows)
659    return; // no need to do anything
660  int numberColumns=model_->numberColumns();
661  double *  COIN_RESTRICT rowScale=model_->rowScale2();
662  double *  COIN_RESTRICT inverseRowScale=model_->inverseRowScale2();
663  double *  COIN_RESTRICT columnScale=model_->columnScale2();
664  double *  COIN_RESTRICT inverseColumnScale=model_->inverseColumnScale2();
665  // we are going to mark bits we are interested in
666  int whichArray=model_->getAvailableArrayPublic();
667  char *  COIN_RESTRICT usefulColumn = reinterpret_cast<char *>(model_->usefulArray(whichArray)->getIndices());
668  memset(usefulColumn,1,numberColumns);
669  const double *  COIN_RESTRICT rowLower = model_->rowLower();
670  const double *  COIN_RESTRICT rowUpper = model_->rowUpper();
671  const double *  COIN_RESTRICT columnLower = model_->columnLower();
672  const double *  COIN_RESTRICT columnUpper = model_->columnUpper();
673  //#define LEAVE_FIXED
674  // mark empty and fixed columns
675  // get matrix data pointers
676  const int *  COIN_RESTRICT row = matrix_->getIndices();
677  const CoinBigIndex *  COIN_RESTRICT columnStart = matrix_->getVectorStarts();
678  double *  COIN_RESTRICT elementByColumn = matrix_->getMutableElements();
679  CoinPackedMatrix *  COIN_RESTRICT rowCopy = reverseOrderedCopy();
680  const int * column = rowCopy->getIndices();
681  const CoinBigIndex *  COIN_RESTRICT rowStart = rowCopy->getVectorStarts();
682  const double *  COIN_RESTRICT element = rowCopy->getElements();
683  assert (numberAlreadyScaled>=0&&numberAlreadyScaled<numberRows);
684#ifndef LEAVE_FIXED
685  for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
686    if (columnUpper[iColumn] == columnLower[iColumn])
687      usefulColumn[iColumn]=0;
688  }
689#endif
690  double overallLargest = -1.0e-20;
691  double overallSmallest = 1.0e20;
692  if (!numberAlreadyScaled) {
693    // may be redundant
694    CoinFillN(rowScale,numberRows,1.0);
695    CoinFillN(columnScale,numberColumns,1.0);
696    CoinBigIndex start = 0;
697    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
698      CoinBigIndex end = columnStart[iColumn+1];
699      if (usefulColumn[iColumn]) {
700        if (end>start) {
701          for (CoinBigIndex j = start; j < end; j++) {
702            double value = fabs(elementByColumn[j]);
703            overallLargest = CoinMax(overallLargest, value);
704            overallSmallest = CoinMin(overallSmallest, value);
705          }
706        } else {
707          usefulColumn[iColumn]=0;
708        }
709      }
710      start=end;
711    }
712  } else {
713    CoinBigIndex start = rowStart[numberAlreadyScaled];
714    for (int iRow = numberAlreadyScaled; iRow < numberRows; iRow++) {
715      rowScale[iRow]=1.0;
716      CoinBigIndex end = rowStart[iRow+1];
717      for (CoinBigIndex j = start; j < end; j++) {
718        int iColumn=column[j];
719        if (usefulColumn[iColumn]) {
720          double value = fabs(elementByColumn[j])*columnScale[iColumn];
721          overallLargest = CoinMax(overallLargest, value);
722          overallSmallest = CoinMin(overallSmallest, value);
723        }
724      }
725    }
726  }
727  if ((overallSmallest >= 0.5 && overallLargest <= 2.0)||!doScaling) {
728    //printf("no scaling\n");
729    delete rowCopy;
730    model_->clearArraysPublic(whichArray);
731    CoinFillN(inverseRowScale+numberAlreadyScaled,numberRows-numberAlreadyScaled,1.0);
732    if (!numberAlreadyScaled) 
733      CoinFillN(inverseColumnScale,numberColumns,1.0);
734    //moveLargestToStart();
735    return ;
736  }
737  // need to scale
738  double largest;
739  double smallest;
740  int scalingMethod=model_->scalingFlag();
741  if (scalingMethod == 4) {
742    // As auto
743    scalingMethod = 3;
744  } else if (scalingMethod == 5) {
745    // As geometric
746    scalingMethod = 2;
747  }
748  double savedOverallRatio = 0.0;
749  double tolerance = 5.0 * model_->primalTolerance();
750  bool finished = false;
751  // if scalingMethod 3 then may change
752  bool extraDetails = (model_->logLevel() > 2);
753  bool secondTime=false;
754  while (!finished) {
755    int numberPass = !numberAlreadyScaled ? 3 : 1;
756    overallLargest = -1.0e-20;
757    overallSmallest = 1.0e20;
758    if (!secondTime) {
759      secondTime=true;
760    } else {
761      CoinFillN ( rowScale, numberRows, 1.0);
762      CoinFillN ( columnScale, numberColumns, 1.0);
763    } 
764    if (scalingMethod == 1 || scalingMethod == 3) {
765      // Maximum in each row
766      for (int iRow = numberAlreadyScaled; iRow < numberRows; iRow++) {
767        largest = 1.0e-10;
768        for (CoinBigIndex j = rowStart[iRow]; j < rowStart[iRow+1]; j++) {
769          int iColumn = column[j];
770          if (usefulColumn[iColumn]) {
771            double value = fabs(element[j]);
772            largest = CoinMax(largest, value);
773            assert (largest < 1.0e40);
774          }
775        }
776        rowScale[iRow] = 1.0 / largest;
777#ifdef COIN_DEVELOP
778        if (extraDetails) {
779          overallLargest = CoinMax(overallLargest, largest);
780          overallSmallest = CoinMin(overallSmallest, largest);
781        }
782#endif
783      }
784    } else {
785      while (numberPass) {
786        overallLargest = 0.0;
787        overallSmallest = 1.0e50;
788        numberPass--;
789        // Geometric mean on row scales
790        for (int iRow = numberAlreadyScaled; iRow < numberRows; iRow++) {
791          largest = 1.0e-50;
792          smallest = 1.0e50;
793          for (CoinBigIndex j = rowStart[iRow]; j < rowStart[iRow+1]; j++) {
794            int iColumn = column[j];
795            if (usefulColumn[iColumn]) {
796              double value = fabs(element[j]);
797              value *= columnScale[iColumn];
798              largest = CoinMax(largest, value);
799              smallest = CoinMin(smallest, value);
800            }
801          }
802#ifdef SQRT_ARRAY
803          rowScale[iRow] = smallest * largest;
804#else
805          rowScale[iRow] = 1.0 / sqrt(smallest * largest);
806#endif
807          if (extraDetails) {
808            overallLargest = CoinMax(largest * rowScale[iRow], overallLargest);
809            overallSmallest = CoinMin(smallest * rowScale[iRow], overallSmallest);
810          }
811        }
812        if (model_->scalingFlag() == 5)
813          break; // just scale rows
814#ifdef SQRT_ARRAY
815        CoinAbcInverseSqrts(rowScale, numberRows);
816#endif
817        if (!inBranchAndBound)
818        model_->messageHandler()->message(CLP_PACKEDSCALE_WHILE, *model_->messagesPointer())
819          << overallSmallest
820          << overallLargest
821          << CoinMessageEol;
822        // skip last column round
823        if (numberPass == 1)
824          break;
825        // Geometric mean on column scales
826        for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
827          if (usefulColumn[iColumn]) {
828            largest = 1.0e-50;
829            smallest = 1.0e50;
830            for (CoinBigIndex j = columnStart[iColumn];
831                 j < columnStart[iColumn+1]; j++) {
832              int iRow = row[j];
833              double value = fabs(elementByColumn[j]);
834              value *= rowScale[iRow];
835              largest = CoinMax(largest, value);
836              smallest = CoinMin(smallest, value);
837            }
838#ifdef SQRT_ARRAY
839            columnScale[iColumn] = smallest * largest;
840#else
841            columnScale[iColumn] = 1.0 / sqrt(smallest * largest);
842#endif
843          }
844        }
845#ifdef SQRT_ARRAY
846        CoinAbcInverseSqrts(columnScale, numberColumns);
847#endif
848      }
849    }
850    // If ranges will make horrid then scale
851    for (int iRow = numberAlreadyScaled; iRow < numberRows; iRow++) {
852      double difference = rowUpper[iRow] - rowLower[iRow];
853      double scaledDifference = difference * rowScale[iRow];
854      if (scaledDifference > tolerance && scaledDifference < 1.0e-4) {
855        // make gap larger
856        rowScale[iRow] *= 1.0e-4 / scaledDifference;
857        rowScale[iRow] = CoinMax(1.0e-10, CoinMin(1.0e10, rowScale[iRow]));
858        //printf("Row %d difference %g scaled diff %g => %g\n",iRow,difference,
859        // scaledDifference,difference*rowScale[iRow]);
860      }
861    }
862    // final pass to scale columns so largest is reasonable
863    // See what smallest will be if largest is 1.0
864    if (model_->scalingFlag() != 5) {
865      overallSmallest = 1.0e50;
866      for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
867        if (usefulColumn[iColumn]) {
868          largest = 1.0e-20;
869          smallest = 1.0e50;
870          for (CoinBigIndex j = columnStart[iColumn];j < columnStart[iColumn+1]; j++) {
871            int iRow = row[j];
872            double value = fabs(elementByColumn[j] * rowScale[iRow]);
873            largest = CoinMax(largest, value);
874            smallest = CoinMin(smallest, value);
875          }
876          if (overallSmallest * largest > smallest)
877            overallSmallest = smallest / largest;
878        }
879      }
880    }
881    if (scalingMethod == 1 || scalingMethod == 2) {
882      finished = true;
883    } else if (savedOverallRatio == 0.0 && scalingMethod != 4) {
884      savedOverallRatio = overallSmallest;
885      scalingMethod = 4;
886    } else {
887      assert (scalingMethod == 4);
888      if (overallSmallest > 2.0 * savedOverallRatio) {
889        finished = true; // geometric was better
890        if (model_->scalingFlag() == 4) {
891          // If in Branch and bound change
892          if ((model_->specialOptions() & 1024) != 0) {
893            model_->scaling(2);
894          }
895        }
896      } else {
897        scalingMethod = 1; // redo equilibrium
898        if (model_->scalingFlag() == 4) {
899          // If in Branch and bound change
900          if ((model_->specialOptions() & 1024) != 0) {
901            model_->scaling(1);
902          }
903        }
904      }
905#if 0
906      if (extraDetails) {
907        if (finished)
908          printf("equilibrium ratio %g, geometric ratio %g , geo chosen\n",
909                 savedOverallRatio, overallSmallest);
910        else
911          printf("equilibrium ratio %g, geometric ratio %g , equi chosen\n",
912                 savedOverallRatio, overallSmallest);
913      }
914#endif
915    }
916  }
917  //#define RANDOMIZE
918#ifdef RANDOMIZE
919  // randomize by up to 10%
920  for (int iRow = numberAlreadyScaled; iRow < numberRows; iRow++) {
921    double value = 0.5 - randomNumberGenerator_.randomDouble(); //between -0.5 to + 0.5
922    rowScale[iRow] *= (1.0 + 0.1 * value);
923  }
924#endif
925  overallLargest = 1.0;
926  if (overallSmallest < 1.0e-1)
927    overallLargest = 1.0 / sqrt(overallSmallest);
928  overallLargest = CoinMin(100.0, overallLargest);
929  overallSmallest = 1.0e50;
930  char *  COIN_RESTRICT usedRow = reinterpret_cast<char *>(inverseRowScale);
931  memset(usedRow, 0, numberRows);
932  //printf("scaling %d\n",model_->scalingFlag());
933  if (model_->scalingFlag() != 5) {
934    for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
935      if (usefulColumn[iColumn]) {
936        largest = 1.0e-20;
937        smallest = 1.0e50;
938        for (CoinBigIndex j = columnStart[iColumn];j < columnStart[iColumn+1]; j++) {
939          int iRow = row[j];
940          usedRow[iRow] = 1;
941          double value = fabs(elementByColumn[j] * rowScale[iRow]);
942          largest = CoinMax(largest, value);
943          smallest = CoinMin(smallest, value);
944        }
945        columnScale[iColumn] = overallLargest / largest;
946        //columnScale[iColumn]=CoinMax(1.0e-10,CoinMin(1.0e10,columnScale[iColumn]));
947#ifdef RANDOMIZE
948        if (!numberAlreadyScaled) {
949          double value = 0.5 - randomNumberGenerator_.randomDouble(); //between -0.5 to + 0.5
950          columnScale[iColumn] *= (1.0 + 0.1 * value);
951        }
952#endif
953        double difference = columnUpper[iColumn] - columnLower[iColumn];
954        if (difference < 1.0e-5 * columnScale[iColumn]) {
955          // make gap larger
956          columnScale[iColumn] = difference / 1.0e-5;
957          //printf("Column %d difference %g scaled diff %g => %g\n",iColumn,difference,
958          // scaledDifference,difference*columnScale[iColumn]);
959        }
960        double value = smallest * columnScale[iColumn];
961        if (overallSmallest > value)
962          overallSmallest = value;
963        //overallSmallest = CoinMin(overallSmallest,smallest*columnScale[iColumn]);
964      } else {
965        assert(columnScale[iColumn] == 1.0);
966        //columnScale[iColumn]=1.0;
967      }
968    }
969    for (int iRow = numberAlreadyScaled; iRow < numberRows; iRow++) {
970      if (!usedRow[iRow]) {
971        rowScale[iRow] = 1.0;
972      }
973    }
974  }
975  if (!inBranchAndBound)
976  model_->messageHandler()->message(CLP_PACKEDSCALE_FINAL, *model_->messagesPointer())
977    << overallSmallest
978    << overallLargest
979    << CoinMessageEol;
980  if (overallSmallest < 1.0e-13) {
981    // Change factorization zero tolerance
982    double newTolerance = CoinMax(1.0e-15 * (overallSmallest / 1.0e-13),
983                                  1.0e-18);
984    if (model_->factorization()->zeroTolerance() > newTolerance)
985      model_->factorization()->zeroTolerance(newTolerance);
986    newTolerance = CoinMax(overallSmallest * 0.5, 1.0e-18);
987    model_->setZeroTolerance(newTolerance);
988#ifndef NDEBUG
989    assert (newTolerance<0.0); // just so we can fix
990#endif
991  }
992  // make copy (could do faster by using previous values)
993  // could just do partial
994  CoinAbcReciprocal(inverseRowScale+numberAlreadyScaled,numberRows-numberAlreadyScaled,
995                rowScale+numberAlreadyScaled);
996  if (!numberAlreadyScaled) 
997    CoinAbcReciprocal(inverseColumnScale,numberColumns,columnScale);
998  // Do scaled copy //NO and move largest to start
999  CoinBigIndex start = 0;
1000  for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
1001    double scale=columnScale[iColumn];
1002    CoinBigIndex end = columnStart[iColumn+1];
1003    for (CoinBigIndex j = start; j < end; j++) {
1004      double value=elementByColumn[j];
1005      int iRow = row[j];
1006      if (iRow>=numberAlreadyScaled) {
1007        value*= scale*rowScale[iRow];
1008        elementByColumn[j]= value;
1009      }
1010    }
1011    start=end;
1012  }
1013  delete rowCopy;
1014#if 0
1015  if (model_->rowCopy()) {
1016    // need to replace row by row
1017    CoinPackedMatrix * rowCopy = NULL;
1018    //static_cast< AbcMatrix*>(model_->rowCopy());
1019    double * element = rowCopy->getMutableElements();
1020    const int * column = rowCopy->getIndices();
1021    const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
1022    // scale row copy
1023    for (iRow = 0; iRow < numberRows; iRow++) {
1024      CoinBigIndex j;
1025      double scale = rowScale[iRow];
1026      double * elementsInThisRow = element + rowStart[iRow];
1027      const int * columnsInThisRow = column + rowStart[iRow];
1028      int number = rowStart[iRow+1] - rowStart[iRow];
1029      assert (number <= numberColumns);
1030      for (j = 0; j < number; j++) {
1031        int iColumn = columnsInThisRow[j];
1032        elementsInThisRow[j] *= scale * columnScale[iColumn];
1033      }
1034    }
1035  }
1036#endif
1037  model_->clearArraysPublic(whichArray);
1038}
1039/* Return <code>y + A * scalar *x</code> in <code>y</code>.
1040   @pre <code>x</code> must be of size <code>numColumns()</code>
1041   @pre <code>y</code> must be of size <code>numRows()</code> */
1042//scaled versions
1043void
1044AbcMatrix::timesModifyExcludingSlacks(double scalar,
1045                                      const double * x, double * y) const
1046{
1047  int numberTotal = model_->numberTotal();
1048  int maximumRows = model_->maximumAbcNumberRows();
1049  // get matrix data pointers
1050  const int *  COIN_RESTRICT row = matrix_->getIndices();
1051  const CoinBigIndex *  COIN_RESTRICT columnStart = matrix_->getVectorStarts()-maximumRows;
1052  const double *  COIN_RESTRICT elementByColumn = matrix_->getElements();
1053  for (int iColumn = maximumRows; iColumn < numberTotal; iColumn++) {
1054    double value = x[iColumn];
1055    if (value) {
1056      CoinBigIndex start = columnStart[iColumn];
1057      CoinBigIndex end = columnStart[iColumn+1];
1058      value *= scalar;
1059      for (CoinBigIndex j = start; j < end; j++) {
1060        int iRow = row[j];
1061        y[iRow] += value * elementByColumn[j];
1062      }
1063    }
1064  }
1065}
1066/* Return <code>y + A * scalar(+-1) *x</code> in <code>y</code>.
1067   @pre <code>x</code> must be of size <code>numColumns()+numRows()</code>
1068   @pre <code>y</code> must be of size <code>numRows()</code> */
1069void 
1070AbcMatrix::timesModifyIncludingSlacks(double scalar,
1071                                      const double * x, double * y) const
1072{
1073  int numberRows=model_->numberRows();
1074  int numberTotal = model_->numberTotal();
1075  int maximumRows = model_->maximumAbcNumberRows();
1076  // makes no sense for x==y??
1077  assert (x!=y);
1078  // For now just by column and assumes already scaled (use reallyScale)
1079  // get matrix data pointers
1080  const int * COIN_RESTRICT row = matrix_->getIndices();
1081  const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts()-maximumRows;
1082  const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
1083  if (scalar==1.0) {
1084    // add
1085    if (x!=y) {
1086      for (int i=0;i<numberRows;i++) 
1087        y[i]+=x[i];
1088    } else {
1089      for (int i=0;i<numberRows;i++) 
1090        y[i]+=y[i];
1091    }
1092    for (int iColumn = maximumRows; iColumn < numberTotal; iColumn++) {
1093      double value = x[iColumn];
1094      if (value) {
1095        CoinBigIndex start = columnStart[iColumn];
1096        CoinBigIndex next = columnStart[iColumn+1];
1097        for (CoinBigIndex j = start; j < next; j++) {
1098          int jRow = row[j];
1099          y[jRow]+=value*elementByColumn[j];
1100        }
1101      }
1102    }
1103  } else {
1104    // subtract
1105    if (x!=y) {
1106      for (int i=0;i<numberRows;i++) 
1107        y[i]-=x[i];
1108    } else {
1109      for (int i=0;i<numberRows;i++) 
1110        y[i]=0.0;
1111    }
1112    for (int iColumn = maximumRows; iColumn < numberTotal; iColumn++) {
1113      double value = x[iColumn];
1114      if (value) {
1115        CoinBigIndex start = columnStart[iColumn];
1116        CoinBigIndex next = columnStart[iColumn+1];
1117        for (CoinBigIndex j = start; j < next; j++) {
1118          int jRow = row[j];
1119          y[jRow]-=value*elementByColumn[j];
1120        }
1121      }
1122    }
1123  }
1124}
1125/* Return A * scalar(+-1) *x</code> in <code>y</code>.
1126   @pre <code>x</code> must be of size <code>numColumns()+numRows()</code>
1127   @pre <code>y</code> must be of size <code>numRows()</code> */
1128void 
1129AbcMatrix::timesIncludingSlacks(double scalar,
1130                                const double * x, double * y) const
1131{
1132  int numberRows=model_->numberRows();
1133  int numberTotal = model_->numberTotal();
1134  int maximumRows = model_->maximumAbcNumberRows();
1135  // For now just by column and assumes already scaled (use reallyScale)
1136  // get matrix data pointers
1137  const int * COIN_RESTRICT row = matrix_->getIndices();
1138  const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts()-maximumRows;
1139  const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
1140  if (scalar==1.0) {
1141    // add
1142    if (x!=y) {
1143      for (int i=0;i<numberRows;i++) 
1144        y[i]=x[i];
1145    } 
1146    for (int iColumn = maximumRows; iColumn < numberTotal; iColumn++) {
1147      double value = x[iColumn];
1148      if (value) {
1149        CoinBigIndex start = columnStart[iColumn];
1150        CoinBigIndex next = columnStart[iColumn+1];
1151        for (CoinBigIndex j = start; j < next; j++) {
1152          int jRow = row[j];
1153          y[jRow]+=value*elementByColumn[j];
1154        }
1155      }
1156    }
1157  } else {
1158    // subtract
1159    if (x!=y) {
1160      for (int i=0;i<numberRows;i++) 
1161        y[i]=-x[i];
1162    } else {
1163      for (int i=0;i<numberRows;i++) 
1164        y[i]=-y[i];
1165    }
1166    for (int iColumn = maximumRows; iColumn < numberTotal; iColumn++) {
1167      double value = x[iColumn];
1168      if (value) {
1169        CoinBigIndex start = columnStart[iColumn];
1170        CoinBigIndex next = columnStart[iColumn+1];
1171        for (CoinBigIndex j = start; j < next; j++) {
1172          int jRow = row[j];
1173          y[jRow]-=value*elementByColumn[j];
1174        }
1175      }
1176    }
1177  }
1178}
1179extern double parallelDual4(AbcSimplexDual * );
1180static double firstPass(AbcSimplex * model,int iBlock,
1181                 double upperTheta, int & freeSequence,
1182                 CoinPartitionedVector & tableauRow,
1183                 CoinPartitionedVector & candidateList) 
1184{
1185  int *  COIN_RESTRICT index = tableauRow.getIndices();
1186  double *  COIN_RESTRICT array = tableauRow.denseVector();
1187  double *  COIN_RESTRICT arrayCandidate=candidateList.denseVector();
1188  int *  COIN_RESTRICT indexCandidate = candidateList.getIndices();
1189  const double *  COIN_RESTRICT abcDj = model->djRegion();
1190  const unsigned char *  COIN_RESTRICT internalStatus = model->internalStatus();
1191  // do first pass to get possibles
1192  double bestPossible = 0.0;
1193  // We can also see if infeasible or pivoting on free
1194  double tentativeTheta = 1.0e25; // try with this much smaller as guess
1195  double acceptablePivot = model->currentAcceptablePivot();
1196  double dualT=-model->currentDualTolerance();
1197  // fixed will have been taken out by now
1198  const double multiplier[] = { 1.0, -1.0};
1199  freeSequence=-1;
1200  int firstIn=model->abcMatrix()->blockStart(iBlock);
1201  int numberNonZero=tableauRow.getNumElements(iBlock)+firstIn;
1202  int numberRemaining=firstIn;
1203  //first=tableauRow.getNumElements();
1204  // could pass in last and if numberNonZero==last-firstIn scan as well
1205  if (model->ordinaryVariables()) {
1206    for (int i = firstIn; i < numberNonZero; i++) {
1207      int iSequence = index[i];
1208      double tableauValue=array[i];
1209      unsigned char iStatus=internalStatus[iSequence]&7;
1210      double mult = multiplier[iStatus];
1211      double alpha = tableauValue * mult;
1212      double oldValue = abcDj[iSequence] * mult;
1213      double value = oldValue - tentativeTheta * alpha;
1214      if (value < dualT) {
1215        bestPossible = CoinMax(bestPossible, alpha);
1216        value = oldValue - upperTheta * alpha;
1217        if (value < dualT && alpha >= acceptablePivot) {
1218          upperTheta = (oldValue - dualT) / alpha;
1219        }
1220        // add to list
1221        arrayCandidate[numberRemaining] = alpha;
1222        indexCandidate[numberRemaining++] = iSequence;
1223      }
1224    }
1225  } else {
1226    double badFree = 0.0;
1227    double freeAlpha = model->currentAcceptablePivot();
1228    int freeSequenceIn=model->freeSequenceIn();
1229    double currentDualTolerance = model->currentDualTolerance();
1230    for (int i = firstIn; i < numberNonZero; i++) {
1231      int iSequence = index[i];
1232      double tableauValue=array[i];
1233      unsigned char iStatus=internalStatus[iSequence]&7;
1234      if ((iStatus&4)==0) {
1235        double mult = multiplier[iStatus];
1236        double alpha = tableauValue * mult;
1237        double oldValue = abcDj[iSequence] * mult;
1238        double value = oldValue - tentativeTheta * alpha;
1239        if (value < dualT) {
1240          bestPossible = CoinMax(bestPossible, alpha);
1241          value = oldValue - upperTheta * alpha;
1242          if (value < dualT && alpha >= acceptablePivot) {
1243            upperTheta = (oldValue - dualT) / alpha;
1244          }
1245          // add to list
1246          arrayCandidate[numberRemaining] = alpha;
1247          indexCandidate[numberRemaining++] = iSequence;
1248        }
1249      } else {
1250        bool keep;
1251        bestPossible = CoinMax(bestPossible, fabs(tableauValue));
1252        double oldValue = abcDj[iSequence];
1253        // If free has to be very large - should come in via dualRow
1254        //if (getInternalStatus(iSequence+addSequence)==isFree&&fabs(tableauValue)<1.0e-3)
1255        //break;
1256        if (oldValue > currentDualTolerance) {
1257          keep = true;
1258        } else if (oldValue < -currentDualTolerance) {
1259          keep = true;
1260        } else {
1261          if (fabs(tableauValue) > CoinMax(10.0 * acceptablePivot, 1.0e-5)) {
1262            keep = true;
1263          } else {
1264            keep = false;
1265            badFree = CoinMax(badFree, fabs(tableauValue));
1266          }
1267        }
1268        if (keep) {
1269#ifdef PAN
1270          if (model->fakeSuperBasic(iSequence)>=0) {
1271#endif
1272            if (iSequence==freeSequenceIn)
1273              tableauValue=COIN_DBL_MAX;
1274            // free - choose largest
1275            if (fabs(tableauValue) > fabs(freeAlpha)) {
1276              freeAlpha = tableauValue;
1277              freeSequence = iSequence;
1278            }
1279#ifdef PAN
1280          }
1281#endif
1282        }
1283      }
1284    }
1285  }
1286  //firstInX=numberNonZero-firstIn;
1287  //lastInX=-1;//numberRemaining-lastInX;
1288  tableauRow.setNumElementsPartition(iBlock,numberNonZero-firstIn);
1289  candidateList.setNumElementsPartition(iBlock,numberRemaining-firstIn);
1290  return upperTheta;
1291}
1292// gets sorted tableau row and a possible value of theta
1293double 
1294AbcMatrix::dualColumn1Row(int iBlock, double upperTheta, int & freeSequence,
1295                          const CoinIndexedVector & update,
1296                          CoinPartitionedVector & tableauRow,
1297                          CoinPartitionedVector & candidateList) const
1298{
1299  int maximumRows = model_->maximumAbcNumberRows();
1300  int number=update.getNumElements();
1301  const double *  COIN_RESTRICT pi=update.denseVector();
1302  const int *  COIN_RESTRICT piIndex = update.getIndices();
1303  const CoinBigIndex * COIN_RESTRICT rowStart = rowStart_;
1304  int numberRows=model_->numberRows();
1305  const CoinBigIndex * COIN_RESTRICT rowEnd = rowStart+numberRows*numberRowBlocks_;
1306  // count down
1307  int nColumns;
1308  int firstIn=blockStart_[iBlock];
1309  int first=firstIn;
1310  if (!first)
1311    first=maximumRows;
1312  int last=blockStart_[iBlock+1];
1313  nColumns=last-first;
1314  int target=nColumns;
1315  rowStart += iBlock*numberRows;
1316  rowEnd = rowStart+numberRows;
1317  for (int i=0;i<number;i++) {
1318    int iRow=piIndex[i];
1319    CoinBigIndex end = rowEnd[iRow];
1320    target -= end-rowStart[iRow];
1321    if (target<0)
1322      break;
1323  }
1324  if (target>0) {
1325    //printf("going to few %d ops %d\n",number,nColumns-target);
1326    return dualColumn1RowFew(iBlock, upperTheta, freeSequence,
1327                          update, tableauRow, candidateList);
1328  }
1329  const unsigned char *  COIN_RESTRICT internalStatus = model_->internalStatus();
1330  int *  COIN_RESTRICT index = tableauRow.getIndices();
1331  double *  COIN_RESTRICT array = tableauRow.denseVector();
1332  const double * COIN_RESTRICT element = element_;
1333  const int * COIN_RESTRICT column = column_;
1334  for (int i=0;i<number;i++) {
1335    int iRow=piIndex[i];
1336    double piValue = pi[iRow];
1337    CoinBigIndex end = rowEnd[iRow];
1338    for (CoinBigIndex j=rowStart[iRow];j<end;j++) {
1339      int iColumn = column[j];
1340      assert (iColumn>=first&&iColumn<last);
1341      double value = element[j];
1342      array[iColumn]+= piValue*value;
1343    }
1344  }
1345  int numberNonZero;
1346  int numberRemaining;
1347  if (iBlock==0) {
1348#if 1
1349    upperTheta=static_cast<AbcSimplexDual *>(model_)->dualColumn1A();
1350    numberNonZero=tableauRow.getNumElements(0);
1351    numberRemaining=candidateList.getNumElements(0);
1352#else
1353    numberNonZero=0;
1354    for (int i=0;i<number;i++) {
1355      int iRow=piIndex[i];
1356      unsigned char type=internalStatus[iRow];
1357      if ((type&7)<6) {
1358        index[numberNonZero]=iRow;
1359        double piValue=pi[iRow];
1360        array[numberNonZero++]=piValue;
1361      }
1362    }
1363    numberRemaining=0;
1364#endif
1365  } else {
1366    numberNonZero=firstIn;
1367    numberRemaining=firstIn;
1368  }
1369  double *  COIN_RESTRICT arrayCandidate=candidateList.denseVector();
1370  int *  COIN_RESTRICT indexCandidate = candidateList.getIndices();
1371  //printf("first %d last %d firstIn %d lastIn %d\n",
1372  //     first,last,firstIn,lastIn);
1373  const double *  COIN_RESTRICT abcDj = model_->djRegion();
1374  // do first pass to get possibles
1375  double bestPossible = 0.0;
1376  // We can also see if infeasible or pivoting on free
1377  double tentativeTheta = 1.0e25; // try with this much smaller as guess
1378  double acceptablePivot = model_->currentAcceptablePivot();
1379  double dualT=-model_->currentDualTolerance();
1380  const double multiplier[] = { 1.0, -1.0};
1381  double zeroTolerance = model_->zeroTolerance();
1382  freeSequence=-1;
1383  if (model_->ordinaryVariables()) {
1384    for (int iSequence = first; iSequence < last; iSequence++) {
1385      double tableauValue=array[iSequence];
1386      if (tableauValue) {
1387        array[iSequence]=0.0;
1388        if (fabs(tableauValue)>zeroTolerance) {
1389          unsigned char iStatus=internalStatus[iSequence]&7;
1390          if (iStatus<4) {
1391            index[numberNonZero]=iSequence;
1392            array[numberNonZero++]=tableauValue;
1393            double mult = multiplier[iStatus];
1394            double alpha = tableauValue * mult;
1395            double oldValue = abcDj[iSequence] * mult;
1396            double value = oldValue - tentativeTheta * alpha;
1397            if (value < dualT) {
1398              bestPossible = CoinMax(bestPossible, alpha);
1399              value = oldValue - upperTheta * alpha;
1400              if (value < dualT && alpha >= acceptablePivot) {
1401                upperTheta = (oldValue - dualT) / alpha;
1402              }
1403              // add to list
1404              arrayCandidate[numberRemaining] = alpha;
1405              indexCandidate[numberRemaining++] = iSequence;
1406            }
1407          }
1408        }
1409      }
1410    }
1411  } else {
1412    double badFree = 0.0;
1413    double freeAlpha = model_->currentAcceptablePivot();
1414    int freeSequenceIn=model_->freeSequenceIn();
1415    //printf("block %d freeSequence %d acceptable %g\n",iBlock,freeSequenceIn,freeAlpha);
1416    double currentDualTolerance = model_->currentDualTolerance();
1417    for (int iSequence = first; iSequence < last; iSequence++) {
1418      double tableauValue=array[iSequence];
1419      if (tableauValue) {
1420        array[iSequence]=0.0;
1421        if (fabs(tableauValue)>zeroTolerance) {
1422          unsigned char iStatus=internalStatus[iSequence]&7;
1423          if (iStatus<6) {
1424            if ((iStatus&4)==0) {
1425              index[numberNonZero]=iSequence;
1426              array[numberNonZero++]=tableauValue;
1427              double mult = multiplier[iStatus];
1428              double alpha = tableauValue * mult;
1429              double oldValue = abcDj[iSequence] * mult;
1430              double value = oldValue - tentativeTheta * alpha;
1431              if (value < dualT) {
1432                bestPossible = CoinMax(bestPossible, alpha);
1433                value = oldValue - upperTheta * alpha;
1434                if (value < dualT && alpha >= acceptablePivot) {
1435                  upperTheta = (oldValue - dualT) / alpha;
1436                }
1437                // add to list
1438                arrayCandidate[numberRemaining] = alpha;
1439                indexCandidate[numberRemaining++] = iSequence;
1440              }
1441            } else {
1442              bool keep;
1443              index[numberNonZero]=iSequence;
1444              array[numberNonZero++]=tableauValue;
1445              bestPossible = CoinMax(bestPossible, fabs(tableauValue));
1446              double oldValue = abcDj[iSequence];
1447              // If free has to be very large - should come in via dualRow
1448              //if (getInternalStatus(iSequence+addSequence)==isFree&&fabs(tableauValue)<1.0e-3)
1449              //break;
1450              // may be fake super basic
1451              if (oldValue > currentDualTolerance) {
1452                keep = true;
1453              } else if (oldValue < -currentDualTolerance) {
1454                keep = true;
1455              } else {
1456                if (fabs(tableauValue) > CoinMax(10.0 * acceptablePivot, 1.0e-5)) {
1457                  keep = true;
1458                } else {
1459                  keep = false;
1460                  badFree = CoinMax(badFree, fabs(tableauValue));
1461                }
1462              }
1463#if 0
1464              if (iSequence==freeSequenceIn)
1465                assert (keep);
1466#endif
1467              if (keep) {
1468#ifdef PAN
1469                if (model_->fakeSuperBasic(iSequence)>=0) {
1470#endif
1471                  if (iSequence==freeSequenceIn)
1472                    tableauValue=COIN_DBL_MAX;
1473                  // free - choose largest
1474                  if (fabs(tableauValue) > fabs(freeAlpha)) {
1475                    freeAlpha = tableauValue;
1476                    freeSequence = iSequence;
1477                  }
1478#ifdef PAN
1479                }
1480#endif
1481              }
1482            }
1483          }
1484        }
1485      }
1486    }
1487  }
1488#if 0
1489  if (model_->freeSequenceIn()>=first&&model_->freeSequenceIn()<last)
1490    assert (freeSequence==model_->freeSequenceIn());
1491  extern int xxInfo[6][8];
1492  xxInfo[0][iBlock]=first;
1493  xxInfo[1][iBlock]=last;
1494  xxInfo[2][iBlock]=firstIn;
1495  xxInfo[3][iBlock]=numberNonZero-firstIn;
1496  xxInfo[4][iBlock]=numberRemaining-firstIn;
1497#endif
1498  tableauRow.setNumElementsPartition(iBlock,numberNonZero-firstIn);
1499  candidateList.setNumElementsPartition(iBlock,numberRemaining-firstIn);
1500  return upperTheta;
1501}
1502// gets sorted tableau row and a possible value of theta
1503double 
1504AbcMatrix::dualColumn1Row2(double upperTheta, int & freeSequence,
1505                           const CoinIndexedVector & update,
1506                           CoinPartitionedVector & tableauRow,
1507                           CoinPartitionedVector & candidateList) const
1508{
1509  //int first=model_->maximumAbcNumberRows();
1510  assert(update.getNumElements()==2);
1511  const double *  COIN_RESTRICT pi=update.denseVector();
1512  const int *  COIN_RESTRICT piIndex = update.getIndices();
1513  int *  COIN_RESTRICT index = tableauRow.getIndices();
1514  double *  COIN_RESTRICT array = tableauRow.denseVector();
1515  const CoinBigIndex * COIN_RESTRICT rowStart = rowStart_;
1516  int numberRows=model_->numberRows();
1517  const CoinBigIndex * COIN_RESTRICT rowEnd = rowStart+numberRows*numberRowBlocks_;
1518  const double * COIN_RESTRICT element = element_;
1519  const int * COIN_RESTRICT column = column_;
1520  int iRow0 = piIndex[0];
1521  int iRow1 = piIndex[1];
1522  CoinBigIndex end0 = rowEnd[iRow0];
1523  CoinBigIndex end1 = rowEnd[iRow1];
1524  if (end0-rowStart[iRow0]>end1-rowStart[iRow1]) {
1525    int temp=iRow0;
1526    iRow0=iRow1;
1527    iRow1=temp;
1528  }
1529  CoinBigIndex start = rowStart[iRow0];
1530  CoinBigIndex end = rowEnd[iRow0];
1531  double piValue = pi[iRow0];
1532  double *  COIN_RESTRICT arrayCandidate=candidateList.denseVector();
1533  int numberNonZero;
1534  numberNonZero=tableauRow.getNumElements(0);
1535  int n=numberNonZero;
1536  for (CoinBigIndex j=start;j<end;j++) {
1537    int iSequence = column[j];
1538    double value = element[j];
1539    arrayCandidate[iSequence]= piValue*value; 
1540    index[n++]=iSequence;
1541  }
1542  start = rowStart[iRow1];
1543  end = rowEnd[iRow1];
1544  piValue = pi[iRow1];
1545  for (CoinBigIndex j=start;j<end;j++) {
1546    int iSequence = column[j];
1547    double value = element[j];
1548    value *= piValue;
1549    if (!arrayCandidate[iSequence]) {
1550      arrayCandidate[iSequence]= value; 
1551      index[n++]=iSequence;
1552    } else {   
1553      arrayCandidate[iSequence] += value; 
1554    }
1555  }
1556  // pack down
1557  const unsigned char *  COIN_RESTRICT internalStatus = model_->internalStatus();
1558  double zeroTolerance = model_->zeroTolerance();
1559  while (numberNonZero<n) {
1560    int iSequence=index[numberNonZero];
1561    double value=arrayCandidate[iSequence];
1562    arrayCandidate[iSequence]=0.0;
1563    unsigned char iStatus=internalStatus[iSequence]&7;
1564    if (fabs(value)>zeroTolerance&&iStatus<6) {
1565      index[numberNonZero]=iSequence;
1566      array[numberNonZero++]=value;
1567    } else {
1568      // kill
1569      n--;
1570      index[numberNonZero]=index[n];
1571    }
1572  }
1573  tableauRow.setNumElementsPartition(0,numberNonZero);
1574  return firstPass(model_,0, upperTheta, freeSequence,
1575                   tableauRow,
1576                   candidateList) ;
1577}
1578//static int ixxxxxx=1;
1579// gets sorted tableau row and a possible value of theta
1580double 
1581AbcMatrix::dualColumn1RowFew(int iBlock, double upperTheta, int & freeSequence,
1582                           const CoinIndexedVector & update,
1583                           CoinPartitionedVector & tableauRow,
1584                           CoinPartitionedVector & candidateList) const
1585{
1586  //int first=model_->maximumAbcNumberRows();
1587  int number = update.getNumElements();
1588  const double *  COIN_RESTRICT pi=update.denseVector();
1589  const int *  COIN_RESTRICT piIndex = update.getIndices();
1590  int *  COIN_RESTRICT index = tableauRow.getIndices();
1591  double *  COIN_RESTRICT array = tableauRow.denseVector();
1592  const CoinBigIndex * COIN_RESTRICT rowStart = rowStart_;
1593  int numberRows=model_->numberRows();
1594  const CoinBigIndex * COIN_RESTRICT rowEnd = rowStart+numberRows*numberRowBlocks_;
1595  const double * COIN_RESTRICT element = element_;
1596  const int * COIN_RESTRICT column = column_;
1597  double *  COIN_RESTRICT arrayCandidate=candidateList.denseVector();
1598  int numberNonZero;
1599  assert (iBlock>=0); 
1600  const unsigned char *  COIN_RESTRICT internalStatus = model_->internalStatus();
1601  int firstIn=blockStart_[iBlock];
1602  if (iBlock==0) {
1603    numberNonZero=0;
1604    for (int i=0;i<number;i++) {
1605      int iRow=piIndex[i];
1606      unsigned char type=internalStatus[iRow];
1607      if ((type&7)<6) {
1608        index[numberNonZero]=iRow;
1609        double piValue=pi[iRow];
1610        array[numberNonZero++]=piValue;
1611      }
1612    }
1613  } else {
1614    numberNonZero=firstIn;
1615  }
1616  int n=numberNonZero;
1617  rowStart += iBlock*numberRows;
1618  rowEnd = rowStart+numberRows;
1619  for (int i=0;i<number;i++) {
1620    int iRow=piIndex[i];
1621    double piValue = pi[iRow];
1622    CoinBigIndex end = rowEnd[iRow];
1623    for (CoinBigIndex j=rowStart[iRow];j<end;j++) {
1624      int iColumn = column[j];
1625      double value = element[j]*piValue;
1626      double oldValue=arrayCandidate[iColumn];
1627      value += oldValue;
1628      if (!oldValue) {
1629        arrayCandidate[iColumn]= value; 
1630        index[n++]=iColumn;
1631      } else if (value) {       
1632        arrayCandidate[iColumn] = value; 
1633      } else {
1634        arrayCandidate[iColumn] = COIN_INDEXED_REALLY_TINY_ELEMENT;
1635      }
1636    }
1637  }
1638  // pack down
1639  double zeroTolerance = model_->zeroTolerance();
1640  while (numberNonZero<n) {
1641    int iSequence=index[numberNonZero];
1642    double value=arrayCandidate[iSequence];
1643    arrayCandidate[iSequence]=0.0;
1644    unsigned char iStatus=internalStatus[iSequence]&7;
1645    if (fabs(value)>zeroTolerance&&iStatus<6) {
1646      index[numberNonZero]=iSequence;
1647      array[numberNonZero++]=value;
1648    } else {
1649      // kill
1650      n--;
1651      index[numberNonZero]=index[n];
1652    }
1653  }
1654  tableauRow.setNumElementsPartition(iBlock,numberNonZero-firstIn);
1655  upperTheta = firstPass(model_,iBlock, upperTheta, freeSequence,
1656                   tableauRow,
1657                   candidateList) ;
1658  return upperTheta;
1659}
1660// gets sorted tableau row and a possible value of theta
1661double 
1662AbcMatrix::dualColumn1Row1(double upperTheta, int & freeSequence,
1663                           const CoinIndexedVector & update,
1664                           CoinPartitionedVector & tableauRow,
1665                           CoinPartitionedVector & candidateList) const
1666{
1667  assert(update.getNumElements()==1);
1668  int iRow = update.getIndices()[0];
1669  double piValue = update.denseVector()[iRow];
1670  int *  COIN_RESTRICT index = tableauRow.getIndices();
1671  double *  COIN_RESTRICT array = tableauRow.denseVector();
1672  const CoinBigIndex * COIN_RESTRICT rowStart = rowStart_;
1673  int numberRows=model_->numberRows();
1674  const CoinBigIndex * COIN_RESTRICT rowEnd = rowStart+numberRows*numberRowBlocks_;
1675  CoinBigIndex start = rowStart[iRow];
1676  CoinBigIndex end = rowEnd[iRow];
1677  const double * COIN_RESTRICT element = element_;
1678  const int * COIN_RESTRICT column = column_;
1679  int numberNonZero;
1680  int numberRemaining;
1681  numberNonZero=tableauRow.getNumElements(0);
1682  numberRemaining = candidateList.getNumElements(0);
1683  double *  COIN_RESTRICT arrayCandidate=candidateList.denseVector();
1684  int *  COIN_RESTRICT indexCandidate = candidateList.getIndices();
1685  const double *  COIN_RESTRICT abcDj = model_->djRegion();
1686  const unsigned char *  COIN_RESTRICT internalStatus = model_->internalStatus();
1687  // do first pass to get possibles
1688  double bestPossible = 0.0;
1689  // We can also see if infeasible or pivoting on free
1690  double tentativeTheta = 1.0e25; // try with this much smaller as guess
1691  double acceptablePivot = model_->currentAcceptablePivot();
1692  double dualT=-model_->currentDualTolerance();
1693  const double multiplier[] = { 1.0, -1.0};
1694  double zeroTolerance = model_->zeroTolerance();
1695  freeSequence=-1;
1696  if (model_->ordinaryVariables()) {
1697    for (CoinBigIndex j=start;j<end;j++) {
1698      int iSequence = column[j];
1699      double value = element[j];
1700      double tableauValue = piValue*value;
1701      if (fabs(tableauValue)>zeroTolerance) {
1702        unsigned char iStatus=internalStatus[iSequence]&7;
1703        if (iStatus<4) {
1704          index[numberNonZero]=iSequence;
1705          array[numberNonZero++]=tableauValue;
1706          double mult = multiplier[iStatus];
1707          double alpha = tableauValue * mult;
1708          double oldValue = abcDj[iSequence] * mult;
1709          double value = oldValue - tentativeTheta * alpha;
1710          if (value < dualT) {
1711            bestPossible = CoinMax(bestPossible, alpha);
1712            value = oldValue - upperTheta * alpha;
1713            if (value < dualT && alpha >= acceptablePivot) {
1714              upperTheta = (oldValue - dualT) / alpha;
1715            }
1716            // add to list
1717            arrayCandidate[numberRemaining] = alpha;
1718            indexCandidate[numberRemaining++] = iSequence;
1719          }
1720        }
1721      }
1722    }
1723  } else {
1724    double badFree = 0.0;
1725    double freeAlpha = model_->currentAcceptablePivot();
1726    int freeSequenceIn=model_->freeSequenceIn();
1727    double currentDualTolerance = model_->currentDualTolerance();
1728    for (CoinBigIndex j=start;j<end;j++) {
1729      int iSequence = column[j];
1730      double value = element[j];
1731      double tableauValue = piValue*value;
1732      if (fabs(tableauValue)>zeroTolerance) {
1733        unsigned char iStatus=internalStatus[iSequence]&7;
1734        if (iStatus<6) {
1735          if ((iStatus&4)==0) {
1736            index[numberNonZero]=iSequence;
1737            array[numberNonZero++]=tableauValue;
1738            double mult = multiplier[iStatus];
1739            double alpha = tableauValue * mult;
1740            double oldValue = abcDj[iSequence] * mult;
1741            double value = oldValue - tentativeTheta * alpha;
1742            if (value < dualT) {
1743              bestPossible = CoinMax(bestPossible, alpha);
1744              value = oldValue - upperTheta * alpha;
1745              if (value < dualT && alpha >= acceptablePivot) {
1746                upperTheta = (oldValue - dualT) / alpha;
1747              }
1748              // add to list
1749              arrayCandidate[numberRemaining] = alpha;
1750              indexCandidate[numberRemaining++] = iSequence;
1751            }
1752          } else {
1753            bool keep;
1754            index[numberNonZero]=iSequence;
1755            array[numberNonZero++]=tableauValue;
1756            bestPossible = CoinMax(bestPossible, fabs(tableauValue));
1757            double oldValue = abcDj[iSequence];
1758            // If free has to be very large - should come in via dualRow
1759            //if (getInternalStatus(iSequence+addSequence)==isFree&&fabs(tableauValue)<1.0e-3)
1760            //break;
1761            if (oldValue > currentDualTolerance) {
1762              keep = true;
1763            } else if (oldValue < -currentDualTolerance) {
1764              keep = true;
1765            } else {
1766              if (fabs(tableauValue) > CoinMax(10.0 * acceptablePivot, 1.0e-5)) {
1767                keep = true;
1768              } else {
1769                keep = false;
1770                badFree = CoinMax(badFree, fabs(tableauValue));
1771              }
1772            }
1773            if (keep) {
1774#ifdef PAN
1775              if (model_->fakeSuperBasic(iSequence)>=0) {
1776#endif
1777                if (iSequence==freeSequenceIn)
1778                  tableauValue=COIN_DBL_MAX;
1779                // free - choose largest
1780                if (fabs(tableauValue) > fabs(freeAlpha)) {
1781                  freeAlpha = tableauValue;
1782                  freeSequence = iSequence;
1783                }
1784#ifdef PAN
1785              }
1786#endif
1787            }
1788          }
1789        }
1790      }
1791    }
1792  }
1793  tableauRow.setNumElementsPartition(0,numberNonZero);
1794  candidateList.setNumElementsPartition(0,numberRemaining);
1795  return upperTheta;
1796}
1797//#define PARALLEL2
1798#ifdef PARALLEL2
1799#undef cilk_for
1800#undef cilk_spawn
1801#undef cilk_sync
1802#include <cilk/cilk.h>
1803#include <cilk/cilk_api.h>
1804#endif
1805#if 0
1806static void compact(int numberBlocks,CoinIndexedVector * vector,const int * starts,const int * lengths)
1807{
1808  int numberNonZeroIn=vector->getNumElements();
1809  int * index = vector->getIndices();
1810  double * array = vector->denseVector();
1811  CoinAbcCompact(numberBlocks,numberNonZeroIn,
1812                 array,starts,lengths);
1813  int numberNonZero = CoinAbcCompact(numberBlocks,numberNonZeroIn,
1814                                     index,starts,lengths);
1815  vector->setNumElements(numberNonZero);
1816}
1817static void compactBoth(int numberBlocks,CoinIndexedVector * vector1,CoinIndexedVector * vector2,
1818                        const int * starts,const int * lengths1, const int * lengths2)
1819{
1820  cilk_spawn compact(numberBlocks,vector1,starts,lengths1);
1821  compact(numberBlocks,vector2,starts,lengths2);
1822  cilk_sync;
1823}
1824#endif
1825void
1826AbcMatrix::rebalance() const
1827{
1828  int maximumRows = model_->maximumAbcNumberRows();
1829  int numberTotal = model_->numberTotal();
1830  /* rebalance
1831     For non-vector version
1832     each basic etc column is 1
1833     each real column is 5+2*nel
1834     each basic slack is 0
1835     each real slack is 3
1836   */
1837#if ABC_PARALLEL
1838  int howOften=CoinMax(model_->factorization()->maximumPivots(),200);
1839  if ((model_->numberIterations()%howOften)==0||!startColumnBlock_[1]) {
1840    int numberCpus=model_->numberCpus();
1841    if (numberCpus>1) {
1842      const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts()-maximumRows;
1843      const unsigned char *  COIN_RESTRICT internalStatus = model_->internalStatus();
1844      int numberRows=model_->numberRows();
1845      int total=0;
1846      for (int iSequence=0;iSequence<numberRows;iSequence++) {
1847        unsigned char iStatus=internalStatus[iSequence]&7;
1848        if (iStatus<4) 
1849          total+=3;
1850      }
1851      int totalSlacks=total;
1852      CoinBigIndex end = columnStart[maximumRows];
1853      for (int iSequence=maximumRows;iSequence<numberTotal;iSequence++) {
1854        CoinBigIndex start=end;
1855        end = columnStart[iSequence+1];
1856        unsigned char iStatus=internalStatus[iSequence]&7;
1857        if (iStatus<4) 
1858          total+=5+2*(end-start);
1859        else
1860          total+=1;
1861      }
1862      int chunk = (total+numberCpus-1)/numberCpus;
1863      total=totalSlacks;
1864      int iCpu=0;
1865      startColumnBlock_[0]=0;
1866      end = columnStart[maximumRows];
1867      for (int iSequence=maximumRows;iSequence<numberTotal;iSequence++) {
1868        CoinBigIndex start=end;
1869        end = columnStart[iSequence+1];
1870        unsigned char iStatus=internalStatus[iSequence]&7;
1871        if (iStatus<4) 
1872          total+=5+2*(end-start);
1873        else
1874          total+=1;
1875        if (total>chunk) {
1876          iCpu++;
1877          total=0;
1878          startColumnBlock_[iCpu]=iSequence;
1879        }
1880      }
1881      assert(iCpu<numberCpus);
1882      iCpu++;
1883      for (;iCpu<=numberCpus;iCpu++)
1884        startColumnBlock_[iCpu]=numberTotal;
1885      numberColumnBlocks_=numberCpus;
1886    } else {
1887      numberColumnBlocks_=1;
1888      startColumnBlock_[0]=0;
1889      startColumnBlock_[1]=numberTotal;
1890    }
1891  }
1892#else
1893  startColumnBlock_[0]=0;
1894  startColumnBlock_[1]=numberTotal;
1895#endif
1896}
1897double 
1898AbcMatrix::dualColumn1(const CoinIndexedVector & update,
1899                                  CoinPartitionedVector & tableauRow,CoinPartitionedVector & candidateList) const
1900{
1901  int numberTotal = model_->numberTotal();
1902  // rebalance
1903  rebalance();
1904  tableauRow.setPackedMode(true);
1905  candidateList.setPackedMode(true);
1906  int number=update.getNumElements();
1907  double upperTheta;
1908  if (rowStart_&&number<3) {
1909#if ABC_INSTRUMENT>1
1910    {
1911      int n=update.getNumElements();
1912      abcPricing[n<19 ? n : 19]++;
1913    }
1914#endif
1915    // always do serially
1916    // do slacks first
1917    int starts[2];
1918    starts[0]=0;
1919    starts[1]=numberTotal;
1920    tableauRow.setPartitions(1,starts);
1921    candidateList.setPartitions(1,starts);
1922    upperTheta = static_cast<AbcSimplexDual *>(model_)->dualColumn1A();
1923    //assert (upperTheta>-100*model_->dualTolerance()||model_->sequenceIn()>=0||model_->lastFirstFree()>=0);
1924    int freeSequence=-1;
1925    // worth using row copy
1926    assert (number);
1927    if (number==2) {
1928      upperTheta=dualColumn1Row2(upperTheta,freeSequence,update,tableauRow,candidateList);
1929    } else {
1930      upperTheta=dualColumn1Row1(upperTheta,freeSequence,update,tableauRow,candidateList);
1931    }
1932    if (freeSequence>=0) {
1933      int numberNonZero=tableauRow.getNumElements(0);
1934      const int *  COIN_RESTRICT index = tableauRow.getIndices();
1935      const double *  COIN_RESTRICT array = tableauRow.denseVector();
1936      // search for free coming in
1937      double freeAlpha=0.0;
1938      int bestSequence=model_->sequenceIn();
1939      if (bestSequence>=0)
1940        freeAlpha=model_->alpha();
1941      index = tableauRow.getIndices();
1942      array = tableauRow.denseVector();
1943      // free variable - search
1944      for (int k=0;k<numberNonZero;k++) {
1945        if (freeSequence==index[k]) {
1946          if (fabs(freeAlpha)<fabs(array[k])) {
1947            freeAlpha=array[k];
1948          }
1949          break;
1950        }
1951      }
1952      if (model_->sequenceIn()<0||fabs(freeAlpha)>fabs(model_->alpha())) {
1953        double oldValue = model_->djRegion()[freeSequence];
1954        model_->setSequenceIn(freeSequence);
1955        model_->setAlpha(freeAlpha);
1956        model_->setTheta(oldValue / freeAlpha);
1957      }
1958    }
1959  } else {
1960    // three or more
1961    // need to do better job on dividing up (but wait until vector or by row)
1962    upperTheta = parallelDual4(static_cast<AbcSimplexDual *>(model_));
1963  }
1964  //tableauRow.compact();
1965  //candidateList.compact();
1966#if 0 //ndef NDEBUG
1967  model_->checkArrays();
1968#endif
1969  candidateList.computeNumberElements();
1970  int numberRemaining=candidateList.getNumElements();
1971  if (!numberRemaining&&model_->sequenceIn()<0) {
1972    return COIN_DBL_MAX; // Looks infeasible
1973  } else {
1974    return upperTheta;
1975  }
1976}
1977#define _mm256_broadcast_sd(x) static_cast<__m256d> (__builtin_ia32_vbroadcastsd256 (x))
1978#define _mm256_load_pd(x) *(__m256d *)(x)
1979#define _mm256_store_pd (s, x)  *((__m256d *)s)=x
1980void
1981AbcMatrix::dualColumn1Part(int iBlock, int & sequenceIn, double & upperThetaResult,
1982                           const CoinIndexedVector & update,
1983                           CoinPartitionedVector & tableauRow,CoinPartitionedVector & candidateList) const
1984{
1985  double upperTheta=upperThetaResult;
1986#if 0
1987  double time0=CoinCpuTime();
1988#endif
1989  int maximumRows = model_->maximumAbcNumberRows();
1990  int firstIn=startColumnBlock_[iBlock];
1991  int last = startColumnBlock_[iBlock+1];
1992  int numberNonZero;
1993  int numberRemaining;
1994  int first;
1995  if (firstIn==0) {
1996    upperTheta=static_cast<AbcSimplexDual *>(model_)->dualColumn1A();
1997    numberNonZero=tableauRow.getNumElements(0);
1998    numberRemaining = candidateList.getNumElements(0);
1999    first=maximumRows;
2000  } else {
2001    first=firstIn;
2002    numberNonZero=first;
2003    numberRemaining=first;
2004  }
2005  sequenceIn=-1;
2006  // get matrix data pointers
2007  const int * COIN_RESTRICT row = matrix_->getIndices();
2008  const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts()-maximumRows;
2009  const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
2010  const double *  COIN_RESTRICT pi=update.denseVector();
2011  //const int *  COIN_RESTRICT piIndex = update.getIndices();
2012  int *  COIN_RESTRICT index = tableauRow.getIndices();
2013  double *  COIN_RESTRICT array = tableauRow.denseVector();
2014  // pivot elements
2015  double *  COIN_RESTRICT arrayCandidate=candidateList.denseVector();
2016  // indices
2017  int *  COIN_RESTRICT indexCandidate = candidateList.getIndices();
2018  const double *  COIN_RESTRICT abcDj = model_->djRegion();
2019  const unsigned char *  COIN_RESTRICT internalStatus = model_->internalStatus();
2020  // do first pass to get possibles
2021  double bestPossible = 0.0;
2022  // We can also see if infeasible or pivoting on free
2023  double tentativeTheta = 1.0e25; // try with this much smaller as guess
2024  double acceptablePivot = model_->currentAcceptablePivot();
2025  double dualT=-model_->currentDualTolerance();
2026  const double multiplier[] = { 1.0, -1.0};
2027  double zeroTolerance = model_->zeroTolerance();
2028  if (model_->ordinaryVariables()) {
2029#ifdef COUNT_COPY
2030    if (first>maximumRows||last<model_->numberTotal()||false) {
2031#endif
2032#if 1
2033    for (int iSequence = first; iSequence < last; iSequence++) {
2034      unsigned char iStatus=internalStatus[iSequence]&7;
2035      if (iStatus<4) {
2036        CoinBigIndex start = columnStart[iSequence];
2037        CoinBigIndex next = columnStart[iSequence+1];
2038        double tableauValue = 0.0;
2039        for (CoinBigIndex j = start; j < next; j++) {
2040          int jRow = row[j];
2041          tableauValue += pi[jRow] * elementByColumn[j];
2042        }
2043        if (fabs(tableauValue)>zeroTolerance) {
2044#else
2045    cilk_for (int iSequence = first; iSequence < last; iSequence++) {
2046      unsigned char iStatus=internalStatus[iSequence]&7;
2047      if (iStatus<4) {
2048        CoinBigIndex start = columnStart[iSequence];
2049        CoinBigIndex next = columnStart[iSequence+1];
2050        double tableauValue = 0.0;
2051        for (CoinBigIndex j = start; j < next; j++) {
2052          int jRow = row[j];
2053          tableauValue += pi[jRow] * elementByColumn[j];
2054        }
2055        array[iSequence]=tableauValue;
2056      }
2057    }
2058    //printf("time %.6g\n",CoinCpuTime()-time0);
2059    for (int iSequence = first; iSequence < last; iSequence++) {
2060      double tableauValue=array[iSequence];
2061      if (tableauValue) {
2062        array[iSequence]=0.0;
2063        if (fabs(tableauValue)>zeroTolerance) {
2064          unsigned char iStatus=internalStatus[iSequence]&7;
2065#endif
2066          index[numberNonZero]=iSequence;
2067          array[numberNonZero++]=tableauValue;
2068          double mult = multiplier[iStatus];
2069          double alpha = tableauValue * mult;
2070          double oldValue = abcDj[iSequence] * mult;
2071          double value = oldValue - tentativeTheta * alpha;
2072          if (value < dualT) {
2073            bestPossible = CoinMax(bestPossible, alpha);
2074            value = oldValue - upperTheta * alpha;
2075            if (value < dualT && alpha >= acceptablePivot) {
2076              upperTheta = (oldValue - dualT) / alpha;
2077            }
2078            // add to list
2079            arrayCandidate[numberRemaining] = alpha;
2080            indexCandidate[numberRemaining++] = iSequence;
2081          }
2082        }
2083      }
2084    }
2085#ifdef COUNT_COPY
2086    } else {
2087      const double * COIN_RESTRICT element = countElement_;
2088      const int * COIN_RESTRICT row = countRow_;
2089      double temp[4] __attribute__ ((aligned (32)));
2090      memset(temp,0,sizeof(temp));
2091      for (int iCount=smallestCount_;iCount<largestCount_;iCount++) {
2092        int iStart=countFirst_[iCount];
2093        int number=countFirst_[iCount+1]-iStart;
2094        if (!number)
2095          continue;
2096        const int * COIN_RESTRICT blockRow = row+countStart_[iCount];
2097        const double * COIN_RESTRICT blockElement = element+countStart_[iCount];
2098        const int * COIN_RESTRICT sequence = countRealColumn_+iStart;
2099        // really need to sort and swap to avoid tests
2100        int numberBlocks=number>>2;
2101        number &= 3;
2102        for (int iBlock=0;iBlock<numberBlocks;iBlock++) {
2103          double tableauValue0=0.0;
2104          double tableauValue1=0.0;
2105          double tableauValue2=0.0;
2106          double tableauValue3=0.0;
2107          for (int j=0;j<iCount;j++) {
2108            tableauValue0 += pi[blockRow[0]]*blockElement[0];
2109            tableauValue1 += pi[blockRow[1]]*blockElement[1];
2110            tableauValue2 += pi[blockRow[2]]*blockElement[2];
2111            tableauValue3 += pi[blockRow[3]]*blockElement[3];
2112            blockRow+=4;
2113            blockElement+=4;
2114          }
2115          if (fabs(tableauValue0)>zeroTolerance) {
2116            int iSequence=sequence[0];
2117            unsigned char iStatus=internalStatus[iSequence]&7;
2118            if (iStatus<4) {
2119              index[numberNonZero]=iSequence;
2120              array[numberNonZero++]=tableauValue0;
2121              double mult = multiplier[iStatus];
2122              double alpha = tableauValue0 * mult;
2123              double oldValue = abcDj[iSequence] * mult;
2124              double value = oldValue - tentativeTheta * alpha;
2125              if (value < dualT) {
2126                bestPossible = CoinMax(bestPossible, alpha);
2127                value = oldValue - upperTheta * alpha;
2128                if (value < dualT && alpha >= acceptablePivot) {
2129                  upperTheta = (oldValue - dualT) / alpha;
2130                }
2131                // add to list
2132                arrayCandidate[numberRemaining] = alpha;
2133                indexCandidate[numberRemaining++] = iSequence;
2134              }
2135            }
2136          }
2137          if (fabs(tableauValue1)>zeroTolerance) {
2138            int iSequence=sequence[1];
2139            unsigned char iStatus=internalStatus[iSequence]&7;
2140            if (iStatus<4) {
2141              index[numberNonZero]=iSequence;
2142              array[numberNonZero++]=tableauValue1;
2143              double mult = multiplier[iStatus];
2144              double alpha = tableauValue1 * mult;
2145              double oldValue = abcDj[iSequence] * mult;
2146              double value = oldValue - tentativeTheta * alpha;
2147              if (value < dualT) {
2148                bestPossible = CoinMax(bestPossible, alpha);
2149                value = oldValue - upperTheta * alpha;
2150                if (value < dualT && alpha >= acceptablePivot) {
2151                  upperTheta = (oldValue - dualT) / alpha;
2152                }
2153                // add to list
2154                arrayCandidate[numberRemaining] = alpha;
2155                indexCandidate[numberRemaining++] = iSequence;
2156              }
2157            }
2158          }
2159          if (fabs(tableauValue2)>zeroTolerance) {
2160            int iSequence=sequence[2];
2161            unsigned char iStatus=internalStatus[iSequence]&7;
2162            if (iStatus<4) {
2163              index[numberNonZero]=iSequence;
2164              array[numberNonZero++]=tableauValue2;
2165              double mult = multiplier[iStatus];
2166              double alpha = tableauValue2 * mult;
2167              double oldValue = abcDj[iSequence] * mult;
2168              double value = oldValue - tentativeTheta * alpha;
2169              if (value < dualT) {
2170                bestPossible = CoinMax(bestPossible, alpha);
2171                value = oldValue - upperTheta * alpha;
2172                if (value < dualT && alpha >= acceptablePivot) {
2173                  upperTheta = (oldValue - dualT) / alpha;
2174                }
2175                // add to list
2176                arrayCandidate[numberRemaining] = alpha;
2177                indexCandidate[numberRemaining++] = iSequence;
2178              }
2179            }
2180          }
2181          if (fabs(tableauValue3)>zeroTolerance) {
2182            int iSequence=sequence[3];
2183            unsigned char iStatus=internalStatus[iSequence]&7;
2184            if (iStatus<4) {
2185              index[numberNonZero]=iSequence;
2186              array[numberNonZero++]=tableauValue3;
2187              double mult = multiplier[iStatus];
2188              double alpha = tableauValue3 * mult;
2189              double oldValue = abcDj[iSequence] * mult;
2190              double value = oldValue - tentativeTheta * alpha;
2191              if (value < dualT) {
2192                bestPossible = CoinMax(bestPossible, alpha);
2193                value = oldValue - upperTheta * alpha;
2194                if (value < dualT && alpha >= acceptablePivot) {
2195                  upperTheta = (oldValue - dualT) / alpha;
2196                }
2197                // add to list
2198                arrayCandidate[numberRemaining] = alpha;
2199                indexCandidate[numberRemaining++] = iSequence;
2200              }
2201            }
2202          }
2203          sequence+=4;
2204        }
2205        for (int i=0;i<number;i++) {
2206          int iSequence=sequence[i];
2207          unsigned char iStatus=internalStatus[iSequence]&7;
2208          if (iStatus<4) {
2209            double tableauValue=0.0;
2210            for (int j=0;j<iCount;j++) {
2211              int iRow=blockRow[4*j];
2212              tableauValue += pi[iRow]*blockElement[4*j];
2213            }
2214            if (fabs(tableauValue)>zeroTolerance) {
2215              index[numberNonZero]=iSequence;
2216              array[numberNonZero++]=tableauValue;
2217              double mult = multiplier[iStatus];
2218              double alpha = tableauValue * mult;
2219              double oldValue = abcDj[iSequence] * mult;
2220              double value = oldValue - tentativeTheta * alpha;
2221              if (value < dualT) {
2222                bestPossible = CoinMax(bestPossible, alpha);
2223                value = oldValue - upperTheta * alpha;
2224                if (value < dualT && alpha >= acceptablePivot) {
2225                  upperTheta = (oldValue - dualT) / alpha;
2226                }
2227                // add to list
2228                arrayCandidate[numberRemaining] = alpha;
2229                indexCandidate[numberRemaining++] = iSequence;
2230              }
2231            }
2232          }
2233          blockRow++;
2234          blockElement++;
2235        }
2236      }
2237      int numberColumns=model_->numberColumns();
2238      // large ones
2239      const CoinBigIndex * COIN_RESTRICT largeStart = countStartLarge_-countFirst_[MAX_COUNT];
2240      for (int i=countFirst_[MAX_COUNT];i<numberColumns;i++) {
2241        int iSequence=countRealColumn_[i];
2242        unsigned char iStatus=internalStatus[iSequence]&7;
2243        if (iStatus<4) {
2244          CoinBigIndex start = largeStart[i];
2245          CoinBigIndex next = largeStart[i+1];
2246          double tableauValue = 0.0;
2247          for (CoinBigIndex j = start; j < next; j++) {
2248            int jRow = row[j];
2249            tableauValue += pi[jRow] * element[j];
2250          }
2251          if (fabs(tableauValue)>zeroTolerance) {
2252            index[numberNonZero]=iSequence;
2253            array[numberNonZero++]=tableauValue;
2254            double mult = multiplier[iStatus];
2255            double alpha = tableauValue * mult;
2256            double oldValue = abcDj[iSequence] * mult;
2257            double value = oldValue - tentativeTheta * alpha;
2258            if (value < dualT) {
2259              bestPossible = CoinMax(bestPossible, alpha);
2260              value = oldValue - upperTheta * alpha;
2261              if (value < dualT && alpha >= acceptablePivot) {
2262                upperTheta = (oldValue - dualT) / alpha;
2263              }
2264              // add to list
2265              arrayCandidate[numberRemaining] = alpha;
2266              indexCandidate[numberRemaining++] = iSequence;
2267            }
2268          }
2269        }
2270      }
2271    }
2272#endif
2273  } else {
2274    double badFree = 0.0;
2275    double freePivot = model_->currentAcceptablePivot();
2276    int freeSequenceIn=model_->freeSequenceIn();
2277    double currentDualTolerance = model_->currentDualTolerance();
2278    for (int iSequence = first; iSequence < last; iSequence++) {
2279      unsigned char iStatus=internalStatus[iSequence]&7;
2280      if (iStatus<6) {
2281        CoinBigIndex start = columnStart[iSequence];
2282        CoinBigIndex next = columnStart[iSequence+1];
2283        double tableauValue = 0.0;
2284        for (CoinBigIndex j = start; j < next; j++) {
2285          int jRow = row[j];
2286          tableauValue += pi[jRow] * elementByColumn[j];
2287        }
2288        if (fabs(tableauValue)>zeroTolerance) {
2289          if ((iStatus&4)==0) {
2290            index[numberNonZero]=iSequence;
2291            array[numberNonZero++]=tableauValue;
2292            double mult = multiplier[iStatus];
2293            double alpha = tableauValue * mult;
2294            double oldValue = abcDj[iSequence] * mult;
2295            double value = oldValue - tentativeTheta * alpha;
2296            if (value < dualT) {
2297              bestPossible = CoinMax(bestPossible, alpha);
2298              value = oldValue - upperTheta * alpha;
2299              if (value < dualT && alpha >= acceptablePivot) {
2300                upperTheta = (oldValue - dualT) / alpha;
2301              }
2302              // add to list
2303              arrayCandidate[numberRemaining] = alpha;
2304              indexCandidate[numberRemaining++] = iSequence;
2305            }
2306          } else {
2307            bool keep;
2308            index[numberNonZero]=iSequence;
2309            array[numberNonZero++]=tableauValue;
2310            bestPossible = CoinMax(bestPossible, fabs(tableauValue));
2311            double oldValue = abcDj[iSequence];
2312            // If free has to be very large - should come in via dualRow
2313            //if (getInternalStatus(iSequence+addSequence)==isFree&&fabs(tableauValue)<1.0e-3)
2314            //break;
2315            if (oldValue > currentDualTolerance) {
2316              keep = true;
2317            } else if (oldValue < -currentDualTolerance) {
2318              keep = true;
2319            } else {
2320              if (fabs(tableauValue) > CoinMax(10.0 * acceptablePivot, 1.0e-5)) {
2321                keep = true;
2322              } else {
2323                keep = false;
2324                badFree = CoinMax(badFree, fabs(tableauValue));
2325              }
2326            }
2327            if (keep) {
2328#ifdef PAN
2329              if (model_->fakeSuperBasic(iSequence)>=0) {
2330#endif
2331                if (iSequence==freeSequenceIn)
2332                  tableauValue=COIN_DBL_MAX;
2333                // free - choose largest
2334                if (fabs(tableauValue) > fabs(freePivot)) {
2335                  freePivot = tableauValue;
2336                  sequenceIn = iSequence;
2337                }
2338#ifdef PAN
2339              }
2340#endif
2341            }
2342          }
2343        }
2344      }
2345    }
2346  } 
2347  // adjust
2348  numberNonZero -= firstIn;
2349  numberRemaining -= firstIn;
2350  tableauRow.setNumElementsPartition(iBlock,numberNonZero);
2351  candidateList.setNumElementsPartition(iBlock,numberRemaining);
2352  if (!numberRemaining&&model_->sequenceIn()<0) {
2353 
2354    upperThetaResult=COIN_DBL_MAX; // Looks infeasible
2355  } else {
2356    upperThetaResult=upperTheta;
2357  }
2358}
2359// gets tableau row - returns number of slacks in block
2360int 
2361AbcMatrix::primalColumnRow(int iBlock,bool doByRow,const CoinIndexedVector & updates,
2362                           CoinPartitionedVector & tableauRow) const
2363{
2364  int maximumRows = model_->maximumAbcNumberRows();
2365  int first=tableauRow.startPartition(iBlock);
2366  int last=tableauRow.startPartition(iBlock+1);
2367  if (doByRow) { 
2368    assert(first==blockStart_[iBlock]);
2369    assert(last==blockStart_[iBlock+1]);
2370  } else {
2371    assert(first==startColumnBlock_[iBlock]);
2372    assert(last==startColumnBlock_[iBlock+1]);
2373  }
2374  int numberSlacks=updates.getNumElements();
2375  int numberNonZero=0;
2376  const double *  COIN_RESTRICT pi=updates.denseVector();
2377  const int *  COIN_RESTRICT piIndex = updates.getIndices();
2378  int *  COIN_RESTRICT index = tableauRow.getIndices()+first;
2379  double *  COIN_RESTRICT array = tableauRow.denseVector()+first;
2380  const unsigned char *  COIN_RESTRICT internalStatus = model_->internalStatus();
2381  if (!iBlock) {
2382    numberNonZero=numberSlacks;
2383    for (int i=0;i<numberSlacks;i++) {
2384      int iRow=piIndex[i];
2385      index[i]=iRow;
2386      array[i]=pi[iRow]; // ? what if small or basic
2387    }
2388    first=maximumRows;
2389  }
2390  double zeroTolerance = model_->zeroTolerance();
2391  if (doByRow) {
2392    int numberRows=model_->numberRows();
2393    const CoinBigIndex * COIN_RESTRICT rowStart = rowStart_+iBlock*numberRows;
2394    const CoinBigIndex * COIN_RESTRICT rowEnd = rowStart+numberRows;
2395    const double * COIN_RESTRICT element = element_;
2396    const int * COIN_RESTRICT column = column_;
2397    if (numberSlacks>1) {
2398      double *  COIN_RESTRICT arrayTemp = tableauRow.denseVector();
2399      for (int i=0;i<numberSlacks;i++) {
2400        int iRow=piIndex[i];
2401        double piValue = pi[iRow];
2402        CoinBigIndex end = rowEnd[iRow];
2403        for (CoinBigIndex j=rowStart[iRow];j<end;j++) {
2404          int iColumn = column[j];
2405          arrayTemp[iColumn] += element[j]*piValue;
2406        }
2407      }
2408      for (int iSequence=first;iSequence<last;iSequence++) {
2409        double tableauValue = arrayTemp[iSequence];
2410        if (tableauValue) {
2411          arrayTemp[iSequence]=0.0;
2412          if (fabs(tableauValue)>zeroTolerance) {
2413            unsigned char iStatus=internalStatus[iSequence]&7;
2414            if (iStatus<6) {
2415              index[numberNonZero]=iSequence;
2416              array[numberNonZero++]=tableauValue;
2417            }
2418          }
2419        }
2420      }
2421    } else {
2422      int iRow=piIndex[0];
2423      double piValue = pi[iRow];
2424      CoinBigIndex end = rowEnd[iRow];
2425      for (CoinBigIndex j=rowStart[iRow];j<end;j++) {
2426        int iSequence = column[j];
2427        double tableauValue = element[j]*piValue;
2428        if (fabs(tableauValue)>zeroTolerance) {
2429          unsigned char iStatus=internalStatus[iSequence]&7;
2430          if (iStatus<6) {
2431            index[numberNonZero]=iSequence;
2432            array[numberNonZero++]=tableauValue;
2433          }
2434        }
2435      }
2436    }
2437  } else {
2438    // get matrix data pointers
2439    const int * COIN_RESTRICT row = matrix_->getIndices();
2440    const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts()-maximumRows;
2441    const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
2442    const double *  COIN_RESTRICT pi=updates.denseVector();
2443    for (int iSequence = first; iSequence < last; iSequence++) {
2444      unsigned char iStatus=internalStatus[iSequence]&7;
2445      if (iStatus<6) {
2446        CoinBigIndex start = columnStart[iSequence];
2447        CoinBigIndex next = columnStart[iSequence+1];
2448        double tableauValue = 0.0;
2449        for (CoinBigIndex j = start; j < next; j++) {
2450          int jRow = row[j];
2451          tableauValue += pi[jRow] * elementByColumn[j];
2452        }
2453        if (fabs(tableauValue)>zeroTolerance) {
2454          index[numberNonZero]=iSequence;
2455          array[numberNonZero++]=tableauValue;
2456        }
2457      }
2458    }
2459  }
2460  tableauRow.setNumElementsPartition(iBlock,numberNonZero);
2461  return numberSlacks;
2462}
2463// Get sequenceIn when Dantzig
2464 int 
2465AbcMatrix::pivotColumnDantzig(const CoinIndexedVector & updates,
2466                              CoinPartitionedVector & spare) const
2467 {
2468   int maximumRows = model_->maximumAbcNumberRows();
2469   // rebalance
2470   rebalance();
2471   spare.setPackedMode(true);
2472   bool useRowCopy=false;
2473   if (rowStart_) {
2474     // see if worth using row copy
2475     int number=updates.getNumElements();
2476     assert (number);
2477     useRowCopy = (number<<2)<maximumRows;
2478   }
2479   const int * starts;
2480   if (useRowCopy)
2481     starts=blockStart();
2482   else
2483     starts=startColumnBlock();
2484#if ABC_PARALLEL
2485#define NUMBER_BLOCKS NUMBER_ROW_BLOCKS
2486  int numberBlocks=CoinMin(NUMBER_BLOCKS,model_->numberCpus());
2487#else
2488#define NUMBER_BLOCKS 1
2489  int numberBlocks=NUMBER_BLOCKS;
2490#endif
2491#if ABC_PARALLEL
2492  if (model_->parallelMode()==0)
2493    numberBlocks=1;
2494#endif
2495   spare.setPartitions(numberBlocks,starts);
2496   int which[NUMBER_BLOCKS];
2497   double best[NUMBER_BLOCKS];
2498   for (int i=0;i<numberBlocks-1;i++) 
2499     which[i]=cilk_spawn pivotColumnDantzig(i,useRowCopy,updates,spare,best[i]); 
2500   which[numberBlocks-1]=pivotColumnDantzig(numberBlocks-1,useRowCopy,updates,
2501                                                       spare,best[numberBlocks-1]);
2502   cilk_sync;
2503   int bestSequence=-1;
2504   double bestValue=model_->dualTolerance();
2505   for (int i=0;i<numberBlocks;i++) {
2506     if (best[i]>bestValue) {
2507       bestValue=best[i];
2508       bestSequence=which[i];
2509     }
2510   }
2511   return bestSequence;
2512}
2513 // Get sequenceIn when Dantzig (One block)
2514 int 
2515   AbcMatrix::pivotColumnDantzig(int iBlock, bool doByRow,const CoinIndexedVector & updates,
2516                                 CoinPartitionedVector & spare,
2517                              double & bestValue) const
2518 {
2519   // could rewrite to do more directly
2520   int numberSlacks=primalColumnRow(iBlock,doByRow,updates,spare);
2521   double * COIN_RESTRICT reducedCost = model_->djRegion();
2522   int first=spare.startPartition(iBlock);
2523   int last=spare.startPartition(iBlock+1);
2524   int *  COIN_RESTRICT index = spare.getIndices()+first;
2525   double *  COIN_RESTRICT array = spare.denseVector()+first;
2526   int numberNonZero=spare.getNumElements(iBlock);
2527   int bestSequence=-1;
2528   double bestDj=model_->dualTolerance();
2529   double bestFreeDj = model_->dualTolerance();
2530   int bestFreeSequence = -1;
2531   // redo LOTS so signs for slacks and columns same
2532   if (!first) {
2533     first = model_->maximumAbcNumberRows();
2534     for (int i=0;i<numberSlacks;i++) {
2535       int iSequence = index[i];
2536       double value = reducedCost[iSequence];
2537       value += array[i];
2538       array[i] = 0.0;
2539       reducedCost[iSequence] = value;
2540     }
2541#ifndef CLP_PRIMAL_SLACK_MULTIPLIER
2542#define CLP_PRIMAL_SLACK_MULTIPLIER 1.0
2543#endif
2544     int numberRows=model_->numberRows();
2545     for (int iSequence=0 ; iSequence < numberRows; iSequence++) {
2546       // check flagged variable
2547       if (!model_->flagged(iSequence)) {
2548         double value = reducedCost[iSequence] * CLP_PRIMAL_SLACK_MULTIPLIER;
2549         AbcSimplex::Status status = model_->getInternalStatus(iSequence);
2550         
2551         switch(status) {
2552           
2553         case AbcSimplex::basic:
2554         case AbcSimplex::isFixed:
2555           break;
2556         case AbcSimplex::isFree:
2557         case AbcSimplex::superBasic:
2558           if (fabs(value) > bestFreeDj) {
2559             bestFreeDj = fabs(value);
2560             bestFreeSequence = iSequence;
2561           }
2562           break;
2563         case AbcSimplex::atUpperBound:
2564           if (value > bestDj) {
2565             bestDj = value;
2566             bestSequence = iSequence;
2567           }
2568           break;
2569         case AbcSimplex::atLowerBound:
2570           if (value < -bestDj) {
2571             bestDj = -value;
2572             bestSequence = iSequence;
2573           }
2574         }
2575       }
2576     }
2577   }
2578   for (int i=numberSlacks;i<numberNonZero;i++) {
2579     int iSequence = index[i];
2580     double value = reducedCost[iSequence];
2581     value += array[i];
2582     array[i] = 0.0;
2583     reducedCost[iSequence] = value;
2584   }
2585   for (int iSequence=first ; iSequence < last; iSequence++) {
2586     // check flagged variable
2587     if (!model_->flagged(iSequence)) {
2588       double value = reducedCost[iSequence];
2589       AbcSimplex::Status status = model_->getInternalStatus(iSequence);
2590       
2591       switch(status) {
2592         
2593       case AbcSimplex::basic:
2594       case AbcSimplex::isFixed:
2595         break;
2596       case AbcSimplex::isFree:
2597       case AbcSimplex::superBasic:
2598         if (fabs(value) > bestFreeDj) {
2599           bestFreeDj = fabs(value);
2600           bestFreeSequence = iSequence;
2601         }
2602         break;
2603       case AbcSimplex::atUpperBound:
2604         if (value > bestDj) {
2605           bestDj = value;
2606           bestSequence = iSequence;
2607         }
2608         break;
2609       case AbcSimplex::atLowerBound:
2610         if (value < -bestDj) {
2611           bestDj = -value;
2612           bestSequence = iSequence;
2613         }
2614       }
2615     }
2616   }
2617   spare.setNumElementsPartition(iBlock,0);
2618   // bias towards free
2619   if (bestFreeSequence >= 0 && bestFreeDj > 0.1 * bestDj) {
2620     bestSequence = bestFreeSequence;
2621     bestDj=10.0*bestFreeDj;
2622   }
2623   bestValue=bestDj;
2624   return bestSequence;
2625 } 
2626// gets tableau row and dj row - returns number of slacks in block
2627 int 
2628   AbcMatrix::primalColumnRowAndDjs(int iBlock,const CoinIndexedVector & updateTableau,
2629                                    const CoinIndexedVector & updateDjs,
2630                                    CoinPartitionedVector & tableauRow) const
2631{
2632  int maximumRows = model_->maximumAbcNumberRows();
2633  int first=tableauRow.startPartition(iBlock);
2634  int last=tableauRow.startPartition(iBlock+1);
2635  assert(first==startColumnBlock_[iBlock]);
2636  assert(last==startColumnBlock_[iBlock+1]);
2637  int numberSlacks=updateTableau.getNumElements();
2638  int numberNonZero=0;
2639  const double *  COIN_RESTRICT piTableau=updateTableau.denseVector();
2640  const double *  COIN_RESTRICT pi=updateDjs.denseVector();
2641  int *  COIN_RESTRICT index = tableauRow.getIndices()+first;
2642  double *  COIN_RESTRICT array = tableauRow.denseVector()+first;
2643  const unsigned char *  COIN_RESTRICT internalStatus = model_->internalStatus();
2644  double * COIN_RESTRICT reducedCost = model_->djRegion();
2645  if (!iBlock) {
2646    const int *  COIN_RESTRICT piIndexTableau = updateTableau.getIndices();
2647    numberNonZero=numberSlacks;
2648    for (int i=0;i<numberSlacks;i++) {
2649      int iRow=piIndexTableau[i];
2650      index[i]=iRow;
2651      array[i]=piTableau[iRow]; // ? what if small or basic
2652    }
2653    const int *  COIN_RESTRICT piIndex = updateDjs.getIndices();
2654    int number=updateDjs.getNumElements();
2655    for (int i=0;i<number;i++) {
2656      int iRow=piIndex[i];
2657      reducedCost[iRow] -= pi[iRow]; // sign?
2658    }
2659    first=maximumRows;
2660  }
2661  double zeroTolerance = model_->zeroTolerance();
2662  // get matrix data pointers
2663  const int * COIN_RESTRICT row = matrix_->getIndices();
2664  const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts()-maximumRows;
2665  const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
2666  for (int iSequence = first; iSequence < last; iSequence++) {
2667    unsigned char iStatus=internalStatus[iSequence]&7;
2668    if (iStatus<6) {
2669      CoinBigIndex start = columnStart[iSequence];
2670      CoinBigIndex next = columnStart[iSequence+1];
2671      double tableauValue = 0.0;
2672      double djValue=reducedCost[iSequence];
2673      for (CoinBigIndex j = start; j < next; j++) {
2674        int jRow = row[j];
2675        tableauValue += piTableau[jRow] * elementByColumn[j];
2676        djValue -= pi[jRow] * elementByColumn[j]; // sign?
2677      }
2678      reducedCost[iSequence]=djValue;
2679      if (fabs(tableauValue)>zeroTolerance) {
2680        index[numberNonZero]=iSequence;
2681        array[numberNonZero++]=tableauValue;
2682      }
2683    }
2684  }
2685  tableauRow.setNumElementsPartition(iBlock,numberNonZero);
2686  return numberSlacks;
2687}
2688/* does steepest edge double or triple update
2689   If scaleFactor!=0 then use with tableau row to update djs
2690   otherwise use updateForDjs
2691*/
2692int
2693AbcMatrix::primalColumnDouble(int iBlock,CoinPartitionedVector & updateForTableauRow,
2694                              CoinPartitionedVector & updateForDjs,
2695                              const CoinIndexedVector & updateForWeights,
2696                              CoinPartitionedVector & spareColumn1,
2697                              double * infeasibilities, 
2698                              double referenceIn, double devex,
2699                              // Array for exact devex to say what is in reference framework
2700                              unsigned int * reference,
2701                              double * weights, double scaleFactor) const
2702{
2703  int maximumRows = model_->maximumAbcNumberRows();
2704  int first=startColumnBlock_[iBlock];
2705  int last=startColumnBlock_[iBlock+1];
2706  const double *  COIN_RESTRICT piTableau=updateForTableauRow.denseVector();
2707  const double *  COIN_RESTRICT pi=updateForDjs.denseVector();
2708  const double *  COIN_RESTRICT piWeights=updateForWeights.denseVector();
2709  const unsigned char *  COIN_RESTRICT internalStatus = model_->internalStatus();
2710  double * COIN_RESTRICT reducedCost = model_->djRegion();
2711  double tolerance = model_->currentDualTolerance();
2712  // use spareColumn to track new infeasibilities
2713  int * COIN_RESTRICT newInfeas = spareColumn1.getIndices()+first;
2714  int numberNewInfeas=0;
2715  // we can't really trust infeasibilities if there is dual error
2716  // this coding has to mimic coding in checkDualSolution
2717  double error = CoinMin(1.0e-2, model_->largestDualError());
2718  // allow tolerance at least slightly bigger than standard
2719  tolerance = tolerance +  error;
2720  double mult[2]={1.0,-1.0};
2721  bool updateWeights=devex!=0.0;
2722#define isReference(i) (((reference[i>>5] >> (i & 31)) & 1) != 0)
2723  // Scale factor is other way round so tableau row is 1.0* while dj update is scaleFactor*
2724  if (!iBlock) {
2725    int numberSlacks=updateForTableauRow.getNumElements();
2726    const int *  COIN_RESTRICT piIndexTableau = updateForTableauRow.getIndices();
2727    if (!scaleFactor) {
2728      const int *  COIN_RESTRICT piIndex = updateForDjs.getIndices();
2729      int number=updateForDjs.getNumElements();
2730      for (int i=0;i<number;i++) {
2731        int iRow=piIndex[i];
2732        int iStatus=internalStatus[iRow]&7;
2733        double value=reducedCost[iRow];
2734        value += pi[iRow];
2735        if (iStatus<4) {
2736          reducedCost[iRow]=value;
2737          value *= mult[iStatus];
2738          if (value<0.0) {
2739            if (!infeasibilities[iRow])
2740              newInfeas[numberNewInfeas++]=iRow;
2741#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
2742            infeasibilities[iRow]=value*value*CLP_PRIMAL_SLACK_MULTIPLIER;
2743#else
2744            infeasibilities[iRow]=value*value;
2745#endif
2746          } else {
2747            if (infeasibilities[iRow])
2748              infeasibilities[iRow]=COIN_INDEXED_REALLY_TINY_ELEMENT;
2749          }
2750        }
2751      }
2752    } else if (scaleFactor!=COIN_DBL_MAX) {
2753      for (int i=0;i<numberSlacks;i++) {
2754        int iRow=piIndexTableau[i];
2755        int iStatus=internalStatus[iRow]&7;
2756        if (iStatus<4) {
2757          double value=reducedCost[iRow];
2758          value += scaleFactor*piTableau[iRow]; // sign?
2759          reducedCost[iRow]=value;
2760          value *= mult[iStatus];
2761          if (value<0.0) {
2762            if (!infeasibilities[iRow])
2763              newInfeas[numberNewInfeas++]=iRow;
2764#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
2765            infeasibilities[iRow]=value*value*CLP_PRIMAL_SLACK_MULTIPLIER;
2766#else
2767            infeasibilities[iRow]=value*value;
2768#endif
2769          } else {
2770            if (infeasibilities[iRow])
2771              infeasibilities[iRow]=COIN_INDEXED_REALLY_TINY_ELEMENT;
2772          }
2773        }
2774      }
2775    }
2776    if (updateWeights) {
2777      for (int i=0;i<numberSlacks;i++) {
2778        int iRow=piIndexTableau[i];
2779        double modification=piWeights[iRow];
2780        double thisWeight = weights[iRow];
2781        double pivot = piTableau[iRow];
2782        double pivotSquared = pivot * pivot;
2783        thisWeight += pivotSquared * devex - pivot * modification;
2784        if (thisWeight < DEVEX_TRY_NORM) {
2785          if (referenceIn < 0.0) {
2786            // steepest
2787            thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
2788          } else {
2789            // exact
2790            thisWeight = referenceIn * pivotSquared;
2791            if (isReference(iRow))
2792              thisWeight += 1.0;
2793            thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
2794          }
2795        }
2796        weights[iRow] = thisWeight;
2797      }
2798    }
2799    first=maximumRows;
2800  }
2801  double zeroTolerance = model_->zeroTolerance();
2802  double freeTolerance = FREE_ACCEPT * tolerance;;
2803  // get matrix data pointers
2804  const int * COIN_RESTRICT row = matrix_->getIndices();
2805  const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts()-maximumRows;
2806  const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
2807  bool updateDjs=scaleFactor!=COIN_DBL_MAX;
2808  for (int iSequence = first; iSequence < last; iSequence++) {
2809    unsigned char iStatus=internalStatus[iSequence]&7;
2810    if (iStatus<6) {
2811      CoinBigIndex start = columnStart[iSequence];
2812      CoinBigIndex next = columnStart[iSequence+1];
2813      double tableauValue = 0.0;
2814      double djValue=reducedCost[iSequence];
2815      if (!scaleFactor) {
2816        for (CoinBigIndex j = start; j < next; j++) {
2817          int jRow = row[j];
2818          tableauValue += piTableau[jRow] * elementByColumn[j];
2819          djValue += pi[jRow] * elementByColumn[j];
2820        }
2821      } else {
2822        for (CoinBigIndex j = start; j < next; j++) {
2823          int jRow = row[j];
2824          tableauValue += piTableau[jRow] * elementByColumn[j];
2825        }
2826        djValue += tableauValue*scaleFactor; // sign?
2827      }
2828      if (updateDjs) {
2829        reducedCost[iSequence]=djValue;
2830        if (iStatus<4) {
2831          djValue *= mult[iStatus];
2832          if (djValue<0.0) {
2833            if (!infeasibilities[iSequence])
2834              newInfeas[numberNewInfeas++]=iSequence;
2835            infeasibilities[iSequence]=djValue*djValue;
2836          } else {
2837            if (infeasibilities[iSequence])
2838              infeasibilities[iSequence]=COIN_INDEXED_REALLY_TINY_ELEMENT;
2839          }
2840        } else {
2841          if (fabs(djValue) > freeTolerance) {
2842            // we are going to bias towards free (but only if reasonable)
2843            djValue *= FREE_BIAS;
2844            if (!infeasibilities[iSequence])
2845              newInfeas[numberNewInfeas++]=iSequence;
2846            // store square in list
2847            infeasibilities[iSequence] = djValue * djValue;
2848          } else {
2849            if (infeasibilities[iSequence])
2850              infeasibilities[iSequence]=COIN_INDEXED_REALLY_TINY_ELEMENT;
2851          }
2852        }
2853      }
2854      if (fabs(tableauValue)>zeroTolerance&&updateWeights) {
2855        double modification=0.0;
2856        for (CoinBigIndex j = start; j < next; j++) {
2857          int jRow = row[j];
2858          modification += piWeights[jRow] * elementByColumn[j];
2859        }
2860        double thisWeight = weights[iSequence];
2861        double pivot = tableauValue;
2862        double pivotSquared = pivot * pivot;
2863        thisWeight += pivotSquared * devex - pivot * modification;
2864        if (thisWeight < DEVEX_TRY_NORM) {
2865          if (referenceIn < 0.0) {
2866            // steepest
2867            thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
2868          } else {
2869            // exact
2870            thisWeight = referenceIn * pivotSquared;
2871            if (isReference(iSequence))
2872              thisWeight += 1.0;
2873            thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
2874          }
2875        }
2876        weights[iSequence] = thisWeight;
2877      }
2878    }
2879  }
2880  spareColumn1.setTempNumElementsPartition(iBlock,numberNewInfeas);
2881  int bestSequence=-1;
2882#if 0
2883  if (!iBlock)
2884    first=0;
2885  // not at present - maybe later
2886  double bestDj=tolerance*tolerance;
2887  for (int iSequence=first;iSequence<last;iSequence++) {
2888    if (infeasibilities[iSequence]>bestDj*weights[iSequence]) {
2889      int iStatus=internalStatus[iSequence]&7;
2890      assert (iStatus<6);
2891      bestSequence=iSequence;
2892      bestDj=infeasibilities[iSequence]/weights[iSequence];
2893    }
2894  }
2895#endif
2896  return bestSequence;
2897}
2898/* does steepest edge double or triple update
2899   If scaleFactor!=0 then use with tableau row to update djs
2900   otherwise use updateForDjs
2901*/
2902int
2903AbcMatrix::primalColumnSparseDouble(int iBlock,CoinPartitionedVector & updateForTableauRow,
2904                              CoinPartitionedVector & updateForDjs,
2905                              const CoinIndexedVector & updateForWeights,
2906                              CoinPartitionedVector & spareColumn1,
2907                              double * infeasibilities, 
2908                              double referenceIn, double devex,
2909                              // Array for exact devex to say what is in reference framework
2910                              unsigned int * reference,
2911                              double * weights, double scaleFactor) const
2912{
2913  int maximumRows = model_->maximumAbcNumberRows();
2914  int first=blockStart_[iBlock];
2915  //int last=blockStart_[iBlock+1];
2916  const double *  COIN_RESTRICT piTableau=updateForTableauRow.denseVector();
2917  //const double *  COIN_RESTRICT pi=updateForDjs.denseVector();
2918  const double *  COIN_RESTRICT piWeights=updateForWeights.denseVector();
2919  const unsigned char *  COIN_RESTRICT internalStatus = model_->internalStatus();
2920  double * COIN_RESTRICT reducedCost = model_->djRegion();
2921  double tolerance = model_->currentDualTolerance();
2922  // use spareColumn to track new infeasibilities
2923  int * COIN_RESTRICT newInfeas = spareColumn1.getIndices()+first;
2924  int numberNewInfeas=0;
2925  // we can't really trust infeasibilities if there is dual error
2926  // this coding has to mimic coding in checkDualSolution
2927  double error = CoinMin(1.0e-2, model_->largestDualError());
2928  // allow tolerance at least slightly bigger than standard
2929  tolerance = tolerance +  error;
2930  double mult[2]={1.0,-1.0};
2931  bool updateWeights=devex!=0.0;
2932  int numberSlacks=updateForTableauRow.getNumElements();
2933  const int *  COIN_RESTRICT piIndexTableau = updateForTableauRow.getIndices();
2934#define isReference(i) (((reference[i>>5] >> (i & 31)) & 1) != 0)
2935  // Scale factor is other way round so tableau row is 1.0* while dj update is scaleFactor*
2936  assert (scaleFactor);
2937  if (!iBlock) {
2938    if (scaleFactor!=COIN_DBL_MAX) {
2939      for (int i=0;i<numberSlacks;i++) {
2940        int iRow=piIndexTableau[i];
2941        int iStatus=internalStatus[iRow]&7;
2942        if (iStatus<4) {
2943          double value=reducedCost[iRow];
2944          value += scaleFactor*piTableau[iRow]; // sign?
2945          reducedCost[iRow]=value;
2946          value *= mult[iStatus];
2947          if (value<0.0) {
2948            if (!infeasibilities[iRow])
2949              newInfeas[numberNewInfeas++]=iRow;
2950#ifdef CLP_PRIMAL_SLACK_MULTIPLIER
2951            infeasibilities[iRow]=value*value*CLP_PRIMAL_SLACK_MULTIPLIER;
2952#else
2953            infeasibilities[iRow]=value*value;
2954#endif
2955          } else {
2956            if (infeasibilities[iRow])
2957              infeasibilities[iRow]=COIN_INDEXED_REALLY_TINY_ELEMENT;
2958          }
2959        }
2960      }
2961    }
2962    if (updateWeights) {
2963      for (int i=0;i<numberSlacks;i++) {
2964        int iRow=piIndexTableau[i];
2965        double modification=piWeights[iRow];
2966        double thisWeight = weights[iRow];
2967        double pivot = piTableau[iRow];
2968        double pivotSquared = pivot * pivot;
2969        thisWeight += pivotSquared * devex - pivot * modification;
2970        if (thisWeight < DEVEX_TRY_NORM) {
2971          if (referenceIn < 0.0) {
2972            // steepest
2973            thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
2974          } else {
2975            // exact
2976            thisWeight = referenceIn * pivotSquared;
2977            if (isReference(iRow))
2978              thisWeight += 1.0;
2979            thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
2980          }
2981        }
2982        weights[iRow] = thisWeight;
2983      }
2984    }
2985    first=maximumRows;
2986  }
2987  double zeroTolerance = model_->zeroTolerance();
2988  double freeTolerance = FREE_ACCEPT * tolerance;;
2989  // get matrix data pointers
2990  const int * COIN_RESTRICT row = matrix_->getIndices();
2991  const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts()-maximumRows;
2992  const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
2993  // get tableau row
2994  int *  COIN_RESTRICT index = updateForTableauRow.getIndices()+first;
2995  double *  COIN_RESTRICT array = updateForTableauRow.denseVector();
2996  int numberRows=model_->numberRows();
2997  const CoinBigIndex * COIN_RESTRICT rowStart = rowStart_+iBlock*numberRows;
2998  const CoinBigIndex * COIN_RESTRICT rowEnd = rowStart+numberRows;
2999  const double * COIN_RESTRICT element = element_;
3000  const int * COIN_RESTRICT column = column_;
3001  int numberInRow=0;
3002  if (numberSlacks>2) {
3003    for (int i=0;i<numberSlacks;i++) {
3004      int iRow=piIndexTableau[i];
3005      double piValue = piTableau[iRow];
3006      CoinBigIndex end = rowEnd[iRow];
3007      for (CoinBigIndex j=rowStart[iRow];j<end;j++) {
3008        int iSequence = column[j];
3009        double value = element[j]*piValue;
3010        double oldValue=array[iSequence];
3011        value += oldValue;
3012        if (!oldValue) {
3013          array[iSequence]= value; 
3014          index[numberInRow++]=iSequence;
3015        } else if (value) {     
3016          array[iSequence] = value; 
3017        } else {
3018          array[iSequence] = COIN_INDEXED_REALLY_TINY_ELEMENT;
3019        }
3020      }
3021    }
3022  } else if (numberSlacks==2) {
3023    int iRow0 = piIndexTableau[0];
3024    int iRow1 = piIndexTableau[1];
3025    CoinBigIndex end0 = rowEnd[iRow0];
3026    CoinBigIndex end1 = rowEnd[iRow1];
3027    if (end0-rowStart[iRow0]>end1-rowStart[iRow1]) {
3028      int temp=iRow0;
3029      iRow0=iRow1;
3030      iRow1=temp;
3031    }
3032    CoinBigIndex start = rowStart[iRow0];
3033    CoinBigIndex end = rowEnd[iRow0];
3034    double piValue = piTableau[iRow0];
3035    for (CoinBigIndex j=start;j<end;j++) {
3036      int iSequence = column[j];
3037      double value = element[j];
3038      array[iSequence]= piValue*value; 
3039      index[numberInRow++]=iSequence;
3040    }
3041    start = rowStart[iRow1];
3042    end = rowEnd[iRow1];
3043    piValue = piTableau[iRow1];
3044    for (CoinBigIndex j=start;j<end;j++) {
3045      int iSequence = column[j];
3046      double value = element[j];
3047      value *= piValue;
3048      if (!array[iSequence]) {
3049        array[iSequence]= value; 
3050        index[numberInRow++]=iSequence;
3051      } else { 
3052        value += array[iSequence];
3053        if (value)
3054          array[iSequence]= value; 
3055        else
3056          array[iSequence] = COIN_INDEXED_REALLY_TINY_ELEMENT;
3057      }
3058    }
3059  } else {
3060    int iRow0 = piIndexTableau[0];
3061    CoinBigIndex start = rowStart[iRow0];
3062    CoinBigIndex end = rowEnd[iRow0];
3063    double piValue = piTableau[iRow0];
3064    for (CoinBigIndex j=start;j<end;j++) {
3065      int iSequence = column[j];
3066      double value = element[j];
3067      array[iSequence]= piValue*value; 
3068      index[numberInRow++]=iSequence;
3069    }
3070  }
3071  bool updateDjs=scaleFactor!=COIN_DBL_MAX;
3072  for (int iLook =0 ;iLook<numberInRow;iLook++) {
3073    int iSequence=index[iLook];
3074    unsigned char iStatus=internalStatus[iSequence]&7;
3075    if (iStatus<6) {
3076      double tableauValue = array[iSequence];
3077      array[iSequence]=0.0;
3078      double djValue=reducedCost[iSequence];
3079      djValue += tableauValue*scaleFactor; // sign?
3080      if (updateDjs) {
3081        reducedCost[iSequence]=djValue;
3082        if (iStatus<4) {
3083          djValue *= mult[iStatus];
3084          if (djValue<0.0) {
3085            if (!infeasibilities[iSequence])
3086              newInfeas[numberNewInfeas++]=iSequence;
3087            infeasibilities[iSequence]=djValue*djValue;
3088          } else {
3089            if (infeasibilities[iSequence])
3090              infeasibilities[iSequence]=COIN_INDEXED_REALLY_TINY_ELEMENT;
3091          }
3092        } else {
3093          if (fabs(djValue) > freeTolerance) {
3094            // we are going to bias towards free (but only if reasonable)
3095            djValue *= FREE_BIAS;
3096            if (!infeasibilities[iSequence])
3097              newInfeas[numberNewInfeas++]=iSequence;
3098            // store square in list
3099            infeasibilities[iSequence] = djValue * djValue;
3100          } else {
3101            if (infeasibilities[iSequence])
3102              infeasibilities[iSequence]=COIN_INDEXED_REALLY_TINY_ELEMENT;
3103          }
3104        }
3105      }
3106      if (fabs(tableauValue)>zeroTolerance&&updateWeights) {
3107        CoinBigIndex start = columnStart[iSequence];
3108        CoinBigIndex next = columnStart[iSequence+1];
3109        double modification=0.0;
3110        for (CoinBigIndex j = start; j < next; j++) {
3111          int jRow = row[j];
3112          modification += piWeights[jRow] * elementByColumn[j];
3113        }
3114        double thisWeight = weights[iSequence];
3115        double pivot = tableauValue;
3116        double pivotSquared = pivot * pivot;
3117        thisWeight += pivotSquared * devex - pivot * modification;
3118        if (thisWeight < DEVEX_TRY_NORM) {
3119          if (referenceIn < 0.0) {
3120            // steepest
3121            thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
3122          } else {
3123            // exact
3124            thisWeight = referenceIn * pivotSquared;
3125            if (isReference(iSequence))
3126              thisWeight += 1.0;
3127            thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
3128          }
3129        }
3130        weights[iSequence] = thisWeight;
3131      }
3132    } else {
3133      array[iSequence]=0.0;
3134    }
3135  }
3136  spareColumn1.setTempNumElementsPartition(iBlock,numberNewInfeas);
3137  int bestSequence=-1;
3138#if 0
3139  if (!iBlock)
3140    first=0;
3141  // not at present - maybe later
3142  double bestDj=tolerance*tolerance;
3143  for (int iSequence=first;iSequence<last;iSequence++) {
3144    if (infeasibilities[iSequence]>bestDj*weights[iSequence]) {
3145      int iStatus=internalStatus[iSequence]&7;
3146      assert (iStatus<6);
3147      bestSequence=iSequence;
3148      bestDj=infeasibilities[iSequence]/weights[iSequence];
3149    }
3150  }
3151#endif
3152  return bestSequence;
3153}
3154#if 0
3155/* Chooses best weighted dj
3156 */
3157int
3158AbcMatrix::chooseBestDj(int iBlock,const CoinIndexedVector & infeasibilities,
3159                        const double * weights) const
3160{
3161  return 0;
3162}
3163#endif
3164/* does steepest edge double or triple update
3165   If scaleFactor!=0 then use with tableau row to update djs
3166   otherwise use updateForDjs
3167   Returns best
3168*/
3169int
3170AbcMatrix::primalColumnDouble(CoinPartitionedVector & updateForTableauRow,
3171                              CoinPartitionedVector & updateForDjs,
3172                              const CoinIndexedVector & updateForWeights,
3173                              CoinPartitionedVector & spareColumn1,
3174                              CoinIndexedVector & infeasible, 
3175                              double referenceIn, double devex,
3176                              // Array for exact devex to say what is in reference framework
3177                              unsigned int * reference,
3178                              double * weights, double scaleFactor) const
3179{
3180  //int maximumRows = model_->maximumAbcNumberRows();
3181   // rebalance
3182   rebalance();
3183#ifdef PRICE_IN_ABC_MATRIX
3184   int which[NUMBER_BLOCKS];
3185#endif
3186   double * infeasibilities=infeasible.denseVector();
3187   int bestSequence=-1;
3188   // see if worth using row copy
3189   int maximumRows=model_->maximumAbcNumberRows();
3190   int number=updateForTableauRow.getNumElements();
3191#ifdef GCC_4_9
3192   assert (number);
3193#else
3194   if (!number) {
3195     printf("Null tableau row!\n");
3196   }
3197#endif
3198   bool useRowCopy = (gotRowCopy()&&(number<<2)<maximumRows);
3199   //useRowCopy=false;
3200   if (!scaleFactor)
3201     useRowCopy=false; // look again later
3202   const int * starts;
3203   int numberBlocks;
3204   if (useRowCopy) {
3205     starts=blockStart_;
3206     numberBlocks=numberRowBlocks_;
3207   } else {
3208     starts=startColumnBlock_;
3209     numberBlocks=numberColumnBlocks_;
3210   }
3211   if (useRowCopy) {
3212     for (int i=0;i<numberBlocks;i++) 
3213#ifdef PRICE_IN_ABC_MATRIX
3214       which[i]=
3215#endif
3216         cilk_spawn primalColumnSparseDouble(i,updateForTableauRow,updateForDjs,updateForWeights,
3217                                       spareColumn1,
3218                                       infeasibilities,referenceIn,devex,reference,weights,scaleFactor);
3219     cilk_sync;
3220   } else {
3221     for (int i=0;i<numberBlocks;i++) 
3222#ifdef PRICE_IN_ABC_MATRIX
3223       which[i]=
3224#endif
3225         cilk_spawn primalColumnDouble(i,updateForTableauRow,updateForDjs,updateForWeights,
3226                                       spareColumn1,
3227                                       infeasibilities,referenceIn,devex,reference,weights,scaleFactor);
3228     cilk_sync;
3229   }
3230#ifdef PRICE_IN_ABC_MATRIX
3231   double bestValue=model_->dualTolerance();
3232   int sequenceIn[8]={-1,-1,-1,-1,-1,-1,-1,-1};
3233#endif
3234   assert (numberColumnBlocks_<=8);
3235#ifdef PRICE_IN_ABC_MATRIX
3236   double weightedDj[8];
3237   int nPossible=0;
3238#endif
3239   //const double * reducedCost = model_->djRegion();
3240   // use spareColumn to track new infeasibilities
3241   int numberInfeas=infeasible.getNumElements();
3242   int * COIN_RESTRICT infeas = infeasible.getIndices();
3243   const int * COIN_RESTRICT newInfeasAll = spareColumn1.getIndices();
3244   for (int i=0;i<numberColumnBlocks_;i++) {
3245     const int * COIN_RESTRICT newInfeas = newInfeasAll+starts[i];
3246     int numberNewInfeas=spareColumn1.getNumElements(i);
3247     spareColumn1.setTempNumElementsPartition(i,0);
3248     for (int j=0;j<numberNewInfeas;j++)
3249       infeas[numberInfeas++]=newInfeas[j];
3250#ifdef PRICE_IN_ABC_MATRIX
3251     if (which[i]>=0) {
3252       int iSequence=which[i];
3253       double value=fabs(infeasibilities[iSequence]/weights[iSequence]);
3254       if (value>bestValue) {
3255         bestValue=value;
3256         bestSequence=which[i];
3257       }
3258       weightedDj[nPossible]=-value;
3259       sequenceIn[nPossible++]=iSequence;
3260     }
3261#endif
3262   }
3263   infeasible.setNumElements(numberInfeas);
3264#ifdef PRICE_IN_ABC_MATRIX
3265   CoinSort_2(weightedDj,weightedDj+nPossible,sequenceIn);
3266   //model_->setMultipleSequenceIn(sequenceIn);
3267#endif
3268   return bestSequence;
3269
3270}
3271// Partial pricing
3272void
3273AbcMatrix::partialPricing(double startFraction, double endFraction,
3274                          int & bestSequence, int & numberWanted)
3275{
3276  numberWanted = currentWanted_;
3277  int maximumRows = model_->maximumAbcNumberRows();
3278  // get matrix data pointers
3279  const int * COIN_RESTRICT row = matrix_->getIndices();
3280  const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts()-maximumRows;
3281  const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
3282  int numberColumns=model_->numberColumns();
3283  int start = static_cast<int> (startFraction * numberColumns);
3284  int end = CoinMin(static_cast<int> (endFraction * numberColumns + 1), numberColumns);
3285  // adjust
3286  start += maximumRows;
3287  end += maximumRows;
3288  double tolerance = model_->currentDualTolerance();
3289  double * reducedCost = model_->djRegion();
3290  const double * duals = model_->dualRowSolution();
3291  const double * cost = model_->costRegion();
3292  double bestDj;
3293  if (bestSequence >= 0)
3294    bestDj = fabs(reducedCost[bestSequence]);
3295  else
3296    bestDj = tolerance;
3297  int sequenceOut = model_->sequenceOut();
3298  int saveSequence = bestSequence;
3299  int lastScan = minimumObjectsScan_ < 0 ? end : start + minimumObjectsScan_;
3300  int minNeg = minimumGoodReducedCosts_ == -1 ? numberWanted : minimumGoodReducedCosts_;
3301  for (int iSequence = start; iSequence < end; iSequence++) {
3302    if (iSequence != sequenceOut) {
3303      double value;
3304      AbcSimplex::Status status = model_->getInternalStatus(iSequence);
3305      switch(status) {
3306      case AbcSimplex::basic:
3307      case AbcSimplex::isFixed:
3308        break;
3309      case AbcSimplex::isFree:
3310      case AbcSimplex::superBasic:
3311        value = cost[iSequence];
3312        for (CoinBigIndex j = columnStart[iSequence];
3313             j < columnStart[iSequence+1]; j++) {
3314          int jRow = row[j];
3315          value -= duals[jRow] * elementByColumn[j];
3316        }
3317        value = fabs(value);
3318        if (value > FREE_ACCEPT * tolerance) {
3319          numberWanted--;
3320          // we are going to bias towards free (but only if reasonable)
3321          value *= FREE_BIAS;
3322          if (value > bestDj) {
3323            // check flagged variable and correct dj
3324            if (!model_->flagged(iSequence)) {
3325              bestDj = value;
3326              bestSequence = iSequence;
3327            } else {
3328              // just to make sure we don't exit before got something
3329              numberWanted++;
3330            }
3331          }
3332        }
3333        break;
3334      case AbcSimplex::atUpperBound:
3335        value = cost[iSequence];
3336        // scaled
3337        for (CoinBigIndex j = columnStart[iSequence];
3338             j < columnStart[iSequence+1]; j++) {
3339          int jRow = row[j];
3340          value -= duals[jRow] * elementByColumn[j];
3341        }
3342        if (value > tolerance) {
3343          numberWanted--;
3344          if (value > bestDj) {
3345            // check flagged variable and correct dj
3346            if (!model_->flagged(iSequence)) {
3347              bestDj = value;
3348              bestSequence = iSequence;
3349            } else {
3350              // just to make sure we don't exit before got something
3351              numberWanted++;
3352            }
3353          }
3354        }
3355        break;
3356      case AbcSimplex::atLowerBound:
3357        value = cost[iSequence];
3358        for (CoinBigIndex j = columnStart[iSequence];
3359             j < columnStart[iSequence+1]; j++) {
3360          int jRow = row[j];
3361          value -= duals[jRow] * elementByColumn[j];
3362        }
3363        value = -value;
3364        if (value > tolerance) {
3365          numberWanted--;
3366          if (value > bestDj) {
3367            // check flagged variable and correct dj
3368            if (!model_->flagged(iSequence)) {
3369              bestDj = value;
3370              bestSequence = iSequence;
3371            } else {
3372              // just to make sure we don't exit before got something
3373              numberWanted++;
3374            }
3375          }
3376        }
3377        break;
3378      }
3379    }
3380    if (numberWanted + minNeg < originalWanted_ && iSequence > lastScan) {
3381      // give up
3382      break;
3383    }
3384    if (!numberWanted)
3385      break;
3386  }
3387  if (bestSequence != saveSequence) {
3388    // recompute dj
3389    double value = cost[bestSequence];
3390    for (CoinBigIndex j = columnStart[bestSequence];
3391         j < columnStart[bestSequence+1]; j++) {
3392      int jRow = row[j];
3393      value -= duals[jRow] * elementByColumn[j];
3394    }
3395    reducedCost[bestSequence] = value;
3396    savedBestSequence_ = bestSequence;
3397    savedBestDj_ = reducedCost[savedBestSequence_];
3398  }
3399  currentWanted_ = numberWanted;
3400}
3401/* Return <code>x *A</code> in <code>z</code> but
3402   just for indices Already in z.
3403   Note - z always packed mode */
3404void 
3405AbcMatrix::subsetTransposeTimes(const CoinIndexedVector & x,
3406                                CoinIndexedVector & z) const
3407{
3408  int maximumRows = model_->maximumAbcNumberRows();
3409  // get matrix data pointers
3410  const int * COIN_RESTRICT row = matrix_->getIndices();
3411  const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts()-maximumRows;
3412  const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
3413  const double *  COIN_RESTRICT other=x.denseVector();
3414  int numberNonZero = z.getNumElements();
3415  int *  COIN_RESTRICT index = z.getIndices();
3416  double *  COIN_RESTRICT array = z.denseVector();
3417  int numberRows=model_->numberRows();
3418  for (int i=0;i<numberNonZero;i++) {
3419    int iSequence=index[i];
3420    if (iSequence>=numberRows) {
3421      CoinBigIndex start = columnStart[iSequence];
3422      CoinBigIndex next = columnStart[iSequence+1];
3423      double value = 0.0;
3424      for (CoinBigIndex j = start; j < next; j++) {
3425        int jRow = row[j];
3426        value -= other[jRow] * elementByColumn[j];
3427      }
3428      array[i]=value;
3429    } else { 
3430      array[i]=-other[iSequence];
3431    }
3432  }
3433}
3434// Return <code>-x *A</code> in <code>z</code>
3435void 
3436AbcMatrix::transposeTimes(const CoinIndexedVector & x,
3437                                CoinIndexedVector & z) const
3438{
3439  int maximumRows = model_->maximumAbcNumberRows();
3440  // get matrix data pointers
3441  const int * COIN_RESTRICT row = matrix_->getIndices();
3442  const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts()-maximumRows;
3443  const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
3444  const double *  COIN_RESTRICT other=x.denseVector();
3445  int numberNonZero = 0;
3446  int *  COIN_RESTRICT index = z.getIndices();
3447  double *  COIN_RESTRICT array = z.denseVector();
3448  int numberColumns=model_->numberColumns();
3449  double zeroTolerance=model_->zeroTolerance();
3450  for (int iSequence=maximumRows;iSequence<maximumRows+numberColumns;iSequence++) {
3451    CoinBigIndex start = columnStart[iSequence];
3452    CoinBigIndex next = columnStart[iSequence+1];
3453    double value = 0.0;
3454    for (CoinBigIndex j = start; j < next; j++) {
3455      int jRow = row[j];
3456      value -= other[jRow] * elementByColumn[j];
3457    }
3458    if (fabs(value)>zeroTolerance) {
3459      // TEMP
3460      array[iSequence-maximumRows]=value;
3461      index[numberNonZero++]=iSequence-maximumRows;
3462    }
3463  }
3464  z.setNumElements(numberNonZero);
3465}
3466/* Return A * scalar(+-1) *x + y</code> in <code>y</code>.
3467   @pre <code>x</code> must be of size <code>numRows()</code>
3468   @pre <code>y</code> must be of size <code>numRows()+numColumns()</code> */
3469void 
3470AbcMatrix::transposeTimesNonBasic(double scalar,
3471                                  const double * pi, double * y) const
3472{
3473  int numberTotal = model_->numberTotal();
3474  int maximumRows = model_->maximumAbcNumberRows();
3475  assert (scalar==-1.0);
3476  // get matrix data pointers
3477  const int * COIN_RESTRICT row = matrix_->getIndices();
3478  const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts()-maximumRows;
3479  const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
3480  int numberRows=model_->numberRows();
3481  const unsigned char *  COIN_RESTRICT status=model_->internalStatus();
3482  for (int i=0;i<numberRows;i++) {
3483    y[i]+=scalar*pi[i];
3484  }
3485  for (int iColumn = maximumRows; iColumn < numberTotal;iColumn++) {
3486    unsigned char type=status[iColumn]&7;
3487    if (type<6) {
3488      CoinBigIndex start = columnStart[iColumn];
3489      CoinBigIndex next = columnStart[iColumn+1];
3490      double value = 0.0;
3491      for (CoinBigIndex j = start; j < next; j++) {
3492        int jRow = row[j];
3493        value += scalar*pi[jRow] * elementByColumn[j];
3494      }
3495      y[iColumn]+=value;
3496    } else {
3497      y[iColumn]=0.0; // may not be but .....
3498    }
3499  }
3500}
3501/* Return y - A * x</code> in <code>y</code>.
3502   @pre <code>x</code> must be of size <code>numRows()</code>
3503   @pre <code>y</code> must be of size <code>numRows()+numColumns()</code> */
3504void 
3505AbcMatrix::transposeTimesAll(const double * pi, double * y) const
3506{
3507  int numberTotal = model_->numberTotal();
3508  int maximumRows = model_->maximumAbcNumberRows();
3509  // get matrix data pointers
3510  const int * COIN_RESTRICT row = matrix_->getIndices();
3511  const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts()-maximumRows;
3512  const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
3513  int numberRows=model_->numberRows();
3514  for (int i=0;i<numberRows;i++) {
3515    y[i]-=pi[i];
3516  }
3517  for (int iColumn = maximumRows; iColumn < numberTotal;iColumn++) {
3518    CoinBigIndex start = columnStart[iColumn];
3519    CoinBigIndex next = columnStart[iColumn+1];
3520    double value = 0.0;
3521    for (CoinBigIndex j = start; j < next; j++) {
3522      int jRow = row[j];
3523      value += pi[jRow] * elementByColumn[j];
3524    }
3525    y[iColumn]-=value;
3526  }
3527}
3528/* Return y  + A * scalar(+-1) *x</code> in <code>y</code>.
3529   @pre <code>x</code> must be of size <code>numRows()</code>
3530   @pre <code>y</code> must be of size <code>numRows()</code> */
3531void 
3532AbcMatrix::transposeTimesBasic(double scalar,
3533                               const double * pi, double * y) const
3534{
3535  assert (scalar==-1.0);
3536  int numberRows=model_->numberRows();
3537  //AbcMemset0(y,numberRows);
3538  // get matrix data pointers
3539  const int * COIN_RESTRICT row = matrix_->getIndices();
3540  const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts()-model_->maximumAbcNumberRows();
3541  const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
3542  const int *  COIN_RESTRICT pivotVariable = model_->pivotVariable();
3543  for (int i=0;i<numberRows;i++) {
3544    int iSequence=pivotVariable[i];
3545    double value;
3546    if (iSequence>=numberRows) {
3547      CoinBigIndex start = columnStart[iSequence];
3548      CoinBigIndex next = columnStart[iSequence+1];
3549      value = 0.0;
3550      for (CoinBigIndex j = start; j < next; j++) {
3551        int jRow = row[j];
3552        value += pi[jRow] * elementByColumn[j];
3553      }
3554    } else {
3555      value=pi[iSequence];
3556    }
3557    y[i]+=scalar*value;
3558  }
3559}
3560// Adds multiple of a column (or slack) into an CoinIndexedvector
3561void 
3562AbcMatrix::add(CoinIndexedVector & rowArray, int iSequence, double multiplier) const 
3563{
3564  int maximumRows=model_->maximumAbcNumberRows();
3565  if (iSequence>=maximumRows) {
3566    const int * COIN_RESTRICT row = matrix_->getIndices();
3567    iSequence -= maximumRows;
3568    const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
3569    const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
3570    CoinBigIndex start = columnStart[iSequence];
3571    CoinBigIndex next = columnStart[iSequence+1];
3572    for (CoinBigIndex j = start; j < next; j++) {
3573      int jRow = row[j];
3574      rowArray.quickAdd(jRow,elementByColumn[j]*multiplier);
3575    }
3576  } else {
3577    rowArray.quickAdd(iSequence,multiplier);
3578  } 
3579}
3580/* Unpacks a column into an CoinIndexedVector
3581 */
3582void 
3583AbcMatrix::unpack(CoinIndexedVector & usefulArray,
3584                  int sequence) const 
3585{
3586  usefulArray.clear();
3587  int maximumRows=model_->maximumAbcNumberRows();
3588  if (sequence < maximumRows) {
3589    //slack
3590    usefulArray.insert(sequence , 1.0);
3591  } else {
3592    // column
3593    const CoinBigIndex *  COIN_RESTRICT columnStart = matrix()->getVectorStarts()-maximumRows;
3594    CoinBigIndex start=columnStart[sequence];
3595    CoinBigIndex end=columnStart[sequence+1];
3596    const int *  COIN_RESTRICT row = matrix()->getIndices()+start;
3597    const double *  COIN_RESTRICT elementByColumn = matrix()->getElements()+start;
3598    int *  COIN_RESTRICT index = usefulArray.getIndices();
3599    double *  COIN_RESTRICT array = usefulArray.denseVector();
3600    int numberNonZero=end-start;
3601    for (int j=0; j<numberNonZero;j++) {
3602      int iRow=row[j];
3603      index[j]=iRow;
3604      array[iRow]=elementByColumn[j];
3605    }
3606    usefulArray.setNumElements(numberNonZero);
3607  }
3608}
3609// Row starts
3610CoinBigIndex * 
3611AbcMatrix::rowStart() const
3612{
3613  return rowStart_;
3614}
3615// Row ends
3616CoinBigIndex * 
3617AbcMatrix::rowEnd() const
3618{
3619  //if (numberRowBlocks_<2) {
3620  //return rowStart_+1;
3621  //} else {
3622    return rowStart_+model_->numberRows()*(numberRowBlocks_+1);
3623    //}
3624}
3625// Row elements
3626double * 
3627AbcMatrix::rowElements() const
3628{
3629  return element_;
3630}
3631// Row columnss
3632CoinSimplexInt * 
3633AbcMatrix::rowColumns() const
3634{
3635  return column_;
3636}
Note: See TracBrowser for help on using the repository browser.