source: trunk/Clp/src/ClpCholeskyWssmp.cpp

Last change on this file was 2385, checked in by unxusr, 5 months ago

formatting

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