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

Last change on this file since 1878 was 1878, checked in by forrest, 7 years ago

minor changes to implement Aboca

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