1+ #include " Solver.h"
2+
3+ // Iterative refinement on the 6x6 system
4+
5+ namespace hipo {
6+
7+ void Solver::refine (NewtonDir& delta) {
8+ // compute the residuals of the linear system in it_->ires
9+ it_->residuals6x6 (delta);
10+
11+ double omega = computeOmega (delta);
12+ }
13+
14+ static void updateOmega (double tau, double & omega1, double & omega2, double num,
15+ double den1, double den2, double den3) {
16+ if (den1 + den2 > tau) {
17+ double omega = num / (den1 + den2);
18+ omega1 = std::max (omega1, omega);
19+ } else {
20+ double omega = num / (den1 + den3);
21+ omega2 = std::max (omega2, omega);
22+ }
23+ }
24+
25+ double Solver::computeOmega (const NewtonDir& delta) const {
26+ // Evaluate iterative refinement progress. Based on "Solving sparse linear
27+ // systems with sparse backward error", Arioli, Demmel, Duff.
28+
29+ // Compute |A| * |dx| and |A^T| * |dy|
30+ std::vector<double > abs_prod_A (m_);
31+ std::vector<double > abs_prod_At (n_);
32+ for (Int col = 0 ; col < n_; ++col) {
33+ for (Int el = model_.A ().start_ [col]; el < model_.A ().start_ [col + 1 ];
34+ ++el) {
35+ Int row = model_.A ().index_ [el];
36+ double val = model_.A ().value_ [el];
37+ abs_prod_A[row] += std::abs (val) * std::abs (delta.x [col]);
38+ abs_prod_At[col] += std::abs (val) * std::abs (delta.y [row]);
39+ }
40+ }
41+
42+ double inf_norm_delta{};
43+ inf_norm_delta = std::max (inf_norm_delta, infNorm (delta.x ));
44+ inf_norm_delta = std::max (inf_norm_delta, infNorm (delta.xl ));
45+ inf_norm_delta = std::max (inf_norm_delta, infNorm (delta.xu ));
46+ inf_norm_delta = std::max (inf_norm_delta, infNorm (delta.y ));
47+ inf_norm_delta = std::max (inf_norm_delta, infNorm (delta.zl ));
48+ inf_norm_delta = std::max (inf_norm_delta, infNorm (delta.zu ));
49+
50+ // rhs is in it_->res
51+ // residual of linear system is in it_->ires
52+
53+ // tau_i =
54+ // tau_const * (inf_norm_rows(big 6x6 matrix)_i * inf_norm_delta + |rhs_i|)
55+ const double tau_const = 1000.0 * (5 * n_ + m_) * 1e-16 ;
56+ double omega_1{}, omega_2{};
57+
58+ // First block
59+ for (Int i = 0 ; i < m_; ++i) {
60+ const double tau = tau_const * (model_.infNormRows (i) * inf_norm_delta +
61+ std::abs (it_->res .r1 [i]));
62+ updateOmega (tau, omega_1, omega_2,
63+ // numerator
64+ std::abs (it_->ires .r1 [i]),
65+ // denominators
66+ abs_prod_A[i], std::abs (it_->res .r1 [i]),
67+ model_.oneNormRows (i) * inf_norm_delta);
68+ }
69+
70+ // Second block
71+ for (Int i = 0 ; i < n_; ++i) {
72+ if (model_.hasLb (i)) {
73+ const double tau =
74+ tau_const * (1.0 * inf_norm_delta + std::abs (it_->res .r2 [i]));
75+ updateOmega (tau, omega_1, omega_2,
76+ // numerator
77+ std::abs (it_->ires .r2 [i]),
78+ // denominators
79+ std::abs (delta.x [i]) + std::abs (delta.xl [i]),
80+ std::abs (it_->res .r2 [i]), 2 * inf_norm_delta
81+
82+ );
83+ }
84+ }
85+
86+ // Third block
87+ for (Int i = 0 ; i < n_; ++i) {
88+ if (model_.hasUb (i)) {
89+ const double tau =
90+ tau_const * (1.0 * inf_norm_delta + std::abs (it_->res .r3 [i]));
91+ updateOmega (tau, omega_1, omega_2,
92+ // numerator
93+ std::abs (it_->ires .r3 [i]),
94+ // denominators
95+ std::abs (delta.x [i]) + std::abs (delta.xu [i]),
96+ std::abs (it_->res .r3 [i]), 2 * inf_norm_delta);
97+ }
98+ }
99+
100+ // Fourth block
101+ for (Int i = 0 ; i < n_; ++i) {
102+ const double tau =
103+ tau_const * (std::max (model_.infNormCols (i), 1.0 ) * inf_norm_delta +
104+ std::abs (it_->res .r4 [i]));
105+
106+ double denom1 = abs_prod_At[i];
107+ if (model_.hasLb (i)) denom1 += std::abs (delta.zl [i]);
108+ if (model_.hasUb (i)) denom1 += std::abs (delta.zu [i]);
109+
110+ updateOmega (tau, omega_1, omega_2,
111+ // numerator
112+ std::abs (it_->ires .r4 [i]),
113+ // denominators
114+ denom1, std::abs (it_->res .r4 [i]),
115+ (model_.oneNormCols (i) + 2 ) * inf_norm_delta);
116+ }
117+
118+ // Fifth block
119+ for (Int i = 0 ; i < n_; ++i) {
120+ if (model_.hasLb (i)) {
121+ const double tau =
122+ tau_const * (std::max (std::abs (it_->xl [i]), std::abs (it_->zl [i])) *
123+ inf_norm_delta +
124+ std::abs (it_->res .r5 [i]));
125+
126+ updateOmega (
127+ tau, omega_1, omega_2,
128+ // numerator
129+ std::abs (it_->ires .r5 [i]),
130+ // denominators
131+ std::abs (it_->zl [i]) * std::abs (delta.xl [i]) +
132+ std::abs (it_->xl [i]) * std::abs (delta.zl [i]),
133+ std::abs (it_->res .r5 [i]),
134+ (std::abs (it_->zl [i]) + std::abs (it_->xl [i])) * inf_norm_delta);
135+ }
136+ }
137+
138+ // Sixth block
139+ for (Int i = 0 ; i < n_; ++i) {
140+ if (model_.hasUb (i)) {
141+ const double tau =
142+ tau_const * (std::max (std::abs (it_->xu [i]), std::abs (it_->zu [i])) *
143+ inf_norm_delta +
144+ std::abs (it_->res .r6 [i]));
145+
146+ updateOmega (
147+ tau, omega_1, omega_2,
148+ // numerator
149+ std::abs (it_->ires .r6 [i]),
150+ // denominators
151+ std::abs (it_->zu [i]) * std::abs (delta.xu [i]) +
152+ std::abs (it_->xu [i]) * std::abs (delta.zu [i]),
153+ std::abs (it_->res .r6 [i]),
154+ (std::abs (it_->zu [i]) + std::abs (it_->xu [i])) * inf_norm_delta);
155+ }
156+ }
157+
158+ return omega_1 + omega_2;
159+ }
160+
161+ } // namespace hipo
0 commit comments