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

Last change on this file since 799 was 799, checked in by andreasw, 14 years ago

undid last commit (patches incorrectly applied)

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