Mercurial > octave

annotate libinterp/corefcn/__qp__.cc @ 31235:b542b88ad3b6

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Use space between function name and '(' in sparse-xpow.cc. * sparse-xpow.cc: Use space between function name and '(' in sparse-xpow.cc.
author Rik <rik@octave.org>
date Tue, 20 Sep 2022 14:43:56 -0700
parents 796f54d4ddbf
children e88a07dec498
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
7a51b27574bc [project @ 2005-04-21 16:03:29 by jwe]
jwe
parents: 5290
diff changeset
28 #endif
7a51b27574bc [project @ 2005-04-21 16:03:29 by jwe]
jwe
parents: 5290
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
7a51b27574bc [project @ 2005-04-21 16:03:29 by jwe]
jwe
parents: 5290
diff changeset
36 #include "mx-m-dm.h"
7a51b27574bc [project @ 2005-04-21 16:03:29 by jwe]
jwe
parents: 5290
diff changeset
37 #include "EIG.h"
5290
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
7a51b27574bc [project @ 2005-04-21 16:03:29 by jwe]
jwe
parents: 5290
diff changeset
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
7a51b27574bc [project @ 2005-04-21 16:03:29 by jwe]
jwe
parents: 5290
diff changeset
43 #include "pr-output.h"
7a51b27574bc [project @ 2005-04-21 16:03:29 by jwe]
jwe
parents: 5290
diff changeset
44 #include "utils.h"
5290
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
60 static Matrix
5295
015b15716cbe [project @ 2005-04-21 16:38:38 by jwe]
jwe
parents: 5293
diff changeset
61 null (const Matrix& A, octave_idx_type& rank)
5290
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
62 {
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
63 Matrix retval;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
64
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
65 rank = 0;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
66
23577
80c42f4cca13 maint: Deprecate is_empty and replace with isempty.
Rik <rik@octave.org>
parents: 23455
diff changeset
67 if (! A.isempty ())
5290
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
70
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
71 DiagMatrix S = A_svd.singular_values ();
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
74
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
75 Matrix V = A_svd.right_singular_matrix ();
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
76
5295
015b15716cbe [project @ 2005-04-21 16:38:38 by jwe]
jwe
parents: 5293
diff changeset
77 octave_idx_type A_nr = A.rows ();
015b15716cbe [project @ 2005-04-21 16:38:38 by jwe]
jwe
parents: 5293
diff changeset
78 octave_idx_type A_nc = A.cols ();
5290
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
85
5295
015b15716cbe [project @ 2005-04-21 16:38:38 by jwe]
jwe
parents: 5293
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
91
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
ff87ad14403f [project @ 2007-03-22 18:20:31 by jwe]
jwe
parents: 5823
diff changeset
96
ff87ad14403f [project @ 2007-03-22 18:20:31 by jwe]
jwe
parents: 5823
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
100 }
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
101
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
102 return retval;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
103 }
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
104
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
105 static int
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
106 qp (const Matrix& H, const ColumnVector& q,
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
107 const Matrix& Aeq, const ColumnVector& beq,
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
110 ColumnVector& x, ColumnVector& lambda, int& iter)
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
111 {
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
112 int info = 0;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
113
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
114 iter = 0;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
115
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
118
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
122
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
123 // Filling the current active set.
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
124
5295
015b15716cbe [project @ 2005-04-21 16:38:38 by jwe]
jwe
parents: 5293
diff changeset
125 octave_idx_type n_act = n_eq;
5290
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
126
5295
015b15716cbe [project @ 2005-04-21 16:38:38 by jwe]
jwe
parents: 5293
diff changeset
127 octave_idx_type n_tot = n_eq + n_in;
5290
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
128
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
129 // Equality constraints come first. We won't check the sign of the
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
130 // Lagrange multiplier for those.
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
131
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
132 Matrix Aact = Aeq;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
133 ColumnVector bact = beq;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
134 ColumnVector Wact;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
135
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
136 if (n_in > 0)
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
137 {
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
138 ColumnVector res = Ain*x - bin;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
139
5295
015b15716cbe [project @ 2005-04-21 16:38:38 by jwe]
jwe
parents: 5293
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
152 }
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
153
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
154 // Computing the ???
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
170
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
171 bool done = false;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
172
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
173 double alpha = 0.0;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
174
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
175 Matrix R;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
176 Matrix Y (n, 0, 0.0);
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
177
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
178 ColumnVector g (n, 0.0);
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
179 ColumnVector p (n, 0.0);
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
180
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
181 ColumnVector lambda_tmp (n_in, 0.0);
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
182
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
183 while (! done)
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
184 {
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
185 iter++;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
186
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
187 // Current Gradient
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
188 // g = q + H * x;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
189
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
190 g = q + H * x;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
191
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
205
10154
40dfc0c99116 DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents: 7729
diff changeset
206 R = cholH.chol_matrix ();
5290
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
212
10154
40dfc0c99116 DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents: 7729
diff changeset
213 p = -Hinv * g;
5290
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
220
10154
40dfc0c99116 DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents: 7729
diff changeset
221 p = eigenvecH.column (indminR);
5290
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
230
10154
40dfc0c99116 DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents: 7729
diff changeset
231 // Multipliers are zero.
5290
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
232 lambda_tmp.fill (0.0);
10154
40dfc0c99116 DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents: 7729
diff changeset
233 }
5290
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
234 else
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
235 {
10154
40dfc0c99116 DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents: 7729
diff changeset
236 // There are active constraints.
5290
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
237
10154
40dfc0c99116 DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents: 7729
diff changeset
238 // Computing the null space.
5290
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
239
10154
40dfc0c99116 DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents: 7729
diff changeset
240 octave_idx_type rank;
5290
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
241
10154
40dfc0c99116 DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents: 7729
diff changeset
242 Matrix Z = null (Aact, rank);
5290
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
a2870fd8ac58 [project @ 2008-01-12 07:50:54 by jwe]
jwe
parents: 7035
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
255
10154
40dfc0c99116 DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents: 7729
diff changeset
256 octave_idx_type pR = 0;
5290
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
279
10154
40dfc0c99116 DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents: 7729
diff changeset
280 ColumnVector pz = -rHinv * Zt * g;
5290
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
311
10154
40dfc0c99116 DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents: 7729
diff changeset
312 ColumnVector eVrH = eigenvecrH.column (indminR);
5290
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
321
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
322 // Checking the step-size.
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
323 ColumnVector abs_p (n);
5295
015b15716cbe [project @ 2005-04-21 16:38:38 by jwe]
jwe
parents: 5293
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
326 double max_p = abs_p.max ();
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
327
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
456
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
463 }
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
464
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
465 lambda_tmp = Y.transpose () * (g + H * p);
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
466
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
467 // Reordering the Lagrange multipliers.
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
468
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
469 lambda.resize (n_tot);
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
470 lambda.fill (0.0);
5295
015b15716cbe [project @ 2005-04-21 16:38:38 by jwe]
jwe
parents: 5293
diff changeset
471 for (octave_idx_type i = 0; i < n_eq; i++)
5290
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
472 lambda(i) = lambda_tmp(i);
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
473
5295
015b15716cbe [project @ 2005-04-21 16:38:38 by jwe]
jwe
parents: 5293
diff changeset
474 for (octave_idx_type i = n_eq; i < n_tot; i++)
5290
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
475 {
5295
015b15716cbe [project @ 2005-04-21 16:38:38 by jwe]
jwe
parents: 5293
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
484 }
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
485
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
486 return info;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
487 }
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
497
20892
c07bee629973 2015 Code Sprint: use ovl ().
Rik <rik@octave.org>
parents: 20853
diff changeset
498 const ColumnVector x0 (args(0).vector_value ());
c07bee629973 2015 Code Sprint: use ovl ().
Rik <rik@octave.org>
parents: 20853
diff changeset
499 const Matrix H (args(1).matrix_value ());
c07bee629973 2015 Code Sprint: use ovl ().
Rik <rik@octave.org>
parents: 20853
diff changeset
500 const ColumnVector q (args(2).vector_value ());
c07bee629973 2015 Code Sprint: use ovl ().
Rik <rik@octave.org>
parents: 20853
diff changeset
501 const Matrix Aeq (args(3).matrix_value ());
c07bee629973 2015 Code Sprint: use ovl ().
Rik <rik@octave.org>
parents: 20853
diff changeset
502 const ColumnVector beq (args(4).vector_value ());
c07bee629973 2015 Code Sprint: use ovl ().
Rik <rik@octave.org>
parents: 20853
diff changeset
503 const Matrix Ain (args(5).matrix_value ());
c07bee629973 2015 Code Sprint: use ovl ().
Rik <rik@octave.org>
parents: 20853
diff changeset
504 const ColumnVector bin (args(6).vector_value ());
c07bee629973 2015 Code Sprint: use ovl ().
Rik <rik@octave.org>
parents: 20853
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
509
20892
c07bee629973 2015 Code Sprint: use ovl ().
Rik <rik@octave.org>
parents: 20853
diff changeset
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
c07bee629973 2015 Code Sprint: use ovl ().
Rik <rik@octave.org>
parents: 20853
diff changeset
518 return ovl (x, lambda, info, iter);
5290
41273fff034d [project @ 2005-04-21 14:46:50 by jwe]
jwe
parents:
diff changeset
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