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

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

minor changes to implement Aboca

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