Context Navigation

source: branches/dev/Algorithm/IpPDFullSpaceSolver.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: 29.5 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: IpPDFullSpaceSolver.cpp 529 2005-09-29 21:12:38Z andreasw $
6	//
7	// Authors: Carl Laird, Andreas Waechter IBM 2004-08-13
8
9	#include "IpPDFullSpaceSolver.hpp"
10	#include "IpDebug.hpp"
11
12	namespace Ipopt
13	{
14
15	#ifdef IP_DEBUG
16	static const Index dbg_verbosity = 0;
17	#endif
18
19	PDFullSpaceSolver::PDFullSpaceSolver(AugSystemSolver& augSysSolver,
20	PDPerturbationHandler& perturbHandler)
21	:
22	PDSystemSolver(),
23	augSysSolver_(&augSysSolver),
24	perturbHandler_(&perturbHandler),
25	dummy_cache_(1)
26	{
27	DBG_START_METH("PDFullSpaceSolver::PDFullSpaceSolver",dbg_verbosity);
28	}
29
30	PDFullSpaceSolver::~PDFullSpaceSolver()
31	{
32	DBG_START_METH("PDFullSpaceSolver::~PDFullSpaceSolver()",dbg_verbosity);
33	}
34
35	void PDFullSpaceSolver::RegisterOptions(SmartPtr<RegisteredOptions> roptions)
36	{
37	roptions->AddLowerBoundedIntegerOption(
38	"min_refinement_steps",
39	"Minimum number of iterative refinement steps per linear system solve.",
40	0, 1,
41	"Iterative refinement (on the full unsymmetric system) is performed for "
42	"each right hand side. This option determines the minimum number "
43	"of iterative refinements (i.e. at least \"min_refinement_steps\" "
44	"iterative refinement steps are enforced per right hand side.)");
45	roptions->AddLowerBoundedIntegerOption(
46	"max_refinement_steps",
47	"Maximum number of iterative refinement steps per linear system solve.",
48	0, 10,
49	"Iterative refinement (on the full unsymmetric system) is performed for "
50	"each right hand side. This option determines the maximum number "
51	"of iterative refinement steps.");
52	roptions->AddLowerBoundedNumberOption(
53	"residual_ratio_max",
54	"Iterative refinement tolerance",
55	0.0, true, 1e-10,
56	"Iterative refinement is performed until the residual test ratio is "
57	"less than this tolerance (or until \"max_refinement_steps\" refinement "
58	"steps are performed).");
59	roptions->AddLowerBoundedNumberOption(
60	"residual_ratio_singular",
61	"Threshold for declaring linear system singular after failed iterative refinement.",
62	0.0, true, 1e-5,
63	"If the residual test ratio is larger than this value after failed "
64	"iterative refinement, the algorithm pretends that the linear system is "
65	"singular.");
66	// ToDo Think about following option - are the correct norms used?
67	roptions->AddLowerBoundedNumberOption(
68	"residual_improvement_factor",
69	"Minimal required reduction of residual test ratio in iterative refinement.",
70	0.0, true, 0.999999999,
71	"If the improvement of the residual test ratio made by one iterative "
72	"refinement step is not better than this factor, iterative refinement "
73	"is aborted.");
74	}
75
76
77	bool PDFullSpaceSolver::InitializeImpl(const OptionsList& options,
78	const std::string& prefix)
79	{
80	// Check for the algorithm options
81	options.GetIntegerValue("min_refinement_steps", min_refinement_steps_, prefix);
82	options.GetIntegerValue("max_refinement_steps", max_refinement_steps_, prefix);
83	ASSERT_EXCEPTION(max_refinement_steps_ >= min_refinement_steps_, OPTION_INVALID,
84	"Option \"max_refinement_steps\": This value must be larger than or equal to min_refinement_steps (default 1)");
85
86	options.GetNumericValue("residual_ratio_max", residual_ratio_max_, prefix);
87	options.GetNumericValue("residual_ratio_singular", residual_ratio_singular_, prefix);
88	ASSERT_EXCEPTION(residual_ratio_singular_ > residual_ratio_max_, OPTION_INVALID,
89	"Option \"residual_ratio_singular\": This value must be larger than residual_ratio_max.");
90	options.GetNumericValue("residual_improvement_factor", residual_improvement_factor_, prefix);
91
92	// Reset internal flags and data
93	augsys_improved_ = false;
94
95	if (!augSysSolver_->Initialize(Jnlst(), IpNLP(), IpData(), IpCq(),
96	options, prefix)) {
97	return false;
98	}
99
100	return perturbHandler_->Initialize(Jnlst(), IpNLP(), IpData(), IpCq(),
101	options, prefix);
102	}
103
104	void PDFullSpaceSolver::Solve(Number alpha,
105	Number beta,
106	const IteratesVector& rhs,
107	IteratesVector& res,
108	bool allow_inexact,
109	bool improve_solution /* = false */)
110	{
111	DBG_START_METH("PDFullSpaceSolver::Solve",dbg_verbosity);
112	DBG_ASSERT(!allow_inexact \|\| !improve_solution);
113	DBG_ASSERT(!improve_solution \|\| beta==0.);
114
115	// Timing of PDSystem solver starts here
116	IpData().TimingStats().PDSystemSolverTotal().Start();
117
118	DBG_PRINT_VECTOR(2, "rhs_x", *rhs.x());
119	DBG_PRINT_VECTOR(2, "rhs_s", *rhs.s());
120	DBG_PRINT_VECTOR(2, "rhs_c", *rhs.y_c());
121	DBG_PRINT_VECTOR(2, "rhs_d", *rhs.y_d());
122	DBG_PRINT_VECTOR(2, "rhs_zL", *rhs.z_L());
123	DBG_PRINT_VECTOR(2, "rhs_zU", *rhs.z_U());
124	DBG_PRINT_VECTOR(2, "rhs_vL", *rhs.v_L());
125	DBG_PRINT_VECTOR(2, "rhs_vU", *rhs.v_U());
126	DBG_PRINT_VECTOR(2, "res_x in", *res.x());
127	DBG_PRINT_VECTOR(2, "res_s in", *res.s());
128	DBG_PRINT_VECTOR(2, "res_c in", *res.y_c());
129	DBG_PRINT_VECTOR(2, "res_d in", *res.y_d());
130	DBG_PRINT_VECTOR(2, "res_zL in", *res.z_L());
131	DBG_PRINT_VECTOR(2, "res_zU in", *res.z_U());
132	DBG_PRINT_VECTOR(2, "res_vL in", *res.v_L());
133	DBG_PRINT_VECTOR(2, "res_vU in", *res.v_U());
134
135	// if beta is nonzero, keep a copy of the incoming values in res_ */
136	SmartPtr<IteratesVector> copy_res;
137	if (beta != 0.) {
138	copy_res = res.MakeNewIteratesVectorCopy();
139	}
140
141	// Receive data about matrix
142	SmartPtr<const Vector> x = IpData().curr()->x();
143	SmartPtr<const Vector> s = IpData().curr()->s();
144	SmartPtr<const SymMatrix> W = IpData().W();
145	SmartPtr<const Matrix> J_c = IpCq().curr_jac_c();
146	SmartPtr<const Matrix> J_d = IpCq().curr_jac_d();
147	SmartPtr<const Matrix> Px_L = IpNLP().Px_L();
148	SmartPtr<const Matrix> Px_U = IpNLP().Px_U();
149	SmartPtr<const Matrix> Pd_L = IpNLP().Pd_L();
150	SmartPtr<const Matrix> Pd_U = IpNLP().Pd_U();
151	SmartPtr<const Vector> z_L = IpData().curr()->z_L();
152	SmartPtr<const Vector> z_U = IpData().curr()->z_U();
153	SmartPtr<const Vector> v_L = IpData().curr()->v_L();
154	SmartPtr<const Vector> v_U = IpData().curr()->v_U();
155	SmartPtr<const Vector> slack_x_L = IpCq().curr_slack_x_L();
156	SmartPtr<const Vector> slack_x_U = IpCq().curr_slack_x_U();
157	SmartPtr<const Vector> slack_s_L = IpCq().curr_slack_s_L();
158	SmartPtr<const Vector> slack_s_U = IpCq().curr_slack_s_U();
159	SmartPtr<const Vector> sigma_x = IpCq().curr_sigma_x();
160	SmartPtr<const Vector> sigma_s = IpCq().curr_sigma_s();
161	DBG_PRINT_VECTOR(2, "Sigma_x", *sigma_x);
162	DBG_PRINT_VECTOR(2, "Sigma_s", *sigma_s);
163
164	bool done = false;
165	// The following flag is set to true, if we asked the linear
166	// solver to improve the quality of the solution in
167	// the next solve
168	bool resolve_with_better_quality = false;
169	// the following flag is set to true, if iterative refinement
170	// failed and we want to try if a modified system is able to
171	// remedy that problem by pretending the matrix is singular
172	bool pretend_singular = false;
173	bool pretend_singular_last_time = false;
174
175	// Beginning of loop for solving the system (including all
176	// modifications for the linear system to ensure good solution
177	// quality)
178	while (!done) {
179
180	// if improve_solution is true, we are given already a solution
181	// from the calling function, so we can skip the first solve
182	bool solve_retval = true;
183	if (!improve_solution) {
184	solve_retval =
185	SolveOnce(resolve_with_better_quality, pretend_singular,
186	W, J_c, J_d, Px_L, Px_U, Pd_L, Pd_U, z_L, *z_U,
187	v_L, v_U, slack_x_L, slack_x_U, slack_s_L, slack_s_U,
188	sigma_x, sigma_s, 1., 0., rhs, res);
189	resolve_with_better_quality = false;
190	pretend_singular = false;
191	}
192	improve_solution = false;
193
194	if (!solve_retval) {
195	// If system seems not to be solvable, we set the search
196	// direction to zero, and hope that the line search will take
197	// care of this (e.g. call the restoration phase). ToDo: We
198	// might want to use a more explicit cue later.
199	res.Set(0.0);
200	IpData().TimingStats().PDSystemSolverTotal().End();
201	return;
202	}
203
204	if (allow_inexact) {
205	// no safety checks required
206	break;
207	}
208
209	// Get space for the residual
210	SmartPtr<IteratesVector> resid = res.MakeNewIteratesVector(true);
211
212	// ToDo don't to that after max refinement?
213	ComputeResiduals(W, J_c, J_d, Px_L, Px_U, Pd_L, *Pd_U,
214	z_L, z_U, v_L, v_U, slack_x_L, slack_x_U,
215	slack_s_L, slack_s_U, sigma_x, sigma_s,
216	alpha, beta, rhs, res, *resid);
217
218	Number residual_ratio =
219	ComputeResidualRatio(rhs, res, *resid);
220	Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
221	"residual_ratio = %e\n", residual_ratio);
222	Number residual_ratio_old = residual_ratio;
223
224	// Beginning of loop for iterative refinement
225	Index num_iter_ref = 0;
226	bool quit_refinement = false;
227	while (!allow_inexact && !quit_refinement &&
228	(num_iter_ref < min_refinement_steps_ \|\|
229	residual_ratio > residual_ratio_max_) ) {
230
231	// To the next back solve
232	solve_retval =
233	SolveOnce(resolve_with_better_quality, false,
234	W, J_c, J_d, Px_L, Px_U, Pd_L, Pd_U, z_L, *z_U,
235	v_L, v_U, slack_x_L, slack_x_U, slack_s_L, slack_s_U,
236	sigma_x, sigma_s, -1., 1., *resid, res);
237	ASSERT_EXCEPTION(solve_retval, INTERNAL_ABORT,
238	"SolveOnce returns false.");
239
240	ComputeResiduals(W, J_c, J_d, Px_L, Px_U, Pd_L, *Pd_U,
241	z_L, z_U, v_L, v_U, slack_x_L, slack_x_U,
242	slack_s_L, slack_s_U, sigma_x, sigma_s,
243	alpha, beta, rhs, res, *resid);
244
245	residual_ratio =
246	ComputeResidualRatio(rhs, res, *resid);
247	Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
248	"residual_ratio = %e\n", residual_ratio);
249
250	num_iter_ref++;
251	// Check if we have to give up on iterative refinement
252	if (num_iter_ref>min_refinement_steps_ &&
253	(num_iter_ref>max_refinement_steps_ \|\|
254	residual_ratio>residual_improvement_factor_*residual_ratio_old)) {
255
256	Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
257	"Iterative refinement failed with residual_ratio = %e\n", residual_ratio);
258	quit_refinement = true;
259
260	// Pretend singularity only once - if it didn't help, we
261	// have to live with what we got so far
262	resolve_with_better_quality = false;
263	DBG_PRINT((1, "pretend_singular = %d\n", pretend_singular));
264	if (!pretend_singular_last_time) {
265	// First try if we can ask the augmented system solver to
266	// improve the quality of the solution (only if that hasn't
267	// been done before for this linear system)
268	if (!augsys_improved_) {
269	Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
270	"Asking augmented system solver to improve quality of its solutions.\n");
271	augsys_improved_ = augSysSolver_->IncreaseQuality();
272	if (augsys_improved_) {
273	IpData().Append_info_string("q");
274	resolve_with_better_quality = true;
275	}
276	else {
277	// solver said it cannot improve quality, so let
278	// possibly conclude that the current modification is
279	// singular
280	pretend_singular = true;
281	}
282	}
283	else {
284	// we had already asked the solver before to improve the
285	// quality of the solution, so let's now pretend that the
286	// modification is possibly singular
287	pretend_singular = true;
288	}
289	pretend_singular_last_time = pretend_singular;
290	if (pretend_singular) {
291	// let's only conclude that the current linear system
292	// including modifications is singular, if the residual is
293	// quite bad
294	if (residual_ratio < residual_ratio_singular_) {
295	pretend_singular = false;
296	IpData().Append_info_string("S");
297	Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
298	"Just accept current solution.\n");
299	}
300	else {
301	IpData().Append_info_string("s");
302	Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
303	"Pretend that the current system (including modifications) is singular.\n");
304	}
305	}
306	}
307	else {
308	pretend_singular = false;
309	DBG_PRINT((1,"Resetting pretend_singular to false.\n"));
310	}
311	}
312
313	residual_ratio_old = residual_ratio;
314	} // End of loop for iterative refinement
315
316	done = !(resolve_with_better_quality) && !(pretend_singular);
317
318	} // End of loop for solving the linear system (incl. modifications)
319
320	// Finally let's assemble the res result vectors
321	if (alpha != 0.) {
322	res.Scal(alpha);
323	}
324
325	if (beta != 0.) {
326	res.Axpy(beta, *copy_res);
327	}
328
329	DBG_PRINT_VECTOR(2, "res_x", *res.x());
330	DBG_PRINT_VECTOR(2, "res_s", *res.s());
331	DBG_PRINT_VECTOR(2, "res_c", *res.y_c());
332	DBG_PRINT_VECTOR(2, "res_d", *res.y_d());
333	DBG_PRINT_VECTOR(2, "res_zL", *res.z_L());
334	DBG_PRINT_VECTOR(2, "res_zU", *res.z_U());
335	DBG_PRINT_VECTOR(2, "res_vL", *res.v_L());
336	DBG_PRINT_VECTOR(2, "res_vU", *res.v_U());
337
338	IpData().TimingStats().PDSystemSolverTotal().End();
339	}
340
341	bool PDFullSpaceSolver::SolveOnce(bool resolve_with_better_quality,
342	bool pretend_singular,
343	const SymMatrix& W,
344	const Matrix& J_c,
345	const Matrix& J_d,
346	const Matrix& Px_L,
347	const Matrix& Px_U,
348	const Matrix& Pd_L,
349	const Matrix& Pd_U,
350	const Vector& z_L,
351	const Vector& z_U,
352	const Vector& v_L,
353	const Vector& v_U,
354	const Vector& slack_x_L,
355	const Vector& slack_x_U,
356	const Vector& slack_s_L,
357	const Vector& slack_s_U,
358	const Vector& sigma_x,
359	const Vector& sigma_s,
360	Number alpha,
361	Number beta,
362	const IteratesVector& rhs,
363	IteratesVector& res)
364	{
365	// TO DO LIST:
366	//
367	// 1. decide for reasonable return codes (e.g. fatal error, too
368	// ill-conditioned...)
369	// 2. Make constants parameters that can be set from the outside
370	// 3. Get Information out of Ipopt structures
371	// 4. add heuristic for structurally singular problems
372	// 5. see if it makes sense to distinguish delta_x and delta_s,
373	// or delta_c and delta_d
374	// 6. increase pivot tolerance if number of get evals so too small
375	DBG_START_METH("PDFullSpaceSolver::SolveOnce",dbg_verbosity);
376
377	IpData().TimingStats().PDSystemSolverSolveOnce().Start();
378
379	// Compute the right hand side for the augmented system formulation
380	SmartPtr<Vector> augRhs_x = rhs.x()->MakeNewCopy();
381	Px_L.AddMSinvZ(1.0, slack_x_L, rhs.z_L(), augRhs_x);
382	Px_U.AddMSinvZ(-1.0, slack_x_U, rhs.z_U(), augRhs_x);
383
384	SmartPtr<Vector> augRhs_s = rhs.s()->MakeNewCopy();
385	Pd_L.AddMSinvZ(1.0, slack_s_L, rhs.v_L(), augRhs_s);
386	Pd_U.AddMSinvZ(-1.0, slack_s_U, rhs.v_U(), augRhs_s);
387
388	// Get space into which we can put the solution of the augmented system
389	SmartPtr<IteratesVector> sol = res.MakeNewIteratesVector(true);
390
391	// Now check whether any data has changed
392	std::vector<const TaggedObject*> deps(13);
393	deps[0] = &W;
394	deps[1] = &J_c;
395	deps[2] = &J_d;
396	deps[3] = &z_L;
397	deps[4] = &z_U;
398	deps[5] = &v_L;
399	deps[6] = &v_U;
400	deps[7] = &slack_x_L;
401	deps[8] = &slack_x_U;
402	deps[9] = &slack_s_L;
403	deps[10] = &slack_s_U;
404	deps[11] = &sigma_x;
405	deps[12] = &sigma_s;
406	void* dummy;
407	bool uptodate = dummy_cache_.GetCachedResult(dummy, deps);
408	if (!uptodate) {
409	dummy_cache_.AddCachedResult(dummy, deps);
410	augsys_improved_ = false;
411	}
412	// improve_current_solution can only be true, if that system has
413	// been solved before
414	DBG_ASSERT((!resolve_with_better_quality && !pretend_singular) \|\| uptodate);
415
416	ESymSolverStatus retval;
417	if (uptodate && !pretend_singular) {
418
419	// Get the perturbation values
420	Number delta_x;
421	Number delta_s;
422	Number delta_c;
423	Number delta_d;
424	perturbHandler_->CurrentPerturbation(delta_x, delta_s, delta_c, delta_d);
425
426	// No need to go throught the pain of finding the appropriate
427	// values for the deltas, because the matrix hasn't changed since
428	// the last call. So, just call the Solve Method
429	//
430	// Note: resolve_with_better_quality is true, then the Solve
431	// method has already asked the augSysSolver to increase the
432	// quality at the end solve, and we are now getting the solution
433	// with that better quality
434	retval = augSysSolver_->Solve(&W, &sigma_x, delta_x,
435	&sigma_s, delta_s, &J_c, NULL,
436	delta_c, &J_d, NULL, delta_d,
437	augRhs_x, augRhs_s, rhs.y_c(), rhs.y_d(),
438	sol->x_NonConst(), sol->s_NonConst(),
439	sol->y_c_NonConst(), sol->y_d_NonConst(),
440	false, 0);
441	if (retval!=SYMSOLVER_SUCCESS) {
442	IpData().TimingStats().PDSystemSolverSolveOnce().End();
443	return false;
444	}
445	}
446	else {
447	Index numberOfEVals=rhs.y_c()->Dim()+rhs.y_d()->Dim();
448	// counter for the number of trial evaluations
449	// (ToDo is not at the correct place)
450	Index count = 0;
451
452	// Get the very first perturbation values from the perturbation
453	// Handler
454	Number delta_x;
455	Number delta_s;
456	Number delta_c;
457	Number delta_d;
458	perturbHandler_->ConsiderNewSystem(delta_x, delta_s, delta_c, delta_d);
459
460	retval = SYMSOLVER_SINGULAR;
461	bool fail = false;
462
463	while (retval!= SYMSOLVER_SUCCESS && !fail) {
464
465	if (pretend_singular) {
466	retval = SYMSOLVER_SINGULAR;
467	pretend_singular = false;
468	}
469	else {
470	count++;
471	Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA,
472	"Solving system with delta_x=%e delta_s=%e\n delta_c=%e delta_d=%e\n",
473	delta_x, delta_s, delta_c, delta_d);
474	retval = augSysSolver_->Solve(&W, &sigma_x, delta_x,
475	&sigma_s, delta_s, &J_c, NULL,
476	delta_c, &J_d, NULL, delta_d,
477	augRhs_x, augRhs_s, rhs.y_c(), rhs.y_d(),
478	sol->x_NonConst(), sol->s_NonConst(),
479	sol->y_c_NonConst(), sol->y_d_NonConst(),
480	true, numberOfEVals);
481	}
482	assert(retval!=SYMSOLVER_FATAL_ERROR); //TODO make return code
483	if (retval==SYMSOLVER_SINGULAR &&
484	(rhs.y_c()->Dim()+rhs.y_d()->Dim() > 0) ) {
485
486	// Get new perturbation factors from the perturbation
487	// handlers for the singular case
488	bool pert_return =
489	perturbHandler_->PerturbForSingularity(delta_x, delta_s,
490	delta_c, delta_d);
491	if (!pert_return) {
492	Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
493	"PerturbForSingularity can't be done\n");
494	IpData().TimingStats().PDSystemSolverSolveOnce().End();
495	return false;
496	}
497	}
498	else if (retval==SYMSOLVER_WRONG_INERTIA &&
499	augSysSolver_->NumberOfNegEVals() < numberOfEVals) {
500	Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
501	"Number of negative eigenvalues too small!\n");
502	// If the number of negative eigenvalues is too small, then
503	// we first try to remedy this by asking for better quality
504	// solution (e.g. increasing pivot tolerance), and if that
505	// doesn't help, we assume that the system is singular
506	bool assume_singular = true;
507	if (!augsys_improved_) {
508	Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
509	"Asking augmented system solver to improve quality of its solutions.\n");
510	augsys_improved_ = augSysSolver_->IncreaseQuality();
511	if (augsys_improved_) {
512	IpData().Append_info_string("q");
513	resolve_with_better_quality = true;
514	assume_singular = false;
515	}
516	else {
517	Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
518	"Quality could not be improved\n");
519	}
520	}
521	if (assume_singular) {
522	bool pert_return =
523	perturbHandler_->PerturbForSingularity(delta_x, delta_s,
524	delta_c, delta_d);
525	if (!pert_return) {
526	Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
527	"PerturbForSingularity can't be done for assume singular.\n");
528	IpData().TimingStats().PDSystemSolverSolveOnce().End();
529	return false;
530	}
531	IpData().Append_info_string("a");
532	}
533	}
534	else if (retval==SYMSOLVER_WRONG_INERTIA \|\|
535	retval==SYMSOLVER_SINGULAR) {
536	// Get new perturbation factors from the perturbation
537	// handlers for the case of wrong inertia
538	bool pert_return =
539	perturbHandler_->PerturbForWrongInertia(delta_x, delta_s,
540	delta_c, delta_d);
541	if (!pert_return) {
542	Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
543	"PerturbForWrongInertia can't be done for assume singular.\n");
544	IpData().TimingStats().PDSystemSolverSolveOnce().End();
545	return false;
546	}
547	}
548	} // while (retval!=SYMSOLVER_SUCCESS && !fail) {
549
550	// Some output
551	Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
552	"Number of trial factorizations performed: %d\n",
553	count);
554	Jnlst().Printf(J_DETAILED, J_LINEAR_ALGEBRA,
555	"Perturbation parameters: delta_x=%e delta_s=%e\n delta_c=%e delta_d=%e\n",
556	delta_x, delta_s, delta_c, delta_d);
557	}
558
559	// Compute the remaining sol Vectors
560	Px_L.SinvBlrmZMTdBr(-1., slack_x_L, rhs.z_L(), z_L, sol->x_NonConst(), *sol->z_L_NonConst());
561	Px_U.SinvBlrmZMTdBr(1., slack_x_U, rhs.z_U(), z_U, sol->x_NonConst(), *sol->z_U_NonConst());
562	Pd_L.SinvBlrmZMTdBr(-1., slack_s_L, rhs.v_L(), v_L, sol->s_NonConst(), *sol->v_L_NonConst());
563	Pd_U.SinvBlrmZMTdBr(1., slack_s_U, rhs.v_U(), v_U, sol->s_NonConst(), *sol->v_U_NonConst());
564
565	// Finally let's assemble the res result vectors
566	res.AddOneVector(alpha, *sol, beta);
567
568	IpData().TimingStats().PDSystemSolverSolveOnce().End();
569
570	return true;
571	}
572
573	void PDFullSpaceSolver::ComputeResiduals(
574	const SymMatrix& W,
575	const Matrix& J_c,
576	const Matrix& J_d,
577	const Matrix& Px_L,
578	const Matrix& Px_U,
579	const Matrix& Pd_L,
580	const Matrix& Pd_U,
581	const Vector& z_L,
582	const Vector& z_U,
583	const Vector& v_L,
584	const Vector& v_U,
585	const Vector& slack_x_L,
586	const Vector& slack_x_U,
587	const Vector& slack_s_L,
588	const Vector& slack_s_U,
589	const Vector& sigma_x,
590	const Vector& sigma_s,
591	Number alpha,
592	Number beta,
593	const IteratesVector& rhs,
594	const IteratesVector& res,
595	IteratesVector& resid)
596	{
597	DBG_START_METH("PDFullSpaceSolver::ComputeResiduals", dbg_verbosity);
598
599	DBG_PRINT_VECTOR(2, "res", res);
600	IpData().TimingStats().ComputeResiduals().Start();
601
602	// Get the current sizes of the perturbation factors
603	Number delta_x;
604	Number delta_s;
605	Number delta_c;
606	Number delta_d;
607	perturbHandler_->CurrentPerturbation(delta_x, delta_s, delta_c, delta_d);
608
609	SmartPtr<Vector> tmp;
610
611	// x
612	W.MultVector(1., res.x(), 0., resid.x_NonConst());
613	J_c.TransMultVector(1., res.y_c(), 1., resid.x_NonConst());
614	J_d.TransMultVector(1., res.y_d(), 1., resid.x_NonConst());
615	Px_L.MultVector(-1., res.z_L(), 1., resid.x_NonConst());
616	Px_U.MultVector(1., res.z_U(), 1., resid.x_NonConst());
617	resid.x_NonConst()->AddTwoVectors(delta_x, res.x(), -1., rhs.x(), 1.);
618
619	// s
620	Pd_U.MultVector(1., res.v_U(), 0., resid.s_NonConst());
621	Pd_L.MultVector(-1., res.v_L(), 1., resid.s_NonConst());
622	resid.s_NonConst()->AddTwoVectors(-1., res.y_d(), -1., rhs.s(), 1.);
623	if (delta_s!=0.) {
624	resid.s_NonConst()->Axpy(delta_s, *res.s());
625	}
626
627	// c
628	J_c.MultVector(1., res.x(), 0., resid.y_c_NonConst());
629	resid.y_c_NonConst()->AddTwoVectors(-delta_c, res.y_c(), -1., rhs.y_c(), 1.);
630
631	// d
632	J_d.MultVector(1., res.x(), 0., resid.y_d_NonConst());
633	resid.y_d_NonConst()->AddTwoVectors(-1., res.s(), -1., rhs.y_d(), 1.);
634	if (delta_d!=0.) {
635	resid.y_d_NonConst()->Axpy(-delta_d, *res.y_d());
636	}
637
638	// zL
639	resid.z_L_NonConst()->Copy(*res.z_L());
640	resid.z_L_NonConst()->ElementWiseMultiply(slack_x_L);
641	tmp = z_L.MakeNew();
642	Px_L.TransMultVector(1., res.x(), 0., tmp);
643	tmp->ElementWiseMultiply(z_L);
644	resid.z_L_NonConst()->AddTwoVectors(1., tmp, -1., rhs.z_L(), 1.);
645
646	// zU
647	resid.z_U_NonConst()->Copy(*res.z_U());
648	resid.z_U_NonConst()->ElementWiseMultiply(slack_x_U);
649	tmp = z_U.MakeNew();
650	Px_U.TransMultVector(1., res.x(), 0., tmp);
651	tmp->ElementWiseMultiply(z_U);
652	resid.z_U_NonConst()->AddTwoVectors(-1., tmp, -1., rhs.z_U(), 1.);
653
654	// vL
655	resid.v_L_NonConst()->Copy(*res.v_L());
656	resid.v_L_NonConst()->ElementWiseMultiply(slack_s_L);
657	tmp = v_L.MakeNew();
658	Pd_L.TransMultVector(1., res.s(), 0., tmp);
659	tmp->ElementWiseMultiply(v_L);
660	resid.v_L_NonConst()->AddTwoVectors(1., tmp, -1., rhs.v_L(), 1.);
661
662	// vU
663	resid.v_U_NonConst()->Copy(*res.v_U());
664	resid.v_U_NonConst()->ElementWiseMultiply(slack_s_U);
665	tmp = v_U.MakeNew();
666	Pd_U.TransMultVector(1., res.s(), 0., tmp);
667	tmp->ElementWiseMultiply(v_U);
668	resid.v_U_NonConst()->AddTwoVectors(-1., tmp, -1., rhs.v_U(), 1.);
669
670	DBG_PRINT_VECTOR(2, "resid", resid);
671
672	if (Jnlst().ProduceOutput(J_MOREVECTOR, J_LINEAR_ALGEBRA)) {
673	resid.Print(Jnlst(), J_MOREVECTOR, J_LINEAR_ALGEBRA, "resid");
674	}
675
676	if (Jnlst().ProduceOutput(J_MOREDETAILED, J_LINEAR_ALGEBRA)) {
677	Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA,
678	"max-norm resid_x %e\n", resid.x()->Amax());
679	Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA,
680	"max-norm resid_s %e\n", resid.s()->Amax());
681	Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA,
682	"max-norm resid_c %e\n", resid.y_c()->Amax());
683	Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA,
684	"max-norm resid_d %e\n", resid.y_d()->Amax());
685	Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA,
686	"max-norm resid_zL %e\n", resid.z_L()->Amax());
687	Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA,
688	"max-norm resid_zU %e\n", resid.z_U()->Amax());
689	Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA,
690	"max-norm resid_vL %e\n", resid.v_L()->Amax());
691	Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA,
692	"max-norm resid_vU %e\n", resid.v_U()->Amax());
693	}
694	IpData().TimingStats().ComputeResiduals().End();
695	}
696
697	Number PDFullSpaceSolver::ComputeResidualRatio(const IteratesVector& rhs,
698	const IteratesVector& res,
699	const IteratesVector& resid)
700	{
701	DBG_START_METH("PDFullSpaceSolver::ComputeResidualRatio", dbg_verbosity);
702
703	Number nrm_rhs = rhs.Amax();
704	Number nrm_res = res.Amax();
705	Number nrm_resid = resid.Amax();
706	Jnlst().Printf(J_MOREDETAILED, J_LINEAR_ALGEBRA,
707	"nrm_rhs = %8.2e nrm_sol = %8.2e nrm_resid = %8.2e\n",
708	nrm_rhs, nrm_res, nrm_resid);
709
710	if (nrm_rhs+nrm_res == 0.) {
711	return nrm_resid; // this should be zero
712	}
713	else {
714	// ToDo: determine how to include norm of matrix, and what
715	// safeguard to use against incredibly large solution vectors
716	Number max_cond = 1e6;
717	return nrm_resid/(Min(nrm_res, max_cond*nrm_rhs)+nrm_rhs);
718	}
719	}
720
721	} // namespace Ipopt

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

Download in other formats: