Changeset 2385 for trunk/Clp/src/ClpCholeskyWssmp.cpp
 Timestamp:
 Jan 6, 2019 2:43:06 PM (3 months ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/Clp/src/ClpCholeskyWssmp.cpp
r1723 r2385 4 4 // This code is licensed under the terms of the Eclipse Public License (EPL). 5 5 6 7 6 #include "CoinPragma.hpp" 8 7 #include "CoinHelperFunctions.hpp" … … 21 20 // Default Constructor 22 21 // 23 ClpCholeskyWssmp::ClpCholeskyWssmp 24 25 { 26 22 ClpCholeskyWssmp::ClpCholeskyWssmp(int denseThreshold) 23 : ClpCholeskyBase(denseThreshold) 24 { 25 type_ = 12; 27 26 } 28 27 … … 30 29 // Copy constructor 31 30 // 32 ClpCholeskyWssmp::ClpCholeskyWssmp (const ClpCholeskyWssmp & rhs) 33 : ClpCholeskyBase(rhs) 34 { 35 } 36 31 ClpCholeskyWssmp::ClpCholeskyWssmp(const ClpCholeskyWssmp &rhs) 32 : ClpCholeskyBase(rhs) 33 { 34 } 37 35 38 36 // 39 37 // Destructor 40 38 // 41 ClpCholeskyWssmp::~ClpCholeskyWssmp 39 ClpCholeskyWssmp::~ClpCholeskyWssmp() 42 40 { 43 41 } … … 47 45 // 48 46 ClpCholeskyWssmp & 49 ClpCholeskyWssmp::operator=(const ClpCholeskyWssmp &rhs)50 { 51 52 53 54 47 ClpCholeskyWssmp::operator=(const ClpCholeskyWssmp &rhs) 48 { 49 if (this != &rhs) { 50 ClpCholeskyBase::operator=(rhs); 51 } 52 return *this; 55 53 } 56 54 // 57 55 // Clone 58 56 // 59 ClpCholeskyBase * 60 { 61 57 ClpCholeskyBase *ClpCholeskyWssmp::clone() const 58 { 59 return new ClpCholeskyWssmp(*this); 62 60 } 63 61 // At present I can't get wssmp to work as my libraries seem to be out of sync … … 65 63 #ifndef USE_EKKWSSMP 66 64 extern "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);65 void F77_FUNC(wsetmaxthrds, WSETMAXTHRDS)(const int *NTHREADS); 66 67 void 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); 75 void F77_FUNC_(wsmp_clear, WSMP_CLEAR)(void); 78 76 } 79 77 #else … … 82 80 typedef struct EKKContext EKKContext; 83 81 84 85 82 extern "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); 83 EKKContext *ekk_initializeContext(); 84 void ekk_endContext(EKKContext *context); 85 EKKModel *ekk_newModel(EKKContext *env, const char *name); 86 int ekk_deleteModel(EKKModel *model); 87 } 88 static EKKModel *model = NULL; 89 static EKKContext *context = NULL; 90 extern "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); 96 static 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); 115 111 } 116 112 #endif 117 113 /* 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; 114 int 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 } 153 203 } 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 } 162 257 } 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; 282 277 #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); 289 285 #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.0e20; 304 doubleParameters_[11] = 1.0e15; 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.0e5); 316 model>setDelta(5.0e5); 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.0e20; 301 doubleParameters_[11] = 1.0e15; 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.0e5); 314 model>setDelta(5.0e5); 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; 334 333 } 335 334 /* Does Symbolic factorization given permutation. … … 337 336 user must provide factorize and solve. Otherwise the default factorization is used 338 337 returns nonzero if not enough memory */ 339 int 340 ClpCholeskyWssmp::symbolic() 341 { 342 return 0; 338 int ClpCholeskyWssmp::symbolic() 339 { 340 return 0; 343 341 } 344 342 /* 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) { 343 int 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) { 371 368 #ifdef COIN_DEVELOP 372 373 369 //if (model_>model()>logLevel()&4) 370 std::cout << "large perturbation " << perturbation << std::endl; 374 371 #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 } 394 418 } 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.0e20; 441 largest = CoinMin(largest2, 1.0e11); 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.0e20, largest); 467 if (doubleParameters_[10] > 1.0e3) 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<<" ( "<<newDroppedoldDropped<<" 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.0e20; 439 largest = CoinMin(largest2, 1.0e11); 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.0e20, largest); 465 if (doubleParameters_[10] > 1.0e3) 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<<" ( "<<newDroppedoldDropped<<" 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; 587 586 } 588 587 /* 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; 588 void ClpCholeskyWssmp::solve(double *region) 589 { 590 int i1 = 1; 591 int i0 = 0; 592 integerParameters_[1] = 4; 593 integerParameters_[2] = 4; 596 594 #if 1 597 598 599 595 integerParameters_[5] = 3; 596 doubleParameters_[5] = 1.0e10; 597 integerParameters_[6] = 6; 600 598 #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 } 639 640 #if 0 640 641 if (integerParameters_[5]) { … … 644 645 #endif 645 646 } 647 648 /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2 649 */
Note: See TracChangeset
for help on using the changeset viewer.