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

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/Clp/src/ClpCholeskyMumps.cpp
r2275 r2385 3 3 // Corporation and others. All Rights Reserved. 4 4 // This code is licensed under the terms of the Eclipse Public License (EPL). 5 6 5 7 6 #include "CoinPragma.hpp" … … 36 35 // Default Constructor 37 36 // 38 ClpCholeskyMumps::ClpCholeskyMumps (int denseThreshold,int logLevel)39 40 { 41 mumps_ = (DMUMPS_STRUC_C*)malloc(sizeof(DMUMPS_STRUC_C));42 43 44 45 46 47 48 mumps_>job = JOB_INIT;//initialize mumps49 mumps_>par = 1;//working host for sequential version50 mumps_>sym = 2;//general symmetric matrix51 52 53 54 37 ClpCholeskyMumps::ClpCholeskyMumps(int denseThreshold, int logLevel) 38 : ClpCholeskyBase(denseThreshold) 39 { 40 mumps_ = (DMUMPS_STRUC_C *)malloc(sizeof(DMUMPS_STRUC_C)); 41 type_ = 16; 42 mumps_>n = 0; 43 mumps_>nz = 0; 44 mumps_>a = NULL; 45 mumps_>jcn = NULL; 46 mumps_>irn = NULL; 47 mumps_>job = JOB_INIT; //initialize mumps 48 mumps_>par = 1; //working host for sequential version 49 mumps_>sym = 2; //general symmetric matrix 50 mumps_>comm_fortran = USE_COMM_WORLD; 51 int myid; 52 int justName; 53 MPI_Init(&justName, NULL); 55 54 #ifndef NDEBUG 56 57 assert(!ierr);55 int ierr = MPI_Comm_rank(MPI_COMM_WORLD, &myid); 56 assert(!ierr); 58 57 #else 59 58 MPI_Comm_rank(MPI_COMM_WORLD, &myid); 60 59 #endif 61 60 dmumps_c(mumps_); 62 61 #define ICNTL(I) icntl[(I)1] /* macro s.t. indices match documentation */ 63 62 #define CNTL(I) cntl[(I)1] /* macro s.t. indices match documentation */ 64 65 66 67 68 69 70 71 72 73 74 63 mumps_>ICNTL(5) = 1; // say compressed format 64 mumps_>ICNTL(4) = 2; // log messages 65 mumps_>ICNTL(24) = 1; // Deal with zeros on diagonal 66 mumps_>CNTL(3) = 1.0e20; // drop if diagonal less than this 67 if (!logLevel) { 68 // output off 69 mumps_>ICNTL(1) = 1; 70 mumps_>ICNTL(2) = 1; 71 mumps_>ICNTL(3) = 1; 72 mumps_>ICNTL(4) = 0; 73 } 75 74 } 76 75 … … 78 77 // Copy constructor 79 78 // 80 ClpCholeskyMumps::ClpCholeskyMumps (const ClpCholeskyMumps & rhs) 81 : ClpCholeskyBase(rhs) 82 { 83 abort(); 84 } 85 79 ClpCholeskyMumps::ClpCholeskyMumps(const ClpCholeskyMumps &rhs) 80 : ClpCholeskyBase(rhs) 81 { 82 abort(); 83 } 86 84 87 85 // 88 86 // Destructor 89 87 // 90 ClpCholeskyMumps::~ClpCholeskyMumps 91 { 92 93 94 95 88 ClpCholeskyMumps::~ClpCholeskyMumps() 89 { 90 mumps_>job = JOB_END; 91 dmumps_c(mumps_); /* Terminate instance */ 92 MPI_Finalize(); 93 free(mumps_); 96 94 } 97 95 … … 100 98 // 101 99 ClpCholeskyMumps & 102 ClpCholeskyMumps::operator=(const ClpCholeskyMumps &rhs)103 { 104 105 106 107 108 100 ClpCholeskyMumps::operator=(const ClpCholeskyMumps &rhs) 101 { 102 if (this != &rhs) { 103 ClpCholeskyBase::operator=(rhs); 104 abort(); 105 } 106 return *this; 109 107 } 110 108 // 111 109 // Clone 112 110 // 113 ClpCholeskyBase * 114 { 115 111 ClpCholeskyBase *ClpCholeskyMumps::clone() const 112 { 113 return new ClpCholeskyMumps(*this); 116 114 } 117 115 /* Orders rows and saves pointer to matrix.and model */ 118 int 119 ClpCholeskyMumps::order(ClpInterior * model) 120 { 121 numberRows_ = model>numberRows(); 122 if (doKKT_) { 123 numberRows_ += numberRows_ + model>numberColumns(); 124 printf("finish coding MUMPS KKT!\n"); 125 abort(); 126 } 127 rowsDropped_ = new char [numberRows_]; 128 memset(rowsDropped_, 0, numberRows_); 129 numberRowsDropped_ = 0; 130 model_ = model; 131 rowCopy_ = model>clpMatrix()>reverseOrderedCopy(); 132 const CoinBigIndex * columnStart = model_>clpMatrix()>getVectorStarts(); 133 const int * columnLength = model_>clpMatrix()>getVectorLengths(); 134 const int * row = model_>clpMatrix()>getIndices(); 135 const CoinBigIndex * rowStart = rowCopy_>getVectorStarts(); 136 const int * rowLength = rowCopy_>getVectorLengths(); 137 const int * column = rowCopy_>getIndices(); 138 // We need two arrays for counts 139 int * which = new int [numberRows_]; 140 int * used = new int[numberRows_+1]; 141 CoinZeroN(used, numberRows_); 142 int iRow; 143 sizeFactor_ = 0; 144 for (iRow = 0; iRow < numberRows_; iRow++) { 145 int number = 1; 146 // make sure diagonal exists 147 which[0] = iRow; 148 used[iRow] = 1; 149 if (!rowsDropped_[iRow]) { 150 CoinBigIndex startRow = rowStart[iRow]; 151 CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow]; 152 for (CoinBigIndex k = startRow; k < endRow; k++) { 153 int iColumn = column[k]; 154 CoinBigIndex start = columnStart[iColumn]; 155 CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn]; 156 for (CoinBigIndex j = start; j < end; j++) { 157 int jRow = row[j]; 158 if (jRow >= iRow && !rowsDropped_[jRow]) { 159 if (!used[jRow]) { 160 used[jRow] = 1; 161 which[number++] = jRow; 162 } 163 } 164 } 165 } 166 sizeFactor_ += number; 167 int j; 168 for (j = 0; j < number; j++) 169 used[which[j]] = 0; 116 int ClpCholeskyMumps::order(ClpInterior *model) 117 { 118 numberRows_ = model>numberRows(); 119 if (doKKT_) { 120 numberRows_ += numberRows_ + model>numberColumns(); 121 printf("finish coding MUMPS KKT!\n"); 122 abort(); 123 } 124 rowsDropped_ = new char[numberRows_]; 125 memset(rowsDropped_, 0, numberRows_); 126 numberRowsDropped_ = 0; 127 model_ = model; 128 rowCopy_ = model>clpMatrix()>reverseOrderedCopy(); 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 for (iRow = 0; iRow < numberRows_; iRow++) { 142 int number = 1; 143 // make sure diagonal exists 144 which[0] = iRow; 145 used[iRow] = 1; 146 if (!rowsDropped_[iRow]) { 147 CoinBigIndex startRow = rowStart[iRow]; 148 CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow]; 149 for (CoinBigIndex k = startRow; k < endRow; k++) { 150 int iColumn = column[k]; 151 CoinBigIndex start = columnStart[iColumn]; 152 CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn]; 153 for (CoinBigIndex j = start; j < end; j++) { 154 int jRow = row[j]; 155 if (jRow >= iRow && !rowsDropped_[jRow]) { 156 if (!used[jRow]) { 157 used[jRow] = 1; 158 which[number++] = jRow; 159 } 170 160 } 171 } 172 delete [] which; 173 // NOT COMPRESSED FOR NOW ???  Space for starts 174 mumps_>ICNTL(5) = 0; // say NOT compressed format 175 try { 176 choleskyStart_ = new int[numberRows_+1+sizeFactor_]; 177 } catch (...) { 178 // no memory 179 return 1; 180 } 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 } 223 } 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; 161 } 162 } 163 sizeFactor_ += number; 164 int j; 165 for (j = 0; j < number; j++) 166 used[which[j]] = 0; 167 } 168 } 169 delete[] which; 170 // NOT COMPRESSED FOR NOW ???  Space for starts 171 mumps_>ICNTL(5) = 0; // say NOT compressed format 172 try { 173 choleskyStart_ = new int[numberRows_ + 1 + sizeFactor_]; 174 } catch (...) { 175 // no memory 176 return 1; 177 } 178 // Now we have size  create arrays and fill in 179 try { 180 choleskyRow_ = new int[sizeFactor_]; 181 } catch (...) { 182 // no memory 183 delete[] choleskyStart_; 184 choleskyStart_ = NULL; 185 return 1; 186 } 187 try { 188 sparseFactor_ = new double[sizeFactor_]; 189 } catch (...) { 190 // no memory 191 delete[] choleskyRow_; 192 choleskyRow_ = NULL; 193 delete[] choleskyStart_; 194 choleskyStart_ = NULL; 195 return 1; 196 } 197 198 sizeFactor_ = 0; 199 which = choleskyRow_; 200 for (iRow = 0; iRow < numberRows_; iRow++) { 201 int number = 1; 202 // make sure diagonal exists 203 which[0] = iRow; 204 used[iRow] = 1; 205 choleskyStart_[iRow] = sizeFactor_; 206 if (!rowsDropped_[iRow]) { 207 CoinBigIndex startRow = rowStart[iRow]; 208 CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow]; 209 for (CoinBigIndex k = startRow; k < endRow; k++) { 210 int iColumn = column[k]; 211 CoinBigIndex start = columnStart[iColumn]; 212 CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn]; 213 for (CoinBigIndex j = start; j < end; j++) { 214 int jRow = row[j]; 215 if (jRow >= iRow && !rowsDropped_[jRow]) { 216 if (!used[jRow]) { 217 used[jRow] = 1; 218 which[number++] = jRow; 219 } 234 220 } 235 } 236 choleskyStart_[numberRows_] = sizeFactor_; 237 delete [] used; 238 permuteInverse_ = new int [numberRows_]; 239 permute_ = new int[numberRows_]; 240 // To Fortran and fake 241 for (iRow = 0; iRow < numberRows_ + 1; iRow++) { 242 int k = choleskyStart_[iRow]; 243 int kEnd = choleskyStart_[iRow+1]; 244 k += numberRows_ + 1; 245 kEnd += numberRows_ + 1; 246 for (; k < kEnd; k++) 247 choleskyStart_[k] = iRow + 1; 248 choleskyStart_[iRow]++; 249 } 250 mumps_>nz = sizeFactor_; 251 mumps_>irn = choleskyStart_ + numberRows_ + 1; 252 mumps_>jcn = choleskyRow_; 253 mumps_>a = NULL; 254 for (CoinBigIndex i = 0; i < sizeFactor_; i++) { 255 choleskyRow_[i]++; 221 } 222 } 223 sizeFactor_ += number; 224 int j; 225 for (j = 0; j < number; j++) 226 used[which[j]] = 0; 227 // Sort 228 std::sort(which, which + number); 229 // move which on 230 which += number; 231 } 232 } 233 choleskyStart_[numberRows_] = sizeFactor_; 234 delete[] used; 235 permuteInverse_ = new int[numberRows_]; 236 permute_ = new int[numberRows_]; 237 // To Fortran and fake 238 for (iRow = 0; iRow < numberRows_ + 1; iRow++) { 239 int k = choleskyStart_[iRow]; 240 int kEnd = choleskyStart_[iRow + 1]; 241 k += numberRows_ + 1; 242 kEnd += numberRows_ + 1; 243 for (; k < kEnd; k++) 244 choleskyStart_[k] = iRow + 1; 245 choleskyStart_[iRow]++; 246 } 247 mumps_>nz = sizeFactor_; 248 mumps_>irn = choleskyStart_ + numberRows_ + 1; 249 mumps_>jcn = choleskyRow_; 250 mumps_>a = NULL; 251 for (CoinBigIndex i = 0; i < sizeFactor_; i++) { 252 choleskyRow_[i]++; 256 253 #ifndef NDEBUG 257 assert(mumps_>irn[i] >= 1 && mumps_>irn[i] <= numberRows_);258 assert(mumps_>jcn[i] >= 1 && mumps_>jcn[i] <= numberRows_);254 assert(mumps_>irn[i] >= 1 && mumps_>irn[i] <= numberRows_); 255 assert(mumps_>jcn[i] >= 1 && mumps_>jcn[i] <= numberRows_); 259 256 #endif 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 257 } 258 // validate 259 //mumps code here 260 mumps_>n = numberRows_; 261 mumps_>nelt = numberRows_; 262 mumps_>eltptr = choleskyStart_; 263 mumps_>eltvar = choleskyRow_; 264 mumps_>a_elt = NULL; 265 mumps_>rhs = NULL; 266 mumps_>job = 1; // order 267 dmumps_c(mumps_); 268 mumps_>a = sparseFactor_; 269 if (mumps_>infog[0]) { 270 COIN_DETAIL_PRINT(printf("MUMPS ordering failed error %d %d\n", 271 mumps_>infog[0], mumps_>infog[1])); 272 return 1; 273 } else { 274 double size = mumps_>infog[19]; 275 if (size < 0) 276 size *= 1000000; 277 COIN_DETAIL_PRINT(printf("%g nonzeros, flop count %g\n", size, mumps_>rinfog[0])); 278 } 279 for (iRow = 0; iRow < numberRows_; iRow++) { 280 permuteInverse_[iRow] = iRow; 281 permute_[iRow] = iRow; 282 } 283 return 0; 287 284 } 288 285 /* Does Symbolic factorization given permutation. … … 290 287 user must provide factorize and solve. Otherwise the default factorization is used 291 288 returns nonzero if not enough memory */ 292 int 293 ClpCholeskyMumps::symbolic() 294 { 295 return 0; 289 int ClpCholeskyMumps::symbolic() 290 { 291 return 0; 296 292 } 297 293 /* Factorize  filling in rowsDropped and returning number dropped */ 298 int 299 ClpCholeskyMumps::factorize(const double * diagonal, int * rowsDropped) 300 { 301 const CoinBigIndex * columnStart = model_>clpMatrix()>getVectorStarts(); 302 const int * columnLength = model_>clpMatrix()>getVectorLengths(); 303 const int * row = model_>clpMatrix()>getIndices(); 304 const double * element = model_>clpMatrix()>getElements(); 305 const CoinBigIndex * rowStart = rowCopy_>getVectorStarts(); 306 const int * rowLength = rowCopy_>getVectorLengths(); 307 const int * column = rowCopy_>getIndices(); 308 const double * elementByRow = rowCopy_>getElements(); 309 int numberColumns = model_>clpMatrix()>getNumCols(); 310 int iRow; 311 double * work = new double[numberRows_]; 312 CoinZeroN(work, numberRows_); 313 const double * diagonalSlack = diagonal + numberColumns; 314 int newDropped = 0; 315 //double smallest; 316 //perturbation 317 double perturbation = model_>diagonalPerturbation() * model_>diagonalNorm(); 318 perturbation = 0.0; 319 perturbation = perturbation * perturbation; 320 if (perturbation > 1.0) { 294 int ClpCholeskyMumps::factorize(const double *diagonal, int *rowsDropped) 295 { 296 const CoinBigIndex *columnStart = model_>clpMatrix()>getVectorStarts(); 297 const int *columnLength = model_>clpMatrix()>getVectorLengths(); 298 const int *row = model_>clpMatrix()>getIndices(); 299 const double *element = model_>clpMatrix()>getElements(); 300 const CoinBigIndex *rowStart = rowCopy_>getVectorStarts(); 301 const int *rowLength = rowCopy_>getVectorLengths(); 302 const int *column = rowCopy_>getIndices(); 303 const double *elementByRow = rowCopy_>getElements(); 304 int numberColumns = model_>clpMatrix()>getNumCols(); 305 int iRow; 306 double *work = new double[numberRows_]; 307 CoinZeroN(work, numberRows_); 308 const double *diagonalSlack = diagonal + numberColumns; 309 int newDropped = 0; 310 //double smallest; 311 //perturbation 312 double perturbation = model_>diagonalPerturbation() * model_>diagonalNorm(); 313 perturbation = 0.0; 314 perturbation = perturbation * perturbation; 315 if (perturbation > 1.0) { 321 316 #ifdef COIN_DEVELOP 322 323 317 //if (model_>model()>logLevel()&4) 318 std::cout << "large perturbation " << perturbation << std::endl; 324 319 #endif 325 perturbation = sqrt(perturbation);; 326 perturbation = 1.0; 327 } 328 double delta2 = model_>delta(); // add delta*delta to diagonal 329 delta2 *= delta2; 330 for (iRow = 0; iRow < numberRows_; iRow++) { 331 double * put = sparseFactor_ + choleskyStart_[iRow]  1; // Fortran 332 int * which = choleskyRow_ + choleskyStart_[iRow]  1; // Fortran 333 int number = choleskyStart_[iRow+1]  choleskyStart_[iRow]; 334 if (!rowLength[iRow]) 335 rowsDropped_[iRow] = 1; 336 if (!rowsDropped_[iRow]) { 337 CoinBigIndex startRow = rowStart[iRow]; 338 CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow]; 339 work[iRow] = diagonalSlack[iRow] + delta2; 340 for (CoinBigIndex k = startRow; k < endRow; k++) { 341 int iColumn = column[k]; 342 if (!whichDense_  !whichDense_[iColumn]) { 343 CoinBigIndex start = columnStart[iColumn]; 344 CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn]; 345 double multiplier = diagonal[iColumn] * elementByRow[k]; 346 for (CoinBigIndex j = start; j < end; j++) { 347 int jRow = row[j]; 348 if (jRow >= iRow && !rowsDropped_[jRow]) { 349 double value = element[j] * multiplier; 350 work[jRow] += value; 351 } 352 } 353 } 354 } 355 int j; 356 for (j = 0; j < number; j++) { 357 int jRow = which[j]  1; // to Fortran 358 put[j] = work[jRow]; 359 work[jRow] = 0.0; 360 } 361 } else { 362 // dropped 363 int j; 364 for (j = 1; j < number; j++) { 365 put[j] = 0.0; 366 } 367 put[0] = 1.0; 320 perturbation = sqrt(perturbation); 321 ; 322 perturbation = 1.0; 323 } 324 double delta2 = model_>delta(); // add delta*delta to diagonal 325 delta2 *= delta2; 326 for (iRow = 0; iRow < numberRows_; iRow++) { 327 double *put = sparseFactor_ + choleskyStart_[iRow]  1; // Fortran 328 int *which = choleskyRow_ + choleskyStart_[iRow]  1; // Fortran 329 int number = choleskyStart_[iRow + 1]  choleskyStart_[iRow]; 330 if (!rowLength[iRow]) 331 rowsDropped_[iRow] = 1; 332 if (!rowsDropped_[iRow]) { 333 CoinBigIndex startRow = rowStart[iRow]; 334 CoinBigIndex endRow = rowStart[iRow] + rowLength[iRow]; 335 work[iRow] = diagonalSlack[iRow] + delta2; 336 for (CoinBigIndex k = startRow; k < endRow; k++) { 337 int iColumn = column[k]; 338 if (!whichDense_  !whichDense_[iColumn]) { 339 CoinBigIndex start = columnStart[iColumn]; 340 CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn]; 341 double multiplier = diagonal[iColumn] * elementByRow[k]; 342 for (CoinBigIndex j = start; j < end; j++) { 343 int jRow = row[j]; 344 if (jRow >= iRow && !rowsDropped_[jRow]) { 345 double value = element[j] * multiplier; 346 work[jRow] += value; 347 } 368 348 } 369 } 370 //check sizes 371 double largest2 = maximumAbsElement(sparseFactor_, sizeFactor_); 372 largest2 *= 1.0e20; 373 int numberDroppedBefore = 0; 374 for (iRow = 0; iRow < numberRows_; iRow++) { 375 int dropped = rowsDropped_[iRow]; 376 // Move to int array 377 rowsDropped[iRow] = dropped; 378 if (!dropped) { 379 int start = choleskyStart_[iRow]  1; // to Fortran 380 double diagonal = sparseFactor_[start]; 381 if (diagonal > largest2) { 382 sparseFactor_[start] = CoinMax(diagonal, 1.0e10); 383 } else { 384 sparseFactor_[start] = CoinMax(diagonal, 1.0e10); 385 rowsDropped[iRow] = 2; 386 numberDroppedBefore++; 387 } 388 } 389 } 390 delete [] work; 391 // code here 392 mumps_>a_elt = sparseFactor_; 393 mumps_>rhs = NULL; 394 mumps_>job = 2; // factorize 395 dmumps_c(mumps_); 396 if (mumps_>infog[0]) { 397 COIN_DETAIL_PRINT(printf("MUMPS factorization failed error %d %d\n", 398 mumps_>infog[0], mumps_>infog[1])); 399 } 400 choleskyCondition_ = 1.0; 401 bool cleanCholesky; 402 if (model_>numberIterations() < 2000) 403 cleanCholesky = true; 404 else 405 cleanCholesky = false; 406 if (cleanCholesky) { 407 //drop fresh makes some formADAT easier 408 //int oldDropped=numberRowsDropped_; 409 if (newDropped  numberRowsDropped_) { 410 //std::cout <<"Rank "<<numberRows_newDropped<<" ( "<< 411 // newDropped<<" dropped)"; 412 //if (newDropped>oldDropped) 413 //std::cout<<" ( "<<newDroppedoldDropped<<" dropped this time)"; 414 //std::cout<<std::endl; 415 newDropped = 0; 416 for (int i = 0; i < numberRows_; i++) { 417 int dropped = rowsDropped[i]; 418 rowsDropped_[i] = (char)dropped; 419 if (dropped == 2) { 420 //dropped this time 421 rowsDropped[newDropped++] = i; 422 rowsDropped_[i] = 0; 423 } 424 } 425 numberRowsDropped_ = newDropped; 426 newDropped = (2 + newDropped); 427 } 428 } else { 429 if (newDropped) { 430 newDropped = 0; 431 for (int i = 0; i < numberRows_; i++) { 432 int dropped = rowsDropped[i]; 433 rowsDropped_[i] = (char)dropped; 434 if (dropped == 2) { 435 //dropped this time 436 rowsDropped[newDropped++] = i; 437 rowsDropped_[i] = 1; 438 } 439 } 440 } 441 numberRowsDropped_ += newDropped; 442 if (numberRowsDropped_ && 0) { 443 std::cout << "Rank " << numberRows_  numberRowsDropped_ << " ( " << 444 numberRowsDropped_ << " dropped)"; 445 if (newDropped) { 446 std::cout << " ( " << newDropped << " dropped this time)"; 447 } 448 std::cout << std::endl; 449 } 450 } 451 status_ = 0; 452 return newDropped; 349 } 350 } 351 int j; 352 for (j = 0; j < number; j++) { 353 int jRow = which[j]  1; // to Fortran 354 put[j] = work[jRow]; 355 work[jRow] = 0.0; 356 } 357 } else { 358 // dropped 359 int j; 360 for (j = 1; j < number; j++) { 361 put[j] = 0.0; 362 } 363 put[0] = 1.0; 364 } 365 } 366 //check sizes 367 double largest2 = maximumAbsElement(sparseFactor_, sizeFactor_); 368 largest2 *= 1.0e20; 369 int numberDroppedBefore = 0; 370 for (iRow = 0; iRow < numberRows_; iRow++) { 371 int dropped = rowsDropped_[iRow]; 372 // Move to int array 373 rowsDropped[iRow] = dropped; 374 if (!dropped) { 375 int start = choleskyStart_[iRow]  1; // to Fortran 376 double diagonal = sparseFactor_[start]; 377 if (diagonal > largest2) { 378 sparseFactor_[start] = CoinMax(diagonal, 1.0e10); 379 } else { 380 sparseFactor_[start] = CoinMax(diagonal, 1.0e10); 381 rowsDropped[iRow] = 2; 382 numberDroppedBefore++; 383 } 384 } 385 } 386 delete[] work; 387 // code here 388 mumps_>a_elt = sparseFactor_; 389 mumps_>rhs = NULL; 390 mumps_>job = 2; // factorize 391 dmumps_c(mumps_); 392 if (mumps_>infog[0]) { 393 COIN_DETAIL_PRINT(printf("MUMPS factorization failed error %d %d\n", 394 mumps_>infog[0], mumps_>infog[1])); 395 } 396 choleskyCondition_ = 1.0; 397 bool cleanCholesky; 398 if (model_>numberIterations() < 2000) 399 cleanCholesky = true; 400 else 401 cleanCholesky = false; 402 if (cleanCholesky) { 403 //drop fresh makes some formADAT easier 404 //int oldDropped=numberRowsDropped_; 405 if (newDropped  numberRowsDropped_) { 406 //std::cout <<"Rank "<<numberRows_newDropped<<" ( "<< 407 // newDropped<<" dropped)"; 408 //if (newDropped>oldDropped) 409 //std::cout<<" ( "<<newDroppedoldDropped<<" dropped this time)"; 410 //std::cout<<std::endl; 411 newDropped = 0; 412 for (int i = 0; i < numberRows_; i++) { 413 int dropped = rowsDropped[i]; 414 rowsDropped_[i] = (char)dropped; 415 if (dropped == 2) { 416 //dropped this time 417 rowsDropped[newDropped++] = i; 418 rowsDropped_[i] = 0; 419 } 420 } 421 numberRowsDropped_ = newDropped; 422 newDropped = (2 + newDropped); 423 } 424 } else { 425 if (newDropped) { 426 newDropped = 0; 427 for (int i = 0; i < numberRows_; i++) { 428 int dropped = rowsDropped[i]; 429 rowsDropped_[i] = (char)dropped; 430 if (dropped == 2) { 431 //dropped this time 432 rowsDropped[newDropped++] = i; 433 rowsDropped_[i] = 1; 434 } 435 } 436 } 437 numberRowsDropped_ += newDropped; 438 if (numberRowsDropped_ && 0) { 439 std::cout << "Rank " << numberRows_  numberRowsDropped_ << " ( " << numberRowsDropped_ << " dropped)"; 440 if (newDropped) { 441 std::cout << " ( " << newDropped << " dropped this time)"; 442 } 443 std::cout << std::endl; 444 } 445 } 446 status_ = 0; 447 return newDropped; 453 448 } 454 449 /* Uses factorization to solve. */ 455 void 456 ClpCholeskyMumps::solve (double * region) 457 { 458 mumps_>rhs = region; 459 mumps_>job = 3; // solve 460 dmumps_c(mumps_); 461 } 450 void ClpCholeskyMumps::solve(double *region) 451 { 452 mumps_>rhs = region; 453 mumps_>job = 3; // solve 454 dmumps_c(mumps_); 455 } 456 457 /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2 458 */
Note: See TracChangeset
for help on using the changeset viewer.