source: trunk/Clp/src/ClpCholeskyWssmp.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 Id
File size: 25.6 KB
Line 
1/* $Id: ClpCholeskyWssmp.cpp 1525 2010-02-26 17:27:59Z mjs $ */
2#ifdef WSSMP_BARRIER
3// Copyright (C) 2003, International Business Machines
4// Corporation and others.  All Rights Reserved.
5
6
7
8#include "CoinPragma.hpp"
9#include "CoinHelperFunctions.hpp"
10#include "ClpHelperFunctions.hpp"
11
12#include "ClpInterior.hpp"
13#include "ClpCholeskyWssmp.hpp"
14#include "ClpCholeskyDense.hpp"
15#include "ClpMessage.hpp"
16
17//#############################################################################
18// Constructors / Destructor / Assignment
19//#############################################################################
20
21//-------------------------------------------------------------------
22// Default Constructor
23//-------------------------------------------------------------------
24ClpCholeskyWssmp::ClpCholeskyWssmp (int denseThreshold)
25     : ClpCholeskyBase(denseThreshold)
26{
27     type_ = 12;
28}
29
30//-------------------------------------------------------------------
31// Copy constructor
32//-------------------------------------------------------------------
33ClpCholeskyWssmp::ClpCholeskyWssmp (const ClpCholeskyWssmp & rhs)
34     : ClpCholeskyBase(rhs)
35{
36}
37
38
39//-------------------------------------------------------------------
40// Destructor
41//-------------------------------------------------------------------
42ClpCholeskyWssmp::~ClpCholeskyWssmp ()
43{
44}
45
46//----------------------------------------------------------------
47// Assignment operator
48//-------------------------------------------------------------------
49ClpCholeskyWssmp &
50ClpCholeskyWssmp::operator=(const ClpCholeskyWssmp& rhs)
51{
52     if (this != &rhs) {
53          ClpCholeskyBase::operator=(rhs);
54     }
55     return *this;
56}
57//-------------------------------------------------------------------
58// Clone
59//-------------------------------------------------------------------
60ClpCholeskyBase * ClpCholeskyWssmp::clone() const
61{
62     return new ClpCholeskyWssmp(*this);
63}
64// At present I can't get wssmp to work as my libraries seem to be out of sync
65// so I have linked in ekkwssmp which is an older version
66#if WSSMP_BARRIER==2
67extern "C" void wssmp(int * n,
68                      int * columnStart , int * rowIndex , double * element,
69                      double * diagonal , int * perm , int * invp ,
70                      double * rhs , int * ldb , int * nrhs ,
71                      double * aux , int * naux ,
72                      int   * mrp , int * iparm , double * dparm);
73extern "C" void wsetmaxthrds(int *);
74#elif WSSMP_BARRIER==3
75extern "C" void wssmp_(int * n,
76                       int * columnStart , int * rowIndex , double * element,
77                       double * diagonal , int * perm , int * invp ,
78                       double * rhs , int * ldb , int * nrhs ,
79                       double * aux , int * naux ,
80                       int   * mrp , int * iparm , double * dparm);
81extern "C" void wsetmaxthrds_(int *);
82static void wssmp( int *n, int *ia, int *ja,
83                   double *avals, double *diag, int *perm, int *invp,
84                   double *b, int *ldb, int *nrhs, double *aux, int *
85                   naux, int *mrp, int *iparm, double *dparm)
86{
87     wssmp_(n, ia, ja,
88            avals, diag, perm, invp,
89            b, ldb, nrhs, aux,
90            naux, mrp, iparm, dparm);
91}
92static void wsetmaxthrds(int * n)
93{
94     wsetmaxthrds_(n);
95}
96#else
97/* minimum needed for user */
98typedef struct EKKModel EKKModel;
99typedef struct EKKContext EKKContext;
100
101
102extern "C" {
103     EKKContext *  ekk_initializeContext();
104     void ekk_endContext(EKKContext * context);
105     EKKModel *  ekk_newModel(EKKContext * env, const char * name);
106     int ekk_deleteModel(EKKModel * model);
107}
108static  EKKModel * model = NULL;
109static  EKKContext * context = NULL;
110extern "C" void ekkwssmp(EKKModel *, int * n,
111                         int * columnStart , int * rowIndex , double * element,
112                         double * diagonal , int * perm , int * invp ,
113                         double * rhs , int * ldb , int * nrhs ,
114                         double * aux , int * naux ,
115                         int   * mrp , int * iparm , double * dparm);
116static void wssmp( int *n, int *ia, int *ja,
117                   double *avals, double *diag, int *perm, int *invp,
118                   double *b, int *ldb, int *nrhs, double *aux, int *
119                   naux, int *mrp, int *iparm, double *dparm)
120{
121     if (!context) {
122          /* initialize OSL environment */
123          context = ekk_initializeContext();
124          model = ekk_newModel(context, "");
125     }
126     ekkwssmp(model, n, ia, ja,
127              avals, diag, perm, invp,
128              b, ldb, nrhs, aux,
129              naux, mrp, iparm, dparm);
130     //ekk_deleteModel(model);
131     //ekk_endContext(context);
132}
133#endif
134/* Orders rows and saves pointer to matrix.and model */
135int
136ClpCholeskyWssmp::order(ClpInterior * model)
137{
138     numberRows_ = model->numberRows();
139     rowsDropped_ = new char [numberRows_];
140     memset(rowsDropped_, 0, numberRows_);
141     numberRowsDropped_ = 0;
142     model_ = model;
143     rowCopy_ = model->clpMatrix()->reverseOrderedCopy();
144     // Space for starts
145     choleskyStart_ = new CoinBigIndex[numberRows_+1];
146     const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
147     const int * columnLength = model_->clpMatrix()->getVectorLengths();
148     const int * row = model_->clpMatrix()->getIndices();
149     const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
150     const int * rowLength = rowCopy_->getVectorLengths();
151     const int * column = rowCopy_->getIndices();
152     // We need two arrays for counts
153     int * which = new int [numberRows_];
154     int * used = new int[numberRows_+1];
155     CoinZeroN(used, numberRows_);
156     int iRow;
157     sizeFactor_ = 0;
158     int numberColumns = model->numberColumns();
159     int numberDense = 0;
160     if (denseThreshold_ > 0) {
161          delete [] whichDense_;
162          delete [] denseColumn_;
163          delete dense_;
164          whichDense_ = new char[numberColumns];
165          int iColumn;
166          used[numberRows_] = 0;
167          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
168               int length = columnLength[iColumn];
169               used[length] += 1;
170          }
171          int nLong = 0;
172          int stop = CoinMax(denseThreshold_ / 2, 100);
173          for (iRow = numberRows_; iRow >= stop; iRow--) {
174               if (used[iRow])
175                    printf("%d columns are of length %d\n", used[iRow], iRow);
176               nLong += used[iRow];
177               if (nLong > 50 || nLong > (numberColumns >> 2))
178                    break;
179          }
180          CoinZeroN(used, numberRows_);
181          for (iColumn = 0; iColumn < numberColumns; iColumn++) {
182               if (columnLength[iColumn] < denseThreshold_) {
183                    whichDense_[iColumn] = 0;
184               } else {
185                    whichDense_[iColumn] = 1;
186                    numberDense++;
187               }
188          }
189          if (!numberDense || numberDense > 100) {
190               // free
191               delete [] whichDense_;
192               whichDense_ = NULL;
193               denseColumn_ = NULL;
194               dense_ = NULL;
195          } else {
196               // space for dense columns
197               denseColumn_ = new double [numberDense*numberRows_];
198               // dense cholesky
199               dense_ = new ClpCholeskyDense();
200               dense_->reserveSpace(NULL, numberDense);
201               printf("Taking %d columns as dense\n", numberDense);
202          }
203     }
204     for (iRow = 0; iRow < numberRows_; iRow++) {
205          int number = 1;
206          // make sure diagonal exists
207          which[0] = iRow;
208          used[iRow] = 1;
209          if (!rowsDropped_[iRow]) {
210               CoinBigIndex startRow = rowStart[iRow];
211               CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
212               for (CoinBigIndex k = startRow; k < endRow; k++) {
213                    int iColumn = column[k];
214                    if (!whichDense_ || !whichDense_[iColumn]) {
215                         CoinBigIndex start = columnStart[iColumn];
216                         CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
217                         for (CoinBigIndex j = start; j < end; j++) {
218                              int jRow = row[j];
219                              if (jRow >= iRow && !rowsDropped_[jRow]) {
220                                   if (!used[jRow]) {
221                                        used[jRow] = 1;
222                                        which[number++] = jRow;
223                                   }
224                              }
225                         }
226                    }
227               }
228               sizeFactor_ += number;
229               int j;
230               for (j = 0; j < number; j++)
231                    used[which[j]] = 0;
232          }
233     }
234     delete [] which;
235     // Now we have size - create arrays and fill in
236     try {
237          choleskyRow_ = new int [sizeFactor_];
238     } catch (...) {
239          // no memory
240          delete [] choleskyStart_;
241          choleskyStart_ = NULL;
242          return -1;
243     }
244     try {
245          sparseFactor_ = new double[sizeFactor_];
246     } catch (...) {
247          // no memory
248          delete [] choleskyRow_;
249          choleskyRow_ = NULL;
250          delete [] choleskyStart_;
251          choleskyStart_ = NULL;
252          return -1;
253     }
254
255     sizeFactor_ = 0;
256     which = choleskyRow_;
257     for (iRow = 0; iRow < numberRows_; iRow++) {
258          int number = 1;
259          // make sure diagonal exists
260          which[0] = iRow;
261          used[iRow] = 1;
262          choleskyStart_[iRow] = sizeFactor_;
263          if (!rowsDropped_[iRow]) {
264               CoinBigIndex startRow = rowStart[iRow];
265               CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
266               for (CoinBigIndex k = startRow; k < endRow; k++) {
267                    int iColumn = column[k];
268                    if (!whichDense_ || !whichDense_[iColumn]) {
269                         CoinBigIndex start = columnStart[iColumn];
270                         CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
271                         for (CoinBigIndex j = start; j < end; j++) {
272                              int jRow = row[j];
273                              if (jRow >= iRow && !rowsDropped_[jRow]) {
274                                   if (!used[jRow]) {
275                                        used[jRow] = 1;
276                                        which[number++] = jRow;
277                                   }
278                              }
279                         }
280                    }
281               }
282               sizeFactor_ += number;
283               int j;
284               for (j = 0; j < number; j++)
285                    used[which[j]] = 0;
286               // Sort
287               std::sort(which, which + number);
288               // move which on
289               which += number;
290          }
291     }
292     choleskyStart_[numberRows_] = sizeFactor_;
293     delete [] used;
294     permuteInverse_ = new int [numberRows_];
295     permute_ = new int[numberRows_];
296     integerParameters_[0] = 0;
297     int i0 = 0;
298     int i1 = 1;
299#if WSSMP_BARRIER>=2
300     int i2 = 1;
301     if (model->numberThreads() <= 0)
302          i2 = 1;
303     else
304          i2 = model->numberThreads();
305     wsetmaxthrds(&i2);
306#endif
307     wssmp(&numberRows_, choleskyStart_, choleskyRow_, sparseFactor_,
308           NULL, permute_, permuteInverse_, 0, &numberRows_, &i1,
309           NULL, &i0, NULL, integerParameters_, doubleParameters_);
310     integerParameters_[1] = 1; //order and symbolic
311     integerParameters_[2] = 2;
312     integerParameters_[3] = 0; //CSR
313     integerParameters_[4] = 0; //C style
314     integerParameters_[13] = 1; //reuse initial factorization space
315     integerParameters_[15+0] = 1; //ordering
316     integerParameters_[15+1] = 0;
317     integerParameters_[15+2] = 1;
318     integerParameters_[15+3] = 0;
319     integerParameters_[15+4] = 1;
320     doubleParameters_[10] = 1.0e-20;
321     doubleParameters_[11] = 1.0e-15;
322     wssmp(&numberRows_, choleskyStart_, choleskyRow_, sparseFactor_,
323           NULL, permute_, permuteInverse_, NULL, &numberRows_, &i1,
324           NULL, &i0, NULL, integerParameters_, doubleParameters_);
325     //std::cout<<"Ordering and symbolic factorization took "<<doubleParameters_[0]<<std::endl;
326     if (integerParameters_[63]) {
327          std::cout << "wssmp returning error code of " << integerParameters_[63] << std::endl;
328          return 1;
329     }
330     // Modify gamma and delta if 0.0
331     if (!model->gamma() && !model->delta()) {
332          model->setGamma(5.0e-5);
333          model->setDelta(5.0e-5);
334     }
335     std::cout << integerParameters_[23] << " elements in sparse Cholesky" << std::endl;
336     if (!integerParameters_[23]) {
337          for (int iRow = 0; iRow < numberRows_; iRow++) {
338               permuteInverse_[iRow] = iRow;
339               permute_[iRow] = iRow;
340          }
341          std::cout << "wssmp says no elements - fully dense? - switching to dense" << std::endl;
342          integerParameters_[1] = 2;
343          integerParameters_[2] = 2;
344          integerParameters_[7] = 1; // no permute
345          wssmp(&numberRows_, choleskyStart_, choleskyRow_, sparseFactor_,
346                NULL, permute_, permuteInverse_, NULL, &numberRows_, &i1,
347                NULL, &i0, NULL, integerParameters_, doubleParameters_);
348          std::cout << integerParameters_[23] << " elements in dense Cholesky" << std::endl;
349     }
350     return 0;
351}
352/* Does Symbolic factorization given permutation.
353   This is called immediately after order.  If user provides this then
354   user must provide factorize and solve.  Otherwise the default factorization is used
355   returns non-zero if not enough memory */
356int
357ClpCholeskyWssmp::symbolic()
358{
359     return 0;
360}
361/* Factorize - filling in rowsDropped and returning number dropped */
362int
363ClpCholeskyWssmp::factorize(const double * diagonal, int * rowsDropped)
364{
365     const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
366     const int * columnLength = model_->clpMatrix()->getVectorLengths();
367     const int * row = model_->clpMatrix()->getIndices();
368     const double * element = model_->clpMatrix()->getElements();
369     const CoinBigIndex * rowStart = rowCopy_->getVectorStarts();
370     const int * rowLength = rowCopy_->getVectorLengths();
371     const int * column = rowCopy_->getIndices();
372     const double * elementByRow = rowCopy_->getElements();
373     int numberColumns = model_->clpMatrix()->getNumCols();
374     int iRow;
375     double * work = new double[numberRows_];
376     CoinZeroN(work, numberRows_);
377     const double * diagonalSlack = diagonal + numberColumns;
378     int newDropped = 0;
379     double largest;
380     double smallest;
381     int numberDense = 0;
382     if (dense_)
383          numberDense = dense_->numberRows();
384     //perturbation
385     double perturbation = model_->diagonalPerturbation() * model_->diagonalNorm();
386     perturbation = perturbation * perturbation;
387     if (perturbation > 1.0) {
388#ifdef COIN_DEVELOP
389          //if (model_->model()->logLevel()&4)
390          std::cout << "large perturbation " << perturbation << std::endl;
391#endif
392          perturbation = sqrt(perturbation);;
393          perturbation = 1.0;
394     }
395     if (whichDense_) {
396          double * denseDiagonal = dense_->diagonal();
397          double * dense = denseColumn_;
398          int iDense = 0;
399          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
400               if (whichDense_[iColumn]) {
401                    denseDiagonal[iDense++] = 1.0 / diagonal[iColumn];
402                    CoinZeroN(dense, numberRows_);
403                    CoinBigIndex start = columnStart[iColumn];
404                    CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
405                    for (CoinBigIndex j = start; j < end; j++) {
406                         int jRow = row[j];
407                         dense[jRow] = element[j];
408                    }
409                    dense += numberRows_;
410               }
411          }
412     }
413     double delta2 = model_->delta(); // add delta*delta to diagonal
414     delta2 *= delta2;
415     for (iRow = 0; iRow < numberRows_; iRow++) {
416          double * put = sparseFactor_ + choleskyStart_[iRow];
417          int * which = choleskyRow_ + choleskyStart_[iRow];
418          int number = choleskyStart_[iRow+1] - choleskyStart_[iRow];
419          if (!rowLength[iRow])
420               rowsDropped_[iRow] = 1;
421          if (!rowsDropped_[iRow]) {
422               CoinBigIndex startRow = rowStart[iRow];
423               CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow];
424               work[iRow] = diagonalSlack[iRow] + delta2;
425               for (CoinBigIndex k = startRow; k < endRow; k++) {
426                    int iColumn = column[k];
427                    if (!whichDense_ || !whichDense_[iColumn]) {
428                         CoinBigIndex start = columnStart[iColumn];
429                         CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
430                         double multiplier = diagonal[iColumn] * elementByRow[k];
431                         for (CoinBigIndex j = start; j < end; j++) {
432                              int jRow = row[j];
433                              if (jRow >= iRow && !rowsDropped_[jRow]) {
434                                   double value = element[j] * multiplier;
435                                   work[jRow] += value;
436                              }
437                         }
438                    }
439               }
440               int j;
441               for (j = 0; j < number; j++) {
442                    int jRow = which[j];
443                    put[j] = work[jRow];
444                    work[jRow] = 0.0;
445               }
446          } else {
447               // dropped
448               int j;
449               for (j = 1; j < number; j++) {
450                    put[j] = 0.0;
451               }
452               put[0] = 1.0;
453          }
454     }
455     //check sizes
456     double largest2 = maximumAbsElement(sparseFactor_, sizeFactor_);
457     largest2 *= 1.0e-20;
458     largest = CoinMin(largest2, 1.0e-11);
459     int numberDroppedBefore = 0;
460     for (iRow = 0; iRow < numberRows_; iRow++) {
461          int dropped = rowsDropped_[iRow];
462          // Move to int array
463          rowsDropped[iRow] = dropped;
464          if (!dropped) {
465               CoinBigIndex start = choleskyStart_[iRow];
466               double diagonal = sparseFactor_[start];
467               if (diagonal > largest2) {
468                    sparseFactor_[start] = diagonal + perturbation;
469               } else {
470                    sparseFactor_[start] = diagonal + perturbation;
471                    rowsDropped[iRow] = 2;
472                    numberDroppedBefore++;
473               }
474          }
475     }
476     int i1 = 1;
477     int i0 = 0;
478     integerParameters_[1] = 3;
479     integerParameters_[2] = 3;
480     integerParameters_[10] = 2;
481     //integerParameters_[11]=1;
482     integerParameters_[12] = 2;
483     doubleParameters_[10] = CoinMax(1.0e-20, largest);
484     if (doubleParameters_[10] > 1.0e-3)
485          integerParameters_[9] = 1;
486     else
487          integerParameters_[9] = 0;
488     wssmp(&numberRows_, choleskyStart_, choleskyRow_, sparseFactor_,
489           NULL, permute_, permuteInverse_, NULL, &numberRows_, &i1,
490           NULL, &i0, rowsDropped, integerParameters_, doubleParameters_);
491     //std::cout<<"factorization took "<<doubleParameters_[0]<<std::endl;
492     if (integerParameters_[9]) {
493          std::cout << "scaling applied" << std::endl;
494     }
495     newDropped = integerParameters_[20] + numberDroppedBefore;
496     //if (integerParameters_[20])
497     //std::cout<<integerParameters_[20]<<" rows dropped"<<std::endl;
498     largest = doubleParameters_[3];
499     smallest = doubleParameters_[4];
500     delete [] work;
501     if (model_->messageHandler()->logLevel() > 1)
502          std::cout << "Cholesky - largest " << largest << " smallest " << smallest << std::endl;
503     choleskyCondition_ = largest / smallest;
504     if (integerParameters_[63] < 0)
505          return -1; // out of memory
506     if (whichDense_) {
507          // Update dense columns (just L)
508          // Zero out dropped rows
509          int i;
510          for (i = 0; i < numberDense; i++) {
511               double * a = denseColumn_ + i * numberRows_;
512               for (int j = 0; j < numberRows_; j++) {
513                    if (rowsDropped[j])
514                         a[j] = 0.0;
515               }
516          }
517          integerParameters_[29] = 1;
518          int i0 = 0;
519          integerParameters_[1] = 4;
520          integerParameters_[2] = 4;
521          wssmp(&numberRows_, choleskyStart_, choleskyRow_, sparseFactor_,
522                NULL, permute_, permuteInverse_, denseColumn_, &numberRows_, &numberDense,
523                NULL, &i0, NULL, integerParameters_, doubleParameters_);
524          integerParameters_[29] = 0;
525          dense_->resetRowsDropped();
526          longDouble * denseBlob = dense_->aMatrix();
527          double * denseDiagonal = dense_->diagonal();
528          // Update dense matrix
529          for (i = 0; i < numberDense; i++) {
530               const double * a = denseColumn_ + i * numberRows_;
531               // do diagonal
532               double value = denseDiagonal[i];
533               const double * b = denseColumn_ + i * numberRows_;
534               for (int k = 0; k < numberRows_; k++)
535                    value += a[k] * b[k];
536               denseDiagonal[i] = value;
537               for (int j = i + 1; j < numberDense; j++) {
538                    double value = 0.0;
539                    const double * b = denseColumn_ + j * numberRows_;
540                    for (int k = 0; k < numberRows_; k++)
541                         value += a[k] * b[k];
542                    *denseBlob = value;
543                    denseBlob++;
544               }
545          }
546          // dense cholesky (? long double)
547          int * dropped = new int [numberDense];
548          dense_->factorizePart2(dropped);
549          delete [] dropped;
550     }
551     bool cleanCholesky;
552     if (model_->numberIterations() < 2000)
553          cleanCholesky = true;
554     else
555          cleanCholesky = false;
556     if (cleanCholesky) {
557          //drop fresh makes some formADAT easier
558          //int oldDropped=numberRowsDropped_;
559          if (newDropped || numberRowsDropped_) {
560               //std::cout <<"Rank "<<numberRows_-newDropped<<" ( "<<
561               //  newDropped<<" dropped)";
562               //if (newDropped>oldDropped)
563               //std::cout<<" ( "<<newDropped-oldDropped<<" dropped this time)";
564               //std::cout<<std::endl;
565               newDropped = 0;
566               for (int i = 0; i < numberRows_; i++) {
567                    char dropped = static_cast<char>(rowsDropped[i]);
568                    rowsDropped_[i] = dropped;
569                    rowsDropped_[i] = 0;
570                    if (dropped == 2) {
571                         //dropped this time
572                         rowsDropped[newDropped++] = i;
573                         rowsDropped_[i] = 0;
574                    }
575               }
576               numberRowsDropped_ = newDropped;
577               newDropped = -(2 + newDropped);
578          }
579     } else {
580          if (newDropped) {
581               newDropped = 0;
582               for (int i = 0; i < numberRows_; i++) {
583                    char dropped = static_cast<char>(rowsDropped[i]);
584                    rowsDropped_[i] = dropped;
585                    if (dropped == 2) {
586                         //dropped this time
587                         rowsDropped[newDropped++] = i;
588                         rowsDropped_[i] = 1;
589                    }
590               }
591          }
592          numberRowsDropped_ += newDropped;
593          if (numberRowsDropped_ && 0) {
594               std::cout << "Rank " << numberRows_ - numberRowsDropped_ << " ( " <<
595                         numberRowsDropped_ << " dropped)";
596               if (newDropped) {
597                    std::cout << " ( " << newDropped << " dropped this time)";
598               }
599               std::cout << std::endl;
600          }
601     }
602     status_ = 0;
603     return newDropped;
604}
605/* Uses factorization to solve. */
606void
607ClpCholeskyWssmp::solve (double * region)
608{
609     int i1 = 1;
610     int i0 = 0;
611     integerParameters_[1] = 4;
612     integerParameters_[2] = 4;
613#if 1
614     integerParameters_[5] = 3;
615     doubleParameters_[5] = 1.0e-10;
616     integerParameters_[6] = 6;
617#endif
618     if (!whichDense_) {
619          wssmp(&numberRows_, choleskyStart_, choleskyRow_, sparseFactor_,
620                NULL, permute_, permuteInverse_, region, &numberRows_, &i1,
621                NULL, &i0, NULL, integerParameters_, doubleParameters_);
622     } else {
623          // dense columns
624          integerParameters_[29] = 1;
625          wssmp(&numberRows_, choleskyStart_, choleskyRow_, sparseFactor_,
626                NULL, permute_, permuteInverse_, region, &numberRows_, &i1,
627                NULL, &i0, NULL, integerParameters_, doubleParameters_);
628          // do change;
629          int numberDense = dense_->numberRows();
630          double * change = new double[numberDense];
631          int i;
632          for (i = 0; i < numberDense; i++) {
633               const double * a = denseColumn_ + i * numberRows_;
634               double value = 0.0;
635               for (int iRow = 0; iRow < numberRows_; iRow++)
636                    value += a[iRow] * region[iRow];
637               change[i] = value;
638          }
639          // solve
640          dense_->solve(change);
641          for (i = 0; i < numberDense; i++) {
642               const double * a = denseColumn_ + i * numberRows_;
643               double value = change[i];
644               for (int iRow = 0; iRow < numberRows_; iRow++)
645                    region[iRow] -= value * a[iRow];
646          }
647          delete [] change;
648          // and finish off
649          integerParameters_[29] = 2;
650          integerParameters_[1] = 4;
651          wssmp(&numberRows_, choleskyStart_, choleskyRow_, sparseFactor_,
652                NULL, permute_, permuteInverse_, region, &numberRows_, &i1,
653                NULL, &i0, NULL, integerParameters_, doubleParameters_);
654          integerParameters_[29] = 0;
655     }
656#if 0
657     if (integerParameters_[5]) {
658          std::cout << integerParameters_[5] << " refinements ";
659     }
660     std::cout << doubleParameters_[6] << std::endl;
661#endif
662}
663#endif
Note: See TracBrowser for help on using the repository browser.