1 // Copyright (C) 2012, The Science and Technology Facilities Council. 2 // Copyright (C) 2009, Jonathan Hogg <jhogg41.at.gmail.com>. 3 // Copyright (C) 2004, 2007 International Business Machines and others. 4 // All Rights Reserved. 5 // This code is published under the Eclipse Public License. 6 // 7 // $Id: IpMa97SolverInterface.cpp 2003 2011-06-03 03:53:44Z andreasw $ 8 // 9 // Authors: Jonathan Hogg STFC 2012-12-21 10 // Jonathan Hogg 2009-07-29 11 // Carl Laird, Andreas Waechter IBM 2004-03-17 12 13 #include "IpoptConfig.h" 14 15 #ifdef COIN_HAS_HSL 16 #include "CoinHslConfig.h" 17 #endif 18 19 // if we have MA97 in HSL or the linear solver loader, then we want to build the MA97 interface 20 #if defined(COINHSL_HAS_MA97) || defined(HAVE_LINEARSOLVERLOADER) 21 22 #include "IpMa97SolverInterface.hpp" 23 #include <iostream> 24 #include <stdio.h> 25 #include <cmath> 26 using namespace std; 27 28 /* Uncomment the following line to enable the ma97_dump_matrix option. 29 * This option requires a version of the coinhsl that supports this function. 30 */ 31 //#define MA97_DUMP_MATRIX 32 33 #ifdef MA97_DUMP_MATRIX 34 extern "C" { 35 extern void F77_FUNC (dump_mat_csc, DUMP_MAT_CSC) ( 36 const ipfint *factidx, 37 const ipfint *n, 38 const ipfint *ptr, 39 const ipfint *row, 40 const double *a); 41 } 42 #endif 43 44 namespace Ipopt 45 { 46 ~Ma97SolverInterface()47 Ma97SolverInterface::~Ma97SolverInterface() 48 { 49 delete [] val_; 50 if(scaling_) delete [] scaling_; 51 52 ma97_finalise(&akeep_, &fkeep_); 53 } 54 RegisterOptions(SmartPtr<RegisteredOptions> roptions)55 void Ma97SolverInterface::RegisterOptions(SmartPtr<RegisteredOptions> roptions) 56 { 57 roptions->AddIntegerOption( 58 "ma97_print_level", 59 "Debug printing level for the linear solver MA97", 60 0, "" 61 /* 62 "<0 no printing.\n" 63 "0 Error and warning messages only.\n" 64 "=1 Limited diagnostic printing.\n" 65 ">1 Additional diagnostic printing."*/); 66 roptions->AddLowerBoundedIntegerOption( 67 "ma97_nemin", 68 "Node Amalgamation parameter", 69 1, 8, 70 "Two nodes in elimination tree are merged if result has fewer than " 71 "ma97_nemin variables."); 72 roptions->AddLowerBoundedNumberOption( 73 "ma97_small", 74 "Zero Pivot Threshold", 75 0.0, false, 1e-20, 76 "Any pivot less than ma97_small is treated as zero."); 77 roptions->AddBoundedNumberOption( 78 "ma97_u", 79 "Pivoting Threshold", 80 0.0, false, 0.5, false, 1e-8, 81 "See MA97 documentation."); 82 roptions->AddBoundedNumberOption( 83 "ma97_umax", 84 "Maximum Pivoting Threshold", 85 0.0, false, 0.5, false, 1e-4, 86 "See MA97 documentation."); 87 roptions->AddStringOption5( 88 "ma97_scaling", 89 "Specifies strategy for scaling in HSL_MA97 linear solver", 90 "dynamic", 91 "none", "Do not scale the linear system matrix", 92 "mc30", "Scale all linear system matrices using MC30", 93 "mc64", "Scale all linear system matrices using MC64", 94 "mc77", "Scale all linear system matrices using MC77 [1,3,0]", 95 "dynamic", "Dynamically select scaling according to rules specified by ma97_scalingX and ma97_switchX options.", 96 ""); 97 roptions->AddStringOption4( 98 "ma97_scaling1", 99 "First scaling.", 100 "mc64", 101 "none", "No scaling", 102 "mc30", "Scale linear system matrix using MC30", 103 "mc64", "Scale linear system matrix using MC64", 104 "mc77", "Scale linear system matrix using MC77 [1,3,0]", 105 "If ma97_scaling=dynamic, this scaling is used according to the trigger " 106 "ma97_switch1. If ma97_switch2 is triggered it is disabled."); 107 roptions->AddStringOption9( 108 "ma97_switch1", 109 "First switch, determine when ma97_scaling1 is enabled.", 110 "od_hd_reuse", 111 "never", "Scaling is never enabled.", 112 "at_start", "Scaling to be used from the very start.", 113 "at_start_reuse", "Scaling to be used on first iteration, then reused thereafter.", 114 "on_demand", "Scaling to be used after Ipopt request improved solution (i.e. iterative refinement has failed).", 115 "on_demand_reuse", "As on_demand, but reuse scaling from previous itr", 116 "high_delay", "Scaling to be used after more than 0.05*n delays are present", 117 "high_delay_reuse", "Scaling to be used only when previous itr created more that 0.05*n additional delays, otherwise reuse scaling from previous itr", 118 "od_hd", "Combination of on_demand and high_delay", 119 "od_hd_reuse", "Combination of on_demand_reuse and high_delay_reuse", 120 "If ma97_scaling=dynamic, ma97_scaling1 is enabled according to this condition. If ma97_switch2 occurs this option is henceforth ignored."); 121 roptions->AddStringOption4( 122 "ma97_scaling2", 123 "Second scaling.", 124 "mc64", 125 "none", "No scaling", 126 "mc30", "Scale linear system matrix using MC30", 127 "mc64", "Scale linear system matrix using MC64", 128 "mc77", "Scale linear system matrix using MC77 [1,3,0]", 129 "If ma97_scaling=dynamic, this scaling is used according to the trigger " 130 "ma97_switch2. If ma97_switch3 is triggered it is disabled."); 131 roptions->AddStringOption9( 132 "ma97_switch2", 133 "Second switch, determine when ma97_scaling2 is enabled.", 134 "never", 135 "never", "Scaling is never enabled.", 136 "at_start", "Scaling to be used from the very start.", 137 "at_start_reuse", "Scaling to be used on first iteration, then reused thereafter.", 138 "on_demand", "Scaling to be used after Ipopt request improved solution (i.e. iterative refinement has failed).", 139 "on_demand_reuse", "As on_demand, but reuse scaling from previous itr", 140 "high_delay", "Scaling to be used after more than 0.05*n delays are present", 141 "high_delay_reuse", "Scaling to be used only when previous itr created more that 0.05*n additional delays, otherwise reuse scaling from previous itr", 142 "od_hd", "Combination of on_demand and high_delay", 143 "od_hd_reuse", "Combination of on_demand_reuse and high_delay_reuse", 144 "If ma97_scaling=dynamic, ma97_scaling2 is enabled according to this condition. If ma97_switch3 occurs this option is henceforth ignored."); 145 roptions->AddStringOption4( 146 "ma97_scaling3", 147 "Third scaling.", 148 "mc64", 149 "none", "No scaling", 150 "mc30", "Scale linear system matrix using MC30", 151 "mc64", "Scale linear system matrix using MC64", 152 "mc77", "Scale linear system matrix using MC77 [1,3,0]", 153 "If ma97_scaling=dynamic, this scaling is used according to the trigger " 154 "ma97_switch3."); 155 roptions->AddStringOption9( 156 "ma97_switch3", 157 "Third switch, determine when ma97_scaling3 is enabled.", 158 "never", 159 "never", "Scaling is never enabled.", 160 "at_start", "Scaling to be used from the very start.", 161 "at_start_reuse", "Scaling to be used on first iteration, then reused thereafter.", 162 "on_demand", "Scaling to be used after Ipopt request improved solution (i.e. iterative refinement has failed).", 163 "on_demand_reuse", "As on_demand, but reuse scaling from previous itr", 164 "high_delay", "Scaling to be used after more than 0.05*n delays are present", 165 "high_delay_reuse", "Scaling to be used only when previous itr created more that 0.05*n additional delays, otherwise reuse scaling from previous itr", 166 "od_hd", "Combination of on_demand and high_delay", 167 "od_hd_reuse", "Combination of on_demand_reuse and high_delay_reuse", 168 "If ma97_scaling=dynamic, ma97_scaling3 is enabled according to this condition."); 169 roptions->AddStringOption7( 170 "ma97_order", 171 "Controls type of ordering used by HSL_MA97", 172 "auto", 173 "auto", "Use HSL_MA97 heuristic to guess best of AMD and METIS", 174 "best", "Try both AMD and MeTiS, pick best", 175 "amd", "Use the HSL_MC68 approximate minimum degree algorithm", 176 "metis", "Use the MeTiS nested dissection algorithm", 177 "matched-auto", "Use the HSL_MC80 matching with heuristic choice of AMD or METIS", 178 "matched-metis", "Use the HSL_MC80 matching based ordering with METIS", 179 "matched-amd", "Use the HSL_MC80 matching based ordering with AMD", 180 ""); 181 #ifdef MA97_DUMP_MATRIX 182 roptions->AddStringOption2( 183 "ma97_dump_matrix", 184 "Controls whether HSL_MA97 dumps each matrix to a file", 185 "no", 186 "no", "Do not dump matrix", 187 "yes", "Do dump matrix", 188 ""); 189 #endif 190 roptions->AddStringOption2( 191 "ma97_solve_blas3", 192 "Controls if blas2 or blas3 routines are used for solve", 193 "no", 194 "no", "Use BLAS2 (faster, some implementations bit incompatible)", 195 "yes", "Use BLAS3 (slower)", 196 ""); 197 } 198 ScaleNameToNum(const std::string & name)199 int Ma97SolverInterface::ScaleNameToNum(const std::string& name) { 200 if(name=="none") return 0; 201 if(name=="mc64") return 1; 202 if(name=="mc77") return 2; 203 if(name=="mc30") return 4; 204 assert(0); 205 return -1; 206 } 207 InitializeImpl(const OptionsList & options,const std::string & prefix)208 bool Ma97SolverInterface::InitializeImpl(const OptionsList& options, 209 const std::string& prefix) 210 { 211 ma97_default_control(&control_); 212 control_.f_arrays = 1; // Use Fortran numbering (faster) 213 control_.action = 0; // false, shuold exit with error on singularity 214 215 options.GetIntegerValue("ma97_print_level", control_.print_level, 216 prefix); 217 options.GetIntegerValue("ma97_nemin", control_.nemin, prefix); 218 options.GetNumericValue("ma97_small", control_.small, prefix); 219 options.GetNumericValue("ma97_u", control_.u, prefix); 220 options.GetNumericValue("ma97_umax", umax_, prefix); 221 std::string order_method, scaling_method, rescale_strategy; 222 options.GetStringValue("ma97_order", order_method, prefix); 223 if(order_method == "metis") { 224 ordering_ = ORDER_METIS; 225 } else if(order_method == "amd") { 226 ordering_ = ORDER_AMD; 227 } else if(order_method == "best") { 228 ordering_ = ORDER_BEST; 229 } else if(order_method == "matched-metis") { 230 ordering_ = ORDER_MATCHED_METIS; 231 } else if(order_method == "matched-amd") { 232 ordering_ = ORDER_MATCHED_AMD; 233 } else if(order_method == "matched-auto") { 234 ordering_ = ORDER_MATCHED_AUTO; 235 } else { 236 ordering_ = ORDER_AUTO; 237 } 238 options.GetStringValue("ma97_scaling", scaling_method, prefix); 239 current_level_ = 0; 240 if(scaling_method == "dynamic") { 241 scaling_type_ = 0; 242 string switch_val[3], scale_val[3]; 243 options.GetStringValue("ma97_switch1", switch_val[0], prefix); 244 options.GetStringValue("ma97_scaling1", scale_val[0], prefix); 245 options.GetStringValue("ma97_switch2", switch_val[1], prefix); 246 options.GetStringValue("ma97_scaling2", scale_val[1], prefix); 247 options.GetStringValue("ma97_switch3", switch_val[2], prefix); 248 options.GetStringValue("ma97_scaling3", scale_val[2], prefix); 249 for(int i=0; i<3; i++) { 250 scaling_val_[i] = ScaleNameToNum(scale_val[i]); 251 if(switch_val[i]=="never") 252 switch_[i] = SWITCH_NEVER; 253 else if(switch_val[i]=="at_start") { 254 switch_[i] = SWITCH_AT_START; 255 scaling_type_ = scaling_val_[i]; 256 current_level_ = i; 257 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 258 "HSL_MA97: Enabled scaling level %d on initialization\n", current_level_); 259 } 260 else if(switch_val[i]=="at_start_reuse") { 261 switch_[i] = SWITCH_AT_START_REUSE; 262 scaling_type_ = scaling_val_[i]; 263 current_level_ = i; 264 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 265 "HSL_MA97: Enabled scaling level %d on initialization\n", current_level_); 266 } 267 else if(switch_val[i]=="on_demand") 268 switch_[i] = SWITCH_ON_DEMAND; 269 else if(switch_val[i]=="on_demand_reuse") 270 switch_[i] = SWITCH_ON_DEMAND_REUSE; 271 else if(switch_val[i]=="high_delay") 272 switch_[i] = SWITCH_NDELAY; 273 else if(switch_val[i]=="high_delay_reuse") 274 switch_[i] = SWITCH_NDELAY_REUSE; 275 else if(switch_val[i]=="od_hd") 276 switch_[i] = SWITCH_OD_ND; 277 else if(switch_val[i]=="od_hd_reuse") 278 switch_[i] = SWITCH_OD_ND_REUSE; 279 } 280 } else { 281 switch_[0] = SWITCH_AT_START; 282 switch_[1] = SWITCH_NEVER; 283 switch_[2] = SWITCH_NEVER; 284 scaling_type_ = ScaleNameToNum(scaling_method); 285 } 286 #ifdef MA97_DUMP_MATRIX 287 options.GetBoolValue("ma97_dump_matrix", dump_, prefix); 288 #endif 289 bool solve_blas3; 290 options.GetBoolValue("ma97_solve_blas3", solve_blas3, prefix); 291 control_.solve_blas3 = solve_blas3 ? 1 : 0; 292 293 // Set whether we scale on first iteration or not 294 switch(switch_[current_level_]) { 295 case SWITCH_NEVER: 296 case SWITCH_ON_DEMAND: 297 case SWITCH_ON_DEMAND_REUSE: 298 case SWITCH_NDELAY: 299 case SWITCH_NDELAY_REUSE: 300 case SWITCH_OD_ND: 301 case SWITCH_OD_ND_REUSE: 302 rescale_ = false; 303 break; 304 case SWITCH_AT_START: 305 case SWITCH_AT_START_REUSE: 306 rescale_ = true; 307 break; 308 } 309 // Set scaling 310 control_.scaling = scaling_type_; 311 312 return true; // All is well 313 } 314 315 /* Method for initializing internal stuctures. Here, ndim gives 316 * the number of rows and columns of the matrix, nonzeros give 317 * the number of nonzero elements, and ia and ja give the 318 * positions of the nonzero elements, given in the matrix format 319 * determined by MatrixFormat. 320 */ InitializeStructure(Index dim,Index nonzeros,const Index * ia,const Index * ja)321 ESymSolverStatus Ma97SolverInterface::InitializeStructure(Index dim, 322 Index nonzeros, const Index* ia, const Index* ja) 323 { 324 struct ma97_info info, info2; 325 void *akeep_amd, *akeep_metis; 326 327 // Store size for later use 328 ndim_ = dim; 329 330 // Setup memory for values 331 if (val_!=NULL) delete[] val_; 332 val_ = new double[nonzeros]; 333 334 // Check if analyse needs to be postponed 335 if(ordering_ == ORDER_MATCHED_AMD || ordering_ == ORDER_MATCHED_METIS) { 336 // Ordering requires values. Just signal success and return 337 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 338 "HSL_MA97: Delaying analyse until values are available\n"); 339 switch(ordering_) 340 { 341 case ORDER_MATCHED_AMD: 342 control_.ordering = 7; // HSL_MC80 with AMD 343 break; 344 case ORDER_MATCHED_METIS: 345 control_.ordering = 8; // HSL_MC80 with METIS 346 break; 347 } 348 return SYMSOLVER_SUCCESS; 349 } 350 351 if (HaveIpData()) { 352 IpData().TimingStats().LinearSystemSymbolicFactorization().Start(); 353 } 354 355 // perform analyse 356 if(ordering_ == ORDER_BEST) { 357 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 358 "HSL_MA97: Use best of AMD or MeTiS:\n"); 359 control_.ordering = 1; // AMD 360 ma97_analyse(0, dim, ia, ja, NULL, &akeep_amd, &control_, &info2, NULL); 361 if (info2.flag<0) return SYMSOLVER_FATAL_ERROR; 362 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 363 "AMD nfactor = %d, nflops = %d:\n", 364 info2.num_factor, info2.num_flops); 365 control_.ordering = 3; // METIS 366 ma97_analyse(0, dim, ia, ja, NULL, &akeep_metis, &control_, &info, NULL); 367 if (info.flag<0) return SYMSOLVER_FATAL_ERROR; 368 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 369 "MeTiS nfactor = %d, nflops = %d:\n", 370 info.num_factor, info.num_flops); 371 if(info.num_flops > info2.num_flops) { 372 // Use AMD 373 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 374 "HSL_MA97: Choose AMD\n"); 375 akeep_ = akeep_amd; 376 ma97_free_akeep(&akeep_metis); 377 info = info2; 378 } else { 379 // Use MeTiS 380 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 381 "HSL_MA97: Choose MeTiS\n"); 382 akeep_ = akeep_metis; 383 ma97_free_akeep(&akeep_amd); 384 } 385 } else { 386 switch(ordering_) 387 { 388 case ORDER_AMD: 389 case ORDER_MATCHED_AMD: 390 control_.ordering = 1; // AMD 391 break; 392 case ORDER_METIS: 393 case ORDER_MATCHED_METIS: 394 control_.ordering = 3; // METIS 395 break; 396 case ORDER_AUTO: 397 case ORDER_MATCHED_AUTO: 398 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 399 "HSL_MA97: Make heuristic choice of AMD or MeTiS\n"); 400 control_.ordering = 5; // Use heuristic to pick which to use 401 } 402 ma97_analyse(0, dim, ia, ja, NULL, &akeep_, &control_, &info, NULL); 403 switch(info.ordering) { 404 case 1: 405 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 406 "HSL_MA97: Used AMD\n"); 407 if(ordering_ == ORDER_MATCHED_AUTO) 408 ordering_ = ORDER_MATCHED_AMD; 409 break; 410 case 3: 411 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 412 "HSL_MA97: Used MeTiS\n"); 413 if(ordering_ == ORDER_MATCHED_AUTO) 414 ordering_ = ORDER_MATCHED_METIS; 415 break; 416 default: 417 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 418 "HSL_MA97: Used ordering %d\n", info.ordering); 419 break; 420 } 421 } 422 423 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 424 "HSL_MA97: PREDICTED nfactor %d, maxfront %d\n", 425 info.num_factor, info.maxfront); 426 427 if (HaveIpData()) { 428 IpData().TimingStats().LinearSystemSymbolicFactorization().End(); 429 } 430 431 if (info.flag>=0) { 432 return SYMSOLVER_SUCCESS; 433 } 434 else { 435 return SYMSOLVER_FATAL_ERROR; 436 } 437 } 438 439 /* Solve operation for multiple right hand sides. Solves the 440 * linear system A * x = b with multiple right hand sides, where 441 * A is the symmtric indefinite matrix. Here, ia and ja give the 442 * positions of the values (in the required matrix data format). 443 * The actual values of the matrix will have been given to this 444 * object by copying them into the array provided by 445 * GetValuesArrayPtr. ia and ja are identical to the ones given 446 * to InitializeStructure. The flag new_matrix is set to true, 447 * if the values of the matrix has changed, and a refactorzation 448 * is required. 449 * 450 * The return code is SYMSOLV_SUCCESS if the factorization and 451 * solves were successful, SYMSOLV_SINGULAR if the linear system 452 * is singular, and SYMSOLV_WRONG_INERTIA if check_NegEVals is 453 * true and the number of negative eigenvalues in the matrix does 454 * not match numberOfNegEVals. If SYMSOLV_CALL_AGAIN is 455 * returned, then the calling function will request the pointer 456 * for the array for storing a again (with GetValuesPtr), write 457 * the values of the nonzero elements into it, and call this 458 * MultiSolve method again with the same right-hand sides. (This 459 * can be done, for example, if the linear solver realized it 460 * does not have sufficient memory and needs to redo the 461 * factorization; e.g., for MA27.) 462 * 463 * The number of right-hand sides is given by nrhs, the values of 464 * the right-hand sides are given in rhs_vals (one full right-hand 465 * side stored immediately after the other), and solutions are 466 * to be returned in the same array. 467 * 468 * check_NegEVals will not be chosen true, if ProvidesInertia() 469 * returns false. 470 */ MultiSolve(bool new_matrix,const Index * ia,const Index * ja,Index nrhs,double * rhs_vals,bool check_NegEVals,Index numberOfNegEVals)471 ESymSolverStatus Ma97SolverInterface::MultiSolve(bool new_matrix, 472 const Index* ia, const Index* ja, Index nrhs, double* rhs_vals, 473 bool check_NegEVals, Index numberOfNegEVals) 474 { 475 struct ma97_info info; 476 Number t1=0,t2; 477 478 if (new_matrix || pivtol_changed_) { 479 480 #ifdef MA97_DUMP_MATRIX 481 if(dump_) { 482 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 483 "Dumping matrix %d\n", fctidx_); 484 F77_FUNC (dump_mat_csc, DUMP_MAT_CSC) 485 (&fctidx_, &ndim_, ia, ja, val_); 486 fctidx_++; 487 } 488 #endif 489 490 // Set scaling option 491 if(rescale_) { 492 control_.scaling = scaling_type_; 493 if(scaling_type_!=0 && scaling_==NULL) 494 scaling_ = new double[ndim_]; // alloc if not already 495 } else { 496 control_.scaling = 0; // None or user (depends if scaling_ is alloc'd) 497 } 498 499 if((ordering_ == ORDER_MATCHED_AMD || ordering_ == ORDER_MATCHED_METIS) && rescale_) { 500 /* 501 * Perform delayed analyse 502 */ 503 if (HaveIpData()) { 504 IpData().TimingStats().LinearSystemSymbolicFactorization().Start(); 505 } 506 507 switch(ordering_) 508 { 509 case ORDER_MATCHED_AMD: 510 control_.ordering = 7; // HSL_MC80 with AMD 511 break; 512 case ORDER_MATCHED_METIS: 513 control_.ordering = 8; // HSL_MC80 with METIS 514 break; 515 } 516 517 ma97_analyse(0, ndim_, ia, ja, val_, &akeep_, &control_, &info, NULL); 518 if(scaling_type_ == 1) 519 control_.scaling = 3; // use mc64 from ordering 520 521 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 522 "HSL_MA97: PREDICTED nfactor %d, maxfront %d\n", 523 info.num_factor, info.maxfront); 524 525 if (HaveIpData()) { 526 IpData().TimingStats().LinearSystemSymbolicFactorization().End(); 527 } 528 529 if (info.flag==6 || info.flag==-7) { 530 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 531 "In Ma97SolverInterface::Factorization: " 532 "Singular system, estimated rank %d of %d\n", 533 info.matrix_rank, ndim_); 534 return SYMSOLVER_SINGULAR; 535 } 536 if (info.flag<0) return SYMSOLVER_FATAL_ERROR; 537 538 } 539 540 if (HaveIpData()) { 541 t1 = IpData().TimingStats().LinearSystemFactorization().TotalWallclockTime(); 542 IpData().TimingStats().LinearSystemFactorization().Start(); 543 } 544 ma97_factor(4, ia, ja, val_, &akeep_, &fkeep_, &control_, &info, scaling_); 545 //ma97_factor_solve(4, ia, ja, val_, nrhs, rhs_vals, ndim_, &akeep_, &fkeep_, 546 // &control_, &info, scaling_); 547 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 548 "HSL_MA97: delays %d, nfactor %d, nflops %ld, maxfront %d\n", 549 info.num_delay, info.num_factor, info.num_flops, info.maxfront); 550 if (HaveIpData()) { 551 IpData().TimingStats().LinearSystemFactorization().End(); 552 t2 = IpData().TimingStats().LinearSystemFactorization().TotalWallclockTime(); 553 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 554 "Ma97SolverInterface::Factorization: " 555 "ma97_factor_solve took %10.3f\n", 556 t2-t1); 557 } 558 if (info.flag==7 || info.flag==-7) { 559 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 560 "In Ma97SolverInterface::Factorization: " 561 "Singular system, estimated rank %d of %d\n", 562 info.matrix_rank, ndim_); 563 return SYMSOLVER_SINGULAR; 564 } 565 for(int i=current_level_; i<3; i++) { 566 switch(switch_[i]) { 567 case SWITCH_NEVER: 568 case SWITCH_AT_START: 569 case SWITCH_ON_DEMAND: 570 // Nothing to do here 571 break; 572 case SWITCH_AT_START_REUSE: 573 rescale_ = false; // Scaled exactly once, never changed again 574 break; 575 case SWITCH_ON_DEMAND_REUSE: 576 if(i==current_level_ && rescale_) rescale_ = false; 577 break; 578 case SWITCH_NDELAY_REUSE: 579 case SWITCH_OD_ND_REUSE: 580 if(rescale_) numdelay_ = info.num_delay; // Need to do this before we reset rescale_ 581 if(i==current_level_ && rescale_) rescale_ = false; 582 // Falls through to: 583 case SWITCH_NDELAY: 584 case SWITCH_OD_ND: 585 if(rescale_) numdelay_ = info.num_delay; 586 if(info.num_delay - numdelay_ > 0.05*ndim_) { 587 // number of delays has signficantly increased, so trigger 588 current_level_ = i; 589 scaling_type_ = scaling_val_[i]; 590 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 591 "HSL_MA97: Enabling scaling %d due to excess delays\n", i); 592 rescale_ = true; 593 } 594 break; 595 } 596 } 597 if (info.flag<0) { 598 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 599 "In Ma97SolverInterface::Factorization: " 600 "Unhandled error. info.flag = %d\n", 601 info.flag); 602 return SYMSOLVER_FATAL_ERROR; 603 } 604 if (check_NegEVals && info.num_neg!=numberOfNegEVals) { 605 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 606 "In Ma97SolverInterface::Factorization: " 607 "info.num_neg = %d, but numberOfNegEVals = %d\n", 608 info.num_neg, numberOfNegEVals); 609 return SYMSOLVER_WRONG_INERTIA; 610 } 611 612 if (HaveIpData()) { 613 IpData().TimingStats().LinearSystemBackSolve().Start(); 614 } 615 ma97_solve(0, nrhs, rhs_vals, ndim_, &akeep_, &fkeep_, &control_, &info); 616 if (HaveIpData()) { 617 IpData().TimingStats().LinearSystemBackSolve().End(); 618 } 619 620 numneg_ = info.num_neg; 621 pivtol_changed_ = false; 622 } 623 else { 624 if (HaveIpData()) { 625 IpData().TimingStats().LinearSystemBackSolve().Start(); 626 } 627 ma97_solve(0, nrhs, rhs_vals, ndim_, &akeep_, &fkeep_, &control_, &info); 628 if (HaveIpData()) { 629 IpData().TimingStats().LinearSystemBackSolve().End(); 630 } 631 } 632 633 return SYMSOLVER_SUCCESS; 634 } 635 636 IncreaseQuality()637 bool Ma97SolverInterface::IncreaseQuality() 638 { 639 for(int i=current_level_; i<3; i++) { 640 switch(switch_[i]) { 641 case SWITCH_ON_DEMAND: 642 case SWITCH_ON_DEMAND_REUSE: 643 case SWITCH_OD_ND: 644 case SWITCH_OD_ND_REUSE: 645 rescale_ = true; 646 current_level_ = i; 647 scaling_type_ = scaling_val_[i]; 648 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 649 "HSL_MA97: Enabling scaling %d due to failure of iterative refinement\n", current_level_); 650 break; 651 default: ; 652 } 653 } 654 655 if (control_.u >= umax_) { 656 return false; 657 } 658 pivtol_changed_ = true; 659 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 660 "Indreasing pivot tolerance for HSL_MA97 from %7.2e ", 661 control_.u); 662 control_.u = Min(umax_, pow(control_.u,0.75)); 663 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 664 "to %7.2e.\n", 665 control_.u); 666 return true; 667 } 668 669 } // namespace Ipopt 670 671 #endif /* COINHSL_HAS_MA97 or HAVE_LINEARSOLVERLOADER */ 672