1 // Copyright (C) 2011, Science and Technology Facilities Council. 2 // Copyright (C) 2009, Jonathan Hogg <jdh41.at.cantab.net>. 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$ 8 // 9 // Authors: Jonathan Hogg STFC 2011-03-15 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 do not have MA86 in HSL or the linear solver loader, then we want to build the MA86 interface 20 #if defined(COINHSL_HAS_MA86) || defined(HAVE_LINEARSOLVERLOADER) 21 22 #include "IpMa86SolverInterface.hpp" 23 #include <iostream> 24 #include <cmath> 25 using namespace std; 26 27 extern "C" 28 { 29 #include "hsl_mc68i.h" 30 } 31 32 namespace Ipopt 33 { 34 ~Ma86SolverInterface()35 Ma86SolverInterface::~Ma86SolverInterface() 36 { 37 delete [] val_; 38 delete [] order_; 39 40 if(keep_) ma86_finalise(&keep_, &control_); 41 } 42 RegisterOptions(SmartPtr<RegisteredOptions> roptions)43 void Ma86SolverInterface::RegisterOptions(SmartPtr<RegisteredOptions> roptions) 44 { 45 roptions->AddIntegerOption( 46 "ma86_print_level", 47 "Debug printing level for the linear solver MA86", 48 -1, 49 ""); 50 /* 51 "<0 no printing.\n" 52 "0 Error and warning messages only.\n" 53 "=1 Limited diagnostic printing.\n" 54 ">1 Additional diagnostic printing."); 55 */ 56 roptions->AddLowerBoundedIntegerOption( 57 "ma86_nemin", 58 "Node Amalgamation parameter", 59 1, 32, 60 "Two nodes in elimination tree are merged if result has fewer than " 61 "ma86_nemin variables."); 62 roptions->AddLowerBoundedNumberOption( 63 "ma86_small", 64 "Zero Pivot Threshold", 65 0.0, false, 1e-20, 66 "Any pivot less than ma86_small is treated as zero."); 67 roptions->AddLowerBoundedNumberOption( 68 "ma86_static", 69 "Static Pivoting Threshold", 70 0.0, false, 0.0, 71 "See MA86 documentation. Either ma86_static=0.0 or " 72 "ma86_static>ma86_small. ma86_static=0.0 disables static pivoting."); 73 roptions->AddBoundedNumberOption( 74 "ma86_u", 75 "Pivoting Threshold", 76 0.0, false, 0.5, false, 1e-8, 77 "See MA86 documentation."); 78 roptions->AddBoundedNumberOption( 79 "ma86_umax", 80 "Maximum Pivoting Threshold", 81 0.0, false, 0.5, false, 1e-4, 82 "Maximum value to which u will be increased to improve quality."); 83 roptions->AddStringOption3( 84 "ma86_scaling", 85 "Controls scaling of matrix", 86 "mc64", 87 "none", "Do not scale the linear system matrix", 88 "mc64", "Scale linear system matrix using MC64", 89 "mc77", "Scale linear system matrix using MC77 [1,3,0]", 90 "This option controls scaling for the solver HSL_MA86."); 91 roptions->AddStringOption3( 92 "ma86_order", 93 "Controls type of ordering used by HSL_MA86", 94 #ifdef COINHSL_HAS_METIS 95 "auto", 96 #else 97 "amd", 98 #endif 99 "auto", "Try both AMD and MeTiS, pick best", 100 "amd", "Use the HSL_MC68 approximate minimum degree algorithm", 101 "metis", "Use the MeTiS nested dissection algorithm (if available)", 102 "This option controls ordering for the solver HSL_MA86."); 103 } 104 InitializeImpl(const OptionsList & options,const std::string & prefix)105 bool Ma86SolverInterface::InitializeImpl(const OptionsList& options, 106 const std::string& prefix) 107 { 108 ma86_default_control(&control_); 109 control_.f_arrays = 1; // Use Fortran numbering (faster) 110 /* Note: we can't set control_.action = false as we need to know the 111 * intertia. (Otherwise we just enter the restoration phase and fail) */ 112 113 options.GetIntegerValue("ma86_print_level", control_.diagnostics_level, 114 prefix); 115 options.GetIntegerValue("ma86_nemin", control_.nemin, prefix); 116 options.GetNumericValue("ma86_small", control_.small_, prefix); 117 options.GetNumericValue("ma86_static", control_.static_, prefix); 118 options.GetNumericValue("ma86_u", control_.u, prefix); 119 options.GetNumericValue("ma86_umax", umax_, prefix); 120 std::string order_method, scaling_method; 121 options.GetStringValue("ma86_order", order_method, prefix); 122 if(order_method == "metis") { 123 ordering_ = ORDER_METIS; 124 } else if(order_method == "amd") { 125 ordering_ = ORDER_AMD; 126 } else { 127 ordering_ = ORDER_AUTO; 128 } 129 options.GetStringValue("ma86_scaling", scaling_method, prefix); 130 if(scaling_method == "mc64") { 131 control_.scaling = 1; 132 } else if(scaling_method == "mc77") { 133 control_.scaling = 2; 134 } else { 135 control_.scaling = 0; 136 } 137 138 return true; // All is well 139 } 140 141 /* Method for initializing internal stuctures. Here, ndim gives 142 * the number of rows and columns of the matrix, nonzeros give 143 * the number of nonzero elements, and ia and ja give the 144 * positions of the nonzero elements, given in the matrix format 145 * determined by MatrixFormat. 146 */ InitializeStructure(Index dim,Index nonzeros,const Index * ia,const Index * ja)147 ESymSolverStatus Ma86SolverInterface::InitializeStructure(Index dim, 148 Index nonzeros, const Index* ia, const Index* ja) 149 { 150 struct ma86_info info, info2; 151 struct mc68_control control68; 152 struct mc68_info info68; 153 Index *order_amd, *order_metis; 154 void *keep_amd, *keep_metis; 155 156 // Store size for later use 157 ndim_ = dim; 158 159 // Determine an ordering 160 mc68_default_control(&control68); 161 control68.f_array_in = 1; // Use Fortran numbering (faster) 162 control68.f_array_out = 1; // Use Fortran numbering (faster) 163 order_amd = NULL; order_metis = NULL; 164 if(ordering_ == ORDER_METIS || ordering_ == ORDER_AUTO) { 165 order_metis = new Index[dim]; 166 mc68_order(3, dim, ia, ja, order_metis, &control68, &info68); /* MeTiS */ 167 if(info68.flag == -5) { 168 // MeTiS not available 169 ordering_ = ORDER_AMD; 170 delete[] order_metis; 171 } else if(info68.flag<0) { 172 return SYMSOLVER_FATAL_ERROR; 173 } 174 } 175 if(ordering_ == ORDER_AMD || ordering_ == ORDER_AUTO) { 176 order_amd = new Index[dim]; 177 mc68_order(1, dim, ia, ja, order_amd, &control68, &info68); /* AMD */ 178 } 179 if(info68.flag<0) return SYMSOLVER_FATAL_ERROR; 180 181 if (HaveIpData()) { 182 IpData().TimingStats().LinearSystemSymbolicFactorization().Start(); 183 } 184 185 // perform analyse 186 if(ordering_ == ORDER_AUTO) { 187 ma86_analyse(dim, ia, ja, order_amd, &keep_amd, &control_, &info2); 188 if (info2.flag<0) return SYMSOLVER_FATAL_ERROR; 189 ma86_analyse(dim, ia, ja, order_metis, &keep_metis, &control_, &info); 190 if (info.flag<0) return SYMSOLVER_FATAL_ERROR; 191 if(info.num_flops > info2.num_flops) { 192 // Use AMD 193 //cout << "Choose AMD\n"; 194 order_ = order_amd; 195 keep_ = keep_amd; 196 delete[] order_metis; 197 ma86_finalise(&keep_metis, &control_); 198 } else { 199 // Use MeTiS 200 //cout << "Choose MeTiS\n"; 201 order_ = order_metis; 202 keep_ = keep_metis; 203 delete[] order_amd; 204 ma86_finalise(&keep_amd, &control_); 205 } 206 } else { 207 if(ordering_ == ORDER_AMD) order_ = order_amd; 208 if(ordering_ == ORDER_METIS) order_ = order_metis; 209 ma86_analyse(dim, ia, ja, order_, &keep_, &control_, &info); 210 } 211 212 if (HaveIpData()) { 213 IpData().TimingStats().LinearSystemSymbolicFactorization().End(); 214 } 215 216 // Setup memory for values 217 if (val_!=NULL) delete[] val_; 218 val_ = new double[nonzeros]; 219 220 if (info.flag>=0) { 221 return SYMSOLVER_SUCCESS; 222 } 223 else { 224 return SYMSOLVER_FATAL_ERROR; 225 } 226 } 227 228 /* Solve operation for multiple right hand sides. Solves the 229 * linear system A * x = b with multiple right hand sides, where 230 * A is the symmtric indefinite matrix. Here, ia and ja give the 231 * positions of the values (in the required matrix data format). 232 * The actual values of the matrix will have been given to this 233 * object by copying them into the array provided by 234 * GetValuesArrayPtr. ia and ja are identical to the ones given 235 * to InitializeStructure. The flag new_matrix is set to true, 236 * if the values of the matrix has changed, and a refactorzation 237 * is required. 238 * 239 * The return code is SYMSOLV_SUCCESS if the factorization and 240 * solves were successful, SYMSOLV_SINGULAR if the linear system 241 * is singular, and SYMSOLV_WRONG_INERTIA if check_NegEVals is 242 * true and the number of negative eigenvalues in the matrix does 243 * not match numberOfNegEVals. If SYMSOLV_CALL_AGAIN is 244 * returned, then the calling function will request the pointer 245 * for the array for storing a again (with GetValuesPtr), write 246 * the values of the nonzero elements into it, and call this 247 * MultiSolve method again with the same right-hand sides. (This 248 * can be done, for example, if the linear solver realized it 249 * does not have sufficient memory and needs to redo the 250 * factorization; e.g., for MA27.) 251 * 252 * The number of right-hand sides is given by nrhs, the values of 253 * the right-hand sides are given in rhs_vals (one full right-hand 254 * side stored immediately after the other), and solutions are 255 * to be returned in the same array. 256 * 257 * check_NegEVals will not be chosen true, if ProvidesInertia() 258 * returns false. 259 */ MultiSolve(bool new_matrix,const Index * ia,const Index * ja,Index nrhs,double * rhs_vals,bool check_NegEVals,Index numberOfNegEVals)260 ESymSolverStatus Ma86SolverInterface::MultiSolve(bool new_matrix, 261 const Index* ia, const Index* ja, Index nrhs, double* rhs_vals, 262 bool check_NegEVals, Index numberOfNegEVals) 263 { 264 struct ma86_info info; 265 266 if (new_matrix || pivtol_changed_) { 267 268 269 if (HaveIpData()) { 270 IpData().TimingStats().LinearSystemFactorization().Start(); 271 } 272 //ma86_factor(ndim_, ia, ja, val_, order_, &keep_, &control_, &info); 273 ma86_factor_solve(ndim_, ia, ja, val_, order_, &keep_, &control_, &info, 274 nrhs, ndim_, rhs_vals, NULL); 275 if (HaveIpData()) { 276 IpData().TimingStats().LinearSystemFactorization().End(); 277 } 278 if (info.flag<0) return SYMSOLVER_FATAL_ERROR; 279 if (info.flag==2 || info.flag==-3) return SYMSOLVER_SINGULAR; 280 if (check_NegEVals && info.num_neg!=numberOfNegEVals) 281 return SYMSOLVER_WRONG_INERTIA; 282 283 /*if (HaveIpData()) { 284 IpData().TimingStats().LinearSystemBackSolve().Start(); 285 } 286 ma86_solve(0, 1, ndim_, rhs_vals, order_, &keep_, &control_, &info); 287 if (HaveIpData()) { 288 IpData().TimingStats().LinearSystemBackSolve().End(); 289 }*/ 290 291 numneg_ = info.num_neg; 292 pivtol_changed_ = false; 293 } 294 else { 295 if (HaveIpData()) { 296 IpData().TimingStats().LinearSystemBackSolve().Start(); 297 } 298 ma86_solve(0, nrhs, ndim_, rhs_vals, order_, &keep_, &control_, &info, 299 NULL); 300 if (HaveIpData()) { 301 IpData().TimingStats().LinearSystemBackSolve().End(); 302 } 303 } 304 305 return SYMSOLVER_SUCCESS; 306 } 307 308 IncreaseQuality()309 bool Ma86SolverInterface::IncreaseQuality() 310 { 311 if (control_.u >= umax_) { 312 return false; 313 } 314 pivtol_changed_ = true; 315 316 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 317 "Increasing pivot tolerance for HSL_MA86 from %7.2e ", 318 control_.u); 319 control_.u = Min(umax_, pow(control_.u,0.75)); 320 Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, 321 "to %7.2e.\n", 322 control_.u); 323 return true; 324 } 325 326 } // namespace Ipopt 327 328 #endif /* COINHSL_HAS_MA86 or HAVE_LINEARSOLVERLOADER */ 329