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

Last change on this file since 1525 was 1525, checked in by mjs, 10 years ago

Formatted .cpp, .hpp, .c, .h files with "astyle -A4 -p". This matches the formatting used in the grand CBC reorganization.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 62.4 KB
Line 
1/* $Id: ClpCholeskyDense.cpp 1525 2010-02-26 17:27:59Z mjs $ */
2/* Copyright (C) 2003, International Business Machines Corporation
3   and others.  All Rights Reserved. */
4#include "CoinPragma.hpp"
5#include "CoinHelperFunctions.hpp"
6#include "ClpHelperFunctions.hpp"
7
8#include "ClpInterior.hpp"
9#include "ClpCholeskyDense.hpp"
10#include "ClpMessage.hpp"
11#include "ClpQuadraticObjective.hpp"
12
13/*#############################################################################*/
14/* Constructors / Destructor / Assignment*/
15/*#############################################################################*/
16
17/*-------------------------------------------------------------------*/
18/* Default Constructor */
19/*-------------------------------------------------------------------*/
20ClpCholeskyDense::ClpCholeskyDense ()
21     : ClpCholeskyBase(),
22       borrowSpace_(false)
23{
24     type_ = 11;;
25}
26
27/*-------------------------------------------------------------------*/
28/* Copy constructor */
29/*-------------------------------------------------------------------*/
30ClpCholeskyDense::ClpCholeskyDense (const ClpCholeskyDense & rhs)
31     : ClpCholeskyBase(rhs),
32       borrowSpace_(rhs.borrowSpace_)
33{
34     assert(!rhs.borrowSpace_ || !rhs.sizeFactor_); /* can't do if borrowing space*/
35}
36
37
38/*-------------------------------------------------------------------*/
39/* Destructor */
40/*-------------------------------------------------------------------*/
41ClpCholeskyDense::~ClpCholeskyDense ()
42{
43     if (borrowSpace_) {
44          /* set NULL*/
45          sparseFactor_ = NULL;
46          workDouble_ = NULL;
47          diagonal_ = NULL;
48     }
49}
50
51/*----------------------------------------------------------------*/
52/* Assignment operator */
53/*-------------------------------------------------------------------*/
54ClpCholeskyDense &
55ClpCholeskyDense::operator=(const ClpCholeskyDense& rhs)
56{
57     if (this != &rhs) {
58          assert(!rhs.borrowSpace_ || !rhs.sizeFactor_); /* can't do if borrowing space*/
59          ClpCholeskyBase::operator=(rhs);
60          borrowSpace_ = rhs.borrowSpace_;
61     }
62     return *this;
63}
64/*-------------------------------------------------------------------*/
65/* Clone*/
66/*-------------------------------------------------------------------*/
67ClpCholeskyBase * ClpCholeskyDense::clone() const
68{
69     return new ClpCholeskyDense(*this);
70}
71/* If not power of 2 then need to redo a bit*/
72#define BLOCK 16
73#define BLOCKSHIFT 4
74/* Block unroll if power of 2 and at least 8*/
75#define BLOCKUNROLL
76
77#define BLOCKSQ ( BLOCK*BLOCK )
78#define BLOCKSQSHIFT ( BLOCKSHIFT+BLOCKSHIFT )
79#define number_blocks(x) (((x)+BLOCK-1)>>BLOCKSHIFT)
80#define number_rows(x) ((x)<<BLOCKSHIFT)
81#define number_entries(x) ((x)<<BLOCKSQSHIFT)
82/* Gets space */
83int
84ClpCholeskyDense::reserveSpace(const ClpCholeskyBase * factor, int numberRows)
85{
86     numberRows_ = numberRows;
87     int numberBlocks = (numberRows_ + BLOCK - 1) >> BLOCKSHIFT;
88     /* allow one stripe extra*/
89     numberBlocks = numberBlocks + ((numberBlocks * (numberBlocks + 1)) / 2);
90     sizeFactor_ = numberBlocks * BLOCKSQ;
91     /*#define CHOL_COMPARE*/
92#ifdef CHOL_COMPARE
93     sizeFactor_ += 95000;
94#endif
95     if (!factor) {
96          sparseFactor_ = new longDouble [sizeFactor_];
97          rowsDropped_ = new char [numberRows_];
98          memset(rowsDropped_, 0, numberRows_);
99          workDouble_ = new longDouble[numberRows_];
100          diagonal_ = new longDouble[numberRows_];
101     } else {
102          borrowSpace_ = true;
103          int numberFull = factor->numberRows();
104          sparseFactor_ = factor->sparseFactor() + (factor->size() - sizeFactor_);
105          workDouble_ = factor->workDouble() + (numberFull - numberRows_);
106          diagonal_ = factor->diagonal() + (numberFull - numberRows_);
107     }
108     numberRowsDropped_ = 0;
109     return 0;
110}
111/* Returns space needed */
112CoinBigIndex
113ClpCholeskyDense::space( int numberRows) const
114{
115     int numberBlocks = (numberRows + BLOCK - 1) >> BLOCKSHIFT;
116     /* allow one stripe extra*/
117     numberBlocks = numberBlocks + ((numberBlocks * (numberBlocks + 1)) / 2);
118     CoinBigIndex sizeFactor = numberBlocks * BLOCKSQ;
119#ifdef CHOL_COMPARE
120     sizeFactor += 95000;
121#endif
122     return sizeFactor;
123}
124/* Orders rows and saves pointer to matrix.and model */
125int
126ClpCholeskyDense::order(ClpInterior * model)
127{
128     model_ = model;
129     int numberRows;
130     int numberRowsModel = model_->numberRows();
131     int numberColumns = model_->numberColumns();
132     if (!doKKT_) {
133          numberRows = numberRowsModel;
134     } else {
135          numberRows = 2 * numberRowsModel + numberColumns;
136     }
137     reserveSpace(NULL, numberRows);
138     rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
139     return 0;
140}
141/* Does Symbolic factorization given permutation.
142   This is called immediately after order.  If user provides this then
143   user must provide factorize and solve.  Otherwise the default factorization is used
144   returns non-zero if not enough memory */
145int
146ClpCholeskyDense::symbolic()
147{
148     return 0;
149}
150/* Factorize - filling in rowsDropped and returning number dropped */
151int
152ClpCholeskyDense::factorize(const CoinWorkDouble * diagonal, int * rowsDropped)
153{
154     assert (!borrowSpace_);
155     const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
156     const int * columnLength = model_->clpMatrix()->getVectorLengths();
157     const int * row = model_->clpMatrix()->getIndices();
158     const double * element = model_->clpMatrix()->getElements();
159     const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
160     const int * rowLength = rowCopy_->getVectorLengths();
161     const int * column = rowCopy_->getIndices();
162     const double * elementByRow = rowCopy_->getElements();
163     int numberColumns = model_->clpMatrix()->getNumCols();
164     CoinZeroN(sparseFactor_, sizeFactor_);
165     /*perturbation*/
166     CoinWorkDouble perturbation = model_->diagonalPerturbation() * model_->diagonalNorm();
167     perturbation = perturbation * perturbation;
168     if (perturbation > 1.0) {
169#ifdef COIN_DEVELOP
170          /*if (model_->model()->logLevel()&4) */
171          std::cout << "large perturbation " << perturbation << std::endl;
172#endif
173          perturbation = CoinSqrt(perturbation);;
174          perturbation = 1.0;
175     }
176     int iRow;
177     int newDropped = 0;
178     CoinWorkDouble largest = 1.0;
179     CoinWorkDouble smallest = COIN_DBL_MAX;
180     CoinWorkDouble delta2 = model_->delta(); /* add delta*delta to diagonal*/
181     delta2 *= delta2;
182     if (!doKKT_) {
183          longDouble * work = sparseFactor_;
184          work--; /* skip diagonal*/
185          int addOffset = numberRows_ - 1;
186          const CoinWorkDouble * diagonalSlack = diagonal + numberColumns;
187          /* largest in initial matrix*/
188          CoinWorkDouble largest2 = 1.0e-20;
189          for (iRow = 0; iRow < numberRows_; iRow++) {
190               if (!rowsDropped_[iRow]) {
191                    CoinBigIndex startRow = rowStart[iRow];
192                    CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
193                    CoinWorkDouble diagonalValue = diagonalSlack[iRow] + delta2;
194                    for (CoinBigIndex k = startRow; k < endRow; k++) {
195                         int iColumn = column[k];
196                         CoinBigIndex start = columnStart[iColumn];
197                         CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
198                         CoinWorkDouble multiplier = diagonal[iColumn] * elementByRow[k];
199                         for (CoinBigIndex j = start; j < end; j++) {
200                              int jRow = row[j];
201                              if (!rowsDropped_[jRow]) {
202                                   if (jRow > iRow) {
203                                        CoinWorkDouble value = element[j] * multiplier;
204                                        work[jRow] += value;
205                                   } else if (jRow == iRow) {
206                                        CoinWorkDouble value = element[j] * multiplier;
207                                        diagonalValue += value;
208                                   }
209                              }
210                         }
211                    }
212                    for (int j = iRow + 1; j < numberRows_; j++)
213                         largest2 = CoinMax(largest2, CoinAbs(work[j]));
214                    diagonal_[iRow] = diagonalValue;
215                    largest2 = CoinMax(largest2, CoinAbs(diagonalValue));
216               } else {
217                    /* dropped*/
218                    diagonal_[iRow] = 1.0;
219               }
220               addOffset--;
221               work += addOffset;
222          }
223          /*check sizes*/
224          largest2 *= 1.0e-20;
225          largest = CoinMin(largest2, CHOL_SMALL_VALUE);
226          int numberDroppedBefore = 0;
227          for (iRow = 0; iRow < numberRows_; iRow++) {
228               int dropped = rowsDropped_[iRow];
229               /* Move to int array*/
230               rowsDropped[iRow] = dropped;
231               if (!dropped) {
232                    CoinWorkDouble diagonal = diagonal_[iRow];
233                    if (diagonal > largest2) {
234                         diagonal_[iRow] = diagonal + perturbation;
235                    } else {
236                         diagonal_[iRow] = diagonal + perturbation;
237                         rowsDropped[iRow] = 2;
238                         numberDroppedBefore++;
239                    }
240               }
241          }
242          doubleParameters_[10] = CoinMax(1.0e-20, largest);
243          integerParameters_[20] = 0;
244          doubleParameters_[3] = 0.0;
245          doubleParameters_[4] = COIN_DBL_MAX;
246          integerParameters_[34] = 0; /* say all must be positive*/
247#ifdef CHOL_COMPARE
248          if (numberRows_ < 200)
249               factorizePart3(rowsDropped);
250#endif
251          factorizePart2(rowsDropped);
252          newDropped = integerParameters_[20] + numberDroppedBefore;
253          largest = doubleParameters_[3];
254          smallest = doubleParameters_[4];
255          if (model_->messageHandler()->logLevel() > 1)
256               std::cout << "Cholesky - largest " << largest << " smallest " << smallest << std::endl;
257          choleskyCondition_ = largest / smallest;
258          /*drop fresh makes some formADAT easier*/
259          if (newDropped || numberRowsDropped_) {
260               newDropped = 0;
261               for (int i = 0; i < numberRows_; i++) {
262                    char dropped = static_cast<char>(rowsDropped[i]);
263                    rowsDropped_[i] = dropped;
264                    if (dropped == 2) {
265                         /*dropped this time*/
266                         rowsDropped[newDropped++] = i;
267                         rowsDropped_[i] = 0;
268                    }
269               }
270               numberRowsDropped_ = newDropped;
271               newDropped = -(2 + newDropped);
272          }
273     } else {
274          /* KKT*/
275          CoinPackedMatrix * quadratic = NULL;
276          ClpQuadraticObjective * quadraticObj =
277               (dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject()));
278          if (quadraticObj)
279               quadratic = quadraticObj->quadraticObjective();
280          /* matrix*/
281          int numberRowsModel = model_->numberRows();
282          int numberColumns = model_->numberColumns();
283          int numberTotal = numberColumns + numberRowsModel;
284          longDouble * work = sparseFactor_;
285          work--; /* skip diagonal*/
286          int addOffset = numberRows_ - 1;
287          int iColumn;
288          if (!quadratic) {
289               for (iColumn = 0; iColumn < numberColumns; iColumn++) {
290                    CoinWorkDouble value = diagonal[iColumn];
291                    if (CoinAbs(value) > 1.0e-100) {
292                         value = 1.0 / value;
293                         largest = CoinMax(largest, CoinAbs(value));
294                         diagonal_[iColumn] = -value;
295                         CoinBigIndex start = columnStart[iColumn];
296                         CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
297                         for (CoinBigIndex j = start; j < end; j++) {
298                              /*choleskyRow_[numberElements]=row[j]+numberTotal;*/
299                              /*sparseFactor_[numberElements++]=element[j];*/
300                              work[row[j] + numberTotal] = element[j];
301                              largest = CoinMax(largest, CoinAbs(element[j]));
302                         }
303                    } else {
304                         diagonal_[iColumn] = -value;
305                    }
306                    addOffset--;
307                    work += addOffset;
308               }
309          } else {
310               /* Quadratic*/
311               const int * columnQuadratic = quadratic->getIndices();
312               const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
313               const int * columnQuadraticLength = quadratic->getVectorLengths();
314               const double * quadraticElement = quadratic->getElements();
315               for (iColumn = 0; iColumn < numberColumns; iColumn++) {
316                    CoinWorkDouble value = diagonal[iColumn];
317                    CoinBigIndex j;
318                    if (CoinAbs(value) > 1.0e-100) {
319                         value = 1.0 / value;
320                         for (j = columnQuadraticStart[iColumn];
321                                   j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
322                              int jColumn = columnQuadratic[j];
323                              if (jColumn > iColumn) {
324                                   work[jColumn] = -quadraticElement[j];
325                              } else if (iColumn == jColumn) {
326                                   value += quadraticElement[j];
327                              }
328                         }
329                         largest = CoinMax(largest, CoinAbs(value));
330                         diagonal_[iColumn] = -value;
331                         CoinBigIndex start = columnStart[iColumn];
332                         CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
333                         for (j = start; j < end; j++) {
334                              work[row[j] + numberTotal] = element[j];
335                              largest = CoinMax(largest, CoinAbs(element[j]));
336                         }
337                    } else {
338                         value = 1.0e100;
339                         diagonal_[iColumn] = -value;
340                    }
341                    addOffset--;
342                    work += addOffset;
343               }
344          }
345          /* slacks*/
346          for (iColumn = numberColumns; iColumn < numberTotal; iColumn++) {
347               CoinWorkDouble value = diagonal[iColumn];
348               if (CoinAbs(value) > 1.0e-100) {
349                    value = 1.0 / value;
350                    largest = CoinMax(largest, CoinAbs(value));
351               } else {
352                    value = 1.0e100;
353               }
354               diagonal_[iColumn] = -value;
355               work[iColumn-numberColumns+numberTotal] = -1.0;
356               addOffset--;
357               work += addOffset;
358          }
359          /* Finish diagonal*/
360          for (iRow = 0; iRow < numberRowsModel; iRow++) {
361               diagonal_[iRow+numberTotal] = delta2;
362          }
363          /*check sizes*/
364          largest *= 1.0e-20;
365          largest = CoinMin(largest, CHOL_SMALL_VALUE);
366          doubleParameters_[10] = CoinMax(1.0e-20, largest);
367          integerParameters_[20] = 0;
368          doubleParameters_[3] = 0.0;
369          doubleParameters_[4] = COIN_DBL_MAX;
370          /* Set up LDL cutoff*/
371          integerParameters_[34] = numberTotal;
372          /* KKT*/
373          int * rowsDropped2 = new int[numberRows_];
374          CoinZeroN(rowsDropped2, numberRows_);
375#ifdef CHOL_COMPARE
376          if (numberRows_ < 200)
377               factorizePart3(rowsDropped2);
378#endif
379          factorizePart2(rowsDropped2);
380          newDropped = integerParameters_[20];
381          largest = doubleParameters_[3];
382          smallest = doubleParameters_[4];
383          if (model_->messageHandler()->logLevel() > 1)
384               std::cout << "Cholesky - largest " << largest << " smallest " << smallest << std::endl;
385          choleskyCondition_ = largest / smallest;
386          /* Should save adjustments in ..R_*/
387          int n1 = 0, n2 = 0;
388          CoinWorkDouble * primalR = model_->primalR();
389          CoinWorkDouble * dualR = model_->dualR();
390          for (iRow = 0; iRow < numberTotal; iRow++) {
391               if (rowsDropped2[iRow]) {
392                    n1++;
393                    /*printf("row region1 %d dropped\n",iRow);*/
394                    /*rowsDropped_[iRow]=1;*/
395                    rowsDropped_[iRow] = 0;
396                    primalR[iRow] = doubleParameters_[20];
397               } else {
398                    rowsDropped_[iRow] = 0;
399                    primalR[iRow] = 0.0;
400               }
401          }
402          for (; iRow < numberRows_; iRow++) {
403               if (rowsDropped2[iRow]) {
404                    n2++;
405                    /*printf("row region2 %d dropped\n",iRow);*/
406                    /*rowsDropped_[iRow]=1;*/
407                    rowsDropped_[iRow] = 0;
408                    dualR[iRow-numberTotal] = doubleParameters_[34];
409               } else {
410                    rowsDropped_[iRow] = 0;
411                    dualR[iRow-numberTotal] = 0.0;
412               }
413          }
414     }
415     return 0;
416}
417/* Factorize - filling in rowsDropped and returning number dropped */
418void
419ClpCholeskyDense::factorizePart3( int * rowsDropped)
420{
421     int iColumn;
422     longDouble * xx = sparseFactor_;
423     longDouble * yy = diagonal_;
424     diagonal_ = sparseFactor_ + 40000;
425     sparseFactor_ = diagonal_ + numberRows_;
426     /*memcpy(sparseFactor_,xx,sizeFactor_*sizeof(double));*/
427     CoinMemcpyN(xx, 40000, sparseFactor_);
428     CoinMemcpyN(yy, numberRows_, diagonal_);
429     int numberDropped = 0;
430     CoinWorkDouble largest = 0.0;
431     CoinWorkDouble smallest = COIN_DBL_MAX;
432     double dropValue = doubleParameters_[10];
433     int firstPositive = integerParameters_[34];
434     longDouble * work = sparseFactor_;
435     /* Allow for triangular*/
436     int addOffset = numberRows_ - 1;
437     work--;
438     for (iColumn = 0; iColumn < numberRows_; iColumn++) {
439          int iRow;
440          int addOffsetNow = numberRows_ - 1;;
441          longDouble * workNow = sparseFactor_ - 1 + iColumn;
442          CoinWorkDouble diagonalValue = diagonal_[iColumn];
443          for (iRow = 0; iRow < iColumn; iRow++) {
444               double aj = *workNow;
445               addOffsetNow--;
446               workNow += addOffsetNow;
447               diagonalValue -= aj * aj * workDouble_[iRow];
448          }
449          bool dropColumn = false;
450          if (iColumn < firstPositive) {
451               /* must be negative*/
452               if (diagonalValue <= -dropValue) {
453                    smallest = CoinMin(smallest, -diagonalValue);
454                    largest = CoinMax(largest, -diagonalValue);
455                    workDouble_[iColumn] = diagonalValue;
456                    diagonalValue = 1.0 / diagonalValue;
457               } else {
458                    dropColumn = true;
459                    workDouble_[iColumn] = -1.0e100;
460                    diagonalValue = 0.0;
461                    integerParameters_[20]++;
462               }
463          } else {
464               /* must be positive*/
465               if (diagonalValue >= dropValue) {
466                    smallest = CoinMin(smallest, diagonalValue);
467                    largest = CoinMax(largest, diagonalValue);
468                    workDouble_[iColumn] = diagonalValue;
469                    diagonalValue = 1.0 / diagonalValue;
470               } else {
471                    dropColumn = true;
472                    workDouble_[iColumn] = 1.0e100;
473                    diagonalValue = 0.0;
474                    integerParameters_[20]++;
475               }
476          }
477          if (!dropColumn) {
478               diagonal_[iColumn] = diagonalValue;
479               for (iRow = iColumn + 1; iRow < numberRows_; iRow++) {
480                    double value = work[iRow];
481                    workNow = sparseFactor_ - 1;
482                    int addOffsetNow = numberRows_ - 1;;
483                    for (int jColumn = 0; jColumn < iColumn; jColumn++) {
484                         double aj = workNow[iColumn];
485                         double multiplier = workDouble_[jColumn];
486                         double ai = workNow[iRow];
487                         addOffsetNow--;
488                         workNow += addOffsetNow;
489                         value -= aj * ai * multiplier;
490                    }
491                    work[iRow] = value * diagonalValue;
492               }
493          } else {
494               /* drop column*/
495               rowsDropped[iColumn] = 2;
496               numberDropped++;
497               diagonal_[iColumn] = 0.0;
498               for (iRow = iColumn + 1; iRow < numberRows_; iRow++) {
499                    work[iRow] = 0.0;
500               }
501          }
502          addOffset--;
503          work += addOffset;
504     }
505     doubleParameters_[3] = largest;
506     doubleParameters_[4] = smallest;
507     integerParameters_[20] = numberDropped;
508     sparseFactor_ = xx;
509     diagonal_ = yy;
510}
511/*#define POS_DEBUG*/
512#ifdef POS_DEBUG
513static int counter = 0;
514int ClpCholeskyDense::bNumber(const longDouble * array, int &iRow, int &iCol)
515{
516     int numberBlocks = (numberRows_ + BLOCK - 1) >> BLOCKSHIFT;
517     longDouble * a = sparseFactor_ + BLOCKSQ * numberBlocks;
518     int k = array - a;
519     assert ((k % BLOCKSQ) == 0);
520     iCol = 0;
521     int kk = k >> BLOCKSQSHIFT;
522     /*printf("%d %d %d %d\n",k,kk,BLOCKSQ,BLOCKSQSHIFT);*/
523     k = kk;
524     while (k >= numberBlocks) {
525          iCol++;
526          k -= numberBlocks;
527          numberBlocks--;
528     }
529     iRow = iCol + k;
530     counter++;
531     if(counter > 100000)
532          exit(77);
533     return kk;
534}
535#endif
536/* Factorize - filling in rowsDropped and returning number dropped */
537void
538ClpCholeskyDense::factorizePart2( int * rowsDropped)
539{
540     int iColumn;
541     /*longDouble * xx = sparseFactor_;*/
542     /*longDouble * yy = diagonal_;*/
543     /*diagonal_ = sparseFactor_ + 40000;*/
544     /*sparseFactor_=diagonal_ + numberRows_;*/
545     /*memcpy(sparseFactor_,xx,sizeFactor_*sizeof(double));*/
546     /*memcpy(sparseFactor_,xx,40000*sizeof(double));*/
547     /*memcpy(diagonal_,yy,numberRows_*sizeof(double));*/
548     int numberBlocks = (numberRows_ + BLOCK - 1) >> BLOCKSHIFT;
549     /* later align on boundary*/
550     longDouble * a = sparseFactor_ + BLOCKSQ * numberBlocks;
551     int n = numberRows_;
552     int nRound = numberRows_ & (~(BLOCK - 1));
553     /* adjust if exact*/
554     if (nRound == n)
555          nRound -= BLOCK;
556     int sizeLastBlock = n - nRound;
557     int get = n * (n - 1) / 2; /* ? as no diagonal*/
558     int block = numberBlocks * (numberBlocks + 1) / 2;
559     int ifOdd;
560     int rowLast;
561     if (sizeLastBlock != BLOCK) {
562          longDouble * aa = &a[(block-1)*BLOCKSQ];
563          rowLast = nRound - 1;
564          ifOdd = 1;
565          int put = BLOCKSQ;
566          /* do last separately*/
567          put -= (BLOCK - sizeLastBlock) * (BLOCK + 1);
568          for (iColumn = numberRows_ - 1; iColumn >= nRound; iColumn--) {
569               int put2 = put;
570               put -= BLOCK;
571               for (int iRow = numberRows_ - 1; iRow > iColumn; iRow--) {
572                    aa[--put2] = sparseFactor_[--get];
573                    assert (aa + put2 >= sparseFactor_ + get);
574               }
575               /* save diagonal as well*/
576               aa[--put2] = diagonal_[iColumn];
577          }
578          n = nRound;
579          block--;
580     } else {
581          /* exact fit*/
582          rowLast = numberRows_ - 1;
583          ifOdd = 0;
584     }
585     /* Now main loop*/
586     int nBlock = 0;
587     for (; n > 0; n -= BLOCK) {
588          longDouble * aa = &a[(block-1)*BLOCKSQ];
589          longDouble * aaLast = NULL;
590          int put = BLOCKSQ;
591          int putLast = 0;
592          /* see if we have small block*/
593          if (ifOdd) {
594               aaLast = &a[(block-1)*BLOCKSQ];
595               aa = aaLast - BLOCKSQ;
596               putLast = BLOCKSQ - BLOCK + sizeLastBlock;
597          }
598          for (iColumn = n - 1; iColumn >= n - BLOCK; iColumn--) {
599               if (aaLast) {
600                    /* last bit*/
601                    for (int iRow = numberRows_ - 1; iRow > rowLast; iRow--) {
602                         aaLast[--putLast] = sparseFactor_[--get];
603                         assert (aaLast + putLast >= sparseFactor_ + get);
604                    }
605                    putLast -= BLOCK - sizeLastBlock;
606               }
607               longDouble * aPut = aa;
608               int j = rowLast;
609               for (int jBlock = 0; jBlock <= nBlock; jBlock++) {
610                    int put2 = put;
611                    int last = CoinMax(j - BLOCK, iColumn);
612                    for (int iRow = j; iRow > last; iRow--) {
613                         aPut[--put2] = sparseFactor_[--get];
614                         assert (aPut + put2 >= sparseFactor_ + get);
615                    }
616                    if (j - BLOCK < iColumn) {
617                         /* save diagonal as well*/
618                         aPut[--put2] = diagonal_[iColumn];
619                    }
620                    j -= BLOCK;
621                    aPut -= BLOCKSQ;
622               }
623               put -= BLOCK;
624          }
625          nBlock++;
626          block -= nBlock + ifOdd;
627     }
628     ClpCholeskyDenseC  info;
629     info.diagonal_ = diagonal_;
630     info.doubleParameters_[0] = doubleParameters_[10];
631     info.integerParameters_[0] = integerParameters_[34];
632#ifndef CLP_CILK
633     ClpCholeskyCfactor(&info, a, numberRows_, numberBlocks,
634                        diagonal_, workDouble_, rowsDropped);
635#else
636     info.a = a;
637     info.n = numberRows_;
638     info.numberBlocks = numberBlocks;
639     info.work = workDouble_;
640     info.rowsDropped = rowsDropped;
641     info.integerParameters_[1] = model_->numberThreads();
642     ClpCholeskySpawn(&info);
643#endif
644     double largest = 0.0;
645     double smallest = COIN_DBL_MAX;
646     int numberDropped = 0;
647     for (int i = 0; i < numberRows_; i++) {
648          if (diagonal_[i]) {
649               largest = CoinMax(largest, CoinAbs(diagonal_[i]));
650               smallest = CoinMin(smallest, CoinAbs(diagonal_[i]));
651          } else {
652               numberDropped++;
653          }
654     }
655     doubleParameters_[3] = CoinMax(doubleParameters_[3], 1.0 / smallest);
656     doubleParameters_[4] = CoinMin(doubleParameters_[4], 1.0 / largest);
657     integerParameters_[20] += numberDropped;
658}
659/* Non leaf recursive factor*/
660void
661ClpCholeskyCfactor(ClpCholeskyDenseC * thisStruct, longDouble * a, int n, int numberBlocks,
662                   longDouble * diagonal, longDouble * work, int * rowsDropped)
663{
664     if (n <= BLOCK) {
665          ClpCholeskyCfactorLeaf(thisStruct, a, n, diagonal, work, rowsDropped);
666     } else {
667          int nb = number_blocks((n + 1) >> 1);
668          int nThis = number_rows(nb);
669          longDouble * aother;
670          int nLeft = n - nThis;
671          int nintri = (nb * (nb + 1)) >> 1;
672          int nbelow = (numberBlocks - nb) * nb;
673          ClpCholeskyCfactor(thisStruct, a, nThis, numberBlocks, diagonal, work, rowsDropped);
674          ClpCholeskyCtriRec(thisStruct, a, nThis, a + number_entries(nb), diagonal, work, nLeft, nb, 0, numberBlocks);
675          aother = a + number_entries(nintri + nbelow);
676          ClpCholeskyCrecTri(thisStruct, a + number_entries(nb), nLeft, nThis, nb, 0, aother, diagonal, work, numberBlocks);
677          ClpCholeskyCfactor(thisStruct, aother, nLeft,
678                             numberBlocks - nb, diagonal + nThis, work + nThis, rowsDropped);
679     }
680}
681/* Non leaf recursive triangle rectangle update*/
682void
683ClpCholeskyCtriRec(ClpCholeskyDenseC * thisStruct, longDouble * aTri, int nThis, longDouble * aUnder,
684                   longDouble * diagonal, longDouble * work,
685                   int nLeft, int iBlock, int jBlock,
686                   int numberBlocks)
687{
688     if (nThis <= BLOCK && nLeft <= BLOCK) {
689          ClpCholeskyCtriRecLeaf(/*thisStruct,*/ aTri, aUnder, diagonal, work, nLeft);
690     } else if (nThis < nLeft) {
691          int nb = number_blocks((nLeft + 1) >> 1);
692          int nLeft2 = number_rows(nb);
693          ClpCholeskyCtriRec(thisStruct, aTri, nThis, aUnder, diagonal, work, nLeft2, iBlock, jBlock, numberBlocks);
694          ClpCholeskyCtriRec(thisStruct, aTri, nThis, aUnder + number_entries(nb), diagonal, work, nLeft - nLeft2,
695                             iBlock + nb, jBlock, numberBlocks);
696     } else {
697          int nb = number_blocks((nThis + 1) >> 1);
698          int nThis2 = number_rows(nb);
699          longDouble * aother;
700          int kBlock = jBlock + nb;
701          int i;
702          int nintri = (nb * (nb + 1)) >> 1;
703          int nbelow = (numberBlocks - nb) * nb;
704          ClpCholeskyCtriRec(thisStruct, aTri, nThis2, aUnder, diagonal, work, nLeft, iBlock, jBlock, numberBlocks);
705          /* and rectangular update */
706          i = ((numberBlocks - jBlock) * (numberBlocks - jBlock - 1) -
707               (numberBlocks - jBlock - nb) * (numberBlocks - jBlock - nb - 1)) >> 1;
708          aother = aUnder + number_entries(i);
709          ClpCholeskyCrecRec(thisStruct, aTri + number_entries(nb), nThis - nThis2, nLeft, nThis2, aUnder, aother,
710                             work, kBlock, jBlock, numberBlocks);
711          ClpCholeskyCtriRec(thisStruct, aTri + number_entries(nintri + nbelow), nThis - nThis2, aother, diagonal + nThis2,
712                             work + nThis2, nLeft,
713                             iBlock - nb, kBlock - nb, numberBlocks - nb);
714     }
715}
716/* Non leaf recursive rectangle triangle update*/
717void
718ClpCholeskyCrecTri(ClpCholeskyDenseC * thisStruct, longDouble * aUnder, int nTri, int nDo,
719                   int iBlock, int jBlock, longDouble * aTri,
720                   longDouble * diagonal, longDouble * work,
721                   int numberBlocks)
722{
723     if (nTri <= BLOCK && nDo <= BLOCK) {
724          ClpCholeskyCrecTriLeaf(/*thisStruct,*/ aUnder, aTri,/*diagonal,*/work, nTri);
725     } else if (nTri < nDo) {
726          int nb = number_blocks((nDo + 1) >> 1);
727          int nDo2 = number_rows(nb);
728          longDouble * aother;
729          int i;
730          ClpCholeskyCrecTri(thisStruct, aUnder, nTri, nDo2, iBlock, jBlock, aTri, diagonal, work, numberBlocks);
731          i = ((numberBlocks - jBlock) * (numberBlocks - jBlock - 1) -
732               (numberBlocks - jBlock - nb) * (numberBlocks - jBlock - nb - 1)) >> 1;
733          aother = aUnder + number_entries(i);
734          ClpCholeskyCrecTri(thisStruct, aother, nTri, nDo - nDo2, iBlock - nb, jBlock, aTri, diagonal + nDo2,
735                             work + nDo2, numberBlocks - nb);
736     } else {
737          int nb = number_blocks((nTri + 1) >> 1);
738          int nTri2 = number_rows(nb);
739          longDouble * aother;
740          int i;
741          ClpCholeskyCrecTri(thisStruct, aUnder, nTri2, nDo, iBlock, jBlock, aTri, diagonal, work, numberBlocks);
742          /* and rectangular update */
743          i = ((numberBlocks - iBlock) * (numberBlocks - iBlock + 1) -
744               (numberBlocks - iBlock - nb) * (numberBlocks - iBlock - nb + 1)) >> 1;
745          aother = aTri + number_entries(nb);
746          ClpCholeskyCrecRec(thisStruct, aUnder, nTri2, nTri - nTri2, nDo, aUnder + number_entries(nb), aother,
747                             work, iBlock, jBlock, numberBlocks);
748          ClpCholeskyCrecTri(thisStruct, aUnder + number_entries(nb), nTri - nTri2, nDo, iBlock + nb, jBlock,
749                             aTri + number_entries(i), diagonal, work, numberBlocks);
750     }
751}
752/* Non leaf recursive rectangle rectangle update,
753   nUnder is number of rows in iBlock,
754   nUnderK is number of rows in kBlock
755*/
756void
757ClpCholeskyCrecRec(ClpCholeskyDenseC * thisStruct, longDouble * above, int nUnder, int nUnderK,
758                   int nDo, longDouble * aUnder, longDouble *aOther,
759                   longDouble * work,
760                   int iBlock, int jBlock,
761                   int numberBlocks)
762{
763     if (nDo <= BLOCK && nUnder <= BLOCK && nUnderK <= BLOCK) {
764          assert (nDo == BLOCK && nUnder == BLOCK);
765          ClpCholeskyCrecRecLeaf(/*thisStruct,*/ above , aUnder ,  aOther, work, nUnderK);
766     } else if (nDo <= nUnderK && nUnder <= nUnderK) {
767          int nb = number_blocks((nUnderK + 1) >> 1);
768          int nUnder2 = number_rows(nb);
769          ClpCholeskyCrecRec(thisStruct, above, nUnder, nUnder2, nDo, aUnder, aOther, work,
770                             iBlock, jBlock, numberBlocks);
771          ClpCholeskyCrecRec(thisStruct, above, nUnder, nUnderK - nUnder2, nDo, aUnder + number_entries(nb),
772                             aOther + number_entries(nb), work, iBlock, jBlock, numberBlocks);
773     } else if (nUnderK <= nDo && nUnder <= nDo) {
774          int nb = number_blocks((nDo + 1) >> 1);
775          int nDo2 = number_rows(nb);
776          int i;
777          ClpCholeskyCrecRec(thisStruct, above, nUnder, nUnderK, nDo2, aUnder, aOther, work,
778                             iBlock, jBlock, numberBlocks);
779          i = ((numberBlocks - jBlock) * (numberBlocks - jBlock - 1) -
780               (numberBlocks - jBlock - nb) * (numberBlocks - jBlock - nb - 1)) >> 1;
781          ClpCholeskyCrecRec(thisStruct, above + number_entries(i), nUnder, nUnderK, nDo - nDo2,
782                             aUnder + number_entries(i),
783                             aOther, work + nDo2,
784                             iBlock - nb, jBlock, numberBlocks - nb);
785     } else {
786          int nb = number_blocks((nUnder + 1) >> 1);
787          int nUnder2 = number_rows(nb);
788          int i;
789          ClpCholeskyCrecRec(thisStruct, above, nUnder2, nUnderK, nDo, aUnder, aOther, work,
790                             iBlock, jBlock, numberBlocks);
791          i = ((numberBlocks - iBlock) * (numberBlocks - iBlock - 1) -
792               (numberBlocks - iBlock - nb) * (numberBlocks - iBlock - nb - 1)) >> 1;
793          ClpCholeskyCrecRec(thisStruct, above + number_entries(nb), nUnder - nUnder2, nUnderK, nDo, aUnder,
794                             aOther + number_entries(i), work, iBlock + nb, jBlock, numberBlocks);
795     }
796}
797/* Leaf recursive factor*/
798void
799ClpCholeskyCfactorLeaf(ClpCholeskyDenseC * thisStruct, longDouble * a, int n,
800                       longDouble * diagonal, longDouble * work, int * rowsDropped)
801{
802     double dropValue = thisStruct->doubleParameters_[0];
803     int firstPositive = thisStruct->integerParameters_[0];
804     int rowOffset = diagonal - thisStruct->diagonal_;
805     int i, j, k;
806     CoinWorkDouble t00, temp1;
807     longDouble * aa;
808     aa = a - BLOCK;
809     for (j = 0; j < n; j ++) {
810          bool dropColumn;
811          CoinWorkDouble useT00;
812          aa += BLOCK;
813          t00 = aa[j];
814          for (k = 0; k < j; ++k) {
815               CoinWorkDouble multiplier = work[k];
816               t00 -= a[j + k * BLOCK] * a[j + k * BLOCK] * multiplier;
817          }
818          dropColumn = false;
819          useT00 = t00;
820          if (j + rowOffset < firstPositive) {
821               /* must be negative*/
822               if (t00 <= -dropValue) {
823                    /*aa[j]=t00;*/
824                    t00 = 1.0 / t00;
825               } else {
826                    dropColumn = true;
827                    /*aa[j]=-1.0e100;*/
828                    useT00 = -1.0e-100;
829                    t00 = 0.0;
830               }
831          } else {
832               /* must be positive*/
833               if (t00 >= dropValue) {
834                    /*aa[j]=t00;*/
835                    t00 = 1.0 / t00;
836               } else {
837                    dropColumn = true;
838                    /*aa[j]=1.0e100;*/
839                    useT00 = 1.0e-100;
840                    t00 = 0.0;
841               }
842          }
843          if (!dropColumn) {
844               diagonal[j] = t00;
845               work[j] = useT00;
846               temp1 = t00;
847               for (i = j + 1; i < n; i++) {
848                    t00 = aa[i];
849                    for (k = 0; k < j; ++k) {
850                         CoinWorkDouble multiplier = work[k];
851                         t00 -= a[i + k * BLOCK] * a[j + k * BLOCK] * multiplier;
852                    }
853                    aa[i] = t00 * temp1;
854               }
855          } else {
856               /* drop column*/
857               rowsDropped[j+rowOffset] = 2;
858               diagonal[j] = 0.0;
859               /*aa[j]=1.0e100;*/
860               work[j] = 1.0e100;
861               for (i = j + 1; i < n; i++) {
862                    aa[i] = 0.0;
863               }
864          }
865     }
866}
867/* Leaf recursive triangle rectangle update*/
868void
869ClpCholeskyCtriRecLeaf(/*ClpCholeskyDenseC * thisStruct,*/ longDouble * aTri, longDouble * aUnder, longDouble * diagonal, longDouble * work,
870          int nUnder)
871{
872#ifdef POS_DEBUG
873     int iru, icu;
874     int iu = bNumber(aUnder, iru, icu);
875     int irt, ict;
876     int it = bNumber(aTri, irt, ict);
877     /*printf("%d %d\n",iu,it);*/
878     printf("trirecleaf  under (%d,%d), tri (%d,%d)\n",
879            iru, icu, irt, ict);
880     assert (diagonal == thisStruct->diagonal_ + ict * BLOCK);
881#endif
882     int j;
883     longDouble * aa;
884#ifdef BLOCKUNROLL
885     if (nUnder == BLOCK) {
886          aa = aTri - 2 * BLOCK;
887          for (j = 0; j < BLOCK; j += 2) {
888               int i;
889               CoinWorkDouble temp0 = diagonal[j];
890               CoinWorkDouble temp1 = diagonal[j+1];
891               aa += 2 * BLOCK;
892               for ( i = 0; i < BLOCK; i += 2) {
893                    CoinWorkDouble at1;
894                    CoinWorkDouble t00 = aUnder[i+j*BLOCK];
895                    CoinWorkDouble t10 = aUnder[i+ BLOCK + j*BLOCK];
896                    CoinWorkDouble t01 = aUnder[i+1+j*BLOCK];
897                    CoinWorkDouble t11 = aUnder[i+1+ BLOCK + j*BLOCK];
898                    int k;
899                    for (k = 0; k < j; ++k) {
900                         CoinWorkDouble multiplier = work[k];
901                         CoinWorkDouble au0 = aUnder[i + k * BLOCK] * multiplier;
902                         CoinWorkDouble au1 = aUnder[i + 1 + k * BLOCK] * multiplier;
903                         CoinWorkDouble at0 = aTri[j + k * BLOCK];
904                         CoinWorkDouble at1 = aTri[j + 1 + k * BLOCK];
905                         t00 -= au0 * at0;
906                         t10 -= au0 * at1;
907                         t01 -= au1 * at0;
908                         t11 -= au1 * at1;
909                    }
910                    t00 *= temp0;
911                    at1 = aTri[j + 1 + j*BLOCK] * work[j];
912                    t10 -= t00 * at1;
913                    t01 *= temp0;
914                    t11 -= t01 * at1;
915                    aUnder[i+j*BLOCK] = t00;
916                    aUnder[i+1+j*BLOCK] = t01;
917                    aUnder[i+ BLOCK + j*BLOCK] = t10 * temp1;
918                    aUnder[i+1+ BLOCK + j*BLOCK] = t11 * temp1;
919               }
920          }
921     } else {
922#endif
923          aa = aTri - BLOCK;
924          for (j = 0; j < BLOCK; j ++) {
925               int i;
926               CoinWorkDouble temp1 = diagonal[j];
927               aa += BLOCK;
928               for (i = 0; i < nUnder; i++) {
929                    int k;
930                    CoinWorkDouble t00 = aUnder[i+j*BLOCK];
931                    for ( k = 0; k < j; ++k) {
932                         CoinWorkDouble multiplier = work[k];
933                         t00 -= aUnder[i + k * BLOCK] * aTri[j + k * BLOCK] * multiplier;
934                    }
935                    aUnder[i+j*BLOCK] = t00 * temp1;
936               }
937          }
938#ifdef BLOCKUNROLL
939     }
940#endif
941}
942/* Leaf recursive rectangle triangle update*/
943void ClpCholeskyCrecTriLeaf(/*ClpCholeskyDenseC * thisStruct,*/ longDouble * aUnder, longDouble * aTri,
944          /*longDouble * diagonal,*/ longDouble * work, int nUnder)
945{
946#ifdef POS_DEBUG
947     int iru, icu;
948     int iu = bNumber(aUnder, iru, icu);
949     int irt, ict;
950     int it = bNumber(aTri, irt, ict);
951     /*printf("%d %d\n",iu,it);*/
952     printf("rectrileaf  under (%d,%d), tri (%d,%d)\n",
953            iru, icu, irt, ict);
954     assert (diagonal == thisStruct->diagonal_ + icu * BLOCK);
955#endif
956     int i, j, k;
957     CoinWorkDouble t00;
958     longDouble * aa;
959#ifdef BLOCKUNROLL
960     if (nUnder == BLOCK) {
961          longDouble * aUnder2 = aUnder - 2;
962          aa = aTri - 2 * BLOCK;
963          for (j = 0; j < BLOCK; j += 2) {
964               CoinWorkDouble t00, t01, t10, t11;
965               aa += 2 * BLOCK;
966               aUnder2 += 2;
967               t00 = aa[j];
968               t01 = aa[j+1];
969               t10 = aa[j+1+BLOCK];
970               for (k = 0; k < BLOCK; ++k) {
971                    CoinWorkDouble multiplier = work[k];
972                    CoinWorkDouble a0 = aUnder2[k * BLOCK];
973                    CoinWorkDouble a1 = aUnder2[1 + k * BLOCK];
974                    CoinWorkDouble x0 = a0 * multiplier;
975                    CoinWorkDouble x1 = a1 * multiplier;
976                    t00 -= a0 * x0;
977                    t01 -= a1 * x0;
978                    t10 -= a1 * x1;
979               }
980               aa[j] = t00;
981               aa[j+1] = t01;
982               aa[j+1+BLOCK] = t10;
983               for (i = j + 2; i < BLOCK; i += 2) {
984                    t00 = aa[i];
985                    t01 = aa[i+BLOCK];
986                    t10 = aa[i+1];
987                    t11 = aa[i+1+BLOCK];
988                    for (k = 0; k < BLOCK; ++k) {
989                         CoinWorkDouble multiplier = work[k];
990                         CoinWorkDouble a0 = aUnder2[k * BLOCK] * multiplier;
991                         CoinWorkDouble a1 = aUnder2[1 + k * BLOCK] * multiplier;
992                         t00 -= aUnder[i + k * BLOCK] * a0;
993                         t01 -= aUnder[i + k * BLOCK] * a1;
994                         t10 -= aUnder[i + 1 + k * BLOCK] * a0;
995                         t11 -= aUnder[i + 1 + k * BLOCK] * a1;
996                    }
997                    aa[i] = t00;
998                    aa[i+BLOCK] = t01;
999                    aa[i+1] = t10;
1000                    aa[i+1+BLOCK] = t11;
1001               }
1002          }
1003     } else {
1004#endif
1005          aa = aTri - BLOCK;
1006          for (j = 0; j < nUnder; j ++) {
1007               aa += BLOCK;
1008               for (i = j; i < nUnder; i++) {
1009                    t00 = aa[i];
1010                    for (k = 0; k < BLOCK; ++k) {
1011                         CoinWorkDouble multiplier = work[k];
1012                         t00 -= aUnder[i + k * BLOCK] * aUnder[j + k * BLOCK] * multiplier;
1013                    }
1014                    aa[i] = t00;
1015               }
1016          }
1017#ifdef BLOCKUNROLL
1018     }
1019#endif
1020}
1021/* Leaf recursive rectangle rectangle update,
1022   nUnder is number of rows in iBlock,
1023   nUnderK is number of rows in kBlock
1024*/
1025void
1026ClpCholeskyCrecRecLeaf(/*ClpCholeskyDenseC * thisStruct,*/
1027     const longDouble * COIN_RESTRICT above,
1028     const longDouble * COIN_RESTRICT aUnder,
1029     longDouble * COIN_RESTRICT aOther,
1030     const longDouble * COIN_RESTRICT work,
1031     int nUnder)
1032{
1033#ifdef POS_DEBUG
1034     int ira, ica;
1035     int ia = bNumber(above, ira, ica);
1036     int iru, icu;
1037     int iu = bNumber(aUnder, iru, icu);
1038     int iro, ico;
1039     int io = bNumber(aOther, iro, ico);
1040     /*printf("%d %d %d\n",ia,iu,io);*/
1041     printf("recrecleaf above (%d,%d), under (%d,%d), other (%d,%d)\n",
1042            ira, ica, iru, icu, iro, ico);
1043#endif
1044     int i, j, k;
1045     longDouble * aa;
1046#ifdef BLOCKUNROLL
1047     aa = aOther - 4 * BLOCK;
1048     if (nUnder == BLOCK) {
1049          /*#define INTEL*/
1050#ifdef INTEL
1051          aa += 2 * BLOCK;
1052          for (j = 0; j < BLOCK; j += 2) {
1053               aa += 2 * BLOCK;
1054               for (i = 0; i < BLOCK; i += 2) {
1055                    CoinWorkDouble t00 = aa[i+0*BLOCK];
1056                    CoinWorkDouble t10 = aa[i+1*BLOCK];
1057                    CoinWorkDouble t01 = aa[i+1+0*BLOCK];
1058                    CoinWorkDouble t11 = aa[i+1+1*BLOCK];
1059                    for (k = 0; k < BLOCK; k++) {
1060                         CoinWorkDouble multiplier = work[k];
1061                         CoinWorkDouble a00 = aUnder[i+k*BLOCK] * multiplier;
1062                         CoinWorkDouble a01 = aUnder[i+1+k*BLOCK] * multiplier;
1063                         t00 -= a00 * above[j + 0 + k * BLOCK];
1064                         t10 -= a00 * above[j + 1 + k * BLOCK];
1065                         t01 -= a01 * above[j + 0 + k * BLOCK];
1066                         t11 -= a01 * above[j + 1 + k * BLOCK];
1067                    }
1068                    aa[i+0*BLOCK] = t00;
1069                    aa[i+1*BLOCK] = t10;
1070                    aa[i+1+0*BLOCK] = t01;
1071                    aa[i+1+1*BLOCK] = t11;
1072               }
1073          }
1074#else
1075          for (j = 0; j < BLOCK; j += 4) {
1076               aa += 4 * BLOCK;
1077               for (i = 0; i < BLOCK; i += 4) {
1078                    CoinWorkDouble t00 = aa[i+0+0*BLOCK];
1079                    CoinWorkDouble t10 = aa[i+0+1*BLOCK];
1080                    CoinWorkDouble t20 = aa[i+0+2*BLOCK];
1081                    CoinWorkDouble t30 = aa[i+0+3*BLOCK];
1082                    CoinWorkDouble t01 = aa[i+1+0*BLOCK];
1083                    CoinWorkDouble t11 = aa[i+1+1*BLOCK];
1084                    CoinWorkDouble t21 = aa[i+1+2*BLOCK];
1085                    CoinWorkDouble t31 = aa[i+1+3*BLOCK];
1086                    CoinWorkDouble t02 = aa[i+2+0*BLOCK];
1087                    CoinWorkDouble t12 = aa[i+2+1*BLOCK];
1088                    CoinWorkDouble t22 = aa[i+2+2*BLOCK];
1089                    CoinWorkDouble t32 = aa[i+2+3*BLOCK];
1090                    CoinWorkDouble t03 = aa[i+3+0*BLOCK];
1091                    CoinWorkDouble t13 = aa[i+3+1*BLOCK];
1092                    CoinWorkDouble t23 = aa[i+3+2*BLOCK];
1093                    CoinWorkDouble t33 = aa[i+3+3*BLOCK];
1094                    const longDouble * COIN_RESTRICT aUnderNow = aUnder + i;
1095                    const longDouble * COIN_RESTRICT aboveNow = above + j;
1096                    for (k = 0; k < BLOCK; k++) {
1097                         CoinWorkDouble multiplier = work[k];
1098                         CoinWorkDouble a00 = aUnderNow[0] * multiplier;
1099                         CoinWorkDouble a01 = aUnderNow[1] * multiplier;
1100                         CoinWorkDouble a02 = aUnderNow[2] * multiplier;
1101                         CoinWorkDouble a03 = aUnderNow[3] * multiplier;
1102                         t00 -= a00 * aboveNow[0];
1103                         t10 -= a00 * aboveNow[1];
1104                         t20 -= a00 * aboveNow[2];
1105                         t30 -= a00 * aboveNow[3];
1106                         t01 -= a01 * aboveNow[0];
1107                         t11 -= a01 * aboveNow[1];
1108                         t21 -= a01 * aboveNow[2];
1109                         t31 -= a01 * aboveNow[3];
1110                         t02 -= a02 * aboveNow[0];
1111                         t12 -= a02 * aboveNow[1];
1112                         t22 -= a02 * aboveNow[2];
1113                         t32 -= a02 * aboveNow[3];
1114                         t03 -= a03 * aboveNow[0];
1115                         t13 -= a03 * aboveNow[1];
1116                         t23 -= a03 * aboveNow[2];
1117                         t33 -= a03 * aboveNow[3];
1118                         aUnderNow += BLOCK;
1119                         aboveNow += BLOCK;
1120                    }
1121                    aa[i+0+0*BLOCK] = t00;
1122                    aa[i+0+1*BLOCK] = t10;
1123                    aa[i+0+2*BLOCK] = t20;
1124                    aa[i+0+3*BLOCK] = t30;
1125                    aa[i+1+0*BLOCK] = t01;
1126                    aa[i+1+1*BLOCK] = t11;
1127                    aa[i+1+2*BLOCK] = t21;
1128                    aa[i+1+3*BLOCK] = t31;
1129                    aa[i+2+0*BLOCK] = t02;
1130                    aa[i+2+1*BLOCK] = t12;
1131                    aa[i+2+2*BLOCK] = t22;
1132                    aa[i+2+3*BLOCK] = t32;
1133                    aa[i+3+0*BLOCK] = t03;
1134                    aa[i+3+1*BLOCK] = t13;
1135                    aa[i+3+2*BLOCK] = t23;
1136                    aa[i+3+3*BLOCK] = t33;
1137               }
1138          }
1139#endif
1140     } else {
1141          int odd = nUnder & 1;
1142          int n = nUnder - odd;
1143          for (j = 0; j < BLOCK; j += 4) {
1144               aa += 4 * BLOCK;
1145               for (i = 0; i < n; i += 2) {
1146                    CoinWorkDouble t00 = aa[i+0*BLOCK];
1147                    CoinWorkDouble t10 = aa[i+1*BLOCK];
1148                    CoinWorkDouble t20 = aa[i+2*BLOCK];
1149                    CoinWorkDouble t30 = aa[i+3*BLOCK];
1150                    CoinWorkDouble t01 = aa[i+1+0*BLOCK];
1151                    CoinWorkDouble t11 = aa[i+1+1*BLOCK];
1152                    CoinWorkDouble t21 = aa[i+1+2*BLOCK];
1153                    CoinWorkDouble t31 = aa[i+1+3*BLOCK];
1154                    const longDouble * COIN_RESTRICT aUnderNow = aUnder + i;
1155                    const longDouble * COIN_RESTRICT aboveNow = above + j;
1156                    for (k = 0; k < BLOCK; k++) {
1157                         CoinWorkDouble multiplier = work[k];
1158                         CoinWorkDouble a00 = aUnderNow[0] * multiplier;
1159                         CoinWorkDouble a01 = aUnderNow[1] * multiplier;
1160                         t00 -= a00 * aboveNow[0];
1161                         t10 -= a00 * aboveNow[1];
1162                         t20 -= a00 * aboveNow[2];
1163                         t30 -= a00 * aboveNow[3];
1164                         t01 -= a01 * aboveNow[0];
1165                         t11 -= a01 * aboveNow[1];
1166                         t21 -= a01 * aboveNow[2];
1167                         t31 -= a01 * aboveNow[3];
1168                         aUnderNow += BLOCK;
1169                         aboveNow += BLOCK;
1170                    }
1171                    aa[i+0*BLOCK] = t00;
1172                    aa[i+1*BLOCK] = t10;
1173                    aa[i+2*BLOCK] = t20;
1174                    aa[i+3*BLOCK] = t30;
1175                    aa[i+1+0*BLOCK] = t01;
1176                    aa[i+1+1*BLOCK] = t11;
1177                    aa[i+1+2*BLOCK] = t21;
1178                    aa[i+1+3*BLOCK] = t31;
1179               }
1180               if (odd) {
1181                    CoinWorkDouble t0 = aa[n+0*BLOCK];
1182                    CoinWorkDouble t1 = aa[n+1*BLOCK];
1183                    CoinWorkDouble t2 = aa[n+2*BLOCK];
1184                    CoinWorkDouble t3 = aa[n+3*BLOCK];
1185                    CoinWorkDouble a0;
1186                    for (k = 0; k < BLOCK; k++) {
1187                         a0 = aUnder[n+k*BLOCK] * work[k];
1188                         t0 -= a0 * above[j + 0 + k * BLOCK];
1189                         t1 -= a0 * above[j + 1 + k * BLOCK];
1190                         t2 -= a0 * above[j + 2 + k * BLOCK];
1191                         t3 -= a0 * above[j + 3 + k * BLOCK];
1192                    }
1193                    aa[n+0*BLOCK] = t0;
1194                    aa[n+1*BLOCK] = t1;
1195                    aa[n+2*BLOCK] = t2;
1196                    aa[n+3*BLOCK] = t3;
1197               }
1198          }
1199     }
1200#else
1201     aa = aOther - BLOCK;
1202     for (j = 0; j < BLOCK; j ++) {
1203          aa += BLOCK;
1204          for (i = 0; i < nUnder; i++) {
1205               CoinWorkDouble t00 = aa[i+0*BLOCK];
1206               for (k = 0; k < BLOCK; k++) {
1207                    CoinWorkDouble a00 = aUnder[i+k*BLOCK] * work[k];
1208                    t00 -= a00 * above[j + k * BLOCK];
1209               }
1210               aa[i] = t00;
1211          }
1212     }
1213#endif
1214}
1215/* Uses factorization to solve. */
1216void
1217ClpCholeskyDense::solve (CoinWorkDouble * region)
1218{
1219#ifdef CHOL_COMPARE
1220     double * region2 = NULL;
1221     if (numberRows_ < 200) {
1222          longDouble * xx = sparseFactor_;
1223          longDouble * yy = diagonal_;
1224          diagonal_ = sparseFactor_ + 40000;
1225          sparseFactor_ = diagonal_ + numberRows_;
1226          region2 = ClpCopyOfArray(region, numberRows_);
1227          int iRow, iColumn;
1228          int addOffset = numberRows_ - 1;
1229          longDouble * work = sparseFactor_ - 1;
1230          for (iColumn = 0; iColumn < numberRows_; iColumn++) {
1231               double value = region2[iColumn];
1232               for (iRow = iColumn + 1; iRow < numberRows_; iRow++)
1233                    region2[iRow] -= value * work[iRow];
1234               addOffset--;
1235               work += addOffset;
1236          }
1237          for (iColumn = numberRows_ - 1; iColumn >= 0; iColumn--) {
1238               double value = region2[iColumn] * diagonal_[iColumn];
1239               work -= addOffset;
1240               addOffset++;
1241               for (iRow = iColumn + 1; iRow < numberRows_; iRow++)
1242                    value -= region2[iRow] * work[iRow];
1243               region2[iColumn] = value;
1244          }
1245          sparseFactor_ = xx;
1246          diagonal_ = yy;
1247     }
1248#endif
1249     /*longDouble * xx = sparseFactor_;*/
1250     /*longDouble * yy = diagonal_;*/
1251     /*diagonal_ = sparseFactor_ + 40000;*/
1252     /*sparseFactor_=diagonal_ + numberRows_;*/
1253     int numberBlocks = (numberRows_ + BLOCK - 1) >> BLOCKSHIFT;
1254     /* later align on boundary*/
1255     longDouble * a = sparseFactor_ + BLOCKSQ * numberBlocks;
1256     int iBlock;
1257     longDouble * aa = a;
1258     int iColumn;
1259     for (iBlock = 0; iBlock < numberBlocks; iBlock++) {
1260          int nChunk;
1261          int jBlock;
1262          int iDo = iBlock * BLOCK;
1263          int base = iDo;
1264          if (iDo + BLOCK > numberRows_) {
1265               nChunk = numberRows_ - iDo;
1266          } else {
1267               nChunk = BLOCK;
1268          }
1269          solveF1(aa, nChunk, region + iDo);
1270          for (jBlock = iBlock + 1; jBlock < numberBlocks; jBlock++) {
1271               base += BLOCK;
1272               aa += BLOCKSQ;
1273               if (base + BLOCK > numberRows_) {
1274                    nChunk = numberRows_ - base;
1275               } else {
1276                    nChunk = BLOCK;
1277               }
1278               solveF2(aa, nChunk, region + iDo, region + base);
1279          }
1280          aa += BLOCKSQ;
1281     }
1282     /* do diagonal outside*/
1283     for (iColumn = 0; iColumn < numberRows_; iColumn++)
1284          region[iColumn] *= diagonal_[iColumn];
1285     int offset = ((numberBlocks * (numberBlocks + 1)) >> 1);
1286     aa = a + number_entries(offset - 1);
1287     int lBase = (numberBlocks - 1) * BLOCK;
1288     for (iBlock = numberBlocks - 1; iBlock >= 0; iBlock--) {
1289          int nChunk;
1290          int jBlock;
1291          int triBase = iBlock * BLOCK;
1292          int iBase = lBase;
1293          for (jBlock = iBlock + 1; jBlock < numberBlocks; jBlock++) {
1294               if (iBase + BLOCK > numberRows_) {
1295                    nChunk = numberRows_ - iBase;
1296               } else {
1297                    nChunk = BLOCK;
1298               }
1299               solveB2(aa, nChunk, region + triBase, region + iBase);
1300               iBase -= BLOCK;
1301               aa -= BLOCKSQ;
1302          }
1303          if (triBase + BLOCK > numberRows_) {
1304               nChunk = numberRows_ - triBase;
1305          } else {
1306               nChunk = BLOCK;
1307          }
1308          solveB1(aa, nChunk, region + triBase);
1309          aa -= BLOCKSQ;
1310     }
1311#ifdef CHOL_COMPARE
1312     if (numberRows_ < 200) {
1313          for (int i = 0; i < numberRows_; i++) {
1314               assert(CoinAbs(region[i] - region2[i]) < 1.0e-3);
1315          }
1316          delete [] region2;
1317     }
1318#endif
1319}
1320/* Forward part of solve 1*/
1321void
1322ClpCholeskyDense::solveF1(longDouble * a, int n, CoinWorkDouble * region)
1323{
1324     int j, k;
1325     CoinWorkDouble t00;
1326     for (j = 0; j < n; j ++) {
1327          t00 = region[j];
1328          for (k = 0; k < j; ++k) {
1329               t00 -= region[k] * a[j + k * BLOCK];
1330          }
1331          /*t00*=a[j + j * BLOCK];*/
1332          region[j] = t00;
1333     }
1334}
1335/* Forward part of solve 2*/
1336void
1337ClpCholeskyDense::solveF2(longDouble * a, int n, CoinWorkDouble * region, CoinWorkDouble * region2)
1338{
1339     int j, k;
1340#ifdef BLOCKUNROLL
1341     if (n == BLOCK) {
1342          for (k = 0; k < BLOCK; k += 4) {
1343               CoinWorkDouble t0 = region2[0];
1344               CoinWorkDouble t1 = region2[1];
1345               CoinWorkDouble t2 = region2[2];
1346               CoinWorkDouble t3 = region2[3];
1347               t0 -= region[0] * a[0 + 0 * BLOCK];
1348               t1 -= region[0] * a[1 + 0 * BLOCK];
1349               t2 -= region[0] * a[2 + 0 * BLOCK];
1350               t3 -= region[0] * a[3 + 0 * BLOCK];
1351
1352               t0 -= region[1] * a[0 + 1 * BLOCK];
1353               t1 -= region[1] * a[1 + 1 * BLOCK];
1354               t2 -= region[1] * a[2 + 1 * BLOCK];
1355               t3 -= region[1] * a[3 + 1 * BLOCK];
1356
1357               t0 -= region[2] * a[0 + 2 * BLOCK];
1358               t1 -= region[2] * a[1 + 2 * BLOCK];
1359               t2 -= region[2] * a[2 + 2 * BLOCK];
1360               t3 -= region[2] * a[3 + 2 * BLOCK];
1361
1362               t0 -= region[3] * a[0 + 3 * BLOCK];
1363               t1 -= region[3] * a[1 + 3 * BLOCK];
1364               t2 -= region[3] * a[2 + 3 * BLOCK];
1365               t3 -= region[3] * a[3 + 3 * BLOCK];
1366
1367               t0 -= region[4] * a[0 + 4 * BLOCK];
1368               t1 -= region[4] * a[1 + 4 * BLOCK];
1369               t2 -= region[4] * a[2 + 4 * BLOCK];
1370               t3 -= region[4] * a[3 + 4 * BLOCK];
1371
1372               t0 -= region[5] * a[0 + 5 * BLOCK];
1373               t1 -= region[5] * a[1 + 5 * BLOCK];
1374               t2 -= region[5] * a[2 + 5 * BLOCK];
1375               t3 -= region[5] * a[3 + 5 * BLOCK];
1376
1377               t0 -= region[6] * a[0 + 6 * BLOCK];
1378               t1 -= region[6] * a[1 + 6 * BLOCK];
1379               t2 -= region[6] * a[2 + 6 * BLOCK];
1380               t3 -= region[6] * a[3 + 6 * BLOCK];
1381
1382               t0 -= region[7] * a[0 + 7 * BLOCK];
1383               t1 -= region[7] * a[1 + 7 * BLOCK];
1384               t2 -= region[7] * a[2 + 7 * BLOCK];
1385               t3 -= region[7] * a[3 + 7 * BLOCK];
1386#if BLOCK>8
1387               t0 -= region[8] * a[0 + 8 * BLOCK];
1388               t1 -= region[8] * a[1 + 8 * BLOCK];
1389               t2 -= region[8] * a[2 + 8 * BLOCK];
1390               t3 -= region[8] * a[3 + 8 * BLOCK];
1391
1392               t0 -= region[9] * a[0 + 9 * BLOCK];
1393               t1 -= region[9] * a[1 + 9 * BLOCK];
1394               t2 -= region[9] * a[2 + 9 * BLOCK];
1395               t3 -= region[9] * a[3 + 9 * BLOCK];
1396
1397               t0 -= region[10] * a[0 + 10 * BLOCK];
1398               t1 -= region[10] * a[1 + 10 * BLOCK];
1399               t2 -= region[10] * a[2 + 10 * BLOCK];
1400               t3 -= region[10] * a[3 + 10 * BLOCK];
1401
1402               t0 -= region[11] * a[0 + 11 * BLOCK];
1403               t1 -= region[11] * a[1 + 11 * BLOCK];
1404               t2 -= region[11] * a[2 + 11 * BLOCK];
1405               t3 -= region[11] * a[3 + 11 * BLOCK];
1406
1407               t0 -= region[12] * a[0 + 12 * BLOCK];
1408               t1 -= region[12] * a[1 + 12 * BLOCK];
1409               t2 -= region[12] * a[2 + 12 * BLOCK];
1410               t3 -= region[12] * a[3 + 12 * BLOCK];
1411
1412               t0 -= region[13] * a[0 + 13 * BLOCK];
1413               t1 -= region[13] * a[1 + 13 * BLOCK];
1414               t2 -= region[13] * a[2 + 13 * BLOCK];
1415               t3 -= region[13] * a[3 + 13 * BLOCK];
1416
1417               t0 -= region[14] * a[0 + 14 * BLOCK];
1418               t1 -= region[14] * a[1 + 14 * BLOCK];
1419               t2 -= region[14] * a[2 + 14 * BLOCK];
1420               t3 -= region[14] * a[3 + 14 * BLOCK];
1421
1422               t0 -= region[15] * a[0 + 15 * BLOCK];
1423               t1 -= region[15] * a[1 + 15 * BLOCK];
1424               t2 -= region[15] * a[2 + 15 * BLOCK];
1425               t3 -= region[15] * a[3 + 15 * BLOCK];
1426#endif
1427               region2[0] = t0;
1428               region2[1] = t1;
1429               region2[2] = t2;
1430               region2[3] = t3;
1431               region2 += 4;
1432               a += 4;
1433          }
1434     } else {
1435#endif
1436          for (k = 0; k < n; ++k) {
1437               CoinWorkDouble t00 = region2[k];
1438               for (j = 0; j < BLOCK; j ++) {
1439                    t00 -= region[j] * a[k + j * BLOCK];
1440               }
1441               region2[k] = t00;
1442          }
1443#ifdef BLOCKUNROLL
1444     }
1445#endif
1446}
1447/* Backward part of solve 1*/
1448void
1449ClpCholeskyDense::solveB1(longDouble * a, int n, CoinWorkDouble * region)
1450{
1451     int j, k;
1452     CoinWorkDouble t00;
1453     for (j = n - 1; j >= 0; j --) {
1454          t00 = region[j];
1455          for (k = j + 1; k < n; ++k) {
1456               t00 -= region[k] * a[k + j * BLOCK];
1457          }
1458          /*t00*=a[j + j * BLOCK];*/
1459          region[j] = t00;
1460     }
1461}
1462/* Backward part of solve 2*/
1463void
1464ClpCholeskyDense::solveB2(longDouble * a, int n, CoinWorkDouble * region, CoinWorkDouble * region2)
1465{
1466     int j, k;
1467#ifdef BLOCKUNROLL
1468     if (n == BLOCK) {
1469          for (j = 0; j < BLOCK; j += 4) {
1470               CoinWorkDouble t0 = region[0];
1471               CoinWorkDouble t1 = region[1];
1472               CoinWorkDouble t2 = region[2];
1473               CoinWorkDouble t3 = region[3];
1474               t0 -= region2[0] * a[0 + 0*BLOCK];
1475               t1 -= region2[0] * a[0 + 1*BLOCK];
1476               t2 -= region2[0] * a[0 + 2*BLOCK];
1477               t3 -= region2[0] * a[0 + 3*BLOCK];
1478
1479               t0 -= region2[1] * a[1 + 0*BLOCK];
1480               t1 -= region2[1] * a[1 + 1*BLOCK];
1481               t2 -= region2[1] * a[1 + 2*BLOCK];
1482               t3 -= region2[1] * a[1 + 3*BLOCK];
1483
1484               t0 -= region2[2] * a[2 + 0*BLOCK];
1485               t1 -= region2[2] * a[2 + 1*BLOCK];
1486               t2 -= region2[2] * a[2 + 2*BLOCK];
1487               t3 -= region2[2] * a[2 + 3*BLOCK];
1488
1489               t0 -= region2[3] * a[3 + 0*BLOCK];
1490               t1 -= region2[3] * a[3 + 1*BLOCK];
1491               t2 -= region2[3] * a[3 + 2*BLOCK];
1492               t3 -= region2[3] * a[3 + 3*BLOCK];
1493
1494               t0 -= region2[4] * a[4 + 0*BLOCK];
1495               t1 -= region2[4] * a[4 + 1*BLOCK];
1496               t2 -= region2[4] * a[4 + 2*BLOCK];
1497               t3 -= region2[4] * a[4 + 3*BLOCK];
1498
1499               t0 -= region2[5] * a[5 + 0*BLOCK];
1500               t1 -= region2[5] * a[5 + 1*BLOCK];
1501               t2 -= region2[5] * a[5 + 2*BLOCK];
1502               t3 -= region2[5] * a[5 + 3*BLOCK];
1503
1504               t0 -= region2[6] * a[6 + 0*BLOCK];
1505               t1 -= region2[6] * a[6 + 1*BLOCK];
1506               t2 -= region2[6] * a[6 + 2*BLOCK];
1507               t3 -= region2[6] * a[6 + 3*BLOCK];
1508
1509               t0 -= region2[7] * a[7 + 0*BLOCK];
1510               t1 -= region2[7] * a[7 + 1*BLOCK];
1511               t2 -= region2[7] * a[7 + 2*BLOCK];
1512               t3 -= region2[7] * a[7 + 3*BLOCK];
1513#if BLOCK>8
1514
1515               t0 -= region2[8] * a[8 + 0*BLOCK];
1516               t1 -= region2[8] * a[8 + 1*BLOCK];
1517               t2 -= region2[8] * a[8 + 2*BLOCK];
1518               t3 -= region2[8] * a[8 + 3*BLOCK];
1519
1520               t0 -= region2[9] * a[9 + 0*BLOCK];
1521               t1 -= region2[9] * a[9 + 1*BLOCK];
1522               t2 -= region2[9] * a[9 + 2*BLOCK];
1523               t3 -= region2[9] * a[9 + 3*BLOCK];
1524
1525               t0 -= region2[10] * a[10 + 0*BLOCK];
1526               t1 -= region2[10] * a[10 + 1*BLOCK];
1527               t2 -= region2[10] * a[10 + 2*BLOCK];
1528               t3 -= region2[10] * a[10 + 3*BLOCK];
1529
1530               t0 -= region2[11] * a[11 + 0*BLOCK];
1531               t1 -= region2[11] * a[11 + 1*BLOCK];
1532               t2 -= region2[11] * a[11 + 2*BLOCK];
1533               t3 -= region2[11] * a[11 + 3*BLOCK];
1534
1535               t0 -= region2[12] * a[12 + 0*BLOCK];
1536               t1 -= region2[12] * a[12 + 1*BLOCK];
1537               t2 -= region2[12] * a[12 + 2*BLOCK];
1538               t3 -= region2[12] * a[12 + 3*BLOCK];
1539
1540               t0 -= region2[13] * a[13 + 0*BLOCK];
1541               t1 -= region2[13] * a[13 + 1*BLOCK];
1542               t2 -= region2[13] * a[13 + 2*BLOCK];
1543               t3 -= region2[13] * a[13 + 3*BLOCK];
1544
1545               t0 -= region2[14] * a[14 + 0*BLOCK];
1546               t1 -= region2[14] * a[14 + 1*BLOCK];
1547               t2 -= region2[14] * a[14 + 2*BLOCK];
1548               t3 -= region2[14] * a[14 + 3*BLOCK];
1549
1550               t0 -= region2[15] * a[15 + 0*BLOCK];
1551               t1 -= region2[15] * a[15 + 1*BLOCK];
1552               t2 -= region2[15] * a[15 + 2*BLOCK];
1553               t3 -= region2[15] * a[15 + 3*BLOCK];
1554#endif
1555               region[0] = t0;
1556               region[1] = t1;
1557               region[2] = t2;
1558               region[3] = t3;
1559               a += 4 * BLOCK;
1560               region += 4;
1561          }
1562     } else {
1563#endif
1564          for (j = 0; j < BLOCK; j ++) {
1565               CoinWorkDouble t00 = region[j];
1566               for (k = 0; k < n; ++k) {
1567                    t00 -= region2[k] * a[k + j * BLOCK];
1568               }
1569               region[j] = t00;
1570          }
1571#ifdef BLOCKUNROLL
1572     }
1573#endif
1574}
Note: See TracBrowser for help on using the repository browser.