1 // Copyright (C) 2012, The Science and Technology Facilities Council (STFC) 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: IpMa97SolverInterface.hpp 2001 2011-06-02 17:43:07Z 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 #ifndef __IPMA97SOLVERINTERFACE_HPP__ 14 #define __IPMA97SOLVERINTERFACE_HPP__ 15 16 #include "IpSparseSymLinearSolverInterface.hpp" 17 extern "C" 18 { 19 #include "hsl_ma97d.h" 20 } 21 22 namespace Ipopt 23 { 24 25 /** Base class for interfaces to symmetric indefinite linear solvers 26 * for sparse matrices. 27 * 28 * This defines the general interface to linear solvers for sparse 29 * symmetric indefinite matrices. The matrices can be provided 30 * either in "triplet format" (like for Harwell's MA27 solver), or 31 * in compressed sparse row (CSR) format for the lower triangular 32 * part of the symmetric matrix. 33 * 34 * The solver should be able to compute the interia of the matrix, 35 * or more specifically, the number of negative eigenvalues in the 36 * factorized matrix. 37 * 38 * This interface is used by the calling objective in the following 39 * way: 40 * 41 * 1. The InitializeImpl method is called at the very beginning 42 * (for every optimization run), which allows the linear solver 43 * object to retrieve options given in the OptionsList (such as 44 * pivot tolerances etc). At this point, some internal data can 45 * also be initialized. 46 * 47 * 2. The calling class calls MatrixFormat to find out which matrix 48 * representation the linear solver requires. The possible options 49 * are Triplet_Format, as well as CSR_Format_0_Offset and 50 * CSR_Format_1_Offset. The difference between the last two is 51 * that for CSR_Format_0_Offset the couning of the element position 52 * in the ia and ja arrays starts are 0 (C-style numbering), 53 * whereas for the other one it starts at 1 (Fortran-style 54 * numbering). 55 * 56 * 3. After this, the InitializeStructure method is called (once). 57 * Here, the structure of the matrix is provided. If the linear 58 * solver requires a symbolic preprocessing phase that can be done 59 * without knowledge of the matrix element values, it can be done 60 * here. 61 * 62 * 4. The calling class will request an array for storing the 63 * actual values for a matrix using the GetValuesArrayPtr method. 64 * This array must be at least as large as the number of nonzeros 65 * in the matrix (as given to this class by the InitializeStructure 66 * method call). After a call of this method, the calling class 67 * will fill this array with the actual values of the matrix. 68 * 69 * 5. Every time lateron, when actual solves of a linear system is 70 * requested, the calling class will call the MultiSolve to request 71 * the solve, possibly for mulitple right-hand sides. The flag 72 * new_matrix then indicates if the values of the matrix have 73 * changed and if a factorization is required, or if an old 74 * factorization can be used to do the solve. 75 * 76 * Note that the GetValuesArrayPtr method will be called before 77 * every call of MultiSolve with new_matrix=true, or before a 78 * renewed call of MultiSolve if the most previous return value was 79 * SYMSOLV_CALL_AGAIN. 80 * 81 * 6. The calling class might request with NumberOfNegEVals the 82 * number of the negative eigenvalues for the original matrix that 83 * were detected during the most recently performed factorization. 84 * 85 * 7. The calling class might ask the linear solver to increase the 86 * quality of the solution. For example, if the linear solver uses 87 * a pivot tolerance, a larger value should be used for the next 88 * solve (which might require a refactorization). 89 * 90 * 8. Finally, when the destructor is called, the internal storage, 91 * also in the linear solver, should be released. 92 * 93 * Note, if the matrix is given in triplet format, entries might be 94 * listed multiple times, in which case the corresponsing elements 95 * have to be added. 96 * 97 * A note for warm starts: If the option 98 * "warm_start_same_structure" is specified with "yes", the 99 * algorithm assumes that a problem with the same sparsity 100 * structure is solved for a repeated time. In that case, the 101 * linear solver might reuse information from the previous 102 * optimization. See Ma27TSolverInterface for an example. 103 */ 104 class Ma97SolverInterface: public SparseSymLinearSolverInterface 105 { 106 private: 107 enum order_opts { 108 ORDER_AUTO, 109 ORDER_BEST, 110 ORDER_AMD, 111 ORDER_METIS, 112 ORDER_MATCHED_AUTO, 113 ORDER_MATCHED_AMD, 114 ORDER_MATCHED_METIS 115 }; 116 enum scale_opts { 117 SWITCH_NEVER, 118 SWITCH_AT_START, 119 SWITCH_AT_START_REUSE, 120 SWITCH_ON_DEMAND, 121 SWITCH_ON_DEMAND_REUSE, 122 SWITCH_NDELAY, 123 SWITCH_NDELAY_REUSE, 124 SWITCH_OD_ND, 125 SWITCH_OD_ND_REUSE 126 }; 127 128 int ndim_; // Number of dimensions 129 double *val_; // Storage for variables 130 int numneg_; // Number of negative pivots in last factorization 131 int numdelay_; // Number of delayed pivots last time we scaled 132 void *akeep_; // Stores pointer to factors (only understood Fortran code!) 133 void *fkeep_; // Stores pointer to factors (only understood Fortran code!) 134 bool pivtol_changed_; // indicates if pivtol has been changed 135 bool rescale_; // Indicates if we shuold rescale next factorization 136 double *scaling_; // Store scaling for reuse if doing dynamic scaling 137 int fctidx_; // Current factorization number to dump to 138 139 /* Options */ 140 struct ma97_control control_; 141 double umax_; 142 int ordering_; 143 int scaling_type_; 144 enum scale_opts switch_[3]; 145 int scaling_val_[3]; 146 int current_level_; 147 bool dump_; 148 149 public: 150 Ma97SolverInterface()151 Ma97SolverInterface() : 152 val_(NULL), numdelay_(0), akeep_(NULL), fkeep_(NULL), pivtol_changed_(false), 153 rescale_(false), scaling_(NULL), fctidx_(0), scaling_type_(0), 154 dump_(false) 155 {} 156 ~Ma97SolverInterface(); 157 158 static void RegisterOptions(SmartPtr<RegisteredOptions> roptions); 159 160 bool InitializeImpl(const OptionsList& options, 161 const std::string& prefix); 162 163 /** @name Methods for requesting solution of the linear system. */ 164 //@{ 165 /** Method for initializing internal stuctures. Here, ndim gives 166 * the number of rows and columns of the matrix, nonzeros give 167 * the number of nonzero elements, and ia and ja give the 168 * positions of the nonzero elements, given in the matrix format 169 * determined by MatrixFormat. 170 */ 171 ESymSolverStatus InitializeStructure(Index dim, Index nonzeros, 172 const Index* ia, 173 const Index* ja); 174 175 /** Method returing an internal array into which the nonzero 176 * elements (in the same order as ja) will be stored by the 177 * calling routine before a call to MultiSolve with a 178 * new_matrix=true (or after a return of MultiSolve with 179 * SYMSOLV_CALL_AGAIN). The returned array must have space for at 180 * least nonzero elements. */ GetValuesArrayPtr()181 double* GetValuesArrayPtr() 182 { 183 return val_; 184 } 185 186 /** Solve operation for multiple right hand sides. Solves the 187 * linear system A * x = b with multiple right hand sides, where 188 * A is the symmtric indefinite matrix. Here, ia and ja give the 189 * positions of the values (in the required matrix data format). 190 * The actual values of the matrix will have been given to this 191 * object by copying them into the array provided by 192 * GetValuesArrayPtr. ia and ja are identical to the ones given 193 * to InitializeStructure. The flag new_matrix is set to true, 194 * if the values of the matrix has changed, and a refactorzation 195 * is required. 196 * 197 * The return code is SYMSOLV_SUCCESS if the factorization and 198 * solves were successful, SYMSOLV_SINGULAR if the linear system 199 * is singular, and SYMSOLV_WRONG_INERTIA if check_NegEVals is 200 * true and the number of negative eigenvalues in the matrix does 201 * not match numberOfNegEVals. If SYMSOLV_CALL_AGAIN is 202 * returned, then the calling function will request the pointer 203 * for the array for storing a again (with GetValuesPtr), write 204 * the values of the nonzero elements into it, and call this 205 * MultiSolve method again with the same right-hand sides. (This 206 * can be done, for example, if the linear solver realized it 207 * does not have sufficient memory and needs to redo the 208 * factorization; e.g., for MA27.) 209 * 210 * The number of right-hand sides is given by nrhs, the values of 211 * the right-hand sides are given in rhs_vals (one full right-hand 212 * side stored immediately after the other), and solutions are 213 * to be returned in the same array. 214 * 215 * check_NegEVals will not be chosen true, if ProvidesInertia() 216 * returns false. 217 */ 218 ESymSolverStatus MultiSolve(bool new_matrix, 219 const Index* ia, 220 const Index* ja, 221 Index nrhs, 222 double* rhs_vals, 223 bool check_NegEVals, 224 Index numberOfNegEVals); 225 226 /** Number of negative eigenvalues detected during last 227 * factorization. Returns the number of negative eigenvalues of 228 * the most recent factorized matrix. This must not be called if 229 * the linear solver does not compute this quantities (see 230 * ProvidesInertia). 231 */ NumberOfNegEVals() const232 Index NumberOfNegEVals() const 233 { 234 return numneg_; 235 } 236 //@} 237 238 //* @name Options of Linear solver */ 239 //@{ 240 /** Request to increase quality of solution for next solve. The 241 * calling class asks linear solver to increase quality of 242 * solution for the next solve (e.g. increase pivot tolerance). 243 * Returns false, if this is not possible (e.g. maximal pivot 244 * tolerance already used.) 245 */ 246 bool IncreaseQuality(); 247 248 /** Query whether inertia is computed by linear solver. Returns 249 * true, if linear solver provides inertia. 250 */ ProvidesInertia() const251 bool ProvidesInertia() const 252 { 253 return true; 254 } 255 256 /** Query of requested matrix type that the linear solver 257 * understands. 258 */ MatrixFormat() const259 EMatrixFormat MatrixFormat() const 260 { 261 return CSR_Format_1_Offset; 262 } 263 //@} 264 265 /** @name Methods related to the detection of linearly dependent 266 * rows in a matrix */ 267 //@{ 268 /** Query whether the indices of linearly dependent rows/columns 269 * can be determined by this linear solver. */ ProvidesDegeneracyDetection() const270 bool ProvidesDegeneracyDetection() const 271 { 272 return false; 273 } 274 /** This method determines the list of row indices of the linearly 275 * dependent rows. */ DetermineDependentRows(const Index * ia,const Index * ja,std::list<Index> & c_deps)276 ESymSolverStatus DetermineDependentRows(const Index* ia, 277 const Index* ja, 278 std::list<Index>& c_deps) 279 { 280 return SYMSOLVER_FATAL_ERROR; 281 } 282 283 /** converts a scalign optoin name to its ma97 option number */ 284 static int ScaleNameToNum(const std::string& name); 285 }; 286 287 } // namespace Ipopt 288 289 #endif 290