source: trunk/Clp/src/ClpCholeskyDense.cpp @ 1367

Last change on this file since 1367 was 1367, checked in by forrest, 11 years ago

try and improve stability - also long double interior point

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 64.1 KB
Line 
1// Copyright (C) 2003, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4
5
6#include "CoinPragma.hpp"
7#include "CoinHelperFunctions.hpp"
8#include "ClpHelperFunctions.hpp"
9
10#include "ClpInterior.hpp"
11#include "ClpCholeskyDense.hpp"
12#include "ClpMessage.hpp"
13#include "ClpQuadraticObjective.hpp"
14
15//#############################################################################
16// Constructors / Destructor / Assignment
17//#############################################################################
18
19//-------------------------------------------------------------------
20// Default Constructor
21//-------------------------------------------------------------------
22ClpCholeskyDense::ClpCholeskyDense () 
23  : ClpCholeskyBase(),
24    borrowSpace_(false)
25{
26  type_=11;;
27}
28
29//-------------------------------------------------------------------
30// Copy constructor
31//-------------------------------------------------------------------
32ClpCholeskyDense::ClpCholeskyDense (const ClpCholeskyDense & rhs) 
33  : ClpCholeskyBase(rhs),
34    borrowSpace_(rhs.borrowSpace_)
35{
36  assert(!rhs.borrowSpace_||!rhs.sizeFactor_); // can't do if borrowing space
37}
38
39
40//-------------------------------------------------------------------
41// Destructor
42//-------------------------------------------------------------------
43ClpCholeskyDense::~ClpCholeskyDense ()
44{
45  if (borrowSpace_) {
46    // set NULL
47    sparseFactor_=NULL;
48    workDouble_=NULL;
49    diagonal_=NULL;
50  }
51}
52
53//----------------------------------------------------------------
54// Assignment operator
55//-------------------------------------------------------------------
56ClpCholeskyDense &
57ClpCholeskyDense::operator=(const ClpCholeskyDense& rhs)
58{
59  if (this != &rhs) {
60    assert(!rhs.borrowSpace_||!rhs.sizeFactor_); // can't do if borrowing space
61    ClpCholeskyBase::operator=(rhs);
62    borrowSpace_=rhs.borrowSpace_;
63  }
64  return *this;
65}
66//-------------------------------------------------------------------
67// Clone
68//-------------------------------------------------------------------
69ClpCholeskyBase * ClpCholeskyDense::clone() const
70{
71  return new ClpCholeskyDense(*this);
72}
73// If not power of 2 then need to redo a bit
74#define BLOCK 16
75#define BLOCKSHIFT 4
76// Block unroll if power of 2 and at least 8
77#define BLOCKUNROLL
78
79#define BLOCKSQ ( BLOCK*BLOCK )
80#define BLOCKSQSHIFT ( BLOCKSHIFT+BLOCKSHIFT )
81#define number_blocks(x) (((x)+BLOCK-1)>>BLOCKSHIFT)
82#define number_rows(x) ((x)<<BLOCKSHIFT)
83#define number_entries(x) ((x)<<BLOCKSQSHIFT)
84/* Gets space */
85int 
86ClpCholeskyDense::reserveSpace(const ClpCholeskyBase * factor, int numberRows) 
87{
88  numberRows_ = numberRows;
89  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
90  // allow one stripe extra
91  numberBlocks = numberBlocks + ((numberBlocks*(numberBlocks+1))/2);
92  sizeFactor_=numberBlocks*BLOCKSQ;
93  //#define CHOL_COMPARE
94#ifdef CHOL_COMPARE 
95  sizeFactor_ += 95000;
96#endif
97  if (!factor) {
98    sparseFactor_ = new longDouble [sizeFactor_];
99    rowsDropped_ = new char [numberRows_];
100    memset(rowsDropped_,0,numberRows_);
101    workDouble_ = new longDouble[numberRows_];
102    diagonal_ = new longDouble[numberRows_];
103  } else {
104    borrowSpace_=true;
105    int numberFull = factor->numberRows();
106    sparseFactor_ = factor->sparseFactor()+(factor->size()-sizeFactor_);
107    workDouble_ = factor->workDouble() + (numberFull-numberRows_);
108    diagonal_ = factor->diagonal() + (numberFull-numberRows_);
109  }
110  numberRowsDropped_=0;
111  return 0;
112}
113/* Returns space needed */
114CoinBigIndex
115ClpCholeskyDense::space( int numberRows) const
116{
117 int numberBlocks = (numberRows+BLOCK-1)>>BLOCKSHIFT;
118  // allow one stripe extra
119  numberBlocks = numberBlocks + ((numberBlocks*(numberBlocks+1))/2);
120  CoinBigIndex sizeFactor=numberBlocks*BLOCKSQ;
121#ifdef CHOL_COMPARE 
122  sizeFactor += 95000;
123#endif
124  return sizeFactor;
125 }
126/* Orders rows and saves pointer to matrix.and model */
127int 
128ClpCholeskyDense::order(ClpInterior * model) 
129{
130  model_=model;
131  int numberRows;
132  int numberRowsModel = model_->numberRows();
133  int numberColumns = model_->numberColumns();
134  if (!doKKT_) {
135    numberRows = numberRowsModel;
136  } else {
137    numberRows = 2*numberRowsModel+numberColumns;
138  }
139  reserveSpace(NULL,numberRows);
140  rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
141  return 0;
142}
143/* Does Symbolic factorization given permutation.
144   This is called immediately after order.  If user provides this then
145   user must provide factorize and solve.  Otherwise the default factorization is used
146   returns non-zero if not enough memory */
147int 
148ClpCholeskyDense::symbolic()
149{
150  return 0;
151}
152/* Factorize - filling in rowsDropped and returning number dropped */
153int 
154ClpCholeskyDense::factorize(const double * diagonal, int * rowsDropped) 
155{
156  assert (!borrowSpace_);
157  const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
158  const int * columnLength = model_->clpMatrix()->getVectorLengths();
159  const int * row = model_->clpMatrix()->getIndices();
160  const double * element = model_->clpMatrix()->getElements();
161  const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
162  const int * rowLength = rowCopy_->getVectorLengths();
163  const int * column = rowCopy_->getIndices();
164  const double * elementByRow = rowCopy_->getElements();
165  int numberColumns=model_->clpMatrix()->getNumCols();
166  CoinZeroN(sparseFactor_,sizeFactor_);
167  //perturbation
168  CoinWorkDouble perturbation=model_->diagonalPerturbation()*model_->diagonalNorm();
169  perturbation=perturbation*perturbation;
170  if (perturbation>1.0) {
171#ifdef COIN_DEVELOP
172    //if (model_->model()->logLevel()&4)
173      std::cout <<"large perturbation "<<perturbation<<std::endl;
174#endif
175    perturbation=CoinSqrt(perturbation);;
176    perturbation=1.0;
177  }
178  int iRow;
179  int newDropped=0;
180  CoinWorkDouble largest=1.0;
181  CoinWorkDouble smallest=COIN_DBL_MAX;
182  CoinWorkDouble delta2 = model_->delta(); // add delta*delta to diagonal
183  delta2 *= delta2;
184  if (!doKKT_) {
185    longDouble * work = sparseFactor_;
186    work--; // skip diagonal
187    int addOffset=numberRows_-1;
188    const double * diagonalSlack = diagonal + numberColumns;
189    // largest in initial matrix
190    CoinWorkDouble largest2=1.0e-20;
191    for (iRow=0;iRow<numberRows_;iRow++) {
192      if (!rowsDropped_[iRow]) {
193        CoinBigIndex startRow=rowStart[iRow];
194        CoinBigIndex endRow=rowStart[iRow]+rowLength[iRow];
195        CoinWorkDouble diagonalValue = diagonalSlack[iRow]+delta2;
196        for (CoinBigIndex k=startRow;k<endRow;k++) {
197          int iColumn=column[k];
198          CoinBigIndex start=columnStart[iColumn];
199          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
200          CoinWorkDouble multiplier = diagonal[iColumn]*elementByRow[k];
201          for (CoinBigIndex j=start;j<end;j++) {
202            int jRow=row[j];
203            if (!rowsDropped_[jRow]) {
204              if (jRow>iRow) {
205                CoinWorkDouble value=element[j]*multiplier;
206                work[jRow] += value;
207              } else if (jRow==iRow) {
208                CoinWorkDouble value=element[j]*multiplier;
209                diagonalValue += value;
210              }
211            }
212          }
213        } 
214        for (int j=iRow+1;j<numberRows_;j++) 
215          largest2 = CoinMax(largest2,CoinAbs(work[j]));
216        diagonal_[iRow]=diagonalValue;
217        largest2 = CoinMax(largest2,CoinAbs(diagonalValue));
218      } else {
219        // dropped
220        diagonal_[iRow]=1.0;
221      }
222      addOffset--;
223      work += addOffset;
224    }
225    //check sizes
226    largest2*=1.0e-20;
227    largest = CoinMin(largest2,CHOL_SMALL_VALUE);
228    int numberDroppedBefore=0;
229    for (iRow=0;iRow<numberRows_;iRow++) {
230      int dropped=rowsDropped_[iRow];
231      // Move to int array
232      rowsDropped[iRow]=dropped;
233      if (!dropped) {
234        CoinWorkDouble diagonal = diagonal_[iRow];
235        if (diagonal>largest2) {
236          diagonal_[iRow]=diagonal+perturbation;
237        } else {
238          diagonal_[iRow]=diagonal+perturbation;
239          rowsDropped[iRow]=2;
240          numberDroppedBefore++;
241        } 
242      }
243    }
244    doubleParameters_[10]=CoinMax(1.0e-20,largest);
245    integerParameters_[20]=0;
246    doubleParameters_[3]=0.0;
247    doubleParameters_[4]=COIN_DBL_MAX;
248    integerParameters_[34]=0; // say all must be positive
249#ifdef CHOL_COMPARE 
250    if (numberRows_<200)
251      factorizePart3(rowsDropped);
252#endif
253    factorizePart2(rowsDropped);
254    newDropped=integerParameters_[20]+numberDroppedBefore;
255    largest=doubleParameters_[3];
256    smallest=doubleParameters_[4];
257    if (model_->messageHandler()->logLevel()>1) 
258      std::cout<<"Cholesky - largest "<<largest<<" smallest "<<smallest<<std::endl;
259    choleskyCondition_=largest/smallest;
260    //drop fresh makes some formADAT easier
261    if (newDropped||numberRowsDropped_) {
262      newDropped=0;
263      for (int i=0;i<numberRows_;i++) {
264        char dropped = static_cast<char>(rowsDropped[i]);
265        rowsDropped_[i]=dropped;
266        if (dropped==2) {
267          //dropped this time
268          rowsDropped[newDropped++]=i;
269          rowsDropped_[i]=0;
270        } 
271      } 
272      numberRowsDropped_=newDropped;
273      newDropped=-(2+newDropped);
274    } 
275  } else {
276    // KKT
277    CoinPackedMatrix * quadratic = NULL;
278    ClpQuadraticObjective * quadraticObj = 
279      (dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject()));
280    if (quadraticObj) 
281      quadratic = quadraticObj->quadraticObjective();
282    // matrix
283    int numberRowsModel = model_->numberRows();
284    int numberColumns = model_->numberColumns();
285    int numberTotal = numberColumns + numberRowsModel;
286    longDouble * work = sparseFactor_;
287    work--; // skip diagonal
288    int addOffset=numberRows_-1;
289    int iColumn;
290    if (!quadratic) {
291      for (iColumn=0;iColumn<numberColumns;iColumn++) {
292        CoinWorkDouble value = diagonal[iColumn];
293        if (CoinAbs(value)>1.0e-100) {
294          value = 1.0/value;
295          largest = CoinMax(largest,CoinAbs(value));
296          diagonal_[iColumn] = -value;
297          CoinBigIndex start=columnStart[iColumn];
298          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
299          for (CoinBigIndex j=start;j<end;j++) {
300            //choleskyRow_[numberElements]=row[j]+numberTotal;
301            //sparseFactor_[numberElements++]=element[j];
302            work[row[j]+numberTotal]=element[j];
303            largest = CoinMax(largest,CoinAbs(element[j]));
304          }
305        } else {
306          diagonal_[iColumn] = -value;
307        }
308        addOffset--;
309        work += addOffset;
310      }
311    } else {
312      // Quadratic
313      const int * columnQuadratic = quadratic->getIndices();
314      const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
315      const int * columnQuadraticLength = quadratic->getVectorLengths();
316      const double * quadraticElement = quadratic->getElements();
317      for (iColumn=0;iColumn<numberColumns;iColumn++) {
318        CoinWorkDouble value = diagonal[iColumn];
319        CoinBigIndex j;
320        if (CoinAbs(value)>1.0e-100) {
321          value = 1.0/value;
322          for (j=columnQuadraticStart[iColumn];
323               j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
324            int jColumn = columnQuadratic[j];
325            if (jColumn>iColumn) {
326              work[jColumn]=-quadraticElement[j];
327            } else if (iColumn==jColumn) {
328              value += quadraticElement[j];
329            }
330          }
331          largest = CoinMax(largest,CoinAbs(value));
332          diagonal_[iColumn] = -value;
333          CoinBigIndex start=columnStart[iColumn];
334          CoinBigIndex end=columnStart[iColumn]+columnLength[iColumn];
335          for (j=start;j<end;j++) {
336            work[row[j]+numberTotal]=element[j];
337            largest = CoinMax(largest,CoinAbs(element[j]));
338          }
339        } else {
340          value = 1.0e100;
341          diagonal_[iColumn] = -value;
342        }
343        addOffset--;
344        work += addOffset;
345      }
346    }
347    // slacks
348    for (iColumn=numberColumns;iColumn<numberTotal;iColumn++) {
349      CoinWorkDouble value = diagonal[iColumn];
350      if (CoinAbs(value)>1.0e-100) {
351        value = 1.0/value;
352        largest = CoinMax(largest,CoinAbs(value));
353      } else {
354        value = 1.0e100;
355      }
356      diagonal_[iColumn] = -value;
357      work[iColumn-numberColumns+numberTotal]=-1.0;
358      addOffset--;
359      work += addOffset;
360    }
361    // Finish diagonal
362    for (iRow=0;iRow<numberRowsModel;iRow++) {
363      diagonal_[iRow+numberTotal]=delta2;
364    }
365    //check sizes
366    largest*=1.0e-20;
367    largest = CoinMin(largest,CHOL_SMALL_VALUE);
368    doubleParameters_[10]=CoinMax(1.0e-20,largest);
369    integerParameters_[20]=0;
370    doubleParameters_[3]=0.0;
371    doubleParameters_[4]=COIN_DBL_MAX;
372    // Set up LDL cutoff
373    integerParameters_[34]=numberTotal;
374    // KKT
375    int * rowsDropped2 = new int[numberRows_];
376    CoinZeroN(rowsDropped2,numberRows_);
377#ifdef CHOL_COMPARE 
378    if (numberRows_<200)
379      factorizePart3(rowsDropped2);
380#endif
381    factorizePart2(rowsDropped2);
382    newDropped=integerParameters_[20];
383    largest=doubleParameters_[3];
384    smallest=doubleParameters_[4];
385    if (model_->messageHandler()->logLevel()>1) 
386      std::cout<<"Cholesky - largest "<<largest<<" smallest "<<smallest<<std::endl;
387    choleskyCondition_=largest/smallest;
388    // Should save adjustments in ..R_
389    int n1=0,n2=0;
390    CoinWorkDouble * primalR = model_->primalR();
391    CoinWorkDouble * dualR = model_->dualR();
392    for (iRow=0;iRow<numberTotal;iRow++) { 
393      if (rowsDropped2[iRow]) {
394        n1++;
395        //printf("row region1 %d dropped\n",iRow);
396        //rowsDropped_[iRow]=1;
397        rowsDropped_[iRow]=0;
398        primalR[iRow]=doubleParameters_[20];
399      } else {
400        rowsDropped_[iRow]=0;
401        primalR[iRow]=0.0;
402      }
403    }
404    for (;iRow<numberRows_;iRow++) {
405      if (rowsDropped2[iRow]) {
406        n2++;
407        //printf("row region2 %d dropped\n",iRow);
408        //rowsDropped_[iRow]=1;
409        rowsDropped_[iRow]=0;
410        dualR[iRow-numberTotal]=doubleParameters_[34];
411      } else {
412        rowsDropped_[iRow]=0;
413        dualR[iRow-numberTotal]=0.0;
414      }
415    }
416  }
417  return 0;
418}
419/* Factorize - filling in rowsDropped and returning number dropped */
420void 
421ClpCholeskyDense::factorizePart3( int * rowsDropped) 
422{
423  int iColumn;
424  longDouble * xx = sparseFactor_;
425  longDouble * yy = diagonal_;
426  diagonal_ = sparseFactor_ + 40000;
427  sparseFactor_=diagonal_ + numberRows_;
428  //memcpy(sparseFactor_,xx,sizeFactor_*sizeof(double));
429  CoinMemcpyN(xx,40000,sparseFactor_);
430  CoinMemcpyN(yy,numberRows_,diagonal_);
431  int numberDropped=0;
432  CoinWorkDouble largest=0.0;
433  CoinWorkDouble smallest=COIN_DBL_MAX;
434  double dropValue = doubleParameters_[10];
435  int firstPositive=integerParameters_[34];
436  longDouble * work = sparseFactor_;
437  // Allow for triangular
438  int addOffset=numberRows_-1;
439  work--;
440  for (iColumn=0;iColumn<numberRows_;iColumn++) {
441    int iRow;
442    int addOffsetNow = numberRows_-1;;
443    longDouble * workNow=sparseFactor_-1+iColumn;
444    CoinWorkDouble diagonalValue = diagonal_[iColumn];
445    for (iRow=0;iRow<iColumn;iRow++) {
446      double aj = *workNow;
447      addOffsetNow--;
448      workNow += addOffsetNow;
449      diagonalValue -= aj*aj*workDouble_[iRow];
450    }
451    bool dropColumn=false;
452    if (iColumn<firstPositive) {
453      // must be negative
454      if (diagonalValue<=-dropValue) {
455        smallest = CoinMin(smallest,-diagonalValue);
456        largest = CoinMax(largest,-diagonalValue);
457        workDouble_[iColumn]=diagonalValue;
458        diagonalValue = 1.0/diagonalValue;
459      } else {
460        dropColumn=true;
461        workDouble_[iColumn]=-1.0e100;
462        diagonalValue=0.0;
463        integerParameters_[20]++;
464      }
465    } else {
466      // must be positive
467      if (diagonalValue>=dropValue) {
468        smallest = CoinMin(smallest,diagonalValue);
469        largest = CoinMax(largest,diagonalValue);
470        workDouble_[iColumn]=diagonalValue;
471        diagonalValue = 1.0/diagonalValue;
472      } else {
473        dropColumn=true;
474        workDouble_[iColumn]=1.0e100;
475        diagonalValue=0.0;
476        integerParameters_[20]++;
477      }
478    }
479    if (!dropColumn) {
480      diagonal_[iColumn]=diagonalValue;
481      for (iRow=iColumn+1;iRow<numberRows_;iRow++) {
482        double value = work[iRow];
483        workNow = sparseFactor_-1; 
484        int addOffsetNow = numberRows_-1;;
485        for (int jColumn=0;jColumn<iColumn;jColumn++) {
486          double aj = workNow[iColumn];
487          double multiplier = workDouble_[jColumn];
488          double ai = workNow[iRow];
489          addOffsetNow--;
490          workNow += addOffsetNow;
491          value -= aj*ai*multiplier;
492        }
493        work[iRow] = value*diagonalValue;
494      }
495    } else {
496      // drop column
497      rowsDropped[iColumn]=2;
498      numberDropped++;
499      diagonal_[iColumn]=0.0;
500      for (iRow=iColumn+1;iRow<numberRows_;iRow++) {
501        work[iRow] = 0.0;
502      }
503    }
504    addOffset--;
505    work += addOffset;
506  }
507  doubleParameters_[3]=largest;
508  doubleParameters_[4]=smallest;
509  integerParameters_[20]=numberDropped;
510  sparseFactor_=xx;
511  diagonal_=yy;
512}
513//#define POS_DEBUG
514#ifdef POS_DEBUG
515static int counter=0;
516int ClpCholeskyDense::bNumber(const longDouble * array, int &iRow, int &iCol)
517{
518  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
519  longDouble * a = sparseFactor_+BLOCKSQ*numberBlocks;
520  int k = array-a;
521  assert ((k%BLOCKSQ)==0);
522  iCol=0;
523  int kk=k>>BLOCKSQSHIFT;
524  //printf("%d %d %d %d\n",k,kk,BLOCKSQ,BLOCKSQSHIFT);
525  k=kk;
526  while (k>=numberBlocks) {
527    iCol++;
528    k -= numberBlocks;
529    numberBlocks--;
530  }
531  iRow = iCol+k;
532  counter++;
533  if(counter>100000)
534    exit(77);
535  return kk;
536 }
537#endif
538/* Factorize - filling in rowsDropped and returning number dropped */
539void 
540ClpCholeskyDense::factorizePart2( int * rowsDropped) 
541{
542  int iColumn;
543  //longDouble * xx = sparseFactor_;
544  //longDouble * yy = diagonal_;
545  //diagonal_ = sparseFactor_ + 40000;
546  //sparseFactor_=diagonal_ + numberRows_;
547  //memcpy(sparseFactor_,xx,sizeFactor_*sizeof(double));
548  //memcpy(sparseFactor_,xx,40000*sizeof(double));
549  //memcpy(diagonal_,yy,numberRows_*sizeof(double));
550  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
551  // later align on boundary
552  longDouble * a = sparseFactor_+BLOCKSQ*numberBlocks;
553  int n=numberRows_;
554  int nRound = numberRows_&(~(BLOCK-1));
555  // adjust if exact
556  if (nRound==n)
557    nRound -= BLOCK;
558  int sizeLastBlock=n-nRound;
559  int get = n*(n-1)/2; // ? as no diagonal
560  int block = numberBlocks*(numberBlocks+1)/2;
561  int ifOdd;
562  int rowLast;
563  if (sizeLastBlock!=BLOCK) {
564    longDouble * aa = &a[(block-1)*BLOCKSQ];
565    rowLast=nRound-1;
566    ifOdd=1;
567    int put = BLOCKSQ;
568    // do last separately
569    put -= (BLOCK-sizeLastBlock)*(BLOCK+1);
570    for (iColumn=numberRows_-1;iColumn>=nRound;iColumn--) {
571      int put2 = put;
572      put -= BLOCK;
573      for (int iRow=numberRows_-1;iRow>iColumn;iRow--) {
574        aa[--put2] = sparseFactor_[--get];
575        assert (aa+put2>=sparseFactor_+get);
576      }
577      // save diagonal as well
578      aa[--put2]=diagonal_[iColumn];
579    }
580    n=nRound;
581    block--;
582  } else {
583    // exact fit
584    rowLast=numberRows_-1;
585    ifOdd=0;
586  }
587  // Now main loop
588  int nBlock=0;
589  for (;n>0;n-=BLOCK) {
590    longDouble * aa = &a[(block-1)*BLOCKSQ];
591    longDouble * aaLast=NULL;
592    int put = BLOCKSQ;
593    int putLast=0;
594    // see if we have small block
595    if (ifOdd) {
596      aaLast = &a[(block-1)*BLOCKSQ];
597      aa = aaLast-BLOCKSQ;
598      putLast = BLOCKSQ-BLOCK+sizeLastBlock;
599    }
600    for (iColumn=n-1;iColumn>=n-BLOCK;iColumn--) {
601      if (aaLast) {
602        // last bit
603        for (int iRow=numberRows_-1;iRow>rowLast;iRow--) {
604          aaLast[--putLast] = sparseFactor_[--get];
605          assert (aaLast+putLast>=sparseFactor_+get);
606        }
607        putLast -= BLOCK-sizeLastBlock;
608      }
609      longDouble * aPut = aa;
610      int j=rowLast;
611      for (int jBlock=0;jBlock<=nBlock;jBlock++) {
612        int put2 = put;
613        int last = CoinMax(j-BLOCK,iColumn);
614        for (int iRow=j;iRow>last;iRow--) {
615          aPut[--put2] = sparseFactor_[--get];
616          assert (aPut+put2>=sparseFactor_+get);
617        }
618        if (j-BLOCK<iColumn) {
619          // save diagonal as well
620          aPut[--put2]=diagonal_[iColumn];
621        }
622        j -= BLOCK;
623        aPut -=BLOCKSQ;
624      }
625      put -= BLOCK;
626    }
627    nBlock++;
628    block -= nBlock+ifOdd;
629  }
630  factor(a,numberRows_,numberBlocks,
631         diagonal_,workDouble_,rowsDropped);
632  //sparseFactor_=xx;
633  //diagonal_=yy;
634}
635// Non leaf recursive factor
636void 
637ClpCholeskyDense::factor(longDouble * a, int n, int numberBlocks,
638            longDouble * diagonal, longDouble * work, int * rowsDropped)
639{
640  if (n<=BLOCK) {
641    factorLeaf(a,n,diagonal,work,rowsDropped);
642  } else {
643    int nb=number_blocks((n+1)>>1);
644    int nThis=number_rows(nb);
645    longDouble * aother;
646    int nLeft=n-nThis;
647    int nintri=(nb*(nb+1))>>1;
648    int nbelow=(numberBlocks-nb)*nb;
649    factor(a,nThis,numberBlocks,diagonal,work,rowsDropped);
650    triRec(a,nThis,a+number_entries(nb),diagonal,work,nLeft,nb,0,numberBlocks);
651    aother=a+number_entries(nintri+nbelow);
652    recTri(a+number_entries(nb),nLeft,nThis,nb,0,aother,diagonal,work,numberBlocks);
653    factor(aother,nLeft,
654        numberBlocks-nb,diagonal+nThis,work+nThis,rowsDropped);
655  }
656}
657// Non leaf recursive triangle rectangle update
658void 
659ClpCholeskyDense::triRec(longDouble * aTri, int nThis, longDouble * aUnder,
660                         longDouble * diagonal, longDouble * work,
661                         int nLeft, int iBlock, int jBlock,
662                         int numberBlocks)
663{
664  if (nThis<=BLOCK&&nLeft<=BLOCK) {
665    triRecLeaf(aTri,aUnder,diagonal,work,nLeft);
666  } else if (nThis<nLeft) {
667    int nb=number_blocks((nLeft+1)>>1);
668    int nLeft2=number_rows(nb);
669    triRec(aTri,nThis,aUnder,diagonal,work,nLeft2,iBlock,jBlock,numberBlocks);
670    triRec(aTri,nThis,aUnder+number_entries(nb),diagonal,work,nLeft-nLeft2,
671          iBlock+nb,jBlock,numberBlocks);
672  } else {
673    int nb=number_blocks((nThis+1)>>1);
674    int nThis2=number_rows(nb);
675    longDouble * aother;
676    int kBlock=jBlock+nb;
677    int i;
678    int nintri=(nb*(nb+1))>>1;
679    int nbelow=(numberBlocks-nb)*nb;
680    triRec(aTri,nThis2,aUnder,diagonal,work,nLeft,iBlock,jBlock,numberBlocks);
681    /* and rectangular update */
682    i=((numberBlocks-jBlock)*(numberBlocks-jBlock-1)-
683      (numberBlocks-jBlock-nb)*(numberBlocks-jBlock-nb-1))>>1;
684    aother=aUnder+number_entries(i);
685    recRec(aTri+number_entries(nb),nThis-nThis2,nLeft,nThis2,aUnder,aother,
686          work,kBlock,jBlock,numberBlocks);
687    triRec(aTri+number_entries(nintri+nbelow),nThis-nThis2,aother,diagonal+nThis2,
688           work+nThis2,nLeft,
689      iBlock-nb,kBlock-nb,numberBlocks-nb);
690  }
691}
692// Non leaf recursive rectangle triangle update
693void 
694ClpCholeskyDense::recTri(longDouble * aUnder, int nTri, int nDo,
695                         int iBlock, int jBlock,longDouble * aTri,
696                         longDouble * diagonal, longDouble * work,
697                         int numberBlocks)
698{
699  if (nTri<=BLOCK&&nDo<=BLOCK) {
700    recTriLeaf(aUnder,aTri,diagonal,work,nTri);
701  } else if (nTri<nDo) {
702    int nb=number_blocks((nDo+1)>>1);
703    int nDo2=number_rows(nb);
704    longDouble * aother;
705    int i;
706    recTri(aUnder,nTri,nDo2,iBlock,jBlock,aTri,diagonal,work,numberBlocks);
707    i=((numberBlocks-jBlock)*(numberBlocks-jBlock-1)-
708      (numberBlocks-jBlock-nb)*(numberBlocks-jBlock-nb-1))>>1;
709    aother=aUnder+number_entries(i);
710    recTri(aother,nTri,nDo-nDo2,iBlock-nb,jBlock,aTri,diagonal+nDo2,
711           work+nDo2,numberBlocks-nb);
712  } else {
713    int nb=number_blocks((nTri+1)>>1);
714    int nTri2=number_rows(nb);
715    longDouble * aother;
716    int i;
717    recTri(aUnder,nTri2,nDo,iBlock,jBlock,aTri,diagonal,work,numberBlocks);
718    /* and rectangular update */
719    i=((numberBlocks-iBlock)*(numberBlocks-iBlock+1)-
720      (numberBlocks-iBlock-nb)*(numberBlocks-iBlock-nb+1))>>1;
721    aother=aTri+number_entries(nb);
722    recRec(aUnder,nTri2,nTri-nTri2,nDo,aUnder+number_entries(nb),aother,
723           work,iBlock,jBlock,numberBlocks);
724    recTri(aUnder+number_entries(nb),nTri-nTri2,nDo,iBlock+nb,jBlock,
725         aTri+number_entries(i),diagonal,work,numberBlocks);
726  }
727}
728/* Non leaf recursive rectangle rectangle update,
729   nUnder is number of rows in iBlock,
730   nUnderK is number of rows in kBlock
731*/
732void 
733ClpCholeskyDense::recRec(longDouble * above, int nUnder, int nUnderK,
734                         int nDo, longDouble * aUnder, longDouble *aOther,
735                         longDouble * work,
736                         int iBlock, int jBlock,
737                         int numberBlocks)
738{
739  if (nDo<=BLOCK&&nUnder<=BLOCK&&nUnderK<=BLOCK) {
740    assert (nDo==BLOCK&&nUnder==BLOCK);
741    recRecLeaf(above , aUnder ,  aOther, work, nUnderK);
742  } else if (nDo<=nUnderK&&nUnder<=nUnderK) {
743    int nb=number_blocks((nUnderK+1)>>1);
744    int nUnder2=number_rows(nb);
745    recRec(above,nUnder,nUnder2,nDo,aUnder,aOther,work,
746               iBlock,jBlock,numberBlocks);
747    recRec(above,nUnder,nUnderK-nUnder2,nDo,aUnder+number_entries(nb),
748        aOther+number_entries(nb),work,iBlock,jBlock,numberBlocks);
749  } else if (nUnderK<=nDo&&nUnder<=nDo) {
750    int nb=number_blocks((nDo+1)>>1);
751    int nDo2=number_rows(nb);
752    int i;
753    recRec(above,nUnder,nUnderK,nDo2,aUnder,aOther,work,
754         iBlock,jBlock,numberBlocks);
755    i=((numberBlocks-jBlock)*(numberBlocks-jBlock-1)-
756      (numberBlocks-jBlock-nb)*(numberBlocks-jBlock-nb-1))>>1;
757    recRec(above+number_entries(i),nUnder,nUnderK,nDo-nDo2,
758           aUnder+number_entries(i),
759           aOther,work+nDo2,
760           iBlock-nb,jBlock,numberBlocks-nb);
761  } else {
762    int nb=number_blocks((nUnder+1)>>1);
763    int nUnder2=number_rows(nb);
764    int i;
765    recRec(above,nUnder2,nUnderK,nDo,aUnder,aOther,work,
766               iBlock,jBlock,numberBlocks);
767    i=((numberBlocks-iBlock)*(numberBlocks-iBlock-1)-
768      (numberBlocks-iBlock-nb)*(numberBlocks-iBlock-nb-1))>>1;
769    recRec(above+number_entries(nb),nUnder-nUnder2,nUnderK,nDo,aUnder,
770        aOther+number_entries(i),work,iBlock+nb,jBlock,numberBlocks);
771  }
772}
773// Leaf recursive factor
774void 
775ClpCholeskyDense::factorLeaf(longDouble * a, int n, 
776                longDouble * diagonal, longDouble * work, int * rowsDropped)
777{
778  CoinWorkDouble largest=doubleParameters_[3];
779  CoinWorkDouble smallest=doubleParameters_[4];
780  double dropValue = doubleParameters_[10];
781  int firstPositive=integerParameters_[34];
782  int rowOffset=diagonal-diagonal_;
783  int numberDropped=0;
784  int i, j, k;
785  CoinWorkDouble t00,temp1;
786  longDouble * aa;
787  aa = a-BLOCK;
788  for (j = 0; j < n; j ++) {
789    aa+=BLOCK;
790    t00 = aa[j];
791    for (k = 0; k < j; ++k) {
792      CoinWorkDouble multiplier = work[k];
793      t00 -= a[j + k * BLOCK] * a[j + k * BLOCK] * multiplier;
794    }
795    bool dropColumn=false;
796    CoinWorkDouble useT00=t00;
797    if (j+rowOffset<firstPositive) {
798      // must be negative
799      if (t00<=-dropValue) {
800        smallest = CoinMin(smallest,-t00);
801        largest = CoinMax(largest,-t00);
802        //aa[j]=t00;
803        t00 = 1.0/t00;
804      } else {
805        dropColumn=true;
806        //aa[j]=-1.0e100;
807        useT00=-1.0e-100;
808        t00=0.0;
809        integerParameters_[20]++;
810      }
811    } else {
812      // must be positive
813      if (t00>=dropValue) {
814        smallest = CoinMin(smallest,t00);
815        largest = CoinMax(largest,t00);
816        //aa[j]=t00;
817        t00 = 1.0/t00;
818      } else {
819        dropColumn=true;
820        //aa[j]=1.0e100;
821        useT00=1.0e-100;
822        t00=0.0;
823        integerParameters_[20]++;
824      }
825    }
826    if (!dropColumn) {
827      diagonal[j]=t00;
828      work[j]=useT00;
829      temp1 = t00;
830      for (i=j+1;i<n;i++) {
831        t00=aa[i];
832        for (k = 0; k < j; ++k) {
833          CoinWorkDouble multiplier = work[k];
834          t00 -= a[i + k * BLOCK] * a[j + k * BLOCK] * multiplier;
835        }
836        aa[i]=t00*temp1;
837      }
838    } else {
839      // drop column
840      rowsDropped[j+rowOffset]=2;
841      numberDropped++;
842      diagonal[j]=0.0;
843      //aa[j]=1.0e100;
844      work[j]=1.0e100;
845      for (i=j+1;i<n;i++) {
846        aa[i]=0.0;
847      }
848    }
849  }
850  doubleParameters_[3]=largest;
851  doubleParameters_[4]=smallest;
852  integerParameters_[20]+=numberDropped;
853}
854// Leaf recursive triangle rectangle update
855void 
856ClpCholeskyDense::triRecLeaf(longDouble * aTri, longDouble * aUnder, longDouble * diagonal, longDouble * work,
857                int nUnder)
858{
859#ifdef POS_DEBUG
860  int iru,icu;
861  int iu=bNumber(aUnder,iru,icu);
862  int irt,ict;
863  int it=bNumber(aTri,irt,ict);
864  //printf("%d %d\n",iu,it);
865  printf("trirecleaf  under (%d,%d), tri (%d,%d)\n",
866         iru,icu,irt,ict);
867  assert (diagonal==diagonal_+ict*BLOCK);
868#endif
869  int j;
870  longDouble * aa;
871#ifdef BLOCKUNROLL
872  if (nUnder==BLOCK) {
873    aa = aTri-2*BLOCK;
874    for (j = 0; j < BLOCK; j +=2) {
875      aa+=2*BLOCK;
876      CoinWorkDouble temp0 = diagonal[j];
877      CoinWorkDouble temp1 = diagonal[j+1];
878      for (int i=0;i<BLOCK;i+=2) {
879        CoinWorkDouble t00=aUnder[i+j*BLOCK];
880        CoinWorkDouble t10=aUnder[i+ BLOCK + j*BLOCK];
881        CoinWorkDouble t01=aUnder[i+1+j*BLOCK];
882        CoinWorkDouble t11=aUnder[i+1+ BLOCK + j*BLOCK];
883        for (int k = 0; k < j; ++k) {
884          CoinWorkDouble multiplier=work[k];
885          CoinWorkDouble au0 = aUnder[i + k * BLOCK]*multiplier;
886          CoinWorkDouble au1 = aUnder[i + 1 + k * BLOCK]*multiplier;
887          CoinWorkDouble at0 = aTri[j + k * BLOCK];
888          CoinWorkDouble at1 = aTri[j + 1 + k * BLOCK];
889          t00 -= au0 * at0;
890          t10 -= au0 * at1;
891          t01 -= au1 * at0;
892          t11 -= au1 * at1;
893        }
894        t00 *= temp0;
895        CoinWorkDouble at1 = aTri[j + 1 + j*BLOCK]*work[j];
896        t10 -= t00 * at1;
897        t01 *= temp0;
898        t11 -= t01 * at1;
899        aUnder[i+j*BLOCK]=t00;
900        aUnder[i+1+j*BLOCK]=t01;
901        aUnder[i+ BLOCK + j*BLOCK]=t10*temp1;
902        aUnder[i+1+ BLOCK + j*BLOCK]=t11*temp1;
903      }
904    }
905  } else {
906#endif
907    aa = aTri-BLOCK;
908    for (j = 0; j < BLOCK; j ++) {
909      aa+=BLOCK;
910      CoinWorkDouble temp1 = diagonal[j];
911      for (int i=0;i<nUnder;i++) {
912        CoinWorkDouble t00=aUnder[i+j*BLOCK];
913        for (int k = 0; k < j; ++k) {
914          CoinWorkDouble multiplier=work[k];
915          t00 -= aUnder[i + k * BLOCK] * aTri[j + k * BLOCK] * multiplier;
916        }
917        aUnder[i+j*BLOCK]=t00*temp1;
918      }
919    }
920#ifdef BLOCKUNROLL
921  }
922#endif
923}
924// Leaf recursive rectangle triangle update
925void ClpCholeskyDense::recTriLeaf(longDouble * aUnder, longDouble * aTri, 
926                                  longDouble * diagonal, longDouble * work, int nUnder)
927{
928#ifdef POS_DEBUG
929  int iru,icu;
930  int iu=bNumber(aUnder,iru,icu);
931  int irt,ict;
932  int it=bNumber(aTri,irt,ict);
933  //printf("%d %d\n",iu,it);
934  printf("rectrileaf  under (%d,%d), tri (%d,%d)\n",
935         iru,icu,irt,ict);
936  assert (diagonal==diagonal_+icu*BLOCK);
937#endif
938  int i, j, k;
939  CoinWorkDouble t00;
940  longDouble * aa;
941#ifdef BLOCKUNROLL
942  if (nUnder==BLOCK) {
943    longDouble * aUnder2 = aUnder-2;
944    aa = aTri-2*BLOCK;
945    for (j = 0; j < BLOCK; j +=2) {
946      CoinWorkDouble t00,t01,t10,t11;
947      aa+=2*BLOCK;
948      aUnder2+=2;
949      t00=aa[j];
950      t01=aa[j+1];
951      t10=aa[j+1+BLOCK];
952      for (k = 0; k < BLOCK; ++k) {
953        CoinWorkDouble multiplier = work[k];
954        CoinWorkDouble a0 = aUnder2[k * BLOCK];
955        CoinWorkDouble a1 = aUnder2[1 + k * BLOCK];
956        CoinWorkDouble x0 = a0 * multiplier;
957        CoinWorkDouble x1 = a1 * multiplier;
958        t00 -= a0 * x0;
959        t01 -= a1 * x0;
960        t10 -= a1 * x1;
961      }
962      aa[j]=t00;
963      aa[j+1]=t01;
964      aa[j+1+BLOCK]=t10;
965      for (i=j+2;i< BLOCK;i+=2) {
966        t00=aa[i];
967        t01=aa[i+BLOCK];
968        t10=aa[i+1];
969        t11=aa[i+1+BLOCK];
970        for (k = 0; k < BLOCK; ++k) {
971          CoinWorkDouble multiplier = work[k];
972          CoinWorkDouble a0 = aUnder2[k * BLOCK]*multiplier;
973          CoinWorkDouble a1 = aUnder2[1 + k * BLOCK]*multiplier;
974          t00 -= aUnder[i + k * BLOCK] * a0;
975          t01 -= aUnder[i + k * BLOCK] * a1;
976          t10 -= aUnder[i + 1 + k * BLOCK] * a0;
977          t11 -= aUnder[i + 1 + k * BLOCK] * a1;
978        }
979        aa[i]=t00;
980        aa[i+BLOCK]=t01;
981        aa[i+1]=t10;
982        aa[i+1+BLOCK]=t11;
983      }
984    }
985  } else {
986#endif
987    aa = aTri-BLOCK;
988    for (j = 0; j < nUnder; j ++) {
989      aa+=BLOCK;
990      for (i=j;i< nUnder;i++) {
991        t00=aa[i];
992        for (k = 0; k < BLOCK; ++k) {
993          CoinWorkDouble multiplier = work[k];
994          t00 -= aUnder[i + k * BLOCK] * aUnder[j + k * BLOCK]*multiplier;
995        }
996        aa[i]=t00;
997      }
998    }
999#ifdef BLOCKUNROLL
1000  }
1001#endif
1002}
1003/* Leaf recursive rectangle rectangle update,
1004   nUnder is number of rows in iBlock,
1005   nUnderK is number of rows in kBlock
1006*/
1007void 
1008ClpCholeskyDense::recRecLeaf(const longDouble * COIN_RESTRICT above, 
1009                             const longDouble * COIN_RESTRICT aUnder, 
1010                             longDouble * COIN_RESTRICT aOther, 
1011                             const longDouble * COIN_RESTRICT work,
1012                             int nUnder)
1013{
1014#ifdef POS_DEBUG
1015  int ira,ica;
1016  int ia=bNumber(above,ira,ica);
1017  int iru,icu;
1018  int iu=bNumber(aUnder,iru,icu);
1019  int iro,ico;
1020  int io=bNumber(aOther,iro,ico);
1021  //printf("%d %d %d\n",ia,iu,io);
1022  printf("recrecleaf above (%d,%d), under (%d,%d), other (%d,%d)\n",
1023         ira,ica,iru,icu,iro,ico);
1024#endif
1025  int i, j, k;
1026  longDouble * aa;
1027#ifdef BLOCKUNROLL
1028  aa = aOther-4*BLOCK;
1029  if (nUnder==BLOCK) {
1030    //#define INTEL
1031#ifdef INTEL
1032    aa+=2*BLOCK;
1033    for (j = 0; j < BLOCK; j +=2) {
1034      aa+=2*BLOCK;
1035      for (i=0;i< BLOCK;i+=2) {
1036        CoinWorkDouble t00=aa[i+0*BLOCK];
1037        CoinWorkDouble t10=aa[i+1*BLOCK];
1038        CoinWorkDouble t01=aa[i+1+0*BLOCK];
1039        CoinWorkDouble t11=aa[i+1+1*BLOCK];
1040        for (k=0;k<BLOCK;k++) {
1041          CoinWorkDouble multiplier = work[k];
1042          CoinWorkDouble a00=aUnder[i+k*BLOCK]*multiplier;
1043          CoinWorkDouble a01=aUnder[i+1+k*BLOCK]*multiplier;
1044          t00 -= a00 * above[j + 0 + k * BLOCK];
1045          t10 -= a00 * above[j + 1 + k * BLOCK];
1046          t01 -= a01 * above[j + 0 + k * BLOCK];
1047          t11 -= a01 * above[j + 1 + k * BLOCK];
1048        }
1049        aa[i+0*BLOCK]=t00;
1050        aa[i+1*BLOCK]=t10;
1051        aa[i+1+0*BLOCK]=t01;
1052        aa[i+1+1*BLOCK]=t11;
1053      }
1054    }
1055#else
1056    for (j = 0; j < BLOCK; j +=4) {
1057      aa+=4*BLOCK;
1058      for (i=0;i< BLOCK;i+=4) {
1059        CoinWorkDouble t00=aa[i+0+0*BLOCK];
1060        CoinWorkDouble t10=aa[i+0+1*BLOCK];
1061        CoinWorkDouble t20=aa[i+0+2*BLOCK];
1062        CoinWorkDouble t30=aa[i+0+3*BLOCK];
1063        CoinWorkDouble t01=aa[i+1+0*BLOCK];
1064        CoinWorkDouble t11=aa[i+1+1*BLOCK];
1065        CoinWorkDouble t21=aa[i+1+2*BLOCK];
1066        CoinWorkDouble t31=aa[i+1+3*BLOCK];
1067        CoinWorkDouble t02=aa[i+2+0*BLOCK];
1068        CoinWorkDouble t12=aa[i+2+1*BLOCK];
1069        CoinWorkDouble t22=aa[i+2+2*BLOCK];
1070        CoinWorkDouble t32=aa[i+2+3*BLOCK];
1071        CoinWorkDouble t03=aa[i+3+0*BLOCK];
1072        CoinWorkDouble t13=aa[i+3+1*BLOCK];
1073        CoinWorkDouble t23=aa[i+3+2*BLOCK];
1074        CoinWorkDouble t33=aa[i+3+3*BLOCK];
1075        const longDouble * COIN_RESTRICT aUnderNow = aUnder+i;
1076        const longDouble * COIN_RESTRICT aboveNow = above+j;
1077        for (k=0;k<BLOCK;k++) {
1078          CoinWorkDouble multiplier = work[k];
1079          CoinWorkDouble a00=aUnderNow[0]*multiplier;
1080          CoinWorkDouble a01=aUnderNow[1]*multiplier;
1081          CoinWorkDouble a02=aUnderNow[2]*multiplier;
1082          CoinWorkDouble a03=aUnderNow[3]*multiplier;
1083          t00 -= a00 * aboveNow[0];
1084          t10 -= a00 * aboveNow[1];
1085          t20 -= a00 * aboveNow[2];
1086          t30 -= a00 * aboveNow[3];
1087          t01 -= a01 * aboveNow[0];
1088          t11 -= a01 * aboveNow[1];
1089          t21 -= a01 * aboveNow[2];
1090          t31 -= a01 * aboveNow[3];
1091          t02 -= a02 * aboveNow[0];
1092          t12 -= a02 * aboveNow[1];
1093          t22 -= a02 * aboveNow[2];
1094          t32 -= a02 * aboveNow[3];
1095          t03 -= a03 * aboveNow[0];
1096          t13 -= a03 * aboveNow[1];
1097          t23 -= a03 * aboveNow[2];
1098          t33 -= a03 * aboveNow[3];
1099          aUnderNow += BLOCK;
1100          aboveNow += BLOCK;
1101        }
1102        aa[i+0+0*BLOCK]=t00;
1103        aa[i+0+1*BLOCK]=t10;
1104        aa[i+0+2*BLOCK]=t20;
1105        aa[i+0+3*BLOCK]=t30;
1106        aa[i+1+0*BLOCK]=t01;
1107        aa[i+1+1*BLOCK]=t11;
1108        aa[i+1+2*BLOCK]=t21;
1109        aa[i+1+3*BLOCK]=t31;
1110        aa[i+2+0*BLOCK]=t02;
1111        aa[i+2+1*BLOCK]=t12;
1112        aa[i+2+2*BLOCK]=t22;
1113        aa[i+2+3*BLOCK]=t32;
1114        aa[i+3+0*BLOCK]=t03;
1115        aa[i+3+1*BLOCK]=t13;
1116        aa[i+3+2*BLOCK]=t23;
1117        aa[i+3+3*BLOCK]=t33;
1118      }
1119    }
1120#endif
1121  } else {
1122    int odd=nUnder&1;
1123    int n=nUnder-odd;
1124    for (j = 0; j < BLOCK; j +=4) {
1125      aa+=4*BLOCK;
1126      for (i=0;i< n;i+=2) {
1127        CoinWorkDouble t00=aa[i+0*BLOCK];
1128        CoinWorkDouble t10=aa[i+1*BLOCK];
1129        CoinWorkDouble t20=aa[i+2*BLOCK];
1130        CoinWorkDouble t30=aa[i+3*BLOCK];
1131        CoinWorkDouble t01=aa[i+1+0*BLOCK];
1132        CoinWorkDouble t11=aa[i+1+1*BLOCK];
1133        CoinWorkDouble t21=aa[i+1+2*BLOCK];
1134        CoinWorkDouble t31=aa[i+1+3*BLOCK];
1135        const longDouble * COIN_RESTRICT aUnderNow = aUnder+i;
1136        const longDouble * COIN_RESTRICT aboveNow = above+j;
1137        for (k=0;k<BLOCK;k++) {
1138          CoinWorkDouble multiplier = work[k];
1139          CoinWorkDouble a00=aUnderNow[0]*multiplier;
1140          CoinWorkDouble a01=aUnderNow[1]*multiplier;
1141          t00 -= a00 * aboveNow[0];
1142          t10 -= a00 * aboveNow[1];
1143          t20 -= a00 * aboveNow[2];
1144          t30 -= a00 * aboveNow[3];
1145          t01 -= a01 * aboveNow[0];
1146          t11 -= a01 * aboveNow[1];
1147          t21 -= a01 * aboveNow[2];
1148          t31 -= a01 * aboveNow[3];
1149          aUnderNow += BLOCK;
1150          aboveNow += BLOCK;
1151        }
1152        aa[i+0*BLOCK]=t00;
1153        aa[i+1*BLOCK]=t10;
1154        aa[i+2*BLOCK]=t20;
1155        aa[i+3*BLOCK]=t30;
1156        aa[i+1+0*BLOCK]=t01;
1157        aa[i+1+1*BLOCK]=t11;
1158        aa[i+1+2*BLOCK]=t21;
1159        aa[i+1+3*BLOCK]=t31;
1160      }
1161      if (odd) {
1162        CoinWorkDouble t0=aa[n+0*BLOCK];
1163        CoinWorkDouble t1=aa[n+1*BLOCK];
1164        CoinWorkDouble t2=aa[n+2*BLOCK];
1165        CoinWorkDouble t3=aa[n+3*BLOCK];
1166        CoinWorkDouble a0;
1167        for (k=0;k<BLOCK;k++) {
1168          a0=aUnder[n+k*BLOCK]*work[k];
1169          t0 -= a0 * above[j + 0 + k * BLOCK];
1170          t1 -= a0 * above[j + 1 + k * BLOCK];
1171          t2 -= a0 * above[j + 2 + k * BLOCK];
1172          t3 -= a0 * above[j + 3 + k * BLOCK];
1173        }
1174        aa[n+0*BLOCK]=t0;
1175        aa[n+1*BLOCK]=t1;
1176        aa[n+2*BLOCK]=t2;
1177        aa[n+3*BLOCK]=t3;
1178      }
1179    }
1180  }
1181#else
1182  aa = aOther-BLOCK;
1183  for (j = 0; j < BLOCK; j ++) {
1184    aa+=BLOCK;
1185    for (i=0;i< nUnder;i++) {
1186      CoinWorkDouble t00=aa[i+0*BLOCK];
1187      for (k=0;k<BLOCK;k++) {
1188        CoinWorkDouble a00=aUnder[i+k*BLOCK]*work[k];
1189        t00 -= a00 * above[j + k * BLOCK];
1190      }
1191      aa[i]=t00;
1192    }
1193  }
1194#endif
1195}
1196/* Uses factorization to solve. */
1197void 
1198ClpCholeskyDense::solve (double * region) 
1199{
1200#ifdef CHOL_COMPARE
1201  double * region2=NULL;
1202  if (numberRows_<200) {
1203    longDouble * xx = sparseFactor_;
1204    longDouble * yy = diagonal_;
1205    diagonal_ = sparseFactor_ + 40000;
1206    sparseFactor_=diagonal_ + numberRows_;
1207    region2 = ClpCopyOfArray(region,numberRows_);
1208    int iRow,iColumn;
1209    int addOffset=numberRows_-1;
1210    longDouble * work = sparseFactor_-1;
1211    for (iColumn=0;iColumn<numberRows_;iColumn++) {
1212      double value = region2[iColumn];
1213      for (iRow=iColumn+1;iRow<numberRows_;iRow++)
1214        region2[iRow] -= value*work[iRow];
1215      addOffset--;
1216      work += addOffset;
1217    }
1218    for (iColumn=numberRows_-1;iColumn>=0;iColumn--) {
1219      double value = region2[iColumn]*diagonal_[iColumn];
1220      work -= addOffset;
1221      addOffset++;
1222      for (iRow=iColumn+1;iRow<numberRows_;iRow++)
1223        value -= region2[iRow]*work[iRow];
1224      region2[iColumn]=value;
1225    }
1226    sparseFactor_=xx;
1227    diagonal_=yy;
1228  }
1229#endif
1230  //longDouble * xx = sparseFactor_;
1231  //longDouble * yy = diagonal_;
1232  //diagonal_ = sparseFactor_ + 40000;
1233  //sparseFactor_=diagonal_ + numberRows_;
1234  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
1235  // later align on boundary
1236  longDouble * a = sparseFactor_+BLOCKSQ*numberBlocks;
1237  int iBlock;
1238  longDouble * aa = a;
1239  int iColumn;
1240  for (iBlock=0;iBlock<numberBlocks;iBlock++) {
1241    int nChunk;
1242    int jBlock;
1243    int iDo = iBlock*BLOCK;
1244    int base=iDo;
1245    if (iDo+BLOCK>numberRows_) {
1246      nChunk=numberRows_-iDo;
1247    } else {
1248      nChunk=BLOCK;
1249    }
1250    solveF1(aa,nChunk,region+iDo);
1251    for (jBlock=iBlock+1;jBlock<numberBlocks;jBlock++) {
1252      base+=BLOCK;
1253      aa+=BLOCKSQ;
1254      if (base+BLOCK>numberRows_) {
1255        nChunk=numberRows_-base;
1256      } else {
1257        nChunk=BLOCK;
1258      }
1259      solveF2(aa,nChunk,region+iDo,region+base);
1260    }
1261    aa+=BLOCKSQ;
1262  }
1263  // do diagonal outside
1264  for (iColumn=0;iColumn<numberRows_;iColumn++) 
1265    region[iColumn] *= diagonal_[iColumn];
1266  int offset=((numberBlocks*(numberBlocks+1))>>1);
1267  aa = a+number_entries(offset-1);
1268  int lBase=(numberBlocks-1)*BLOCK;
1269  for (iBlock=numberBlocks-1;iBlock>=0;iBlock--) {
1270    int nChunk;
1271    int jBlock;
1272    int triBase=iBlock*BLOCK;
1273    int iBase=lBase;
1274    for (jBlock=iBlock+1;jBlock<numberBlocks;jBlock++) {
1275      if (iBase+BLOCK>numberRows_) {
1276        nChunk=numberRows_-iBase;
1277      } else {
1278        nChunk=BLOCK;
1279      }
1280      solveB2(aa,nChunk,region+triBase,region+iBase);
1281      iBase-=BLOCK;
1282      aa-=BLOCKSQ;
1283    }
1284    if (triBase+BLOCK>numberRows_) {
1285      nChunk=numberRows_-triBase;
1286    } else {
1287      nChunk=BLOCK;
1288    }
1289    solveB1(aa,nChunk,region+triBase);
1290    aa-=BLOCKSQ;
1291  }
1292#ifdef CHOL_COMPARE
1293  if (numberRows_<200) {
1294    for (int i=0;i<numberRows_;i++) {
1295      assert(CoinAbs(region[i]-region2[i])<1.0e-3);
1296    }
1297    delete [] region2;
1298  }
1299#endif
1300}
1301// Forward part of solve 1
1302void 
1303ClpCholeskyDense::solveF1(longDouble * a,int n,double * region)
1304{
1305  int j, k;
1306  CoinWorkDouble t00;
1307  for (j = 0; j < n; j ++) {
1308    t00 = region[j];
1309    for (k = 0; k < j; ++k) {
1310      t00 -= region[k] * a[j + k * BLOCK];
1311    }
1312    //t00*=a[j + j * BLOCK];
1313    region[j] = t00;
1314  }
1315}
1316// Forward part of solve 2
1317void 
1318ClpCholeskyDense::solveF2(longDouble * a,int n,double * region, double * region2)
1319{
1320  int j, k;
1321#ifdef BLOCKUNROLL
1322  if (n==BLOCK) {
1323    for (k = 0; k < BLOCK; k+=4) {
1324      CoinWorkDouble t0 = region2[0];
1325      CoinWorkDouble t1 = region2[1];
1326      CoinWorkDouble t2 = region2[2];
1327      CoinWorkDouble t3 = region2[3];
1328      t0 -= region[0] * a[0 + 0 * BLOCK];
1329      t1 -= region[0] * a[1 + 0 * BLOCK];
1330      t2 -= region[0] * a[2 + 0 * BLOCK];
1331      t3 -= region[0] * a[3 + 0 * BLOCK];
1332
1333      t0 -= region[1] * a[0 + 1 * BLOCK];
1334      t1 -= region[1] * a[1 + 1 * BLOCK];
1335      t2 -= region[1] * a[2 + 1 * BLOCK];
1336      t3 -= region[1] * a[3 + 1 * BLOCK];
1337
1338      t0 -= region[2] * a[0 + 2 * BLOCK];
1339      t1 -= region[2] * a[1 + 2 * BLOCK];
1340      t2 -= region[2] * a[2 + 2 * BLOCK];
1341      t3 -= region[2] * a[3 + 2 * BLOCK];
1342
1343      t0 -= region[3] * a[0 + 3 * BLOCK];
1344      t1 -= region[3] * a[1 + 3 * BLOCK];
1345      t2 -= region[3] * a[2 + 3 * BLOCK];
1346      t3 -= region[3] * a[3 + 3 * BLOCK];
1347
1348      t0 -= region[4] * a[0 + 4 * BLOCK];
1349      t1 -= region[4] * a[1 + 4 * BLOCK];
1350      t2 -= region[4] * a[2 + 4 * BLOCK];
1351      t3 -= region[4] * a[3 + 4 * BLOCK];
1352
1353      t0 -= region[5] * a[0 + 5 * BLOCK];
1354      t1 -= region[5] * a[1 + 5 * BLOCK];
1355      t2 -= region[5] * a[2 + 5 * BLOCK];
1356      t3 -= region[5] * a[3 + 5 * BLOCK];
1357
1358      t0 -= region[6] * a[0 + 6 * BLOCK];
1359      t1 -= region[6] * a[1 + 6 * BLOCK];
1360      t2 -= region[6] * a[2 + 6 * BLOCK];
1361      t3 -= region[6] * a[3 + 6 * BLOCK];
1362
1363      t0 -= region[7] * a[0 + 7 * BLOCK];
1364      t1 -= region[7] * a[1 + 7 * BLOCK];
1365      t2 -= region[7] * a[2 + 7 * BLOCK];
1366      t3 -= region[7] * a[3 + 7 * BLOCK];
1367#if BLOCK>8
1368      t0 -= region[8] * a[0 + 8 * BLOCK];
1369      t1 -= region[8] * a[1 + 8 * BLOCK];
1370      t2 -= region[8] * a[2 + 8 * BLOCK];
1371      t3 -= region[8] * a[3 + 8 * BLOCK];
1372
1373      t0 -= region[9] * a[0 + 9 * BLOCK];
1374      t1 -= region[9] * a[1 + 9 * BLOCK];
1375      t2 -= region[9] * a[2 + 9 * BLOCK];
1376      t3 -= region[9] * a[3 + 9 * BLOCK];
1377
1378      t0 -= region[10] * a[0 + 10 * BLOCK];
1379      t1 -= region[10] * a[1 + 10 * BLOCK];
1380      t2 -= region[10] * a[2 + 10 * BLOCK];
1381      t3 -= region[10] * a[3 + 10 * BLOCK];
1382
1383      t0 -= region[11] * a[0 + 11 * BLOCK];
1384      t1 -= region[11] * a[1 + 11 * BLOCK];
1385      t2 -= region[11] * a[2 + 11 * BLOCK];
1386      t3 -= region[11] * a[3 + 11 * BLOCK];
1387
1388      t0 -= region[12] * a[0 + 12 * BLOCK];
1389      t1 -= region[12] * a[1 + 12 * BLOCK];
1390      t2 -= region[12] * a[2 + 12 * BLOCK];
1391      t3 -= region[12] * a[3 + 12 * BLOCK];
1392
1393      t0 -= region[13] * a[0 + 13 * BLOCK];
1394      t1 -= region[13] * a[1 + 13 * BLOCK];
1395      t2 -= region[13] * a[2 + 13 * BLOCK];
1396      t3 -= region[13] * a[3 + 13 * BLOCK];
1397
1398      t0 -= region[14] * a[0 + 14 * BLOCK];
1399      t1 -= region[14] * a[1 + 14 * BLOCK];
1400      t2 -= region[14] * a[2 + 14 * BLOCK];
1401      t3 -= region[14] * a[3 + 14 * BLOCK];
1402
1403      t0 -= region[15] * a[0 + 15 * BLOCK];
1404      t1 -= region[15] * a[1 + 15 * BLOCK];
1405      t2 -= region[15] * a[2 + 15 * BLOCK];
1406      t3 -= region[15] * a[3 + 15 * BLOCK];
1407#endif
1408      region2[0] = t0;
1409      region2[1] = t1;
1410      region2[2] = t2;
1411      region2[3] = t3;
1412      region2+=4;
1413      a+=4;
1414    }
1415  } else {
1416#endif
1417    for (k = 0; k < n; ++k) {
1418      CoinWorkDouble t00 = region2[k];
1419      for (j = 0; j < BLOCK; j ++) {
1420        t00 -= region[j] * a[k + j * BLOCK];
1421      }
1422      region2[k] = t00;
1423    }
1424#ifdef BLOCKUNROLL
1425  }
1426#endif
1427}
1428// Backward part of solve 1
1429void 
1430ClpCholeskyDense::solveB1(longDouble * a,int n,double * region)
1431{
1432  int j, k;
1433  CoinWorkDouble t00;
1434  for (j = n-1; j >=0; j --) {
1435    t00 = region[j];
1436    for (k = j+1; k < n; ++k) {
1437      t00 -= region[k] * a[k + j * BLOCK];
1438    }
1439    //t00*=a[j + j * BLOCK];
1440    region[j] = t00;
1441  }
1442}
1443// Backward part of solve 2
1444void 
1445ClpCholeskyDense::solveB2(longDouble * a,int n,double * region, double * region2)
1446{
1447  int j, k;
1448#ifdef BLOCKUNROLL
1449  if (n==BLOCK) {
1450    for (j = 0; j < BLOCK; j +=4) {
1451      CoinWorkDouble t0 = region[0];
1452      CoinWorkDouble t1 = region[1];
1453      CoinWorkDouble t2 = region[2];
1454      CoinWorkDouble t3 = region[3];
1455      t0 -= region2[0] * a[0 + 0*BLOCK];
1456      t1 -= region2[0] * a[0 + 1*BLOCK];
1457      t2 -= region2[0] * a[0 + 2*BLOCK];
1458      t3 -= region2[0] * a[0 + 3*BLOCK];
1459
1460      t0 -= region2[1] * a[1 + 0*BLOCK];
1461      t1 -= region2[1] * a[1 + 1*BLOCK];
1462      t2 -= region2[1] * a[1 + 2*BLOCK];
1463      t3 -= region2[1] * a[1 + 3*BLOCK];
1464
1465      t0 -= region2[2] * a[2 + 0*BLOCK];
1466      t1 -= region2[2] * a[2 + 1*BLOCK];
1467      t2 -= region2[2] * a[2 + 2*BLOCK];
1468      t3 -= region2[2] * a[2 + 3*BLOCK];
1469
1470      t0 -= region2[3] * a[3 + 0*BLOCK];
1471      t1 -= region2[3] * a[3 + 1*BLOCK];
1472      t2 -= region2[3] * a[3 + 2*BLOCK];
1473      t3 -= region2[3] * a[3 + 3*BLOCK];
1474
1475      t0 -= region2[4] * a[4 + 0*BLOCK];
1476      t1 -= region2[4] * a[4 + 1*BLOCK];
1477      t2 -= region2[4] * a[4 + 2*BLOCK];
1478      t3 -= region2[4] * a[4 + 3*BLOCK];
1479
1480      t0 -= region2[5] * a[5 + 0*BLOCK];
1481      t1 -= region2[5] * a[5 + 1*BLOCK];
1482      t2 -= region2[5] * a[5 + 2*BLOCK];
1483      t3 -= region2[5] * a[5 + 3*BLOCK];
1484
1485      t0 -= region2[6] * a[6 + 0*BLOCK];
1486      t1 -= region2[6] * a[6 + 1*BLOCK];
1487      t2 -= region2[6] * a[6 + 2*BLOCK];
1488      t3 -= region2[6] * a[6 + 3*BLOCK];
1489
1490      t0 -= region2[7] * a[7 + 0*BLOCK];
1491      t1 -= region2[7] * a[7 + 1*BLOCK];
1492      t2 -= region2[7] * a[7 + 2*BLOCK];
1493      t3 -= region2[7] * a[7 + 3*BLOCK];
1494#if BLOCK>8
1495
1496      t0 -= region2[8] * a[8 + 0*BLOCK];
1497      t1 -= region2[8] * a[8 + 1*BLOCK];
1498      t2 -= region2[8] * a[8 + 2*BLOCK];
1499      t3 -= region2[8] * a[8 + 3*BLOCK];
1500
1501      t0 -= region2[9] * a[9 + 0*BLOCK];
1502      t1 -= region2[9] * a[9 + 1*BLOCK];
1503      t2 -= region2[9] * a[9 + 2*BLOCK];
1504      t3 -= region2[9] * a[9 + 3*BLOCK];
1505
1506      t0 -= region2[10] * a[10 + 0*BLOCK];
1507      t1 -= region2[10] * a[10 + 1*BLOCK];
1508      t2 -= region2[10] * a[10 + 2*BLOCK];
1509      t3 -= region2[10] * a[10 + 3*BLOCK];
1510
1511      t0 -= region2[11] * a[11 + 0*BLOCK];
1512      t1 -= region2[11] * a[11 + 1*BLOCK];
1513      t2 -= region2[11] * a[11 + 2*BLOCK];
1514      t3 -= region2[11] * a[11 + 3*BLOCK];
1515
1516      t0 -= region2[12] * a[12 + 0*BLOCK];
1517      t1 -= region2[12] * a[12 + 1*BLOCK];
1518      t2 -= region2[12] * a[12 + 2*BLOCK];
1519      t3 -= region2[12] * a[12 + 3*BLOCK];
1520
1521      t0 -= region2[13] * a[13 + 0*BLOCK];
1522      t1 -= region2[13] * a[13 + 1*BLOCK];
1523      t2 -= region2[13] * a[13 + 2*BLOCK];
1524      t3 -= region2[13] * a[13 + 3*BLOCK];
1525
1526      t0 -= region2[14] * a[14 + 0*BLOCK];
1527      t1 -= region2[14] * a[14 + 1*BLOCK];
1528      t2 -= region2[14] * a[14 + 2*BLOCK];
1529      t3 -= region2[14] * a[14 + 3*BLOCK];
1530
1531      t0 -= region2[15] * a[15 + 0*BLOCK];
1532      t1 -= region2[15] * a[15 + 1*BLOCK];
1533      t2 -= region2[15] * a[15 + 2*BLOCK];
1534      t3 -= region2[15] * a[15 + 3*BLOCK];
1535#endif
1536      region[0] = t0;
1537      region[1] = t1;
1538      region[2] = t2;
1539      region[3] = t3;
1540      a+=4*BLOCK;
1541      region+=4;
1542    }
1543  } else {
1544#endif
1545    for (j = 0; j < BLOCK; j ++) {
1546      CoinWorkDouble t00 = region[j];
1547      for (k = 0; k < n; ++k) {
1548        t00 -= region2[k] * a[k + j * BLOCK];
1549      }
1550      region[j] = t00;
1551    }
1552#ifdef BLOCKUNROLL
1553  }
1554#endif
1555}
1556/* Uses factorization to solve. */
1557void 
1558ClpCholeskyDense::solveLong (CoinWorkDouble * region) 
1559{
1560  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
1561  // later align on boundary
1562  longDouble * a = sparseFactor_+BLOCKSQ*numberBlocks;
1563  int iBlock;
1564  longDouble * aa = a;
1565  int iColumn;
1566  for (iBlock=0;iBlock<numberBlocks;iBlock++) {
1567    int nChunk;
1568    int jBlock;
1569    int iDo = iBlock*BLOCK;
1570    int base=iDo;
1571    if (iDo+BLOCK>numberRows_) {
1572      nChunk=numberRows_-iDo;
1573    } else {
1574      nChunk=BLOCK;
1575    }
1576    solveF1Long(aa,nChunk,region+iDo);
1577    for (jBlock=iBlock+1;jBlock<numberBlocks;jBlock++) {
1578      base+=BLOCK;
1579      aa+=BLOCKSQ;
1580      if (base+BLOCK>numberRows_) {
1581        nChunk=numberRows_-base;
1582      } else {
1583        nChunk=BLOCK;
1584      }
1585      solveF2Long(aa,nChunk,region+iDo,region+base);
1586    }
1587    aa+=BLOCKSQ;
1588  }
1589  // do diagonal outside
1590  for (iColumn=0;iColumn<numberRows_;iColumn++) 
1591    region[iColumn] *= diagonal_[iColumn];
1592  int offset=((numberBlocks*(numberBlocks+1))>>1);
1593  aa = a+number_entries(offset-1);
1594  int lBase=(numberBlocks-1)*BLOCK;
1595  for (iBlock=numberBlocks-1;iBlock>=0;iBlock--) {
1596    int nChunk;
1597    int jBlock;
1598    int triBase=iBlock*BLOCK;
1599    int iBase=lBase;
1600    for (jBlock=iBlock+1;jBlock<numberBlocks;jBlock++) {
1601      if (iBase+BLOCK>numberRows_) {
1602        nChunk=numberRows_-iBase;
1603      } else {
1604        nChunk=BLOCK;
1605      }
1606      solveB2Long(aa,nChunk,region+triBase,region+iBase);
1607      iBase-=BLOCK;
1608      aa-=BLOCKSQ;
1609    }
1610    if (triBase+BLOCK>numberRows_) {
1611      nChunk=numberRows_-triBase;
1612    } else {
1613      nChunk=BLOCK;
1614    }
1615    solveB1Long(aa,nChunk,region+triBase);
1616    aa-=BLOCKSQ;
1617  }
1618}
1619// Forward part of solve 1
1620void 
1621ClpCholeskyDense::solveF1Long(longDouble * a,int n,CoinWorkDouble * region)
1622{
1623  int j, k;
1624  CoinWorkDouble t00;
1625  for (j = 0; j < n; j ++) {
1626    t00 = region[j];
1627    for (k = 0; k < j; ++k) {
1628      t00 -= region[k] * a[j + k * BLOCK];
1629    }
1630    //t00*=a[j + j * BLOCK];
1631    region[j] = t00;
1632  }
1633}
1634// Forward part of solve 2
1635void 
1636ClpCholeskyDense::solveF2Long(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2)
1637{
1638  int j, k;
1639#ifdef BLOCKUNROLL
1640  if (n==BLOCK) {
1641    for (k = 0; k < BLOCK; k+=4) {
1642      CoinWorkDouble t0 = region2[0];
1643      CoinWorkDouble t1 = region2[1];
1644      CoinWorkDouble t2 = region2[2];
1645      CoinWorkDouble t3 = region2[3];
1646      t0 -= region[0] * a[0 + 0 * BLOCK];
1647      t1 -= region[0] * a[1 + 0 * BLOCK];
1648      t2 -= region[0] * a[2 + 0 * BLOCK];
1649      t3 -= region[0] * a[3 + 0 * BLOCK];
1650
1651      t0 -= region[1] * a[0 + 1 * BLOCK];
1652      t1 -= region[1] * a[1 + 1 * BLOCK];
1653      t2 -= region[1] * a[2 + 1 * BLOCK];
1654      t3 -= region[1] * a[3 + 1 * BLOCK];
1655
1656      t0 -= region[2] * a[0 + 2 * BLOCK];
1657      t1 -= region[2] * a[1 + 2 * BLOCK];
1658      t2 -= region[2] * a[2 + 2 * BLOCK];
1659      t3 -= region[2] * a[3 + 2 * BLOCK];
1660
1661      t0 -= region[3] * a[0 + 3 * BLOCK];
1662      t1 -= region[3] * a[1 + 3 * BLOCK];
1663      t2 -= region[3] * a[2 + 3 * BLOCK];
1664      t3 -= region[3] * a[3 + 3 * BLOCK];
1665
1666      t0 -= region[4] * a[0 + 4 * BLOCK];
1667      t1 -= region[4] * a[1 + 4 * BLOCK];
1668      t2 -= region[4] * a[2 + 4 * BLOCK];
1669      t3 -= region[4] * a[3 + 4 * BLOCK];
1670
1671      t0 -= region[5] * a[0 + 5 * BLOCK];
1672      t1 -= region[5] * a[1 + 5 * BLOCK];
1673      t2 -= region[5] * a[2 + 5 * BLOCK];
1674      t3 -= region[5] * a[3 + 5 * BLOCK];
1675
1676      t0 -= region[6] * a[0 + 6 * BLOCK];
1677      t1 -= region[6] * a[1 + 6 * BLOCK];
1678      t2 -= region[6] * a[2 + 6 * BLOCK];
1679      t3 -= region[6] * a[3 + 6 * BLOCK];
1680
1681      t0 -= region[7] * a[0 + 7 * BLOCK];
1682      t1 -= region[7] * a[1 + 7 * BLOCK];
1683      t2 -= region[7] * a[2 + 7 * BLOCK];
1684      t3 -= region[7] * a[3 + 7 * BLOCK];
1685#if BLOCK>8
1686      t0 -= region[8] * a[0 + 8 * BLOCK];
1687      t1 -= region[8] * a[1 + 8 * BLOCK];
1688      t2 -= region[8] * a[2 + 8 * BLOCK];
1689      t3 -= region[8] * a[3 + 8 * BLOCK];
1690
1691      t0 -= region[9] * a[0 + 9 * BLOCK];
1692      t1 -= region[9] * a[1 + 9 * BLOCK];
1693      t2 -= region[9] * a[2 + 9 * BLOCK];
1694      t3 -= region[9] * a[3 + 9 * BLOCK];
1695
1696      t0 -= region[10] * a[0 + 10 * BLOCK];
1697      t1 -= region[10] * a[1 + 10 * BLOCK];
1698      t2 -= region[10] * a[2 + 10 * BLOCK];
1699      t3 -= region[10] * a[3 + 10 * BLOCK];
1700
1701      t0 -= region[11] * a[0 + 11 * BLOCK];
1702      t1 -= region[11] * a[1 + 11 * BLOCK];
1703      t2 -= region[11] * a[2 + 11 * BLOCK];
1704      t3 -= region[11] * a[3 + 11 * BLOCK];
1705
1706      t0 -= region[12] * a[0 + 12 * BLOCK];
1707      t1 -= region[12] * a[1 + 12 * BLOCK];
1708      t2 -= region[12] * a[2 + 12 * BLOCK];
1709      t3 -= region[12] * a[3 + 12 * BLOCK];
1710
1711      t0 -= region[13] * a[0 + 13 * BLOCK];
1712      t1 -= region[13] * a[1 + 13 * BLOCK];
1713      t2 -= region[13] * a[2 + 13 * BLOCK];
1714      t3 -= region[13] * a[3 + 13 * BLOCK];
1715
1716      t0 -= region[14] * a[0 + 14 * BLOCK];
1717      t1 -= region[14] * a[1 + 14 * BLOCK];
1718      t2 -= region[14] * a[2 + 14 * BLOCK];
1719      t3 -= region[14] * a[3 + 14 * BLOCK];
1720
1721      t0 -= region[15] * a[0 + 15 * BLOCK];
1722      t1 -= region[15] * a[1 + 15 * BLOCK];
1723      t2 -= region[15] * a[2 + 15 * BLOCK];
1724      t3 -= region[15] * a[3 + 15 * BLOCK];
1725#endif
1726      region2[0] = t0;
1727      region2[1] = t1;
1728      region2[2] = t2;
1729      region2[3] = t3;
1730      region2+=4;
1731      a+=4;
1732    }
1733  } else {
1734#endif
1735    for (k = 0; k < n; ++k) {
1736      CoinWorkDouble t00 = region2[k];
1737      for (j = 0; j < BLOCK; j ++) {
1738        t00 -= region[j] * a[k + j * BLOCK];
1739      }
1740      region2[k] = t00;
1741    }
1742#ifdef BLOCKUNROLL
1743  }
1744#endif
1745}
1746// Backward part of solve 1
1747void 
1748ClpCholeskyDense::solveB1Long(longDouble * a,int n,CoinWorkDouble * region)
1749{
1750  int j, k;
1751  CoinWorkDouble t00;
1752  for (j = n-1; j >=0; j --) {
1753    t00 = region[j];
1754    for (k = j+1; k < n; ++k) {
1755      t00 -= region[k] * a[k + j * BLOCK];
1756    }
1757    //t00*=a[j + j * BLOCK];
1758    region[j] = t00;
1759  }
1760}
1761// Backward part of solve 2
1762void 
1763ClpCholeskyDense::solveB2Long(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2)
1764{
1765  int j, k;
1766#ifdef BLOCKUNROLL
1767  if (n==BLOCK) {
1768    for (j = 0; j < BLOCK; j +=4) {
1769      CoinWorkDouble t0 = region[0];
1770      CoinWorkDouble t1 = region[1];
1771      CoinWorkDouble t2 = region[2];
1772      CoinWorkDouble t3 = region[3];
1773      t0 -= region2[0] * a[0 + 0*BLOCK];
1774      t1 -= region2[0] * a[0 + 1*BLOCK];
1775      t2 -= region2[0] * a[0 + 2*BLOCK];
1776      t3 -= region2[0] * a[0 + 3*BLOCK];
1777
1778      t0 -= region2[1] * a[1 + 0*BLOCK];
1779      t1 -= region2[1] * a[1 + 1*BLOCK];
1780      t2 -= region2[1] * a[1 + 2*BLOCK];
1781      t3 -= region2[1] * a[1 + 3*BLOCK];
1782
1783      t0 -= region2[2] * a[2 + 0*BLOCK];
1784      t1 -= region2[2] * a[2 + 1*BLOCK];
1785      t2 -= region2[2] * a[2 + 2*BLOCK];
1786      t3 -= region2[2] * a[2 + 3*BLOCK];
1787
1788      t0 -= region2[3] * a[3 + 0*BLOCK];
1789      t1 -= region2[3] * a[3 + 1*BLOCK];
1790      t2 -= region2[3] * a[3 + 2*BLOCK];
1791      t3 -= region2[3] * a[3 + 3*BLOCK];
1792
1793      t0 -= region2[4] * a[4 + 0*BLOCK];
1794      t1 -= region2[4] * a[4 + 1*BLOCK];
1795      t2 -= region2[4] * a[4 + 2*BLOCK];
1796      t3 -= region2[4] * a[4 + 3*BLOCK];
1797
1798      t0 -= region2[5] * a[5 + 0*BLOCK];
1799      t1 -= region2[5] * a[5 + 1*BLOCK];
1800      t2 -= region2[5] * a[5 + 2*BLOCK];
1801      t3 -= region2[5] * a[5 + 3*BLOCK];
1802
1803      t0 -= region2[6] * a[6 + 0*BLOCK];
1804      t1 -= region2[6] * a[6 + 1*BLOCK];
1805      t2 -= region2[6] * a[6 + 2*BLOCK];
1806      t3 -= region2[6] * a[6 + 3*BLOCK];
1807
1808      t0 -= region2[7] * a[7 + 0*BLOCK];
1809      t1 -= region2[7] * a[7 + 1*BLOCK];
1810      t2 -= region2[7] * a[7 + 2*BLOCK];
1811      t3 -= region2[7] * a[7 + 3*BLOCK];
1812#if BLOCK>8
1813
1814      t0 -= region2[8] * a[8 + 0*BLOCK];
1815      t1 -= region2[8] * a[8 + 1*BLOCK];
1816      t2 -= region2[8] * a[8 + 2*BLOCK];
1817      t3 -= region2[8] * a[8 + 3*BLOCK];
1818
1819      t0 -= region2[9] * a[9 + 0*BLOCK];
1820      t1 -= region2[9] * a[9 + 1*BLOCK];
1821      t2 -= region2[9] * a[9 + 2*BLOCK];
1822      t3 -= region2[9] * a[9 + 3*BLOCK];
1823
1824      t0 -= region2[10] * a[10 + 0*BLOCK];
1825      t1 -= region2[10] * a[10 + 1*BLOCK];
1826      t2 -= region2[10] * a[10 + 2*BLOCK];
1827      t3 -= region2[10] * a[10 + 3*BLOCK];
1828
1829      t0 -= region2[11] * a[11 + 0*BLOCK];
1830      t1 -= region2[11] * a[11 + 1*BLOCK];
1831      t2 -= region2[11] * a[11 + 2*BLOCK];
1832      t3 -= region2[11] * a[11 + 3*BLOCK];
1833
1834      t0 -= region2[12] * a[12 + 0*BLOCK];
1835      t1 -= region2[12] * a[12 + 1*BLOCK];
1836      t2 -= region2[12] * a[12 + 2*BLOCK];
1837      t3 -= region2[12] * a[12 + 3*BLOCK];
1838
1839      t0 -= region2[13] * a[13 + 0*BLOCK];
1840      t1 -= region2[13] * a[13 + 1*BLOCK];
1841      t2 -= region2[13] * a[13 + 2*BLOCK];
1842      t3 -= region2[13] * a[13 + 3*BLOCK];
1843
1844      t0 -= region2[14] * a[14 + 0*BLOCK];
1845      t1 -= region2[14] * a[14 + 1*BLOCK];
1846      t2 -= region2[14] * a[14 + 2*BLOCK];
1847      t3 -= region2[14] * a[14 + 3*BLOCK];
1848
1849      t0 -= region2[15] * a[15 + 0*BLOCK];
1850      t1 -= region2[15] * a[15 + 1*BLOCK];
1851      t2 -= region2[15] * a[15 + 2*BLOCK];
1852      t3 -= region2[15] * a[15 + 3*BLOCK];
1853#endif
1854      region[0] = t0;
1855      region[1] = t1;
1856      region[2] = t2;
1857      region[3] = t3;
1858      a+=4*BLOCK;
1859      region+=4;
1860    }
1861  } else {
1862#endif
1863    for (j = 0; j < BLOCK; j ++) {
1864      CoinWorkDouble t00 = region[j];
1865      for (k = 0; k < n; ++k) {
1866        t00 -= region2[k] * a[k + j * BLOCK];
1867      }
1868      region[j] = t00;
1869    }
1870#ifdef BLOCKUNROLL
1871  }
1872#endif
1873}
1874/* Uses factorization to solve. */
1875void 
1876ClpCholeskyDense::solveLongWork (CoinWorkDouble * region) 
1877{
1878  int numberBlocks = (numberRows_+BLOCK-1)>>BLOCKSHIFT;
1879  // later align on boundary
1880  longDouble * a = sparseFactor_+BLOCKSQ*numberBlocks;
1881  int iBlock;
1882  longDouble * aa = a;
1883  int iColumn;
1884  for (iBlock=0;iBlock<numberBlocks;iBlock++) {
1885    int nChunk;
1886    int jBlock;
1887    int iDo = iBlock*BLOCK;
1888    int base=iDo;
1889    if (iDo+BLOCK>numberRows_) {
1890      nChunk=numberRows_-iDo;
1891    } else {
1892      nChunk=BLOCK;
1893    }
1894    solveF1LongWork(aa,nChunk,region+iDo);
1895    for (jBlock=iBlock+1;jBlock<numberBlocks;jBlock++) {
1896      base+=BLOCK;
1897      aa+=BLOCKSQ;
1898      if (base+BLOCK>numberRows_) {
1899        nChunk=numberRows_-base;
1900      } else {
1901        nChunk=BLOCK;
1902      }
1903      solveF2LongWork(aa,nChunk,region+iDo,region+base);
1904    }
1905    aa+=BLOCKSQ;
1906  }
1907  // do diagonal outside
1908  for (iColumn=0;iColumn<numberRows_;iColumn++) 
1909    region[iColumn] *= diagonal_[iColumn];
1910  int offset=((numberBlocks*(numberBlocks+1))>>1);
1911  aa = a+number_entries(offset-1);
1912  int lBase=(numberBlocks-1)*BLOCK;
1913  for (iBlock=numberBlocks-1;iBlock>=0;iBlock--) {
1914    int nChunk;
1915    int jBlock;
1916    int triBase=iBlock*BLOCK;
1917    int iBase=lBase;
1918    for (jBlock=iBlock+1;jBlock<numberBlocks;jBlock++) {
1919      if (iBase+BLOCK>numberRows_) {
1920        nChunk=numberRows_-iBase;
1921      } else {
1922        nChunk=BLOCK;
1923      }
1924      solveB2LongWork(aa,nChunk,region+triBase,region+iBase);
1925      iBase-=BLOCK;
1926      aa-=BLOCKSQ;
1927    }
1928    if (triBase+BLOCK>numberRows_) {
1929      nChunk=numberRows_-triBase;
1930    } else {
1931      nChunk=BLOCK;
1932    }
1933    solveB1LongWork(aa,nChunk,region+triBase);
1934    aa-=BLOCKSQ;
1935  }
1936}
1937// Forward part of solve 1
1938void 
1939ClpCholeskyDense::solveF1LongWork(longDouble * a,int n,CoinWorkDouble * region)
1940{
1941  int j, k;
1942  CoinWorkDouble t00;
1943  for (j = 0; j < n; j ++) {
1944    t00 = region[j];
1945    for (k = 0; k < j; ++k) {
1946      t00 -= region[k] * a[j + k * BLOCK];
1947    }
1948    //t00*=a[j + j * BLOCK];
1949    region[j] = t00;
1950  }
1951}
1952// Forward part of solve 2
1953void 
1954ClpCholeskyDense::solveF2LongWork(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2)
1955{
1956  int j, k;
1957#ifdef BLOCKUNROLL
1958  if (n==BLOCK) {
1959    for (k = 0; k < BLOCK; k+=4) {
1960      CoinWorkDouble t0 = region2[0];
1961      CoinWorkDouble t1 = region2[1];
1962      CoinWorkDouble t2 = region2[2];
1963      CoinWorkDouble t3 = region2[3];
1964      t0 -= region[0] * a[0 + 0 * BLOCK];
1965      t1 -= region[0] * a[1 + 0 * BLOCK];
1966      t2 -= region[0] * a[2 + 0 * BLOCK];
1967      t3 -= region[0] * a[3 + 0 * BLOCK];
1968
1969      t0 -= region[1] * a[0 + 1 * BLOCK];
1970      t1 -= region[1] * a[1 + 1 * BLOCK];
1971      t2 -= region[1] * a[2 + 1 * BLOCK];
1972      t3 -= region[1] * a[3 + 1 * BLOCK];
1973
1974      t0 -= region[2] * a[0 + 2 * BLOCK];
1975      t1 -= region[2] * a[1 + 2 * BLOCK];
1976      t2 -= region[2] * a[2 + 2 * BLOCK];
1977      t3 -= region[2] * a[3 + 2 * BLOCK];
1978
1979      t0 -= region[3] * a[0 + 3 * BLOCK];
1980      t1 -= region[3] * a[1 + 3 * BLOCK];
1981      t2 -= region[3] * a[2 + 3 * BLOCK];
1982      t3 -= region[3] * a[3 + 3 * BLOCK];
1983
1984      t0 -= region[4] * a[0 + 4 * BLOCK];
1985      t1 -= region[4] * a[1 + 4 * BLOCK];
1986      t2 -= region[4] * a[2 + 4 * BLOCK];
1987      t3 -= region[4] * a[3 + 4 * BLOCK];
1988
1989      t0 -= region[5] * a[0 + 5 * BLOCK];
1990      t1 -= region[5] * a[1 + 5 * BLOCK];
1991      t2 -= region[5] * a[2 + 5 * BLOCK];
1992      t3 -= region[5] * a[3 + 5 * BLOCK];
1993
1994      t0 -= region[6] * a[0 + 6 * BLOCK];
1995      t1 -= region[6] * a[1 + 6 * BLOCK];
1996      t2 -= region[6] * a[2 + 6 * BLOCK];
1997      t3 -= region[6] * a[3 + 6 * BLOCK];
1998
1999      t0 -= region[7] * a[0 + 7 * BLOCK];
2000      t1 -= region[7] * a[1 + 7 * BLOCK];
2001      t2 -= region[7] * a[2 + 7 * BLOCK];
2002      t3 -= region[7] * a[3 + 7 * BLOCK];
2003#if BLOCK>8
2004      t0 -= region[8] * a[0 + 8 * BLOCK];
2005      t1 -= region[8] * a[1 + 8 * BLOCK];
2006      t2 -= region[8] * a[2 + 8 * BLOCK];
2007      t3 -= region[8] * a[3 + 8 * BLOCK];
2008
2009      t0 -= region[9] * a[0 + 9 * BLOCK];
2010      t1 -= region[9] * a[1 + 9 * BLOCK];
2011      t2 -= region[9] * a[2 + 9 * BLOCK];
2012      t3 -= region[9] * a[3 + 9 * BLOCK];
2013
2014      t0 -= region[10] * a[0 + 10 * BLOCK];
2015      t1 -= region[10] * a[1 + 10 * BLOCK];
2016      t2 -= region[10] * a[2 + 10 * BLOCK];
2017      t3 -= region[10] * a[3 + 10 * BLOCK];
2018
2019      t0 -= region[11] * a[0 + 11 * BLOCK];
2020      t1 -= region[11] * a[1 + 11 * BLOCK];
2021      t2 -= region[11] * a[2 + 11 * BLOCK];
2022      t3 -= region[11] * a[3 + 11 * BLOCK];
2023
2024      t0 -= region[12] * a[0 + 12 * BLOCK];
2025      t1 -= region[12] * a[1 + 12 * BLOCK];
2026      t2 -= region[12] * a[2 + 12 * BLOCK];
2027      t3 -= region[12] * a[3 + 12 * BLOCK];
2028
2029      t0 -= region[13] * a[0 + 13 * BLOCK];
2030      t1 -= region[13] * a[1 + 13 * BLOCK];
2031      t2 -= region[13] * a[2 + 13 * BLOCK];
2032      t3 -= region[13] * a[3 + 13 * BLOCK];
2033
2034      t0 -= region[14] * a[0 + 14 * BLOCK];
2035      t1 -= region[14] * a[1 + 14 * BLOCK];
2036      t2 -= region[14] * a[2 + 14 * BLOCK];
2037      t3 -= region[14] * a[3 + 14 * BLOCK];
2038
2039      t0 -= region[15] * a[0 + 15 * BLOCK];
2040      t1 -= region[15] * a[1 + 15 * BLOCK];
2041      t2 -= region[15] * a[2 + 15 * BLOCK];
2042      t3 -= region[15] * a[3 + 15 * BLOCK];
2043#endif
2044      region2[0] = t0;
2045      region2[1] = t1;
2046      region2[2] = t2;
2047      region2[3] = t3;
2048      region2+=4;
2049      a+=4;
2050    }
2051  } else {
2052#endif
2053    for (k = 0; k < n; ++k) {
2054      CoinWorkDouble t00 = region2[k];
2055      for (j = 0; j < BLOCK; j ++) {
2056        t00 -= region[j] * a[k + j * BLOCK];
2057      }
2058      region2[k] = t00;
2059    }
2060#ifdef BLOCKUNROLL
2061  }
2062#endif
2063}
2064// Backward part of solve 1
2065void 
2066ClpCholeskyDense::solveB1LongWork(longDouble * a,int n,CoinWorkDouble * region)
2067{
2068  int j, k;
2069  CoinWorkDouble t00;
2070  for (j = n-1; j >=0; j --) {
2071    t00 = region[j];
2072    for (k = j+1; k < n; ++k) {
2073      t00 -= region[k] * a[k + j * BLOCK];
2074    }
2075    //t00*=a[j + j * BLOCK];
2076    region[j] = t00;
2077  }
2078}
2079// Backward part of solve 2
2080void 
2081ClpCholeskyDense::solveB2LongWork(longDouble * a,int n,CoinWorkDouble * region, CoinWorkDouble * region2)
2082{
2083  int j, k;
2084#ifdef BLOCKUNROLL
2085  if (n==BLOCK) {
2086    for (j = 0; j < BLOCK; j +=4) {
2087      CoinWorkDouble t0 = region[0];
2088      CoinWorkDouble t1 = region[1];
2089      CoinWorkDouble t2 = region[2];
2090      CoinWorkDouble t3 = region[3];
2091      t0 -= region2[0] * a[0 + 0*BLOCK];
2092      t1 -= region2[0] * a[0 + 1*BLOCK];
2093      t2 -= region2[0] * a[0 + 2*BLOCK];
2094      t3 -= region2[0] * a[0 + 3*BLOCK];
2095
2096      t0 -= region2[1] * a[1 + 0*BLOCK];
2097      t1 -= region2[1] * a[1 + 1*BLOCK];
2098      t2 -= region2[1] * a[1 + 2*BLOCK];
2099      t3 -= region2[1] * a[1 + 3*BLOCK];
2100
2101      t0 -= region2[2] * a[2 + 0*BLOCK];
2102      t1 -= region2[2] * a[2 + 1*BLOCK];
2103      t2 -= region2[2] * a[2 + 2*BLOCK];
2104      t3 -= region2[2] * a[2 + 3*BLOCK];
2105
2106      t0 -= region2[3] * a[3 + 0*BLOCK];
2107      t1 -= region2[3] * a[3 + 1*BLOCK];
2108      t2 -= region2[3] * a[3 + 2*BLOCK];
2109      t3 -= region2[3] * a[3 + 3*BLOCK];
2110
2111      t0 -= region2[4] * a[4 + 0*BLOCK];
2112      t1 -= region2[4] * a[4 + 1*BLOCK];
2113      t2 -= region2[4] * a[4 + 2*BLOCK];
2114      t3 -= region2[4] * a[4 + 3*BLOCK];
2115
2116      t0 -= region2[5] * a[5 + 0*BLOCK];
2117      t1 -= region2[5] * a[5 + 1*BLOCK];
2118      t2 -= region2[5] * a[5 + 2*BLOCK];
2119      t3 -= region2[5] * a[5 + 3*BLOCK];
2120
2121      t0 -= region2[6] * a[6 + 0*BLOCK];
2122      t1 -= region2[6] * a[6 + 1*BLOCK];
2123      t2 -= region2[6] * a[6 + 2*BLOCK];
2124      t3 -= region2[6] * a[6 + 3*BLOCK];
2125
2126      t0 -= region2[7] * a[7 + 0*BLOCK];
2127      t1 -= region2[7] * a[7 + 1*BLOCK];
2128      t2 -= region2[7] * a[7 + 2*BLOCK];
2129      t3 -= region2[7] * a[7 + 3*BLOCK];
2130#if BLOCK>8
2131
2132      t0 -= region2[8] * a[8 + 0*BLOCK];
2133      t1 -= region2[8] * a[8 + 1*BLOCK];
2134      t2 -= region2[8] * a[8 + 2*BLOCK];
2135      t3 -= region2[8] * a[8 + 3*BLOCK];
2136
2137      t0 -= region2[9] * a[9 + 0*BLOCK];
2138      t1 -= region2[9] * a[9 + 1*BLOCK];
2139      t2 -= region2[9] * a[9 + 2*BLOCK];
2140      t3 -= region2[9] * a[9 + 3*BLOCK];
2141
2142      t0 -= region2[10] * a[10 + 0*BLOCK];
2143      t1 -= region2[10] * a[10 + 1*BLOCK];
2144      t2 -= region2[10] * a[10 + 2*BLOCK];
2145      t3 -= region2[10] * a[10 + 3*BLOCK];
2146
2147      t0 -= region2[11] * a[11 + 0*BLOCK];
2148      t1 -= region2[11] * a[11 + 1*BLOCK];
2149      t2 -= region2[11] * a[11 + 2*BLOCK];
2150      t3 -= region2[11] * a[11 + 3*BLOCK];
2151
2152      t0 -= region2[12] * a[12 + 0*BLOCK];
2153      t1 -= region2[12] * a[12 + 1*BLOCK];
2154      t2 -= region2[12] * a[12 + 2*BLOCK];
2155      t3 -= region2[12] * a[12 + 3*BLOCK];
2156
2157      t0 -= region2[13] * a[13 + 0*BLOCK];
2158      t1 -= region2[13] * a[13 + 1*BLOCK];
2159      t2 -= region2[13] * a[13 + 2*BLOCK];
2160      t3 -= region2[13] * a[13 + 3*BLOCK];
2161
2162      t0 -= region2[14] * a[14 + 0*BLOCK];
2163      t1 -= region2[14] * a[14 + 1*BLOCK];
2164      t2 -= region2[14] * a[14 + 2*BLOCK];
2165      t3 -= region2[14] * a[14 + 3*BLOCK];
2166
2167      t0 -= region2[15] * a[15 + 0*BLOCK];
2168      t1 -= region2[15] * a[15 + 1*BLOCK];
2169      t2 -= region2[15] * a[15 + 2*BLOCK];
2170      t3 -= region2[15] * a[15 + 3*BLOCK];
2171#endif
2172      region[0] = t0;
2173      region[1] = t1;
2174      region[2] = t2;
2175      region[3] = t3;
2176      a+=4*BLOCK;
2177      region+=4;
2178    }
2179  } else {
2180#endif
2181    for (j = 0; j < BLOCK; j ++) {
2182      CoinWorkDouble t00 = region[j];
2183      for (k = 0; k < n; ++k) {
2184        t00 -= region2[k] * a[k + j * BLOCK];
2185      }
2186      region[j] = t00;
2187    }
2188#ifdef BLOCKUNROLL
2189  }
2190#endif
2191}
Note: See TracBrowser for help on using the repository browser.