1 | // Copyright (C) 2005 International Business Machines and others. |
---|
2 | // All Rights Reserved. |
---|
3 | // This code is published under the Common Public License. |
---|
4 | // |
---|
5 | // $Id: IpMa57TSolverInterface.cpp 568 2005-11-01 23:28:48Z andreasw $ |
---|
6 | // |
---|
7 | // Authors: Michael Hagemann Univ of Basel 2005-10-28 |
---|
8 | // original version (based on MA27TSolverInterface.cpp) |
---|
9 | |
---|
10 | #include "IpMa57TSolverInterface.hpp" |
---|
11 | |
---|
12 | #ifdef HAVE_CMATH |
---|
13 | # include <cmath> |
---|
14 | #else |
---|
15 | # ifdef HAVE_MATH_H |
---|
16 | # include <math.h> |
---|
17 | # else |
---|
18 | # error "don't have header file for math" |
---|
19 | # endif |
---|
20 | #endif |
---|
21 | |
---|
22 | #include <iostream> |
---|
23 | |
---|
24 | |
---|
25 | /** Prototypes for MA57's Fortran subroutines */ |
---|
26 | extern "C" |
---|
27 | { |
---|
28 | /* |
---|
29 | * MA57ID -- Initialize solver. |
---|
30 | */ |
---|
31 | extern void F77_FUNC (ma57id, MA57ID) ( |
---|
32 | double *cntl, |
---|
33 | ipfint *icntl); |
---|
34 | |
---|
35 | /* |
---|
36 | * MA57AD -- Symbolic Factorization. |
---|
37 | */ |
---|
38 | extern void F77_FUNC (ma57ad, MA57AD) ( |
---|
39 | ipfint *n, /* Order of matrix. */ |
---|
40 | ipfint *ne, /* Number of entries. */ |
---|
41 | |
---|
42 | const ipfint *irn, /* Matrix nonzero row structure */ |
---|
43 | const ipfint *jcn, /* Matrix nonzero column structure */ |
---|
44 | |
---|
45 | ipfint *lkeep, /* Workspace for the pivot order of lenght 3*n */ |
---|
46 | ipfint *keep, /* Workspace for the pivot order of lenght 3*n */ |
---|
47 | /* Automatically iflag = 0; ikeep pivot order iflag = 1 */ |
---|
48 | ipfint *iwork, /* Integer work space. */ |
---|
49 | ipfint *icntl, /* Integer Control parameter of length 30*/ |
---|
50 | ipfint *info, /* Statistical Information; Integer array of length 20 */ |
---|
51 | double *rinfo); /* Double Control parameter of length 5 */ |
---|
52 | |
---|
53 | /* |
---|
54 | * MA57BD -- Numerical Factorization. |
---|
55 | */ |
---|
56 | extern void F77_FUNC (ma57bd, MA57BD) ( |
---|
57 | ipfint *n, /* Order of matrix. */ |
---|
58 | ipfint *ne, /* Number of entries. */ |
---|
59 | |
---|
60 | double *a, /* Numerical values. */ |
---|
61 | double *fact, /* Entries of factors. */ |
---|
62 | ipfint *lfact, /* Length of array `fact'. */ |
---|
63 | ipfint *ifact, /* Indexing info for factors. */ |
---|
64 | ipfint *lifact, /* Length of array `ifact'. */ |
---|
65 | |
---|
66 | ipfint *lkeep, /* Length of array `keep'. */ |
---|
67 | ipfint *keep, /* Integer array. */ |
---|
68 | |
---|
69 | ipfint *iwork, /* Workspace of length `n'. */ |
---|
70 | |
---|
71 | ipfint *icntl, /* Integer Control parameter of length 20. */ |
---|
72 | double *cntl, /* Double Control parameter of length 5. */ |
---|
73 | ipfint *info, /* Statistical Information; Integer array of length 40. */ |
---|
74 | double *rinfo); /* Statistical Information; Real array of length 20. */ |
---|
75 | |
---|
76 | /* |
---|
77 | * MA57CD -- Solution. |
---|
78 | */ |
---|
79 | extern void F77_FUNC (ma57cd, MA57CD) ( |
---|
80 | ipfint *job, /* Solution job. Solve for... |
---|
81 | JOB <= 1: A |
---|
82 | JOB == 2: PLP^t |
---|
83 | JOB == 3: PDP^t |
---|
84 | JOB >= 4: PL^t P^t */ |
---|
85 | |
---|
86 | ipfint *n, /* Order of matrix. */ |
---|
87 | |
---|
88 | double *fact, /* Entries of factors. */ |
---|
89 | ipfint *lfact, /* Length of array `fact'. */ |
---|
90 | ipfint *ifact, /* Indexing info for factors. */ |
---|
91 | ipfint *lifact, /* Length of array `ifact'. */ |
---|
92 | |
---|
93 | ipfint *nrhs, /* Number of right hand sides. */ |
---|
94 | double *rhs, /* Numerical Values. */ |
---|
95 | ipfint *lrhs, /* Leading dimensions of `rhs'. */ |
---|
96 | |
---|
97 | double *work, /* Real workspace. */ |
---|
98 | ipfint *lwork, /* Length of `work', >= N*NRHS. */ |
---|
99 | ipfint *iwork, /* Integer array of length `n'. */ |
---|
100 | |
---|
101 | ipfint *icntl, /* Integer Control parameter array of length 20. */ |
---|
102 | ipfint *info); /* Statistical Information; Integer array of length 40. */ |
---|
103 | |
---|
104 | /* |
---|
105 | * MC57ED -- Copy arrays. |
---|
106 | */ |
---|
107 | extern void F77_FUNC (ma57ed, MA57ED) ( |
---|
108 | ipfint *n, |
---|
109 | ipfint *ic, /* 0: copy real array. >=1: copy integer array. */ |
---|
110 | ipfint *keep, |
---|
111 | |
---|
112 | double *fact, |
---|
113 | ipfint *lfact, |
---|
114 | double *newfac, |
---|
115 | ipfint *lnew, |
---|
116 | |
---|
117 | ipfint *ifact, |
---|
118 | ipfint *lifact, |
---|
119 | ipfint *newifc, |
---|
120 | ipfint *linew, |
---|
121 | |
---|
122 | ipfint *info); |
---|
123 | } |
---|
124 | |
---|
125 | namespace Ipopt |
---|
126 | { |
---|
127 | #ifdef IP_DEBUG |
---|
128 | static const Index dbg_verbosity = 0; |
---|
129 | #endif |
---|
130 | |
---|
131 | std::string ma57_err_msg[] = { |
---|
132 | "Operation successful.\n", |
---|
133 | |
---|
134 | "Value of N is out of range on a call to MA57A/AD, MA57B/BD, MA57C/CD, or\n" |
---|
135 | "MA57D/DD. Value given is held in INFO(2).\n", |
---|
136 | |
---|
137 | "Value of NE is out of range on a call to MA57A/AD, MA57B/BD, or\n" |
---|
138 | "MA57D/DD. Value given is held in INFO(2).\n", |
---|
139 | |
---|
140 | "Failure due to insufficient REAL space on a call to MA57B/BD. INFO(17)\n" |
---|
141 | "is set to a value that may suffice. INFO(2) is set to value of\n" |
---|
142 | "LFACT. The user can allocate a larger array and copy the contents of\n" |
---|
143 | "FACT into it using MA57E/ED, and recall MA57B/BD.\n", |
---|
144 | |
---|
145 | "Failure due to insufficient INTEGER space on a call to\n" |
---|
146 | "MA57B/BD. INFO(18) is set to a value that may suffice. INFO(2) is set to\n" |
---|
147 | "value of LIFACT. The user can allocate a larger array and copy the\n" |
---|
148 | "contents of IFACT into it using MA57E/ED, and recall MA57B/BD.\n", |
---|
149 | |
---|
150 | "A pivot with magnitude less than or equal to CNTL(2) was found at pivot\n" |
---|
151 | "step INFO(2) when calling MA57B/BD with ICNTL(7) = 2 or 3, or the\n" |
---|
152 | "correction obtained when using matrix modification does not give a pivot\n" |
---|
153 | "greater than CNTL(2) when ICNTL(7) = 4.\n", |
---|
154 | |
---|
155 | "A change of sign of pivots has been detected when ICNTL(7) = 2. INFO(2)\n" |
---|
156 | "is set to the pivot step at which the change was detected on a call to\n" |
---|
157 | "MA57B/BD.\n", |
---|
158 | |
---|
159 | "Either LNEW < LFACT or LINEW < LIFACT on a call to MA57E/ED. INFO(2) is\n" |
---|
160 | "set to LNEW or LINEW as appropriate.\n", |
---|
161 | |
---|
162 | "Iterative refinement fails to converge in specified number of iterations\n" |
---|
163 | "on a call to MA57D/DD.\n", |
---|
164 | |
---|
165 | "Error in permutation array when ICNTL(6)=1 on a call to\n" |
---|
166 | "MA57A/AD. INFO(2) holds first component at which error was detected.\n", |
---|
167 | |
---|
168 | "Value of ICNTL(7) out of range on a call to MA57B/BD. Value given held\n" |
---|
169 | "in INFO(2).\n", |
---|
170 | |
---|
171 | "LRHS < N on a call to MA57C/CD. INFO(2) holds value of LRHS.\n", |
---|
172 | |
---|
173 | "Invalid value for JOB on a call to MA57D/DD. Value given held in\n" |
---|
174 | "INFO(2).\n", |
---|
175 | |
---|
176 | "Invalid value of ICNTL(9) on a call to MA57D/DD. Value given held in\n" |
---|
177 | "INFO(2).\n", |
---|
178 | |
---|
179 | "Failure of MC71A/AD on a call to MA57D/DD with ICNTL(10)> 0.\n", |
---|
180 | |
---|
181 | "LKEEP less than 5*N+NE+MAX(N,NE) +42 on a call to MA57A/AD or\n" |
---|
182 | "MA57B/BD. INFO(2) holds value of LKEEP.\n", |
---|
183 | |
---|
184 | "NRHS less than 1 on call to MA57C/CD. INFO(2) holds value of NRHS.\n", |
---|
185 | |
---|
186 | "LWORK too small on entry to MA57C/CD. INFO(2) holds minimum length\n" |
---|
187 | "required. A positive value of INFO(1) is associated with a warning\n" |
---|
188 | "message that will be output on unit ICNTL(2).\n" |
---|
189 | }; |
---|
190 | |
---|
191 | std::string ma57_wrn_msg[] = { |
---|
192 | "Operation successful.\n", |
---|
193 | |
---|
194 | "Index (in IRN or JCN) out of range on call to MA57A/AD or\n" |
---|
195 | "MA57D/DD. Action taken by subroutine is to ignore any such entries and\n" |
---|
196 | "continue. INFO(3) is set to the number of faulty entries. Details of the\n" |
---|
197 | "first ten are printed on unit ICNTL(2).\n", |
---|
198 | |
---|
199 | "Duplicate indices on call to MA57A/AD or MA57D/DD. Action taken by\n" |
---|
200 | "subroutine is to keep the duplicates and then to sum corresponding reals\n" |
---|
201 | "when MA57B/BD is called. INFO(4) is set to the number of faulty\n" |
---|
202 | "entries. Details of the first ten are printed on unit ICNTL(2).\n", |
---|
203 | |
---|
204 | "Both out-of-range indices and duplicates exist.\n", |
---|
205 | |
---|
206 | "Matrix is rank deficient on exit from MA57B/BD. In this case, a\n" |
---|
207 | "decomposition will still have been produced that will enable the\n" |
---|
208 | "subsequent solution of consistent equations. INFO(25) will be set to the\n" |
---|
209 | "rank of the factorized matrix.\n", |
---|
210 | |
---|
211 | "Pivots have different signs when factorizing a supposedly definite\n" |
---|
212 | "matrix (ICNTL(7) = 3) on call to MA57B/BD. INFO(26) is set to the number\n" |
---|
213 | "of sign changes.\n", |
---|
214 | |
---|
215 | "-", |
---|
216 | "-", |
---|
217 | |
---|
218 | "During error analysis the infinity norm of the computed solution was\n" |
---|
219 | "found to be zero.\n", |
---|
220 | |
---|
221 | "Insufficient real space to complete factorization when MA57B/BD called\n" |
---|
222 | "with ICNTL(8) != 0. User can copy real values to a longer array using\n" |
---|
223 | "MA57E/ED and recall MA57B/BD using this longer array to continue the\n" |
---|
224 | "factorization.\n", |
---|
225 | |
---|
226 | "Insufficient integer space to complete factorization when MA57B/BD\n" |
---|
227 | "called with ICNTL(8) != 0. User can copy integer values to a longer\n" |
---|
228 | "array using MA57E/ED and recall MA57B/BD using this longer array to\n" |
---|
229 | "continue the factorization.\n" |
---|
230 | }; |
---|
231 | |
---|
232 | Ma57TSolverInterface::Ma57TSolverInterface() |
---|
233 | : |
---|
234 | dim_(0), |
---|
235 | nonzeros_(0), |
---|
236 | initialized_(false), |
---|
237 | pivtol_changed_(false), |
---|
238 | refactorize_(false), |
---|
239 | |
---|
240 | wd_keep_(NULL), |
---|
241 | wd_iwork_(NULL), |
---|
242 | wd_fact_(NULL), |
---|
243 | wd_ifact_(NULL), |
---|
244 | a_(NULL) |
---|
245 | { |
---|
246 | DBG_START_METH("Ma57TSolverInterface::Ma57TSolverInterface()", |
---|
247 | dbg_verbosity); |
---|
248 | } |
---|
249 | |
---|
250 | Ma57TSolverInterface::~Ma57TSolverInterface() |
---|
251 | { |
---|
252 | DBG_START_METH("Ma57TSolverInterface::~Ma57TSolverInterface()", |
---|
253 | dbg_verbosity); |
---|
254 | delete [] a_; |
---|
255 | |
---|
256 | delete [] wd_fact_; |
---|
257 | delete [] wd_ifact_; |
---|
258 | |
---|
259 | delete [] wd_iwork_; |
---|
260 | delete [] wd_keep_; |
---|
261 | } |
---|
262 | |
---|
263 | void Ma57TSolverInterface::RegisterOptions(SmartPtr<RegisteredOptions> roptions) |
---|
264 | { |
---|
265 | roptions->AddBoundedNumberOption( |
---|
266 | "ma57_pivtol", |
---|
267 | "Pivot tolerance for the linear solver MA57.", |
---|
268 | 0.0, true, 1.0, true, 1e-8, |
---|
269 | "A smaller number pivots for sparsity, " |
---|
270 | "a larger number pivots for stability."); |
---|
271 | roptions->AddBoundedNumberOption( |
---|
272 | "ma57_pivtolmax", |
---|
273 | "Maximum pivot tolerance for the linear solver MA57.", |
---|
274 | 0.0, true, 1.0, true, 1e-4, |
---|
275 | "Ipopt may increase pivtol as high as pivtolmax " |
---|
276 | "to get a more accurate solution to the linear system."); |
---|
277 | roptions->AddLowerBoundedNumberOption( |
---|
278 | "ma57_pre_alloc", |
---|
279 | "Safety factor for work space memory allocation for the linear solver MA57.", |
---|
280 | 1., false, 3., |
---|
281 | "If 1 is chosen, the suggested amount of work space is used. However, " |
---|
282 | "choosen a larger number might avoid reallocation if the suggest values " |
---|
283 | "do not suffice."); |
---|
284 | } |
---|
285 | |
---|
286 | bool Ma57TSolverInterface::InitializeImpl(const OptionsList& options, |
---|
287 | const std::string& prefix) |
---|
288 | { |
---|
289 | // Obtain the options settings |
---|
290 | options.GetNumericValue("ma57_pivtol", pivtol_, prefix); |
---|
291 | if(options.GetNumericValue("ma57_pivtolmax", pivtolmax_, prefix)) { |
---|
292 | ASSERT_EXCEPTION(pivtolmax_>=pivtol_, OPTION_INVALID, |
---|
293 | "Option \"pivtolmax\": This value must be between " |
---|
294 | "pivtol and 1."); |
---|
295 | } |
---|
296 | else { |
---|
297 | pivtolmax_ = Max(pivtolmax_, pivtol_); |
---|
298 | } |
---|
299 | |
---|
300 | options.GetNumericValue("ma57_pre_alloc", ma57_pre_alloc_, prefix); |
---|
301 | |
---|
302 | // The following option is registered by OrigIpoptNLP |
---|
303 | options.GetBoolValue("warm_start_same_structure", |
---|
304 | warm_start_same_structure_, prefix); |
---|
305 | DBG_ASSERT(!warm_start_same_structure_ && "warm_start_same_structure not yet implemented"); |
---|
306 | |
---|
307 | |
---|
308 | /* Initialize. */ |
---|
309 | F77_FUNC (ma57id, MA57ID) (wd_cntl_, wd_icntl_); |
---|
310 | |
---|
311 | /* Custom settings for MA57. */ |
---|
312 | wd_icntl_[1-1] = 0; /* Error stream */ |
---|
313 | wd_icntl_[2-1] = 0; /* Warning stream. */ |
---|
314 | |
---|
315 | wd_icntl_[4-1] = 1; /* Print statistics. NOT Used. */ |
---|
316 | wd_icntl_[5-1] = 0; /* Print error. */ |
---|
317 | |
---|
318 | // wd_icntl[6-1] = 0; /* Pivoting order. */ |
---|
319 | |
---|
320 | wd_cntl_[1-1] = pivtol_; /* Pivot threshold. */ |
---|
321 | wd_icntl_[7-1] = 1; /* Pivoting strategy. */ |
---|
322 | |
---|
323 | // wd_icntl[8-1] = 0; /* Retry factorization. */ |
---|
324 | |
---|
325 | if (!warm_start_same_structure_) { |
---|
326 | dim_=0; |
---|
327 | nonzeros_=0; |
---|
328 | delete [] a_; |
---|
329 | delete [] wd_fact_; |
---|
330 | delete [] wd_ifact_; |
---|
331 | delete [] wd_iwork_; |
---|
332 | delete [] wd_keep_; |
---|
333 | } |
---|
334 | else { |
---|
335 | ASSERT_EXCEPTION(dim_>0 && nonzeros_>0, INVALID_WARMSTART, |
---|
336 | "Ma57TSolverInterface called with warm_start_same_structure, " |
---|
337 | "but the problem is solved for the first time."); |
---|
338 | } |
---|
339 | |
---|
340 | return true; |
---|
341 | } |
---|
342 | |
---|
343 | ESymSolverStatus |
---|
344 | Ma57TSolverInterface::MultiSolve(bool new_matrix, |
---|
345 | const Index* airn, |
---|
346 | const Index* ajcn, |
---|
347 | Index nrhs, |
---|
348 | double* rhs_vals, |
---|
349 | bool check_NegEVals, |
---|
350 | Index numberOfNegEVals) |
---|
351 | { |
---|
352 | DBG_START_METH("Ma57TSolverInterface::MultiSolve",dbg_verbosity); |
---|
353 | |
---|
354 | // DBG_ASSERT(!check_NegEVals || ProvidesInertia()); |
---|
355 | // DBG_ASSERT(initialized_); |
---|
356 | // DBG_ASSERT(la_!=0); |
---|
357 | |
---|
358 | if (pivtol_changed_) { |
---|
359 | DBG_PRINT((1,"Pivot tolerance has changed.\n")); |
---|
360 | pivtol_changed_ = false; |
---|
361 | // If the pivot tolerance has been changed but the matrix is not |
---|
362 | // new, we have to request the values for the matrix again to do |
---|
363 | // the factorization again. |
---|
364 | if (!new_matrix) { |
---|
365 | DBG_PRINT((1,"Ask caller to call again.\n")); |
---|
366 | refactorize_ = true; |
---|
367 | return SYMSOLVER_CALL_AGAIN; |
---|
368 | } |
---|
369 | } |
---|
370 | |
---|
371 | // check if a factorization has to be done |
---|
372 | DBG_PRINT((1, "new_matrix = %d\n", new_matrix)); |
---|
373 | if (new_matrix || refactorize_) { |
---|
374 | // perform the factorization |
---|
375 | ESymSolverStatus retval; |
---|
376 | retval = Factorization(airn, ajcn, check_NegEVals, numberOfNegEVals); |
---|
377 | if (retval!=SYMSOLVER_SUCCESS) { |
---|
378 | DBG_PRINT((1, "FACTORIZATION FAILED!\n")); |
---|
379 | return retval; // Matrix singular or error occurred |
---|
380 | } |
---|
381 | refactorize_ = false; |
---|
382 | } |
---|
383 | |
---|
384 | // do the backsolve |
---|
385 | return Backsolve(nrhs, rhs_vals); |
---|
386 | } |
---|
387 | |
---|
388 | double* Ma57TSolverInterface::GetValuesArrayPtr() |
---|
389 | { |
---|
390 | DBG_START_METH("Ma57TSolverInterface::GetValuesArrayPtr",dbg_verbosity); |
---|
391 | DBG_ASSERT(initialized_); |
---|
392 | |
---|
393 | return a_; |
---|
394 | } |
---|
395 | |
---|
396 | /** Initialize the local copy of the positions of the nonzero |
---|
397 | elements */ |
---|
398 | ESymSolverStatus |
---|
399 | Ma57TSolverInterface::InitializeStructure( |
---|
400 | Index dim, |
---|
401 | Index nonzeros, |
---|
402 | const Index* airn, |
---|
403 | const Index* ajcn) |
---|
404 | { |
---|
405 | DBG_START_METH("Ma57TSolverInterface::InitializeStructure",dbg_verbosity); |
---|
406 | |
---|
407 | ESymSolverStatus retval = SYMSOLVER_SUCCESS; |
---|
408 | if (!warm_start_same_structure_) { |
---|
409 | dim_ = dim; |
---|
410 | nonzeros_ = nonzeros; |
---|
411 | // for MA57, a_ only has to be as long as the number of nonzero |
---|
412 | // elements |
---|
413 | delete [] a_; |
---|
414 | a_ = new double [nonzeros_]; |
---|
415 | |
---|
416 | // Do the symbolic facotrization |
---|
417 | retval = SymbolicFactorization(airn, ajcn); |
---|
418 | if (retval != SYMSOLVER_SUCCESS ) { |
---|
419 | return retval; |
---|
420 | } |
---|
421 | } |
---|
422 | else { |
---|
423 | ASSERT_EXCEPTION(dim_==dim && nonzeros_==nonzeros, INVALID_WARMSTART, |
---|
424 | "Ma57TSolverInterface called with warm_start_same_structure, " |
---|
425 | "but the problem size has changed."); |
---|
426 | } |
---|
427 | |
---|
428 | initialized_ = true; |
---|
429 | |
---|
430 | return retval; |
---|
431 | } |
---|
432 | |
---|
433 | ESymSolverStatus |
---|
434 | Ma57TSolverInterface::SymbolicFactorization(const Index* airn, |
---|
435 | const Index* ajcn) |
---|
436 | { |
---|
437 | DBG_START_METH("Ma57TSolverInterface::SymbolicFactorization",dbg_verbosity); |
---|
438 | IpData().TimingStats().LinearSystemSymbolicFactorization().Start(); |
---|
439 | |
---|
440 | ipfint n = dim_; |
---|
441 | ipfint ne = nonzeros_; |
---|
442 | |
---|
443 | wd_lkeep_ = 5*n + ne + Max (n,ne) + 42; |
---|
444 | |
---|
445 | wd_iwork_ = new ipfint[5*n]; |
---|
446 | wd_keep_ = new ipfint[wd_lkeep_]; |
---|
447 | |
---|
448 | F77_FUNC (ma57ad, MA57AD) |
---|
449 | (&n, &ne, airn, ajcn, &wd_lkeep_, wd_keep_, wd_iwork_, |
---|
450 | wd_icntl_, wd_info_, wd_rinfo_); |
---|
451 | |
---|
452 | if (wd_info_[0] < 0) { |
---|
453 | Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA, |
---|
454 | "*** Error from MA57AD *** INFO(0) = %d\n", wd_info_[0]); |
---|
455 | } |
---|
456 | |
---|
457 | wd_lfact_ = (ipfint)((Number)wd_info_[8] * ma57_pre_alloc_); |
---|
458 | wd_lifact_ = (ipfint)((Number)wd_info_[9] * ma57_pre_alloc_); |
---|
459 | |
---|
460 | // XXX MH: Why is this necessary? Is `::Factorization' called more |
---|
461 | // than once per object lifetime? Where should allocation take |
---|
462 | // place, then? |
---|
463 | |
---|
464 | // AW: I moved this here now - my understanding is that wd_info[8] |
---|
465 | // and wd_info[9] are computed here, so we can just allocate the |
---|
466 | // amount of memory here. I don't think there is any need to |
---|
467 | // reallocate it later for every factorization |
---|
468 | delete [] wd_fact_; |
---|
469 | delete [] wd_ifact_; |
---|
470 | |
---|
471 | wd_fact_ = new double[wd_lfact_]; |
---|
472 | wd_ifact_ = new int[wd_lifact_]; |
---|
473 | |
---|
474 | Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, |
---|
475 | "Suggested lfact (*%e): %d\n", ma57_pre_alloc_, wd_lfact_); |
---|
476 | Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, |
---|
477 | "Suggested lifact (*%e): %d\n", ma57_pre_alloc_, wd_lifact_); |
---|
478 | |
---|
479 | IpData().TimingStats().LinearSystemSymbolicFactorization().End(); |
---|
480 | return SYMSOLVER_SUCCESS; |
---|
481 | } |
---|
482 | |
---|
483 | ESymSolverStatus |
---|
484 | Ma57TSolverInterface::Factorization(const Index* airn, |
---|
485 | const Index* ajcn, |
---|
486 | bool check_NegEVals, |
---|
487 | Index numberOfNegEVals) |
---|
488 | { |
---|
489 | DBG_START_METH("Ma57TSolverInterface::Factorization",dbg_verbosity); |
---|
490 | IpData().TimingStats().LinearSystemFactorization().Start(); |
---|
491 | |
---|
492 | int fact_error = 1; |
---|
493 | |
---|
494 | ipfint n = dim_; |
---|
495 | ipfint ne = nonzeros_; |
---|
496 | |
---|
497 | while (fact_error > 0) { |
---|
498 | F77_FUNC (ma57bd, MA57BD) |
---|
499 | (&n, &ne, a_, wd_fact_, &wd_lfact_, wd_ifact_, &wd_lifact_, |
---|
500 | &wd_lkeep_, wd_keep_, wd_iwork_, |
---|
501 | wd_icntl_, wd_cntl_, wd_info_, wd_rinfo_); |
---|
502 | |
---|
503 | negevals_ = wd_info_[24-1]; // Number of negative eigenvalues |
---|
504 | |
---|
505 | if (wd_info_[0] == 0) { |
---|
506 | fact_error = 0; |
---|
507 | } |
---|
508 | else if (wd_info_[0] == -3) { |
---|
509 | /* Failure due to insufficient REAL space on a call to MA57B/BD. |
---|
510 | * INFO(17) is set to a value that may suffice. INFO(2) is set |
---|
511 | * to value of LFACT. The user can allocate a larger array and |
---|
512 | * copy the contents of FACT into it using MA57E/ED, and recall |
---|
513 | * MA57B/BD. |
---|
514 | */ |
---|
515 | double *temp; |
---|
516 | ipfint ic = 0; |
---|
517 | |
---|
518 | wd_lfact_ = wd_info_[16]; |
---|
519 | temp = new double[wd_lfact_]; |
---|
520 | |
---|
521 | Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, |
---|
522 | "Reallocating lfact (%d)\n", wd_lfact_); |
---|
523 | |
---|
524 | F77_FUNC (ma57ed, MA57ED) |
---|
525 | (&n, &ic, wd_keep_, |
---|
526 | wd_fact_, &wd_info_[1], temp, &wd_lfact_, |
---|
527 | wd_ifact_, &wd_info_[1], NULL, &wd_lfact_, |
---|
528 | wd_info_); |
---|
529 | |
---|
530 | delete [] wd_fact_; |
---|
531 | wd_fact_ = temp; |
---|
532 | } |
---|
533 | else if (wd_info_[0] == -4) { |
---|
534 | /* Failure due to insufficient INTEGER space on a call to |
---|
535 | * MA57B/BD. INFO(18) is set to a value that may suffice. |
---|
536 | * INFO(2) is set to value of LIFACT. The user can allocate a |
---|
537 | * larger array and copy the contents of IFACT into it using |
---|
538 | * MA57E/ED, and recall MA57B/BD. |
---|
539 | */ |
---|
540 | |
---|
541 | ipfint *temp; |
---|
542 | ipfint ic = 1; |
---|
543 | |
---|
544 | wd_lifact_ = wd_info_[17]; |
---|
545 | temp = new ipfint[wd_lifact_]; |
---|
546 | |
---|
547 | Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, |
---|
548 | "Reallocating lifact (%d)\n", wd_lifact_); |
---|
549 | |
---|
550 | F77_FUNC (ma57ed, MA57ED) |
---|
551 | (&n, &ic, wd_keep_, |
---|
552 | wd_fact_, &wd_info_[1], NULL, &wd_lifact_, |
---|
553 | wd_ifact_, &wd_info_[1], temp, &wd_lifact_, |
---|
554 | wd_info_); |
---|
555 | |
---|
556 | delete [] wd_ifact_; |
---|
557 | wd_ifact_ = temp; |
---|
558 | } |
---|
559 | else if (wd_info_[0] < 0) { |
---|
560 | Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA, |
---|
561 | "Error in MA57BD: %d\n", wd_info_[0]); |
---|
562 | Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, |
---|
563 | "MA57 Error message: %s\n", |
---|
564 | ma57_err_msg[-wd_info_[1-1]].c_str()); |
---|
565 | return SYMSOLVER_FATAL_ERROR; |
---|
566 | } |
---|
567 | // Check if the system is singular. |
---|
568 | else if (wd_info_[0] == 4) { |
---|
569 | IpData().TimingStats().LinearSystemFactorization().End(); |
---|
570 | Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, |
---|
571 | "System singular, rank = %d\n", wd_info_[25-1]); |
---|
572 | return SYMSOLVER_SINGULAR; |
---|
573 | } |
---|
574 | else if (wd_info_[0] > 0) { |
---|
575 | Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA, |
---|
576 | "Warning in MA57BD: %d\n", wd_info_[0]); |
---|
577 | Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, |
---|
578 | "MA57 Warning message: %s\n", |
---|
579 | ma57_wrn_msg[wd_info_[1-1]].c_str()); |
---|
580 | // For now, abort the process so that we don't miss any problems |
---|
581 | return SYMSOLVER_FATAL_ERROR; |
---|
582 | } |
---|
583 | } |
---|
584 | |
---|
585 | // double peak_mem = 1.0e-3 * (wd_lfact*8.0 + wd_lifact*4.0 + wd_lkeep*4.0); |
---|
586 | |
---|
587 | |
---|
588 | // Check whether the number of negative eigenvalues matches the |
---|
589 | // requested count. |
---|
590 | IpData().TimingStats().LinearSystemFactorization().End(); |
---|
591 | if (check_NegEVals && (numberOfNegEVals!=negevals_)) { |
---|
592 | Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, |
---|
593 | "In Ma57TSolverInterface::Factorization: " |
---|
594 | "negevals_ = %d, but numberOfNegEVals = %d\n", |
---|
595 | negevals_, numberOfNegEVals); |
---|
596 | return SYMSOLVER_WRONG_INERTIA; |
---|
597 | } |
---|
598 | |
---|
599 | return SYMSOLVER_SUCCESS; |
---|
600 | } |
---|
601 | |
---|
602 | ESymSolverStatus Ma57TSolverInterface::Backsolve( |
---|
603 | Index nrhs, |
---|
604 | double *rhs_vals) |
---|
605 | { |
---|
606 | DBG_START_METH("Ma27TSolverInterface::Backsolve",dbg_verbosity); |
---|
607 | IpData().TimingStats().LinearSystemBackSolve().Start(); |
---|
608 | |
---|
609 | ipfint n = dim_; |
---|
610 | ipfint job = 1; |
---|
611 | |
---|
612 | ipfint nrhs_X = nrhs; |
---|
613 | ipfint lrhs = n; |
---|
614 | |
---|
615 | ipfint lwork; |
---|
616 | double* work; |
---|
617 | |
---|
618 | lwork = n * nrhs; |
---|
619 | work = new double[lwork]; |
---|
620 | |
---|
621 | // For each right hand side, call MA57CD |
---|
622 | // XXX MH: MA57 can do several RHSs; just do one solve... |
---|
623 | // AW: Ok is the following correct? |
---|
624 | if (DBG_VERBOSITY()>=2) { |
---|
625 | for(Index irhs=0; irhs<nrhs; irhs++) { |
---|
626 | for (Index i=0; i<dim_; i++) { |
---|
627 | DBG_PRINT((2, "rhs[%2d,%5d] = %23.15e\n", irhs, i, rhs_vals[irhs*dim_+i])); |
---|
628 | } |
---|
629 | } |
---|
630 | } |
---|
631 | |
---|
632 | F77_FUNC (ma57cd, MA57CD) |
---|
633 | (&job, &n, wd_fact_, &wd_lfact_, wd_ifact_, &wd_lifact_, |
---|
634 | &nrhs_X, rhs_vals, &lrhs, |
---|
635 | work, &lwork, wd_iwork_, |
---|
636 | wd_icntl_, wd_info_); |
---|
637 | |
---|
638 | if (wd_info_[0] != 0) |
---|
639 | Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA, |
---|
640 | "Error in MA57CD: %d.\n", wd_info_[0]); |
---|
641 | |
---|
642 | if (DBG_VERBOSITY()>=2) { |
---|
643 | for(Index irhs=0; irhs<nrhs; irhs++) { |
---|
644 | for (Index i=0; i<dim_; i++) { |
---|
645 | DBG_PRINT((2, "sol[%2d,%5d] = %23.15e\n", irhs, i, rhs_vals[irhs*dim_+i])); |
---|
646 | } |
---|
647 | } |
---|
648 | } |
---|
649 | |
---|
650 | delete [] work; |
---|
651 | |
---|
652 | IpData().TimingStats().LinearSystemBackSolve().End(); |
---|
653 | return SYMSOLVER_SUCCESS; |
---|
654 | } |
---|
655 | |
---|
656 | Index Ma57TSolverInterface::NumberOfNegEVals() const |
---|
657 | { |
---|
658 | DBG_START_METH("Ma57TSolverInterface::NumberOfNegEVals",dbg_verbosity); |
---|
659 | DBG_ASSERT(ProvidesInertia()); |
---|
660 | DBG_ASSERT(initialized_); |
---|
661 | return negevals_; |
---|
662 | } |
---|
663 | |
---|
664 | bool Ma57TSolverInterface::IncreaseQuality() |
---|
665 | { |
---|
666 | DBG_START_METH("Ma57TSolverInterface::IncreaseQuality",dbg_verbosity); |
---|
667 | if (pivtol_ == pivtolmax_) { |
---|
668 | return false; |
---|
669 | } |
---|
670 | pivtol_changed_ = true; |
---|
671 | |
---|
672 | Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, |
---|
673 | "Indreasing pivot tolerance for MA57 from %7.2e ", |
---|
674 | pivtol_); |
---|
675 | pivtol_ = Min(pivtolmax_, pow(pivtol_,0.75)); |
---|
676 | Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA, |
---|
677 | "to %7.2e.\n", |
---|
678 | pivtol_); |
---|
679 | return true; |
---|
680 | } |
---|
681 | |
---|
682 | } // namespace Ipopt |
---|