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

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/Clp/src/ClpCholeskyUfl.cpp
r1723 r2385 31 31 // Default Constructor 32 32 // 33 ClpCholeskyUfl::ClpCholeskyUfl 34 35 { 36 37 38 39 33 ClpCholeskyUfl::ClpCholeskyUfl(int denseThreshold) 34 : ClpCholeskyBase(denseThreshold) 35 { 36 type_ = 14; 37 L_ = NULL; 38 c_ = NULL; 39 40 40 #ifdef COIN_HAS_CHOLMOD 41 c_ = (cholmod_common*)malloc(sizeof(cholmod_common));42 cholmod_start (c_);43 44 41 c_ = (cholmod_common *)malloc(sizeof(cholmod_common)); 42 cholmod_start(c_); 43 // Can't use supernodal as may not be positive definite 44 c_>supernodal = 0; 45 45 #endif 46 46 } … … 49 49 // Copy constructor 50 50 // 51 ClpCholeskyUfl::ClpCholeskyUfl (const ClpCholeskyUfl & rhs) 52 : ClpCholeskyBase(rhs) 53 { 54 abort(); 55 } 56 51 ClpCholeskyUfl::ClpCholeskyUfl(const ClpCholeskyUfl &rhs) 52 : ClpCholeskyBase(rhs) 53 { 54 abort(); 55 } 57 56 58 57 // 59 58 // Destructor 60 59 // 61 ClpCholeskyUfl::~ClpCholeskyUfl 60 ClpCholeskyUfl::~ClpCholeskyUfl() 62 61 { 63 62 #ifdef COIN_HAS_CHOLMOD 64 cholmod_free_factor (&L_, c_);65 cholmod_finish (c_);66 63 cholmod_free_factor(&L_, c_); 64 cholmod_finish(c_); 65 free(c_); 67 66 #endif 68 67 } … … 72 71 // 73 72 ClpCholeskyUfl & 74 ClpCholeskyUfl::operator=(const ClpCholeskyUfl &rhs)75 { 76 77 78 79 80 73 ClpCholeskyUfl::operator=(const ClpCholeskyUfl &rhs) 74 { 75 if (this != &rhs) { 76 ClpCholeskyBase::operator=(rhs); 77 abort(); 78 } 79 return *this; 81 80 } 82 81 … … 84 83 // Clone 85 84 // 86 ClpCholeskyBase * 87 { 88 85 ClpCholeskyBase *ClpCholeskyUfl::clone() const 86 { 87 return new ClpCholeskyUfl(*this); 89 88 } 90 89 91 90 #ifndef COIN_HAS_CHOLMOD 92 91 /* Orders rows and saves pointer to matrix.and model */ 93 int 94 ClpCholeskyUfl::order(ClpInterior * model) 95 { 96 int iRow; 97 model_ = model; 98 if (preOrder(false, true, doKKT_)) 99 return 1; 100 permuteInverse_ = new int [numberRows_]; 101 permute_ = new int[numberRows_]; 102 double Control[AMD_CONTROL]; 103 double Info[AMD_INFO]; 104 105 amd_defaults(Control); 106 //amd_control(Control); 107 108 int returnCode = amd_order(numberRows_, choleskyStart_, choleskyRow_, 109 permute_, Control, Info); 110 delete [] choleskyRow_; 111 choleskyRow_ = NULL; 112 delete [] choleskyStart_; 113 choleskyStart_ = NULL; 114 //amd_info(Info); 115 116 if (returnCode != AMD_OK) { 117 std::cout << "AMD ordering failed" << std::endl; 118 return 1; 119 } 120 for (iRow = 0; iRow < numberRows_; iRow++) { 121 permuteInverse_[permute_[iRow]] = iRow; 122 } 123 return 0; 92 int ClpCholeskyUfl::order(ClpInterior *model) 93 { 94 int iRow; 95 model_ = model; 96 if (preOrder(false, true, doKKT_)) 97 return 1; 98 permuteInverse_ = new int[numberRows_]; 99 permute_ = new int[numberRows_]; 100 double Control[AMD_CONTROL]; 101 double Info[AMD_INFO]; 102 103 amd_defaults(Control); 104 //amd_control(Control); 105 106 int returnCode = amd_order(numberRows_, choleskyStart_, choleskyRow_, 107 permute_, Control, Info); 108 delete[] choleskyRow_; 109 choleskyRow_ = NULL; 110 delete[] choleskyStart_; 111 choleskyStart_ = NULL; 112 //amd_info(Info); 113 114 if (returnCode != AMD_OK) { 115 std::cout << "AMD ordering failed" << std::endl; 116 return 1; 117 } 118 for (iRow = 0; iRow < numberRows_; iRow++) { 119 permuteInverse_[permute_[iRow]] = iRow; 120 } 121 return 0; 124 122 } 125 123 #else 126 124 /* Orders rows and saves pointer to matrix.and model */ 127 int 128 ClpCholeskyUfl::order(ClpInterior * model) 129 { 130 numberRows_ = model>numberRows(); 131 if (doKKT_) { 132 numberRows_ += numberRows_ + model>numberColumns(); 133 printf("finish coding UFL KKT!\n"); 134 abort(); 135 } 136 rowsDropped_ = new char [numberRows_]; 137 memset(rowsDropped_, 0, numberRows_); 138 numberRowsDropped_ = 0; 139 model_ = model; 140 rowCopy_ = model>clpMatrix()>reverseOrderedCopy(); 141 // Space for starts 142 choleskyStart_ = new CoinBigIndex[numberRows_+1]; 143 const CoinBigIndex * columnStart = model_>clpMatrix()>getVectorStarts(); 144 const int * columnLength = model_>clpMatrix()>getVectorLengths(); 145 const int * row = model_>clpMatrix()>getIndices(); 146 const CoinBigIndex * rowStart = rowCopy_>getVectorStarts(); 147 const int * rowLength = rowCopy_>getVectorLengths(); 148 const int * column = rowCopy_>getIndices(); 149 // We need two arrays for counts 150 int * which = new int [numberRows_]; 151 int * used = new int[numberRows_+1]; 152 CoinZeroN(used, numberRows_); 153 int iRow; 154 sizeFactor_ = 0; 155 for (iRow = 0; iRow < numberRows_; iRow++) { 156 int number = 1; 157 // make sure diagonal exists 158 which[0] = iRow; 159 used[iRow] = 1; 160 if (!rowsDropped_[iRow]) { 161 CoinBigIndex startRow = rowStart[iRow]; 162 CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow]; 163 for (CoinBigIndex k = startRow; k < endRow; k++) { 164 int iColumn = column[k]; 165 CoinBigIndex start = columnStart[iColumn]; 166 CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn]; 167 for (CoinBigIndex j = start; j < end; j++) { 168 int jRow = row[j]; 169 if (jRow >= iRow && !rowsDropped_[jRow]) { 170 if (!used[jRow]) { 171 used[jRow] = 1; 172 which[number++] = jRow; 173 } 174 } 175 } 176 } 177 sizeFactor_ += number; 178 int j; 179 for (j = 0; j < number; j++) 180 used[which[j]] = 0; 125 int ClpCholeskyUfl::order(ClpInterior *model) 126 { 127 numberRows_ = model>numberRows(); 128 if (doKKT_) { 129 numberRows_ += numberRows_ + model>numberColumns(); 130 printf("finish coding UFL KKT!\n"); 131 abort(); 132 } 133 rowsDropped_ = new char[numberRows_]; 134 memset(rowsDropped_, 0, numberRows_); 135 numberRowsDropped_ = 0; 136 model_ = model; 137 rowCopy_ = model>clpMatrix()>reverseOrderedCopy(); 138 // Space for starts 139 choleskyStart_ = new CoinBigIndex[numberRows_ + 1]; 140 const CoinBigIndex *columnStart = model_>clpMatrix()>getVectorStarts(); 141 const int *columnLength = model_>clpMatrix()>getVectorLengths(); 142 const int *row = model_>clpMatrix()>getIndices(); 143 const CoinBigIndex *rowStart = rowCopy_>getVectorStarts(); 144 const int *rowLength = rowCopy_>getVectorLengths(); 145 const int *column = rowCopy_>getIndices(); 146 // We need two arrays for counts 147 int *which = new int[numberRows_]; 148 int *used = new int[numberRows_ + 1]; 149 CoinZeroN(used, numberRows_); 150 int iRow; 151 sizeFactor_ = 0; 152 for (iRow = 0; iRow < numberRows_; iRow++) { 153 int number = 1; 154 // make sure diagonal exists 155 which[0] = iRow; 156 used[iRow] = 1; 157 if (!rowsDropped_[iRow]) { 158 CoinBigIndex startRow = rowStart[iRow]; 159 CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow]; 160 for (CoinBigIndex k = startRow; k < endRow; k++) { 161 int iColumn = column[k]; 162 CoinBigIndex start = columnStart[iColumn]; 163 CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn]; 164 for (CoinBigIndex j = start; j < end; j++) { 165 int jRow = row[j]; 166 if (jRow >= iRow && !rowsDropped_[jRow]) { 167 if (!used[jRow]) { 168 used[jRow] = 1; 169 which[number++] = jRow; 170 } 181 171 } 182 } 183 delete [] which; 184 // Now we have size  create arrays and fill in 185 try { 186 choleskyRow_ = new int [sizeFactor_]; 187 } catch (...) { 188 // no memory 189 delete [] choleskyStart_; 190 choleskyStart_ = NULL; 191 return 1; 192 } 193 try { 194 sparseFactor_ = new double[sizeFactor_]; 195 } catch (...) { 196 // no memory 197 delete [] choleskyRow_; 198 choleskyRow_ = NULL; 199 delete [] choleskyStart_; 200 choleskyStart_ = NULL; 201 return 1; 202 } 203 204 sizeFactor_ = 0; 205 which = choleskyRow_; 206 for (iRow = 0; iRow < numberRows_; iRow++) { 207 int number = 1; 208 // make sure diagonal exists 209 which[0] = iRow; 210 used[iRow] = 1; 211 choleskyStart_[iRow] = sizeFactor_; 212 if (!rowsDropped_[iRow]) { 213 CoinBigIndex startRow = rowStart[iRow]; 214 CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow]; 215 for (CoinBigIndex k = startRow; k < endRow; k++) { 216 int iColumn = column[k]; 217 CoinBigIndex start = columnStart[iColumn]; 218 CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn]; 219 for (CoinBigIndex j = start; j < end; j++) { 220 int jRow = row[j]; 221 if (jRow >= iRow && !rowsDropped_[jRow]) { 222 if (!used[jRow]) { 223 used[jRow] = 1; 224 which[number++] = jRow; 225 } 226 } 227 } 228 } 229 sizeFactor_ += number; 230 int j; 231 for (j = 0; j < number; j++) 232 used[which[j]] = 0; 233 // Sort 234 std::sort(which, which + number); 235 // move which on 236 which += number; 172 } 173 } 174 sizeFactor_ += number; 175 int j; 176 for (j = 0; j < number; j++) 177 used[which[j]] = 0; 178 } 179 } 180 delete[] which; 181 // Now we have size  create arrays and fill in 182 try { 183 choleskyRow_ = new int[sizeFactor_]; 184 } catch (...) { 185 // no memory 186 delete[] choleskyStart_; 187 choleskyStart_ = NULL; 188 return 1; 189 } 190 try { 191 sparseFactor_ = new double[sizeFactor_]; 192 } catch (...) { 193 // no memory 194 delete[] choleskyRow_; 195 choleskyRow_ = NULL; 196 delete[] choleskyStart_; 197 choleskyStart_ = NULL; 198 return 1; 199 } 200 201 sizeFactor_ = 0; 202 which = choleskyRow_; 203 for (iRow = 0; iRow < numberRows_; iRow++) { 204 int number = 1; 205 // make sure diagonal exists 206 which[0] = iRow; 207 used[iRow] = 1; 208 choleskyStart_[iRow] = sizeFactor_; 209 if (!rowsDropped_[iRow]) { 210 CoinBigIndex startRow = rowStart[iRow]; 211 CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow]; 212 for (CoinBigIndex k = startRow; k < endRow; k++) { 213 int iColumn = column[k]; 214 CoinBigIndex start = columnStart[iColumn]; 215 CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn]; 216 for (CoinBigIndex j = start; j < end; j++) { 217 int jRow = row[j]; 218 if (jRow >= iRow && !rowsDropped_[jRow]) { 219 if (!used[jRow]) { 220 used[jRow] = 1; 221 which[number++] = jRow; 222 } 237 223 } 238 } 239 choleskyStart_[numberRows_] = sizeFactor_; 240 delete [] used; 241 permuteInverse_ = new int [numberRows_]; 242 permute_ = new int[numberRows_]; 243 cholmod_sparse A ; 244 A.nrow = numberRows_; 245 A.ncol = numberRows_; 246 A.nzmax = choleskyStart_[numberRows_]; 247 A.p = choleskyStart_; 248 A.i = choleskyRow_; 249 A.x = NULL; 250 A.stype = 1; 251 A.itype = CHOLMOD_INT; 252 A.xtype = CHOLMOD_PATTERN; 253 A.dtype = CHOLMOD_DOUBLE; 254 A.sorted = 1; 255 A.packed = 1; 256 c_>nmethods = 9; 257 c_>postorder = true; 258 //c_>dbound=1.0e20; 259 L_ = cholmod_analyze (&A, c_) ; 260 if (c_>status) { 261 COIN_DETAIL_PRINT(std::cout << "CHOLMOD ordering failed" << std::endl); 262 return 1; 263 } else { 264 COIN_DETAIL_PRINT(printf("%g nonzeros, flop count %g\n", c_>lnz, c_>fl)); 265 } 266 for (iRow = 0; iRow < numberRows_; iRow++) { 267 permuteInverse_[iRow] = iRow; 268 permute_[iRow] = iRow; 269 } 270 return 0; 224 } 225 } 226 sizeFactor_ += number; 227 int j; 228 for (j = 0; j < number; j++) 229 used[which[j]] = 0; 230 // Sort 231 std::sort(which, which + number); 232 // move which on 233 which += number; 234 } 235 } 236 choleskyStart_[numberRows_] = sizeFactor_; 237 delete[] used; 238 permuteInverse_ = new int[numberRows_]; 239 permute_ = new int[numberRows_]; 240 cholmod_sparse A; 241 A.nrow = numberRows_; 242 A.ncol = numberRows_; 243 A.nzmax = choleskyStart_[numberRows_]; 244 A.p = choleskyStart_; 245 A.i = choleskyRow_; 246 A.x = NULL; 247 A.stype = 1; 248 A.itype = CHOLMOD_INT; 249 A.xtype = CHOLMOD_PATTERN; 250 A.dtype = CHOLMOD_DOUBLE; 251 A.sorted = 1; 252 A.packed = 1; 253 c_>nmethods = 9; 254 c_>postorder = true; 255 //c_>dbound=1.0e20; 256 L_ = cholmod_analyze(&A, c_); 257 if (c_>status) { 258 COIN_DETAIL_PRINT(std::cout << "CHOLMOD ordering failed" << std::endl); 259 return 1; 260 } else { 261 COIN_DETAIL_PRINT(printf("%g nonzeros, flop count %g\n", c_>lnz, c_>fl)); 262 } 263 for (iRow = 0; iRow < numberRows_; iRow++) { 264 permuteInverse_[iRow] = iRow; 265 permute_[iRow] = iRow; 266 } 267 return 0; 271 268 } 272 269 #endif … … 276 273 user must provide factorize and solve. Otherwise the default factorization is used 277 274 returns nonzero if not enough memory */ 278 int 279 ClpCholeskyUfl::symbolic() 275 int ClpCholeskyUfl::symbolic() 280 276 { 281 277 #ifdef COIN_HAS_CHOLMOD 282 283 #else 284 278 return 0; 279 #else 280 return ClpCholeskyBase::symbolic(); 285 281 #endif 286 282 } … … 288 284 #ifdef COIN_HAS_CHOLMOD 289 285 /* Factorize  filling in rowsDropped and returning number dropped */ 290 int 291 ClpCholeskyUfl::factorize(const double * diagonal, int * rowsDropped) 292 { 293 const CoinBigIndex * columnStart = model_>clpMatrix()>getVectorStarts(); 294 const int * columnLength = model_>clpMatrix()>getVectorLengths(); 295 const int * row = model_>clpMatrix()>getIndices(); 296 const double * element = model_>clpMatrix()>getElements(); 297 const CoinBigIndex * rowStart = rowCopy_>getVectorStarts(); 298 const int * rowLength = rowCopy_>getVectorLengths(); 299 const int * column = rowCopy_>getIndices(); 300 const double * elementByRow = rowCopy_>getElements(); 301 int numberColumns = model_>clpMatrix()>getNumCols(); 302 int iRow; 303 double * work = new double[numberRows_]; 304 CoinZeroN(work, numberRows_); 305 const double * diagonalSlack = diagonal + numberColumns; 306 int newDropped = 0; 307 double largest; 308 //double smallest; 309 //perturbation 310 double perturbation = model_>diagonalPerturbation() * model_>diagonalNorm(); 311 perturbation = 0.0; 312 perturbation = perturbation * perturbation; 313 if (perturbation > 1.0) { 286 int ClpCholeskyUfl::factorize(const double *diagonal, int *rowsDropped) 287 { 288 const CoinBigIndex *columnStart = model_>clpMatrix()>getVectorStarts(); 289 const int *columnLength = model_>clpMatrix()>getVectorLengths(); 290 const int *row = model_>clpMatrix()>getIndices(); 291 const double *element = model_>clpMatrix()>getElements(); 292 const CoinBigIndex *rowStart = rowCopy_>getVectorStarts(); 293 const int *rowLength = rowCopy_>getVectorLengths(); 294 const int *column = rowCopy_>getIndices(); 295 const double *elementByRow = rowCopy_>getElements(); 296 int numberColumns = model_>clpMatrix()>getNumCols(); 297 int iRow; 298 double *work = new double[numberRows_]; 299 CoinZeroN(work, numberRows_); 300 const double *diagonalSlack = diagonal + numberColumns; 301 int newDropped = 0; 302 double largest; 303 //double smallest; 304 //perturbation 305 double perturbation = model_>diagonalPerturbation() * model_>diagonalNorm(); 306 perturbation = 0.0; 307 perturbation = perturbation * perturbation; 308 if (perturbation > 1.0) { 314 309 #ifdef COIN_DEVELOP 315 //if (model_>model()>logLevel()&4) 316 std::cout << "large perturbation " << perturbation << std::endl; 317 #endif 318 perturbation = sqrt(perturbation);; 319 perturbation = 1.0; 320 } 321 double delta2 = model_>delta(); // add delta*delta to diagonal 322 delta2 *= delta2; 323 for (iRow = 0; iRow < numberRows_; iRow++) { 324 double * put = sparseFactor_ + choleskyStart_[iRow]; 325 int * which = choleskyRow_ + choleskyStart_[iRow]; 326 int number = choleskyStart_[iRow+1]  choleskyStart_[iRow]; 327 if (!rowLength[iRow]) 328 rowsDropped_[iRow] = 1; 329 if (!rowsDropped_[iRow]) { 330 CoinBigIndex startRow = rowStart[iRow]; 331 CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow]; 332 work[iRow] = diagonalSlack[iRow] + delta2; 333 for (CoinBigIndex k = startRow; k < endRow; k++) { 334 int iColumn = column[k]; 335 if (!whichDense_  !whichDense_[iColumn]) { 336 CoinBigIndex start = columnStart[iColumn]; 337 CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn]; 338 double multiplier = diagonal[iColumn] * elementByRow[k]; 339 for (CoinBigIndex j = start; j < end; j++) { 340 int jRow = row[j]; 341 if (jRow >= iRow && !rowsDropped_[jRow]) { 342 double value = element[j] * multiplier; 343 work[jRow] += value; 344 } 345 } 346 } 347 } 348 int j; 349 for (j = 0; j < number; j++) { 350 int jRow = which[j]; 351 put[j] = work[jRow]; 352 work[jRow] = 0.0; 353 } 354 } else { 355 // dropped 356 int j; 357 for (j = 1; j < number; j++) { 358 put[j] = 0.0; 359 } 360 put[0] = 1.0; 310 //if (model_>model()>logLevel()&4) 311 std::cout << "large perturbation " << perturbation << std::endl; 312 #endif 313 perturbation = sqrt(perturbation); 314 ; 315 perturbation = 1.0; 316 } 317 double delta2 = model_>delta(); // add delta*delta to diagonal 318 delta2 *= delta2; 319 for (iRow = 0; iRow < numberRows_; iRow++) { 320 double *put = sparseFactor_ + choleskyStart_[iRow]; 321 int *which = choleskyRow_ + choleskyStart_[iRow]; 322 int number = choleskyStart_[iRow + 1]  choleskyStart_[iRow]; 323 if (!rowLength[iRow]) 324 rowsDropped_[iRow] = 1; 325 if (!rowsDropped_[iRow]) { 326 CoinBigIndex startRow = rowStart[iRow]; 327 CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow]; 328 work[iRow] = diagonalSlack[iRow] + delta2; 329 for (CoinBigIndex k = startRow; k < endRow; k++) { 330 int iColumn = column[k]; 331 if (!whichDense_  !whichDense_[iColumn]) { 332 CoinBigIndex start = columnStart[iColumn]; 333 CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn]; 334 double multiplier = diagonal[iColumn] * elementByRow[k]; 335 for (CoinBigIndex j = start; j < end; j++) { 336 int jRow = row[j]; 337 if (jRow >= iRow && !rowsDropped_[jRow]) { 338 double value = element[j] * multiplier; 339 work[jRow] += value; 340 } 361 341 } 362 } 363 //check sizes 364 double largest2 = maximumAbsElement(sparseFactor_, sizeFactor_); 365 largest2 *= 1.0e20; 366 largest = CoinMin(largest2, 1.0e11); 367 int numberDroppedBefore = 0; 368 for (iRow = 0; iRow < numberRows_; iRow++) { 369 int dropped = rowsDropped_[iRow]; 370 // Move to int array 371 rowsDropped[iRow] = dropped; 372 if (!dropped) { 373 CoinBigIndex start = choleskyStart_[iRow]; 374 double diagonal = sparseFactor_[start]; 375 if (diagonal > largest2) { 376 sparseFactor_[start] = CoinMax(diagonal, 1.0e10); 377 } else { 378 sparseFactor_[start] = CoinMax(diagonal, 1.0e10); 379 rowsDropped[iRow] = 2; 380 numberDroppedBefore++; 381 } 382 } 383 } 384 delete [] work; 385 cholmod_sparse A ; 386 A.nrow = numberRows_; 387 A.ncol = numberRows_; 388 A.nzmax = choleskyStart_[numberRows_]; 389 A.p = choleskyStart_; 390 A.i = choleskyRow_; 391 A.x = sparseFactor_; 392 A.stype = 1; 393 A.itype = CHOLMOD_INT; 394 A.xtype = CHOLMOD_REAL; 395 A.dtype = CHOLMOD_DOUBLE; 396 A.sorted = 1; 397 A.packed = 1; 398 cholmod_factorize (&A, L_, c_) ; /* factorize */ 399 choleskyCondition_ = 1.0; 400 bool cleanCholesky; 401 if (model_>numberIterations() < 2000) 402 cleanCholesky = true; 403 else 404 cleanCholesky = false; 405 if (cleanCholesky) { 406 //drop fresh makes some formADAT easier 407 //int oldDropped=numberRowsDropped_; 408 if (newDropped  numberRowsDropped_) { 409 //std::cout <<"Rank "<<numberRows_newDropped<<" ( "<< 410 // newDropped<<" dropped)"; 411 //if (newDropped>oldDropped) 412 //std::cout<<" ( "<<newDroppedoldDropped<<" dropped this time)"; 413 //std::cout<<std::endl; 414 newDropped = 0; 415 for (int i = 0; i < numberRows_; i++) { 416 int dropped = rowsDropped[i]; 417 rowsDropped_[i] = (char)dropped; 418 if (dropped == 2) { 419 //dropped this time 420 rowsDropped[newDropped++] = i; 421 rowsDropped_[i] = 0; 422 } 423 } 424 numberRowsDropped_ = newDropped; 425 newDropped = (2 + newDropped); 426 } 427 } else { 428 if (newDropped) { 429 newDropped = 0; 430 for (int i = 0; i < numberRows_; i++) { 431 int dropped = rowsDropped[i]; 432 rowsDropped_[i] = (char)dropped; 433 if (dropped == 2) { 434 //dropped this time 435 rowsDropped[newDropped++] = i; 436 rowsDropped_[i] = 1; 437 } 438 } 439 } 440 numberRowsDropped_ += newDropped; 441 if (numberRowsDropped_ && 0) { 442 std::cout << "Rank " << numberRows_  numberRowsDropped_ << " ( " << 443 numberRowsDropped_ << " dropped)"; 444 if (newDropped) { 445 std::cout << " ( " << newDropped << " dropped this time)"; 446 } 447 std::cout << std::endl; 448 } 449 } 450 status_ = 0; 451 return newDropped; 342 } 343 } 344 int j; 345 for (j = 0; j < number; j++) { 346 int jRow = which[j]; 347 put[j] = work[jRow]; 348 work[jRow] = 0.0; 349 } 350 } else { 351 // dropped 352 int j; 353 for (j = 1; j < number; j++) { 354 put[j] = 0.0; 355 } 356 put[0] = 1.0; 357 } 358 } 359 //check sizes 360 double largest2 = maximumAbsElement(sparseFactor_, sizeFactor_); 361 largest2 *= 1.0e20; 362 largest = CoinMin(largest2, 1.0e11); 363 int numberDroppedBefore = 0; 364 for (iRow = 0; iRow < numberRows_; iRow++) { 365 int dropped = rowsDropped_[iRow]; 366 // Move to int array 367 rowsDropped[iRow] = dropped; 368 if (!dropped) { 369 CoinBigIndex start = choleskyStart_[iRow]; 370 double diagonal = sparseFactor_[start]; 371 if (diagonal > largest2) { 372 sparseFactor_[start] = CoinMax(diagonal, 1.0e10); 373 } else { 374 sparseFactor_[start] = CoinMax(diagonal, 1.0e10); 375 rowsDropped[iRow] = 2; 376 numberDroppedBefore++; 377 } 378 } 379 } 380 delete[] work; 381 cholmod_sparse A; 382 A.nrow = numberRows_; 383 A.ncol = numberRows_; 384 A.nzmax = choleskyStart_[numberRows_]; 385 A.p = choleskyStart_; 386 A.i = choleskyRow_; 387 A.x = sparseFactor_; 388 A.stype = 1; 389 A.itype = CHOLMOD_INT; 390 A.xtype = CHOLMOD_REAL; 391 A.dtype = CHOLMOD_DOUBLE; 392 A.sorted = 1; 393 A.packed = 1; 394 cholmod_factorize(&A, L_, c_); /* factorize */ 395 choleskyCondition_ = 1.0; 396 bool cleanCholesky; 397 if (model_>numberIterations() < 2000) 398 cleanCholesky = true; 399 else 400 cleanCholesky = false; 401 if (cleanCholesky) { 402 //drop fresh makes some formADAT easier 403 //int oldDropped=numberRowsDropped_; 404 if (newDropped  numberRowsDropped_) { 405 //std::cout <<"Rank "<<numberRows_newDropped<<" ( "<< 406 // newDropped<<" dropped)"; 407 //if (newDropped>oldDropped) 408 //std::cout<<" ( "<<newDroppedoldDropped<<" dropped this time)"; 409 //std::cout<<std::endl; 410 newDropped = 0; 411 for (int i = 0; i < numberRows_; i++) { 412 int dropped = rowsDropped[i]; 413 rowsDropped_[i] = (char)dropped; 414 if (dropped == 2) { 415 //dropped this time 416 rowsDropped[newDropped++] = i; 417 rowsDropped_[i] = 0; 418 } 419 } 420 numberRowsDropped_ = newDropped; 421 newDropped = (2 + newDropped); 422 } 423 } else { 424 if (newDropped) { 425 newDropped = 0; 426 for (int i = 0; i < numberRows_; i++) { 427 int dropped = rowsDropped[i]; 428 rowsDropped_[i] = (char)dropped; 429 if (dropped == 2) { 430 //dropped this time 431 rowsDropped[newDropped++] = i; 432 rowsDropped_[i] = 1; 433 } 434 } 435 } 436 numberRowsDropped_ += newDropped; 437 if (numberRowsDropped_ && 0) { 438 std::cout << "Rank " << numberRows_  numberRowsDropped_ << " ( " << numberRowsDropped_ << " dropped)"; 439 if (newDropped) { 440 std::cout << " ( " << newDropped << " dropped this time)"; 441 } 442 std::cout << std::endl; 443 } 444 } 445 status_ = 0; 446 return newDropped; 452 447 } 453 448 #else 454 449 /* Factorize  filling in rowsDropped and returning number dropped */ 455 int 456 ClpCholeskyUfl::factorize(const double * diagonal, int * rowsDropped) 450 int ClpCholeskyUfl::factorize(const double *diagonal, int *rowsDropped) 457 451 { 458 452 return ClpCholeskyBase::factorize(diagonal, rowsDropped); … … 462 456 #ifdef COIN_HAS_CHOLMOD 463 457 /* Uses factorization to solve. */ 464 void 465 ClpCholeskyUfl::solve (double * region) 466 { 467 cholmod_dense *x, *b; 468 b = cholmod_allocate_dense (numberRows_, 1, numberRows_, CHOLMOD_REAL, c_) ; 469 CoinMemcpyN(region, numberRows_, (double *) b>x); 470 x = cholmod_solve (CHOLMOD_A, L_, b, c_) ; 471 CoinMemcpyN((double *) x>x, numberRows_, region); 472 cholmod_free_dense (&x, c_) ; 473 cholmod_free_dense (&b, c_) ; 474 } 475 #else 476 void 477 ClpCholeskyUfl::solve (double * region) 458 void ClpCholeskyUfl::solve(double *region) 459 { 460 cholmod_dense *x, *b; 461 b = cholmod_allocate_dense(numberRows_, 1, numberRows_, CHOLMOD_REAL, c_); 462 CoinMemcpyN(region, numberRows_, (double *)b>x); 463 x = cholmod_solve(CHOLMOD_A, L_, b, c_); 464 CoinMemcpyN((double *)x>x, numberRows_, region); 465 cholmod_free_dense(&x, c_); 466 cholmod_free_dense(&b, c_); 467 } 468 #else 469 void ClpCholeskyUfl::solve(double *region) 478 470 { 479 471 ClpCholeskyBase::solve(region); 480 472 } 481 473 #endif 474 475 /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2 476 */
Note: See TracChangeset
for help on using the changeset viewer.