Mercurial > octave
annotate libinterp/corefcn/__qp__.cc @ 31225:3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
author | Arun Giridhar <arungiridhar@gmail.com> |
---|---|
date | Sun, 11 Sep 2022 13:53:38 -0400 |
parents | 796f54d4ddbf |
children | e88a07dec498 |
rev | line source |
---|---|
27923
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
1 //////////////////////////////////////////////////////////////////////// |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
2 // |
30564
796f54d4ddbf
update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents:
30390
diff
changeset
|
3 // Copyright (C) 2000-2022 The Octave Project Developers |
27923
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
4 // |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
5 // See the file COPYRIGHT.md in the top-level directory of this |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
6 // distribution or <https://octave.org/copyright/>. |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
7 // |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
8 // This file is part of Octave. |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
9 // |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
10 // Octave is free software: you can redistribute it and/or modify it |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
11 // under the terms of the GNU General Public License as published by |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
12 // the Free Software Foundation, either version 3 of the License, or |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
13 // (at your option) any later version. |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
14 // |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
15 // Octave is distributed in the hope that it will be useful, but |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
16 // WITHOUT ANY WARRANTY; without even the implied warranty of |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
18 // GNU General Public License for more details. |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
19 // |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
20 // You should have received a copy of the GNU General Public License |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
21 // along with Octave; see the file COPYING. If not, see |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
22 // <https://www.gnu.org/licenses/>. |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
23 // |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
24 //////////////////////////////////////////////////////////////////////// |
5290 | 25 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21301
diff
changeset
|
26 #if defined (HAVE_CONFIG_H) |
21301
40de9f8f23a6
Use '#include "config.h"' rather than <config.h>.
Rik <rik@octave.org>
parents:
21273
diff
changeset
|
27 # include "config.h" |
5293 | 28 #endif |
29 | |
23455
73ff72d3d603
maint: Eliminate <cfloat.h> header from libinterp files
Rik <rik@octave.org>
parents:
23450
diff
changeset
|
30 #include <cmath> |
73ff72d3d603
maint: Eliminate <cfloat.h> header from libinterp files
Rik <rik@octave.org>
parents:
23450
diff
changeset
|
31 |
73ff72d3d603
maint: Eliminate <cfloat.h> header from libinterp files
Rik <rik@octave.org>
parents:
23450
diff
changeset
|
32 #include <limits> |
5290 | 33 |
21269
3c8a3d35661a
better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents:
21200
diff
changeset
|
34 #include "chol.h" |
21273
cbced1c09916
better use of templates for svd classes
John W. Eaton <jwe@octave.org>
parents:
21269
diff
changeset
|
35 #include "svd.h" |
5293 | 36 #include "mx-m-dm.h" |
37 #include "EIG.h" | |
5290 | 38 |
15039
e753177cde93
maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory
Rik <rik@octave.org>
parents:
14501
diff
changeset
|
39 #include "defun.h" |
5293 | 40 #include "error.h" |
21100
e39e05d90788
Switch gripe_XXX to either err_XXX or warn_XXX naming scheme.
Rik <rik@octave.org>
parents:
20940
diff
changeset
|
41 #include "errwarn.h" |
20940
48b2ad5ee801
maint: Rename oct-obj.[cc|h] to ovl.[cc|h] for clarity.
Rik <rik@octave.org>
parents:
20892
diff
changeset
|
42 #include "ovl.h" |
5293 | 43 #include "pr-output.h" |
44 #include "utils.h" | |
5290 | 45 |
29958
32c3a5805893
move DEFUN and DEFMETHOD functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
29359
diff
changeset
|
46 OCTAVE_NAMESPACE_BEGIN |
32c3a5805893
move DEFUN and DEFMETHOD functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
29359
diff
changeset
|
47 |
27030
4a2cb3392014
refactor minimum eigenvalue index search in qp (bug #56037)
Mike Miller <mtmiller@octave.org>
parents:
27027
diff
changeset
|
48 static octave_idx_type |
4a2cb3392014
refactor minimum eigenvalue index search in qp (bug #56037)
Mike Miller <mtmiller@octave.org>
parents:
27027
diff
changeset
|
49 min_index (const ColumnVector& x) |
4a2cb3392014
refactor minimum eigenvalue index search in qp (bug #56037)
Mike Miller <mtmiller@octave.org>
parents:
27027
diff
changeset
|
50 { |
4a2cb3392014
refactor minimum eigenvalue index search in qp (bug #56037)
Mike Miller <mtmiller@octave.org>
parents:
27027
diff
changeset
|
51 double min_val = x.min (); |
4a2cb3392014
refactor minimum eigenvalue index search in qp (bug #56037)
Mike Miller <mtmiller@octave.org>
parents:
27027
diff
changeset
|
52 |
4a2cb3392014
refactor minimum eigenvalue index search in qp (bug #56037)
Mike Miller <mtmiller@octave.org>
parents:
27027
diff
changeset
|
53 for (octave_idx_type i = 0; i < x.numel (); i++) |
4a2cb3392014
refactor minimum eigenvalue index search in qp (bug #56037)
Mike Miller <mtmiller@octave.org>
parents:
27027
diff
changeset
|
54 if (min_val == x.xelem (i)) |
4a2cb3392014
refactor minimum eigenvalue index search in qp (bug #56037)
Mike Miller <mtmiller@octave.org>
parents:
27027
diff
changeset
|
55 return i; |
4a2cb3392014
refactor minimum eigenvalue index search in qp (bug #56037)
Mike Miller <mtmiller@octave.org>
parents:
27027
diff
changeset
|
56 |
4a2cb3392014
refactor minimum eigenvalue index search in qp (bug #56037)
Mike Miller <mtmiller@octave.org>
parents:
27027
diff
changeset
|
57 return 0; |
4a2cb3392014
refactor minimum eigenvalue index search in qp (bug #56037)
Mike Miller <mtmiller@octave.org>
parents:
27027
diff
changeset
|
58 } |
4a2cb3392014
refactor minimum eigenvalue index search in qp (bug #56037)
Mike Miller <mtmiller@octave.org>
parents:
27027
diff
changeset
|
59 |
5290 | 60 static Matrix |
5295 | 61 null (const Matrix& A, octave_idx_type& rank) |
5290 | 62 { |
63 Matrix retval; | |
64 | |
65 rank = 0; | |
66 | |
23577
80c42f4cca13
maint: Deprecate is_empty and replace with isempty.
Rik <rik@octave.org>
parents:
23455
diff
changeset
|
67 if (! A.isempty ()) |
5290 | 68 { |
29961
7d6709900da7
eliminate octave:: namespace tags in DEFUN and DEFMETHOD and more
John W. Eaton <jwe@octave.org>
parents:
29958
diff
changeset
|
69 math::svd<Matrix> A_svd (A); |
5290 | 70 |
71 DiagMatrix S = A_svd.singular_values (); | |
72 | |
15448
0a0912a9ab6e
Replace deprecated DiagArray2<T>::diag calls with DiagArray2<T>::extract_diag
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
15220
diff
changeset
|
73 ColumnVector s = S.extract_diag (); |
5290 | 74 |
75 Matrix V = A_svd.right_singular_matrix (); | |
76 | |
5295 | 77 octave_idx_type A_nr = A.rows (); |
78 octave_idx_type A_nc = A.cols (); | |
5290 | 79 |
23450
855122b993da
maint: Wrap tertiary operator in parentheses "(COND ? x : y)".
Rik <rik@octave.org>
parents:
23220
diff
changeset
|
80 octave_idx_type tmp = (A_nr > A_nc ? A_nr : A_nc); |
5290 | 81 |
15220
61822c866ba1
use std::numeric_limits<T>::epsilon in C++ code
John W. Eaton <jwe@octave.org>
parents:
15195
diff
changeset
|
82 double tol = tmp * s(0) * std::numeric_limits<double>::epsilon (); |
5290 | 83 |
20232
a9574e3c6e9e
Deprecate Array::length() and Sparse::length() in favour of ::numel().
Carnë Draug <carandraug@octave.org>
parents:
19697
diff
changeset
|
84 octave_idx_type n = s.numel (); |
5290 | 85 |
5295 | 86 for (octave_idx_type i = 0; i < n; i++) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
87 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
88 if (s(i) > tol) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
89 rank++; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
90 } |
5290 | 91 |
92 if (rank < A_nc) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
93 retval = V.extract (0, rank, A_nc-1, A_nc-1); |
5290 | 94 else |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
95 retval.resize (A_nc, 0); |
6431 | 96 |
97 for (octave_idx_type i = 0; i < retval.numel (); i++) | |
15220
61822c866ba1
use std::numeric_limits<T>::epsilon in C++ code
John W. Eaton <jwe@octave.org>
parents:
15195
diff
changeset
|
98 if (std::abs (retval(i)) < std::numeric_limits<double>::epsilon ()) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
99 retval(i) = 0; |
5290 | 100 } |
101 | |
102 return retval; | |
103 } | |
104 | |
105 static int | |
106 qp (const Matrix& H, const ColumnVector& q, | |
107 const Matrix& Aeq, const ColumnVector& beq, | |
108 const Matrix& Ain, const ColumnVector& bin, | |
26332
d4211810202a
qp and sqp: Non-fixed tolerance for qp (bug #53506).
Maor Shutman <maorus12@gmail.com>
parents:
25054
diff
changeset
|
109 int maxit, double rtol, |
5290 | 110 ColumnVector& x, ColumnVector& lambda, int& iter) |
111 { | |
112 int info = 0; | |
113 | |
114 iter = 0; | |
115 | |
116 // Problem dimension. | |
20232
a9574e3c6e9e
Deprecate Array::length() and Sparse::length() in favour of ::numel().
Carnë Draug <carandraug@octave.org>
parents:
19697
diff
changeset
|
117 octave_idx_type n = x.numel (); |
5290 | 118 |
119 // Dimension of constraints. | |
20232
a9574e3c6e9e
Deprecate Array::length() and Sparse::length() in favour of ::numel().
Carnë Draug <carandraug@octave.org>
parents:
19697
diff
changeset
|
120 octave_idx_type n_eq = beq.numel (); |
a9574e3c6e9e
Deprecate Array::length() and Sparse::length() in favour of ::numel().
Carnë Draug <carandraug@octave.org>
parents:
19697
diff
changeset
|
121 octave_idx_type n_in = bin.numel (); |
5290 | 122 |
123 // Filling the current active set. | |
124 | |
5295 | 125 octave_idx_type n_act = n_eq; |
5290 | 126 |
5295 | 127 octave_idx_type n_tot = n_eq + n_in; |
5290 | 128 |
129 // Equality constraints come first. We won't check the sign of the | |
130 // Lagrange multiplier for those. | |
131 | |
132 Matrix Aact = Aeq; | |
133 ColumnVector bact = beq; | |
134 ColumnVector Wact; | |
135 | |
136 if (n_in > 0) | |
137 { | |
138 ColumnVector res = Ain*x - bin; | |
139 | |
5295 | 140 for (octave_idx_type i = 0; i < n_in; i++) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
141 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
142 res(i) /= (1.0 + std::abs (bin(i))); |
5290 | 143 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
144 if (res(i) < rtol) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
145 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
146 n_act++; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
147 Aact = Aact.stack (Ain.row (i)); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
148 bact.resize (n_act, bin(i)); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
149 Wact.resize (n_act-n_eq, i); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
150 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
151 } |
5290 | 152 } |
153 | |
154 // Computing the ??? | |
155 | |
20741
a5ab31b52ae8
eliminate more uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20678
diff
changeset
|
156 EIG eigH; |
7578
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7361
diff
changeset
|
157 |
20741
a5ab31b52ae8
eliminate more uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20678
diff
changeset
|
158 try |
a5ab31b52ae8
eliminate more uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20678
diff
changeset
|
159 { |
a5ab31b52ae8
eliminate more uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20678
diff
changeset
|
160 eigH = EIG (H); |
a5ab31b52ae8
eliminate more uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20678
diff
changeset
|
161 } |
29961
7d6709900da7
eliminate octave:: namespace tags in DEFUN and DEFMETHOD and more
John W. Eaton <jwe@octave.org>
parents:
29958
diff
changeset
|
162 catch (execution_exception& ee) |
7578
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7361
diff
changeset
|
163 { |
29163
8f67ad8b3103
maint: Updating naming conventions for exceptions and use const where possible.
Rik <rik@octave.org>
parents:
28662
diff
changeset
|
164 error (ee, "qp: failed to compute eigenvalues of H"); |
7578
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7361
diff
changeset
|
165 } |
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7361
diff
changeset
|
166 |
5290 | 167 ColumnVector eigenvalH = real (eigH.eigenvalues ()); |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
21966
diff
changeset
|
168 Matrix eigenvecH = real (eigH.right_eigenvectors ()); |
27030
4a2cb3392014
refactor minimum eigenvalue index search in qp (bug #56037)
Mike Miller <mtmiller@octave.org>
parents:
27027
diff
changeset
|
169 octave_idx_type indminR = min_index (eigenvalH); |
5290 | 170 |
171 bool done = false; | |
172 | |
173 double alpha = 0.0; | |
174 | |
175 Matrix R; | |
176 Matrix Y (n, 0, 0.0); | |
177 | |
178 ColumnVector g (n, 0.0); | |
179 ColumnVector p (n, 0.0); | |
180 | |
181 ColumnVector lambda_tmp (n_in, 0.0); | |
182 | |
183 while (! done) | |
184 { | |
185 iter++; | |
186 | |
187 // Current Gradient | |
188 // g = q + H * x; | |
189 | |
190 g = q + H * x; | |
191 | |
192 if (n_act == 0) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
193 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
194 // There are no active constraints. |
5290 | 195 |
27030
4a2cb3392014
refactor minimum eigenvalue index search in qp (bug #56037)
Mike Miller <mtmiller@octave.org>
parents:
27027
diff
changeset
|
196 double minReal = eigenvalH.xelem (indminR); |
4a2cb3392014
refactor minimum eigenvalue index search in qp (bug #56037)
Mike Miller <mtmiller@octave.org>
parents:
27027
diff
changeset
|
197 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
198 if (minReal > 0.0) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
199 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
200 // Inverting the Hessian. Using the Cholesky |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
201 // factorization since the Hessian is positive |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
202 // definite. |
5290 | 203 |
29961
7d6709900da7
eliminate octave:: namespace tags in DEFUN and DEFMETHOD and more
John W. Eaton <jwe@octave.org>
parents:
29958
diff
changeset
|
204 math::chol<Matrix> cholH (H); |
5290 | 205 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
206 R = cholH.chol_matrix (); |
5290 | 207 |
29961
7d6709900da7
eliminate octave:: namespace tags in DEFUN and DEFMETHOD and more
John W. Eaton <jwe@octave.org>
parents:
29958
diff
changeset
|
208 Matrix Hinv = math::chol2inv (R); |
5290 | 209 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
210 // Computing the unconstrained step. |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
211 // p = -Hinv * g; |
5290 | 212 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
213 p = -Hinv * g; |
5290 | 214 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
215 info = 0; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
216 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
217 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
218 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
219 // Finding the negative curvature of H. |
5290 | 220 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
221 p = eigenvecH.column (indminR); |
5290 | 222 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
223 // Following the negative curvature of H. |
5290 | 224 |
15220
61822c866ba1
use std::numeric_limits<T>::epsilon in C++ code
John W. Eaton <jwe@octave.org>
parents:
15195
diff
changeset
|
225 if (p.transpose () * g > std::numeric_limits<double>::epsilon ()) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
226 p = -p; |
5290 | 227 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
228 info = 1; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
229 } |
5290 | 230 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
231 // Multipliers are zero. |
5290 | 232 lambda_tmp.fill (0.0); |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
233 } |
5290 | 234 else |
235 { | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
236 // There are active constraints. |
5290 | 237 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
238 // Computing the null space. |
5290 | 239 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
240 octave_idx_type rank; |
5290 | 241 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
242 Matrix Z = null (Aact, rank); |
5290 | 243 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
244 octave_idx_type dimZ = n - rank; |
5290 | 245 |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
246 // FIXME: still remain to handle the case of |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
247 // non-full rank active set matrix. |
5290 | 248 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
249 // Computing the Y matrix (orthogonal to Z) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
250 Y = Aact.pseudo_inverse (); |
7361 | 251 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
252 // Reduced Hessian |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
253 Matrix Zt = Z.transpose (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
254 Matrix rH = Zt * H * Z; |
5290 | 255 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
256 octave_idx_type pR = 0; |
5290 | 257 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
258 if (dimZ > 0) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
259 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
260 // Computing the Cholesky factorization (pR = 0 means |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
261 // that the reduced Hessian was positive definite). |
5290 | 262 |
29961
7d6709900da7
eliminate octave:: namespace tags in DEFUN and DEFMETHOD and more
John W. Eaton <jwe@octave.org>
parents:
29958
diff
changeset
|
263 math::chol<Matrix> cholrH (rH, pR); |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
264 Matrix tR = cholrH.chol_matrix (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
265 if (pR == 0) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
266 R = tR; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
267 } |
5290 | 268 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
269 if (pR == 0) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
270 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
271 info = 0; |
5290 | 272 |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
273 // Computing the step pz. |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
274 if (dimZ > 0) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
275 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
276 // Using the Cholesky factorization to invert rH |
5290 | 277 |
29961
7d6709900da7
eliminate octave:: namespace tags in DEFUN and DEFMETHOD and more
John W. Eaton <jwe@octave.org>
parents:
29958
diff
changeset
|
278 Matrix rHinv = math::chol2inv (R); |
5290 | 279 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
280 ColumnVector pz = -rHinv * Zt * g; |
5290 | 281 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
282 // Global step. |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
283 p = Z * pz; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
284 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
285 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
286 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
287 // Global step. |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
288 p.fill (0.0); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
289 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
290 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
291 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
292 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
293 info = 1; |
5290 | 294 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
295 // Searching for the most negative curvature. |
5290 | 296 |
20741
a5ab31b52ae8
eliminate more uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20678
diff
changeset
|
297 EIG eigrH; |
7578
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7361
diff
changeset
|
298 |
20741
a5ab31b52ae8
eliminate more uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20678
diff
changeset
|
299 try |
a5ab31b52ae8
eliminate more uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20678
diff
changeset
|
300 { |
a5ab31b52ae8
eliminate more uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20678
diff
changeset
|
301 eigrH = EIG (rH); |
a5ab31b52ae8
eliminate more uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20678
diff
changeset
|
302 } |
29961
7d6709900da7
eliminate octave:: namespace tags in DEFUN and DEFMETHOD and more
John W. Eaton <jwe@octave.org>
parents:
29958
diff
changeset
|
303 catch (execution_exception& ee) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
304 { |
29163
8f67ad8b3103
maint: Updating naming conventions for exceptions and use const where possible.
Rik <rik@octave.org>
parents:
28662
diff
changeset
|
305 error (ee, "qp: failed to compute eigenvalues of rH"); |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
306 } |
7578
91f8446ce4ae
handle possible error from EIG
John W. Eaton <jwe@octave.org>
parents:
7361
diff
changeset
|
307 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
308 ColumnVector eigenvalrH = real (eigrH.eigenvalues ()); |
22305
510886d03ef2
eig: new options for choice of algorithm, balancing, and output (patch #8960)
Barbara Locsi <locsi.barbara@gmail.com>
parents:
21966
diff
changeset
|
309 Matrix eigenvecrH = real (eigrH.right_eigenvectors ()); |
27030
4a2cb3392014
refactor minimum eigenvalue index search in qp (bug #56037)
Mike Miller <mtmiller@octave.org>
parents:
27027
diff
changeset
|
310 indminR = min_index (eigenvalrH); |
5290 | 311 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
312 ColumnVector eVrH = eigenvecrH.column (indminR); |
5290 | 313 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
314 // Computing the step pz. |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
315 p = Z * eVrH; |
5290 | 316 |
15220
61822c866ba1
use std::numeric_limits<T>::epsilon in C++ code
John W. Eaton <jwe@octave.org>
parents:
15195
diff
changeset
|
317 if (p.transpose () * g > std::numeric_limits<double>::epsilon ()) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
318 p = -p; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
319 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
320 } |
5290 | 321 |
322 // Checking the step-size. | |
323 ColumnVector abs_p (n); | |
5295 | 324 for (octave_idx_type i = 0; i < n; i++) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
325 abs_p(i) = std::abs (p(i)); |
5290 | 326 double max_p = abs_p.max (); |
327 | |
328 if (max_p < rtol) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
329 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
330 // The step is null. Checking constraints. |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
331 if (n_act - n_eq == 0) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
332 // Solution is found because no inequality |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
333 // constraints are active. |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
334 done = true; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
335 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
336 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
337 // Computing the multipliers only for the inequality |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
338 // constraints that are active. We do NOT compute |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
339 // multipliers for the equality constraints. |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
340 Matrix Yt = Y.transpose (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
341 Yt = Yt.extract_n (n_eq, 0, n_act-n_eq, n); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
342 lambda_tmp = Yt * (g + H * p); |
5290 | 343 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
344 // Checking the multipliers. We remove the most |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
345 // negative from the set (if any). |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
346 double min_lambda = lambda_tmp.min (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
347 if (min_lambda >= 0) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
348 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
349 // Solution is found. |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
350 done = true; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
351 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
352 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
353 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
354 octave_idx_type which_eig = 0; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
355 for (octave_idx_type i = 0; i < n_act; i++) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
356 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
357 if (lambda_tmp(i) == min_lambda) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
358 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
359 which_eig = i; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
360 break; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
361 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
362 } |
5290 | 363 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
364 // At least one multiplier is negative, we |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
365 // remove it from the set. |
5290 | 366 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
367 n_act--; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
368 for (octave_idx_type i = which_eig; i < n_act - n_eq; i++) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
369 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
370 Wact(i) = Wact(i+1); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
371 for (octave_idx_type j = 0; j < n; j++) |
30390
a61e1a0f6024
maint: style check C++ files in libinterp/ ahead of 7.1 release.
Rik <rik@octave.org>
parents:
29961
diff
changeset
|
372 Aact(n_eq+i, j) = Aact(n_eq+i+1, j); |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
373 bact(n_eq+i) = bact(n_eq+i+1); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
374 } |
5290 | 375 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
376 // Resizing the active set. |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
377 Wact.resize (n_act-n_eq); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
378 bact.resize (n_act); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
379 Aact.resize (n_act, n); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
380 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
381 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
382 } |
5290 | 383 else |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
384 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
385 // The step is not null. |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
386 if (n_act - n_eq == n_in) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
387 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
388 // All inequality constraints were active. We can |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
389 // add the whole step. |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
390 x += p; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
391 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
392 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
393 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
394 // Some constraints were not active. Checking if |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
395 // there is a blocking constraint. |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
396 alpha = 1.0; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
397 octave_idx_type is_block = -1; |
5290 | 398 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
399 for (octave_idx_type i = 0; i < n_in; i++) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
400 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
401 bool found = false; |
5290 | 402 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
403 for (octave_idx_type j = 0; j < n_act-n_eq; j++) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
404 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
405 if (Wact(j) == i) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
406 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
407 found = true; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
408 break; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
409 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
410 } |
5290 | 411 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
412 if (! found) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
413 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
414 // The i-th constraint was not in the set. Is it a |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
415 // blocking constraint? |
5290 | 416 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
417 RowVector tmp_row = Ain.row (i); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
418 double tmp = tmp_row * p; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
419 double res = tmp_row * x; |
5290 | 420 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
421 if (tmp < 0.0) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
422 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
423 double alpha_tmp = (bin(i) - res) / tmp; |
5290 | 424 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
425 if (alpha_tmp < alpha) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
426 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
427 alpha = alpha_tmp; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
428 is_block = i; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
429 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
430 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
431 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
432 } |
5290 | 433 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
434 // In is_block there is the index of the blocking |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
435 // constraint (if any). |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
436 if (is_block >= 0) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
437 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
438 // There is a blocking constraint (index in |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
439 // is_block) which is added to the active set. |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
440 n_act++; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
441 Aact = Aact.stack (Ain.row (is_block)); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
442 bact.resize (n_act, bin(is_block)); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
443 Wact.resize (n_act-n_eq, is_block); |
5290 | 444 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
445 // Adding the reduced step |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
446 x += alpha * p; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
447 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
448 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
449 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
450 // There are no blocking constraints. Adding the |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
451 // whole step. |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
452 x += alpha * p; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
453 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
454 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
455 } |
5290 | 456 |
457 if (iter == maxit) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
458 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
459 done = true; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
460 // warning ("qp_main: maximum number of iteration reached"); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
461 info = 3; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
462 } |
5290 | 463 } |
464 | |
465 lambda_tmp = Y.transpose () * (g + H * p); | |
466 | |
467 // Reordering the Lagrange multipliers. | |
468 | |
469 lambda.resize (n_tot); | |
470 lambda.fill (0.0); | |
5295 | 471 for (octave_idx_type i = 0; i < n_eq; i++) |
5290 | 472 lambda(i) = lambda_tmp(i); |
473 | |
5295 | 474 for (octave_idx_type i = n_eq; i < n_tot; i++) |
5290 | 475 { |
5295 | 476 for (octave_idx_type j = 0; j < n_act-n_eq; j++) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
477 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
478 if (Wact(j) == i - n_eq) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
479 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
480 lambda(i) = lambda_tmp(n_eq+j); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
481 break; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
482 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
7729
diff
changeset
|
483 } |
5290 | 484 } |
485 | |
486 return info; | |
487 } | |
488 | |
15039
e753177cde93
maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory
Rik <rik@octave.org>
parents:
14501
diff
changeset
|
489 DEFUN (__qp__, args, , |
21966
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21724
diff
changeset
|
490 doc: /* -*- texinfo -*- |
28662
bf941bb59f72
* __qp__.cc (F__qp__): Mention all input arguments in stub description.
Markus Mützel <markus.muetzel@gmx.de>
parents:
27923
diff
changeset
|
491 @deftypefn {} {[@var{x}, @var{lambda}, @var{info}, @var{iter}] =} __qp__ (@var{x0}, @var{H}, @var{q}, @var{Aeq}, @var{beq}, @var{Ain}, @var{bin}, @var{maxit}, @var{rtol}) |
21966
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21724
diff
changeset
|
492 Undocumented internal function. |
112b20240c87
move docstrings in C++ files out of C strings and into comments
John W. Eaton <jwe@octave.org>
parents:
21724
diff
changeset
|
493 @end deftypefn */) |
5290 | 494 { |
26332
d4211810202a
qp and sqp: Non-fixed tolerance for qp (bug #53506).
Maor Shutman <maorus12@gmail.com>
parents:
25054
diff
changeset
|
495 if (args.length () != 9) |
20802
8bb38ba1bad6
eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents:
20785
diff
changeset
|
496 print_usage (); |
5290 | 497 |
20892 | 498 const ColumnVector x0 (args(0).vector_value ()); |
499 const Matrix H (args(1).matrix_value ()); | |
500 const ColumnVector q (args(2).vector_value ()); | |
501 const Matrix Aeq (args(3).matrix_value ()); | |
502 const ColumnVector beq (args(4).vector_value ()); | |
503 const Matrix Ain (args(5).matrix_value ()); | |
504 const ColumnVector bin (args(6).vector_value ()); | |
505 const int maxit (args(7).int_value ()); | |
26332
d4211810202a
qp and sqp: Non-fixed tolerance for qp (bug #53506).
Maor Shutman <maorus12@gmail.com>
parents:
25054
diff
changeset
|
506 const double rtol (args(8).double_value()); |
5290 | 507 |
20802
8bb38ba1bad6
eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents:
20785
diff
changeset
|
508 int iter = 0; |
5290 | 509 |
20892 | 510 // Copy the initial guess into the working variable |
20802
8bb38ba1bad6
eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents:
20785
diff
changeset
|
511 ColumnVector x = x0; |
20678
4b00afb5e9c3
eliminate more uses of error_state
John W. Eaton <jwe@octave.org>
parents:
20232
diff
changeset
|
512 |
20802
8bb38ba1bad6
eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents:
20785
diff
changeset
|
513 // Reordering the Lagrange multipliers |
8bb38ba1bad6
eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents:
20785
diff
changeset
|
514 ColumnVector lambda; |
8bb38ba1bad6
eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents:
20785
diff
changeset
|
515 |
26332
d4211810202a
qp and sqp: Non-fixed tolerance for qp (bug #53506).
Maor Shutman <maorus12@gmail.com>
parents:
25054
diff
changeset
|
516 int info = qp (H, q, Aeq, beq, Ain, bin, maxit, rtol, x, lambda, iter); |
20802
8bb38ba1bad6
eliminate return statements after calls to print_usage
John W. Eaton <jwe@octave.org>
parents:
20785
diff
changeset
|
517 |
20892 | 518 return ovl (x, lambda, info, iter); |
5290 | 519 } |
12805
3641167e5b75
codesprint: *.cc helper functions do not need tests
Rik <octave@nomad.inbox5.com>
parents:
11586
diff
changeset
|
520 |
3641167e5b75
codesprint: *.cc helper functions do not need tests
Rik <octave@nomad.inbox5.com>
parents:
11586
diff
changeset
|
521 /* |
3641167e5b75
codesprint: *.cc helper functions do not need tests
Rik <octave@nomad.inbox5.com>
parents:
11586
diff
changeset
|
522 ## No test needed for internal helper function. |
3641167e5b75
codesprint: *.cc helper functions do not need tests
Rik <octave@nomad.inbox5.com>
parents:
11586
diff
changeset
|
523 %!assert (1) |
3641167e5b75
codesprint: *.cc helper functions do not need tests
Rik <octave@nomad.inbox5.com>
parents:
11586
diff
changeset
|
524 */ |
29958
32c3a5805893
move DEFUN and DEFMETHOD functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
29359
diff
changeset
|
525 |
32c3a5805893
move DEFUN and DEFMETHOD functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
29359
diff
changeset
|
526 OCTAVE_NAMESPACE_END |