source: stable/1.15/Clp/src/ClpCholeskyWssmp.cpp @ 2018

Last change on this file since 2018 was 1723, checked in by forrest, 9 years ago

out some printf statements

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