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

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

changes for cholesky including code from Anshul Gupta

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