Context Navigation

source: branches/dev/Algorithm/LinearSolvers/IpMa27TSolverInterface.cpp @ 529

Last change on this file since 529 was 529, checked in by andreasw, 15 years ago
in Vector class, copy cached values in Copy and update them in Scal perform iterative refinement only once for adaptive strategy fix bugs in PDPerturbationHandler avoid some overhead in CalculateQualityFunction? minor changes in timing
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 16.8 KB

Line
1	// Copyright (C) 2004, 2005 International Business Machines and others.
2	// All Rights Reserved.
3	// This code is published under the Common Public License.
4	//
5	// $Id: IpMa27TSolverInterface.cpp 529 2005-09-29 21:12:38Z andreasw $
6	//
7	// Authors: Carl Laird, Andreas Waechter IBM 2005-03-17
8
9	#include "IpMa27TSolverInterface.hpp"
10
11	#ifdef HAVE_CMATH
12	# include <cmath>
13	#else
14	# ifdef HAVE_MATH_H
15	# include <math.h>
16	# else
17	# error "don't have header file for math"
18	# endif
19	#endif
20
21	/** Prototypes for MA27's Fortran subroutines */
22	extern "C"
23	{
24	void F77_FUNC(ma27id,MA27ID)(ipfint* ICNTL, double* CNTL);
25	void F77_FUNC(ma27ad,MA27AD)(ipfint N, ipfint NZ, const ipfint IRN, const ipfint ICN,
26	ipfint IW, ipfint LIW, ipfint* IKEEP, ipfint *IW1,
27	ipfint* NSTEPS, ipfint* IFLAG, ipfint* ICNTL,
28	double* CNTL, ipfint INFO, double OPS);
29	void F77_FUNC(ma27bd,MA27BD)(ipfint N, ipfint NZ, const ipfint IRN, const ipfint ICN,
30	double* A, ipfint* LA, ipfint* IW, ipfint* LIW,
31	ipfint* IKEEP, ipfint* NSTEPS, ipfint* MAXFRT,
32	ipfint* IW1, ipfint* ICNTL, double* CNTL,
33	ipfint* INFO);
34	void F77_FUNC(ma27cd,MA27CD)(ipfint N, double A, ipfint* LA, ipfint* IW,
35	ipfint* LIW, double* W, ipfint* MAXFRT,
36	double* RHS, ipfint* IW1, ipfint* NSTEPS,
37	ipfint* ICNTL, double* CNTL);
38	}
39
40	namespace Ipopt
41	{
42	#ifdef IP_DEBUG
43	static const Index dbg_verbosity = 0;
44	#endif
45
46	Ma27TSolverInterface::Ma27TSolverInterface()
47	:
48	dim_(0),
49	nonzeros_(0),
50	initialized_(false),
51	pivtol_changed_(false),
52	refactorize_(false),
53
54	liw_(0),
55	iw_(NULL),
56	ikeep_(NULL),
57	la_(0),
58	a_(NULL),
59
60	la_increase_(false),
61	liw_increase_(false)
62	{
63	DBG_START_METH("Ma27TSolverInterface::Ma27TSolverInterface()",dbg_verbosity);
64	}
65
66	Ma27TSolverInterface::~Ma27TSolverInterface()
67	{
68	DBG_START_METH("Ma27TSolverInterface::~Ma27TSolverInterface()",
69	dbg_verbosity);
70	delete [] iw_;
71	delete [] ikeep_;
72	delete [] a_;
73	}
74
75	void Ma27TSolverInterface::RegisterOptions(SmartPtr<RegisteredOptions> roptions)
76	{
77	roptions->AddBoundedNumberOption(
78	"pivtol",
79	"Pivot tolerance for the linear solver MA27.",
80	0.0, true, 1.0, true, 1e-8,
81	"A smaller number pivots for sparsity, "
82	"a larger number pivots for stability.");
83	roptions->AddBoundedNumberOption(
84	"pivtolmax",
85	"Maximum pivot tolerance.",
86	0.0, true, 1.0, true, 1e-4,
87	"Ipopt may increase pivtol as high as pivtolmax "
88	"to get a more accurate solution to the linear system.");
89	roptions->AddLowerBoundedNumberOption(
90	"liw_init_factor",
91	"Integer workspace memory for MA27.",
92	1.0, false, 5.0,
93	"The initial integer workspace memory = liw_init_factor * memory "
94	"required by unfactored system. Ipopt will increase the workspace "
95	"size by meminc_factor if required.");
96	roptions->AddLowerBoundedNumberOption(
97	"la_init_factor",
98	"Real workspace memory for MA27.",
99	1.0, false, 5.0,
100	"The initial real workspace memory = la_init_factor * memory "
101	"required by unfactored system. Ipopt will increase the workspace"
102	" size by meminc_factor if required.");
103	roptions->AddLowerBoundedNumberOption(
104	"meminc_factor",
105	"Increment factor for workspace size for MA27.",
106	1.0, false, 10.0,
107	"If the integer or real workspace is not large enough, "
108	"Ipopt will increase its size by this factor.");
109	}
110
111	bool Ma27TSolverInterface::InitializeImpl(const OptionsList& options,
112	const std::string& prefix)
113	{
114	options.GetNumericValue("pivtol", pivtol_, prefix);
115	if(options.GetNumericValue("pivtolmax", pivtolmax_, prefix)) {
116	ASSERT_EXCEPTION(pivtolmax_>=pivtol_, OPTION_INVALID,
117	"Option \"pivtolmax\": This value must be between "
118	"pivtol and 1.");
119	}
120	else {
121	pivtolmax_ = Max(pivtolmax_, pivtol_);
122	}
123
124	options.GetNumericValue("liw_init_factor", liw_init_factor_, prefix);
125	options.GetNumericValue("la_init_factor", la_init_factor_, prefix);
126	options.GetNumericValue("meminc_factor", meminc_factor_, prefix);
127	// The following option is registered by OrigIpoptNLP
128	options.GetBoolValue("warm_start_same_structure",
129	warm_start_same_structure_, prefix);
130
131	/* Set the default options for MA27 */
132	F77_FUNC(ma27id,MA27ID)(icntl_, cntl_);
133	#ifndef IP_DEBUG
134
135	icntl_[0] = 0; // Suppress error messages
136	icntl_[1] = 0; // Suppress diagnostic messages
137	#endif
138
139	// Reset all private data
140	initialized_=false;
141	pivtol_changed_ = false;
142	refactorize_ = false;
143
144	la_increase_=false;
145	liw_increase_=false;
146
147	if (!warm_start_same_structure_) {
148	dim_=0;
149	nonzeros_=0;
150	}
151	else {
152	ASSERT_EXCEPTION(dim_>0 && nonzeros_>0, INVALID_WARMSTART,
153	"Ma27TSolverInterface called with warm_start_same_structure, but the problem is solved for the first time.");
154	}
155
156	return true;
157	}
158
159	ESymSolverStatus Ma27TSolverInterface::MultiSolve(bool new_matrix,
160	const Index* airn,
161	const Index* ajcn,
162	Index nrhs,
163	double* rhs_vals,
164	bool check_NegEVals,
165	Index numberOfNegEVals)
166	{
167	DBG_START_METH("Ma27TSolverInterface::MultiSolve",dbg_verbosity);
168	DBG_ASSERT(!check_NegEVals \|\| ProvidesInertia());
169	DBG_ASSERT(initialized_);
170	DBG_ASSERT(la_!=0);
171
172	if (pivtol_changed_) {
173	DBG_PRINT((1,"Pivot tolerance has changed.\n"));
174	pivtol_changed_ = false;
175	// If the pivot tolerance has been changed but the matrix is not
176	// new, we have to request the values for the matrix again to do
177	// the factorization again.
178	if (!new_matrix) {
179	DBG_PRINT((1,"Ask caller to call again.\n"));
180	refactorize_ = true;
181	return SYMSOLVER_CALL_AGAIN;
182	}
183	}
184
185	// check if a factorization has to be done
186	DBG_PRINT((1, "new_matrix = %d\n", new_matrix));
187	if (new_matrix \|\| refactorize_) {
188	// perform the factorization
189	ESymSolverStatus retval;
190	retval = Factorization(airn, ajcn, check_NegEVals, numberOfNegEVals);
191	if (retval!=SYMSOLVER_SUCCESS) {
192	DBG_PRINT((1, "FACTORIZATION FAILED!\n"));
193	return retval; // Matrix singular or error occurred
194	}
195	refactorize_ = false;
196	}
197
198	// do the backsolve
199	return Backsolve(nrhs, rhs_vals);
200	}
201
202	double* Ma27TSolverInterface::GetValuesArrayPtr()
203	{
204	DBG_START_METH("Ma27TSolverInterface::GetValuesArrayPtr",dbg_verbosity);
205	DBG_ASSERT(initialized_);
206
207	// If the size of a is to be increase for the next factorization
208	// anyway, delete the current large array and just return enough
209	// to store the values
210
211	if (la_increase_) {
212	delete [] a_;
213	a_ = new double [nonzeros_];
214	}
215
216	return a_;
217	}
218
219	/** Initialize the local copy of the positions of the nonzero
220	elements */
221	ESymSolverStatus Ma27TSolverInterface::InitializeStructure(Index dim, Index nonzeros,
222	const Index* airn,
223	const Index* ajcn)
224	{
225	DBG_START_METH("Ma27TSolverInterface::InitializeStructure",dbg_verbosity);
226
227	ESymSolverStatus retval = SYMSOLVER_SUCCESS;
228	if (!warm_start_same_structure_) {
229	dim_ = dim;
230	nonzeros_ = nonzeros;
231
232	// Do the symbolic facotrization
233	retval = SymbolicFactorization(airn, ajcn);
234	if (retval != SYMSOLVER_SUCCESS ) {
235	return retval;
236	}
237	}
238	else {
239	ASSERT_EXCEPTION(dim_==dim && nonzeros_==nonzeros, INVALID_WARMSTART,
240	"Ma27TSolverInterface called with warm_start_same_structure, but the problem size has changed.");
241	}
242
243	initialized_ = true;
244
245	return retval;
246	}
247
248	ESymSolverStatus Ma27TSolverInterface::SymbolicFactorization(const Index* airn,
249	const Index* ajcn)
250	{
251	DBG_START_METH("Ma27TSolverInterface::SymbolicFactorization",dbg_verbosity);
252
253	IpData().TimingStats().LinearSystemSymbolicFactorization().Start();
254
255	// Get memory for the IW workspace
256	delete [] iw_;
257
258	// Overstimation factor for LIW (20% recommended in MA27 documentation)
259	const double LiwFact = 2.0; // This is 100% overestimation
260	Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
261	"In Ma27TSolverInterface::InitializeStructure: Using overestimation factor LiwFact = %e\n",
262	LiwFact);
263	liw_ = (ipfint)(LiwFact(double(2nonzeros_+3*dim_+1)));
264	iw_ = new ipfint[liw_];
265
266	// Get memory for IKEEP
267	delete [] ikeep_;
268	ikeep_ = new ipfint[3*dim_];
269
270	// Call MA27AD (cast to ipfint for Index types)
271	ipfint N = dim_;
272	ipfint NZ = nonzeros_;
273	ipfint IFLAG = 0;
274	double OPS;
275	ipfint INFO[20];
276	ipfint* IW1 = new ipfint[2*dim_]; // Get memory for IW1 (only local)
277	F77_FUNC(ma27ad,MA27AD)(&N, &NZ, airn, ajcn, iw_, &liw_, ikeep_,
278	IW1, &nsteps_, &IFLAG, icntl_, cntl_,
279	INFO, &OPS);
280	delete [] IW1; // No longer required
281
282	// Receive several information
283	ipfint iflag = INFO[0]; // Information flag
284	ipfint ierror = INFO[1]; // Error flag
285	ipfint nrlnec = INFO[4]; // recommended value for la
286	ipfint nirnec = INFO[5]; // recommended value for liw
287
288	// Check if error occurred
289	if (iflag!=0) {
290	Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,
291	"* Error from MA27AD * IFLAG = %d IERROR = %d\n", iflag, ierror);
292	if (iflag==1) {
293	Jnlst().Printf(J_ERROR, J_LINEAR_ALGEBRA,
294	"The index a matrix is out of range.\nPlease check your implementation of the Jabobian and Hessian matrices.");
295	}
296	return SYMSOLVER_FATAL_ERROR;
297	}
298
299	// ToDo: try and catch
300	// Reserve memory for iw_ for later calls, based on suggested size
301	delete [] iw_;
302	liw_ = (ipfint)(liw_init_factor_ * (double)(nirnec));
303	iw_ = new ipfint[liw_];
304
305	// Reserve memory for a_
306	delete [] a_;
307	la_ = Max(nonzeros_,(ipfint)(la_init_factor_ * (double)(nrlnec)));
308	a_ = new double[la_];
309
310	IpData().TimingStats().LinearSystemSymbolicFactorization().End();
311
312	return SYMSOLVER_SUCCESS;
313	}
314
315	ESymSolverStatus
316	Ma27TSolverInterface::Factorization(const Index* airn,
317	const Index* ajcn,
318	bool check_NegEVals,
319	Index numberOfNegEVals)
320	{
321	DBG_START_METH("Ma27TSolverInterface::Factorization",dbg_verbosity);
322	// Check if la should be increased
323	IpData().TimingStats().LinearSystemFactorization().Start();
324	if (la_increase_) {
325	double* a_old = a_;
326	ipfint la_old = la_;
327	la_ = (ipfint)(meminc_factor_ * (double)(la_));
328	a_ = new double[la_];
329	for (Index i=0; i<nonzeros_; i++) {
330	a_[i] = a_old[i];
331	}
332	delete [] a_old;
333	la_increase_ = false;
334	Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
335	"In Ma27TSolverInterface::Factorization: Increasing la from %d to %d\n",
336	la_old, la_);
337	}
338
339	// Check if liw should be increased
340	if (liw_increase_) {
341	delete [] iw_;
342	ipfint liw_old = liw_;
343	liw_ = (ipfint)(meminc_factor_ * (double)(liw_));
344	iw_ = new ipfint[liw_];
345	liw_increase_ = false;
346	Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
347	"In Ma27TSolverInterface::Factorization: Increasing liw from %d to %d\n",
348	liw_old, liw_);
349	}
350
351	ipfint iflag; // Information flag
352	ipfint ncmpbr; // Number of double precision compressions
353	ipfint ncmpbi; // Number of integer compressions
354
355	// Call MA27BD; possibly repeatedly if workspaces are too small
356	ipfint N=dim_;
357	ipfint NZ=nonzeros_;
358	ipfint* IW1 = new ipfint[2*dim_];
359	ipfint INFO[20];
360	cntl_[0] = pivtol_; // Set pivot tolerance
361
362	F77_FUNC(ma27bd,MA27BD)(&N, &NZ, airn, ajcn, a_,
363	&la_, iw_, &liw_, ikeep_, &nsteps_,
364	&maxfrt_, IW1, icntl_, cntl_, INFO);
365	delete [] IW1;
366
367	// Receive information about the factorization
368	iflag = INFO[0]; // Information flag
369	ipfint ierror = INFO[1]; // Error flag
370	ncmpbr = INFO[11]; // Number of double compressions
371	ncmpbi = INFO[12]; // Number of integer compressions
372	negevals_ = INFO[14]; // Number of negative eigenvalues
373
374	DBG_PRINT((1,"Return from MA27 iflag = %d and ierror = %d\n",
375	iflag, ierror));
376
377	// Check if factorization failed due to insufficient memory space
378	// iflag==-3 if LIW too small (recommended value in ierror)
379	// iflag==-4 if LA too small (recommended value in ierror)
380	if (iflag==-3 \|\| iflag==-4) {
381	// Increase size of both LIW and LA
382	delete [] iw_;
383	delete [] a_;
384	ipfint liw_old = liw_;
385	ipfint la_old = la_;
386	if(iflag==-3) {
387	liw_ = (ipfint)(meminc_factor_ * (double)(ierror));
388	la_ = (ipfint)(meminc_factor_ * (double)(la_));
389	}
390	else {
391	liw_ = (ipfint)(meminc_factor_ * (double)(liw_));
392	la_ = (ipfint)(meminc_factor_ * (double)(ierror));
393	}
394	iw_ = new ipfint[liw_];
395	a_ = new double[la_];
396	Jnlst().Printf(J_WARNING, J_LINEAR_ALGEBRA,
397	"MA27BD returned iflag=%d.\n Increase liw from %d to %d and la from %d to %d and factorize again.\n",
398	iflag, liw_old, liw_, la_old, la_);
399	IpData().TimingStats().LinearSystemFactorization().End();
400	return SYMSOLVER_CALL_AGAIN;
401	}
402
403	// Check if the system is singular, and if some other error occurred
404	if (iflag==-5 \|\| iflag==3) {
405	IpData().TimingStats().LinearSystemFactorization().End();
406	return SYMSOLVER_SINGULAR;
407	}
408	else if (iflag != 0) {
409	// There is some error
410	IpData().TimingStats().LinearSystemFactorization().End();
411	return SYMSOLVER_FATAL_ERROR;
412	}
413
414	// Check if it might be more efficient to use more memory next time
415	// (if there were too many compressions for this factorization)
416	if (ncmpbr>=10) {
417	la_increase_ = true;
418	Jnlst().Printf(J_WARNING, J_LINEAR_ALGEBRA,
419	"MA27BD returned ncmpbr=%d. Increase la before the next factorization.\n",
420	ncmpbr);
421	}
422	if (ncmpbi>=10) {
423	liw_increase_ = true;
424	Jnlst().Printf(J_WARNING, J_LINEAR_ALGEBRA,
425	"MA27BD returned ncmpbi=%d. Increase liw before the next factorization.\n",
426	ncmpbr);
427	}
428
429	// Check whether the number of negative eigenvalues matches the requested
430	// count
431	if (check_NegEVals && (numberOfNegEVals!=negevals_)) {
432	Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
433	"In Ma27TSolverInterface::Factorization: negevals_ = %d, but numberOfNegEVals = %d\n",
434	negevals_, numberOfNegEVals);
435	IpData().TimingStats().LinearSystemFactorization().End();
436	return SYMSOLVER_WRONG_INERTIA;
437	}
438
439	IpData().TimingStats().LinearSystemFactorization().End();
440	return SYMSOLVER_SUCCESS;
441	}
442
443	ESymSolverStatus Ma27TSolverInterface::Backsolve(Index nrhs,
444	double *rhs_vals)
445	{
446	DBG_START_METH("Ma27TSolverInterface::Backsolve",dbg_verbosity);
447	IpData().TimingStats().LinearSystemBackSolve().Start();
448
449	ipfint N=dim_;
450	double* W = new double[maxfrt_];
451	ipfint* IW1 = new ipfint[nsteps_];
452
453	// For each right hand side, call MA27CD
454	for(Index irhs=0; irhs<nrhs; irhs++) {
455	if (DBG_VERBOSITY()>=2) {
456	for (Index i=0; i<dim_; i++) {
457	DBG_PRINT((2, "rhs[%5d] = %23.15e\n", i, rhs_vals[irhs*dim_+i]));
458	}
459	}
460
461	F77_FUNC(ma27cd,MA27CD)(&N, a_, &la_, iw_, &liw_, W, &maxfrt_,
462	&rhs_vals[irhs*dim_], IW1, &nsteps_,
463	icntl_, cntl_);
464
465	if (DBG_VERBOSITY()>=2) {
466	for (Index i=0; i<dim_; i++) {
467	DBG_PRINT((2, "sol[%5d] = %23.15e\n", i, rhs_vals[irhs*dim_+i]));
468	}
469	}
470	}
471	delete [] W;
472	delete [] IW1;
473
474	IpData().TimingStats().LinearSystemBackSolve().End();
475	return SYMSOLVER_SUCCESS;
476	}
477
478	Index Ma27TSolverInterface::NumberOfNegEVals() const
479	{
480	DBG_START_METH("Ma27TSolverInterface::NumberOfNegEVals",dbg_verbosity);
481	DBG_ASSERT(ProvidesInertia());
482	DBG_ASSERT(initialized_);
483	return negevals_;
484	}
485
486	bool Ma27TSolverInterface::IncreaseQuality()
487	{
488	DBG_START_METH("Ma27TSolverInterface::IncreaseQuality",dbg_verbosity);
489	if (pivtol_ == pivtolmax_) {
490	return false;
491	}
492	pivtol_changed_ = true;
493
494	Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
495	"Indreasing pivot tolerance from %7.2e ",
496	pivtol_);
497	pivtol_ = Min(pivtolmax_, pow(pivtol_,0.75));
498	Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
499	"to %7.2e.\n",
500	pivtol_);
501	return true;
502	}
503
504	} // namespace Ipopt

Note: See TracBrowser for help on using the repository browser.

Download in other formats: