Ignore:
Timestamp:
Jan 6, 2019 2:43:06 PM (3 months ago)
Author:
unxusr
Message:

formatting

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Clp/src/ClpCholeskyWssmp.cpp

    r1723 r2385  
    44// This code is licensed under the terms of the Eclipse Public License (EPL).
    55
    6 
    76#include "CoinPragma.hpp"
    87#include "CoinHelperFunctions.hpp"
     
    2120// Default Constructor
    2221//-------------------------------------------------------------------
    23 ClpCholeskyWssmp::ClpCholeskyWssmp (int denseThreshold)
    24      : ClpCholeskyBase(denseThreshold)
    25 {
    26      type_ = 12;
     22ClpCholeskyWssmp::ClpCholeskyWssmp(int denseThreshold)
     23  : ClpCholeskyBase(denseThreshold)
     24{
     25  type_ = 12;
    2726}
    2827
     
    3029// Copy constructor
    3130//-------------------------------------------------------------------
    32 ClpCholeskyWssmp::ClpCholeskyWssmp (const ClpCholeskyWssmp & rhs)
    33      : ClpCholeskyBase(rhs)
    34 {
    35 }
    36 
     31ClpCholeskyWssmp::ClpCholeskyWssmp(const ClpCholeskyWssmp &rhs)
     32  : ClpCholeskyBase(rhs)
     33{
     34}
    3735
    3836//-------------------------------------------------------------------
    3937// Destructor
    4038//-------------------------------------------------------------------
    41 ClpCholeskyWssmp::~ClpCholeskyWssmp ()
     39ClpCholeskyWssmp::~ClpCholeskyWssmp()
    4240{
    4341}
     
    4745//-------------------------------------------------------------------
    4846ClpCholeskyWssmp &
    49 ClpCholeskyWssmp::operator=(const ClpCholeskyWssmp& rhs)
    50 {
    51      if (this != &rhs) {
    52           ClpCholeskyBase::operator=(rhs);
    53      }
    54      return *this;
     47ClpCholeskyWssmp::operator=(const ClpCholeskyWssmp &rhs)
     48{
     49  if (this != &rhs) {
     50    ClpCholeskyBase::operator=(rhs);
     51  }
     52  return *this;
    5553}
    5654//-------------------------------------------------------------------
    5755// Clone
    5856//-------------------------------------------------------------------
    59 ClpCholeskyBase * ClpCholeskyWssmp::clone() const
    60 {
    61      return new ClpCholeskyWssmp(*this);
     57ClpCholeskyBase *ClpCholeskyWssmp::clone() const
     58{
     59  return new ClpCholeskyWssmp(*this);
    6260}
    6361// At present I can't get wssmp to work as my libraries seem to be out of sync
     
    6563#ifndef USE_EKKWSSMP
    6664extern "C" {
    67 void F77_FUNC(wsetmaxthrds,WSETMAXTHRDS)(const int* NTHREADS);
    68 
    69 void 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);
    77 void F77_FUNC_(wsmp_clear,WSMP_CLEAR)(void);
     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);
    7876}
    7977#else
     
    8280typedef struct EKKContext EKKContext;
    8381
    84 
    8582extern "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 }
    91 static  EKKModel * model = NULL;
    92 static  EKKContext * context = NULL;
    93 extern "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);
    99 static 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);
     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);
    115111}
    116112#endif
    117113/* Orders rows and saves pointer to matrix.and model */
    118 int
    119 ClpCholeskyWssmp::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;
     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            }
    153203          }
    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;
     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            }
    162257          }
    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;
     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;
    282277#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);
     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);
    289285#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;
     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;
    334333}
    335334/* Does Symbolic factorization given permutation.
     
    337336   user must provide factorize and solve.  Otherwise the default factorization is used
    338337   returns non-zero if not enough memory */
    339 int
    340 ClpCholeskyWssmp::symbolic()
    341 {
    342      return 0;
     338int ClpCholeskyWssmp::symbolic()
     339{
     340  return 0;
    343341}
    344342/* Factorize - filling in rowsDropped and returning number dropped */
    345 int
    346 ClpCholeskyWssmp::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) {
     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) {
    371368#ifdef COIN_DEVELOP
    372           //if (model_->model()->logLevel()&4)
    373           std::cout << "large perturbation " << perturbation << std::endl;
     369    //if (model_->model()->logLevel()&4)
     370    std::cout << "large perturbation " << perturbation << std::endl;
    374371#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                }
     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            }
    394418          }
    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;
     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;
    587586}
    588587/* Uses factorization to solve. */
    589 void
    590 ClpCholeskyWssmp::solve (double * region)
    591 {
    592      int i1 = 1;
    593      int i0 = 0;
    594      integerParameters_[1] = 4;
    595      integerParameters_[2] = 4;
     588void ClpCholeskyWssmp::solve(double *region)
     589{
     590  int i1 = 1;
     591  int i0 = 0;
     592  integerParameters_[1] = 4;
     593  integerParameters_[2] = 4;
    596594#if 1
    597      integerParameters_[5] = 3;
    598      doubleParameters_[5] = 1.0e-10;
    599      integerParameters_[6] = 6;
     595  integerParameters_[5] = 3;
     596  doubleParameters_[5] = 1.0e-10;
     597  integerParameters_[6] = 6;
    600598#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      }
     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  }
    639640#if 0
    640641     if (integerParameters_[5]) {
     
    644645#endif
    645646}
     647
     648/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
     649*/
Note: See TracChangeset for help on using the changeset viewer.