octave: libinterp/corefcn/__qp_

annotate libinterp/corefcn/qp.cc @ 21966:112b20240c87

move docstrings in C++ files out of C strings and into comments * __contourc__.cc, __dispatch__.cc, __dsearchn__.cc, __ichol__.cc, __ilu__.cc, __lin_interpn__.cc, __luinc__.cc, __magick_read__.cc, __pchip_deriv__.cc, __qp__.cc, balance.cc, besselj.cc, betainc.cc, bitfcns.cc, bsxfun.cc, cellfun.cc, colloc.cc, conv2.cc, daspk.cc, dasrt.cc, dassl.cc, data.cc, debug.cc, defaults.cc, det.cc, dirfns.cc, dlmread.cc, dot.cc, eig.cc, ellipj.cc, error.cc, fft.cc, fft2.cc, fftn.cc, file-io.cc, filter.cc, find.cc, gammainc.cc, gcd.cc, getgrent.cc, getpwent.cc, getrusage.cc, givens.cc, graphics.cc, hash.cc, help.cc, hess.cc, hex2num.cc, input.cc, inv.cc, kron.cc, load-path.cc, load-save.cc, lookup.cc, ls-oct-text.cc, lsode.cc, lu.cc, mappers.cc, matrix_type.cc, max.cc, mgorth.cc, nproc.cc, oct-hist.cc, octave-link.cc, ordschur.cc, pager.cc, pinv.cc, pr-output.cc, profiler.cc, psi.cc, pt-jit.cc, quad.cc, quadcc.cc, qz.cc, rand.cc, rcond.cc, regexp.cc, schur.cc, sighandlers.cc, sparse.cc, spparms.cc, sqrtm.cc, str2double.cc, strfind.cc, strfns.cc, sub2ind.cc, svd.cc, sylvester.cc, symtab.cc, syscalls.cc, sysdep.cc, time.cc, toplev.cc, tril.cc, tsearch.cc, typecast.cc, urlwrite.cc, utils.cc, variables.cc, __delaunayn__.cc, __eigs__.cc, __fltk_uigetfile__.cc, __glpk__.cc, __init_fltk__.cc, __init_gnuplot__.cc, __osmesa_print__.cc, __voronoi__.cc, amd.cc, audiodevinfo.cc, audioread.cc, ccolamd.cc, chol.cc, colamd.cc, convhulln.cc, dmperm.cc, fftw.cc, qr.cc, symbfact.cc, symrcm.cc, ov-base.cc, ov-bool-mat.cc, ov-cell.cc, ov-class.cc, ov-classdef.cc, ov-fcn-handle.cc, ov-fcn-inline.cc, ov-flt-re-mat.cc, ov-int16.cc, ov-int32.cc, ov-int64.cc, ov-int8.cc, ov-java.cc, ov-null-mat.cc, ov-oncleanup.cc, ov-range.cc, ov-re-mat.cc, ov-struct.cc, ov-typeinfo.cc, ov-uint16.cc, ov-uint32.cc, ov-uint64.cc, ov-uint8.cc, ov-usr-fcn.cc, ov.cc, octave.cc, pt-arg-list.cc, pt-binop.cc, pt-eval.cc, pt-mat.cc, lex.ll, oct-parse.in.yy: Docstrings are now comments instead of C strings. * build-aux/mk-opts.pl: Emit docstrings as comments instead of C strings. * DASPK-opts.in, LSODE-opts.in: Don't quote " in docstring fragments. * builtins.h: Include builtin-defun-decls.h unconditionally. * defun.h (DEFUN, DEFUNX, DEFCONSTFUN): Simply emit declaration. (DEFALIAS): Always expand to nothing. * defun-dld.h: No special macro expansions for MAKE_BUILTINS. (DEFUN_DLD): Use FORWARD_DECLARE_FUN. (DEFUNX_DLD): Use FORWARD_DECLARE_FUNX. * defun-int.h: No special macro expansions for MAKE_BUILTINS. (FORWARD_DECLARE_FUN, FORWARD_DECLARE_FUNX): New macros. (DEFINE_FUN_INSTALLER_FUN): If compiling an Octave source file, pass "external-doc" to DEFINE_FUNX_INSTALLER_FUN. (DEFUN_INTERNAL, DEFCONSTFUN_INTERNAL, DEFUNX_INTERNAL, DEFALIAS_INTERNAL): Delete. * common.mk (move_if_change_rule): New macro. (simple_move_if_change_rule): Define using move_if_change_rule. * find-defun-files.sh (DEFUN_PATTERN): Update. Don't transform file name extension to ".df". * libinterp/mk-pkg-add, gendoc.pl: Operate directly on source files. * mkbuiltins: New argument, SRCDIR. Operate directly on source files. * mkdefs: Delete. * libinterp/module.mk (BUILT_SOURCES): Update list to contain only files included in other source files. (GENERATED_MAKE_BUILTINS_INCS, DEF_FILES): Delete. (LIBINTERP_BUILT_DISTFILES): Include $(OPT_HANDLERS) here. (LIBINTERP_BUILT_NODISTFILES): Not here. Remove $(ALL_DEF_FILES from the list. (libinterp_EXTRA_DIST): Remove mkdefs from the list. (FOUND_DEFUN_FILES): Rename from SRC_DEF_FILES. (DLDFCN_DEFUN_FILES): Rename from DLDFCN_DEF_FILES. (SRC_DEFUN_FILES): Rename from SRC_DEF_FILES. (ALL_DEFUN_FILES): Rename from ALL_DEF_FILES. (%.df: %.cc): Delete pattern rule. (libinterp/build-env-features.cc, libinterp/builtins.cc, libinterp/dldfcn/PKG_ADD): Use mv instead of move-if-change. (libinterp/builtins.cc, libinterp/builtin-defun-decls.h): Update mkbuiltins command. ($(srcdir)/libinterp/DOCSTRINGS): Update gendoc.pl command. * liboctave/module.mk (BUILT_SOURCES): Don't include liboctave-build-info.cc in the list.

author	John W. Eaton <jwe@octave.org>
date	Tue, 21 Jun 2016 16:07:51 -0400
parents	aba2e6293dd8
children	510886d03ef2

rev	line source
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	1 /*
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	2
19697 4197fc428c7d maint: Update copyright notices for 2015. John W. Eaton <jwe@octave.org> parents: 17787 diff changeset	3 Copyright (C) 2000-2015 Gabriele Pannocchia
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	4
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	5 This file is part of Octave.
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	6
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	7 Octave is free software; you can redistribute it and/or modify it
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	8 under the terms of the GNU General Public License as published by the
7016 93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6945 diff changeset	9 Free Software Foundation; either version 3 of the License, or (at your
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6945 diff changeset	10 option) any later version.
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	11
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	12 Octave is distributed in the hope that it will be useful, but WITHOUT
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	15 for more details.
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	16
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	17 You should have received a copy of the GNU General Public License
7016 93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6945 diff changeset	18 along with Octave; see the file COPYING. If not, see
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 6945 diff changeset	19 <http://www.gnu.org/licenses/>.
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	20
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	21 */
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	22
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21301 diff changeset	23 #if defined (HAVE_CONFIG_H)
21301 40de9f8f23a6 Use '#include "config.h"' rather than <config.h>. Rik <rik@octave.org> parents: 21273 diff changeset	24 # include "config.h"
5293 7a51b27574bc [project @ 2005-04-21 16:03:29 by jwe] jwe parents: 5290 diff changeset	25 #endif
7a51b27574bc [project @ 2005-04-21 16:03:29 by jwe] jwe parents: 5290 diff changeset	26
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	27 #include <cfloat>
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	28
21269 3c8a3d35661a better use of templates for Cholesky factorization John W. Eaton <jwe@octave.org> parents: 21200 diff changeset	29 #include "chol.h"
21273 cbced1c09916 better use of templates for svd classes John W. Eaton <jwe@octave.org> parents: 21269 diff changeset	30 #include "svd.h"
5293 7a51b27574bc [project @ 2005-04-21 16:03:29 by jwe] jwe parents: 5290 diff changeset	31 #include "mx-m-dm.h"
7a51b27574bc [project @ 2005-04-21 16:03:29 by jwe] jwe parents: 5290 diff changeset	32 #include "EIG.h"
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	33
15039 e753177cde93 maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory Rik <rik@octave.org> parents: 14501 diff changeset	34 #include "defun.h"
5293 7a51b27574bc [project @ 2005-04-21 16:03:29 by jwe] jwe parents: 5290 diff changeset	35 #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	36 #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	37 #include "ovl.h"
5293 7a51b27574bc [project @ 2005-04-21 16:03:29 by jwe] jwe parents: 5290 diff changeset	38 #include "pr-output.h"
7a51b27574bc [project @ 2005-04-21 16:03:29 by jwe] jwe parents: 5290 diff changeset	39 #include "utils.h"
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	40
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	41 static Matrix
5295 015b15716cbe [project @ 2005-04-21 16:38:38 by jwe] jwe parents: 5293 diff changeset	42 null (const Matrix& A, octave_idx_type& rank)
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	43 {
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	44 Matrix retval;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	45
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	46 rank = 0;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	47
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	48 if (! A.is_empty ())
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	49 {
21273 cbced1c09916 better use of templates for svd classes John W. Eaton <jwe@octave.org> parents: 21269 diff changeset	50 svd<Matrix> A_svd (A);
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	51
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	52 DiagMatrix S = A_svd.singular_values ();
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	53
15448 0a0912a9ab6e Replace deprecated DiagArray2<T>::diag calls with DiagArray2<T>::extract_diag Jordi Gutiérrez Hermoso <jordigh@octave.org> parents: 15220 diff changeset	54 ColumnVector s = S.extract_diag ();
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	55
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	56 Matrix V = A_svd.right_singular_matrix ();
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	57
5295 015b15716cbe [project @ 2005-04-21 16:38:38 by jwe] jwe parents: 5293 diff changeset	58 octave_idx_type A_nr = A.rows ();
015b15716cbe [project @ 2005-04-21 16:38:38 by jwe] jwe parents: 5293 diff changeset	59 octave_idx_type A_nc = A.cols ();
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	60
5295 015b15716cbe [project @ 2005-04-21 16:38:38 by jwe] jwe parents: 5293 diff changeset	61 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	62
15220 61822c866ba1 use std::numeric_limits<T>::epsilon in C++ code John W. Eaton <jwe@octave.org> parents: 15195 diff changeset	63 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	64
20232 a9574e3c6e9e Deprecate Array::length() and Sparse::length() in favour of ::numel(). Carnë Draug <carandraug@octave.org> parents: 19697 diff changeset	65 octave_idx_type n = s.numel ();
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	66
5295 015b15716cbe [project @ 2005-04-21 16:38:38 by jwe] jwe parents: 5293 diff changeset	67 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	68 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	69 if (s(i) > tol)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	70 rank++;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	71 }
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	72
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	73 if (rank < A_nc)
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	74 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	75 else
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	76 retval.resize (A_nc, 0);
6431 ff87ad14403f [project @ 2007-03-22 18:20:31 by jwe] jwe parents: 5823 diff changeset	77
ff87ad14403f [project @ 2007-03-22 18:20:31 by jwe] jwe parents: 5823 diff changeset	78 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	79 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	80 retval(i) = 0;
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	81 }
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	82
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	83 return retval;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	84 }
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	85
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	86 static int
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	87 qp (const Matrix& H, const ColumnVector& q,
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	88 const Matrix& Aeq, const ColumnVector& beq,
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	89 const Matrix& Ain, const ColumnVector& bin,
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	90 int maxit,
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	91 ColumnVector& x, ColumnVector& lambda, int& iter)
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	92 {
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	93 int info = 0;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	94
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	95 iter = 0;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	96
15220 61822c866ba1 use std::numeric_limits<T>::epsilon in C++ code John W. Eaton <jwe@octave.org> parents: 15195 diff changeset	97 double rtol = sqrt (std::numeric_limits<double>::epsilon ());
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	98
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	99 // Problem dimension.
20232 a9574e3c6e9e Deprecate Array::length() and Sparse::length() in favour of ::numel(). Carnë Draug <carandraug@octave.org> parents: 19697 diff changeset	100 octave_idx_type n = x.numel ();
5290 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 // Dimension of constraints.
20232 a9574e3c6e9e Deprecate Array::length() and Sparse::length() in favour of ::numel(). Carnë Draug <carandraug@octave.org> parents: 19697 diff changeset	103 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	104 octave_idx_type n_in = bin.numel ();
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	105
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	106 // Filling the current active set.
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	107
5295 015b15716cbe [project @ 2005-04-21 16:38:38 by jwe] jwe parents: 5293 diff changeset	108 octave_idx_type n_act = n_eq;
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	109
5295 015b15716cbe [project @ 2005-04-21 16:38:38 by jwe] jwe parents: 5293 diff changeset	110 octave_idx_type n_tot = n_eq + n_in;
5290 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 // 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	113 // Lagrange multiplier for those.
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	114
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	115 Matrix Aact = Aeq;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	116 ColumnVector bact = beq;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	117 ColumnVector Wact;
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 if (n_in > 0)
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	120 {
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	121 ColumnVector res = Ain*x - bin;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	122
5295 015b15716cbe [project @ 2005-04-21 16:38:38 by jwe] jwe parents: 5293 diff changeset	123 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	124 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	125 res(i) /= (1.0 + std::abs (bin(i)));
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	126
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	127 if (res(i) < rtol)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	128 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	129 n_act++;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	130 Aact = Aact.stack (Ain.row (i));
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	131 bact.resize (n_act, bin(i));
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	132 Wact.resize (n_act-n_eq, i);
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	133 }
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	134 }
5290 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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	137 // Computing the ???
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	138
20741 a5ab31b52ae8 eliminate more uses of error_state John W. Eaton <jwe@octave.org> parents: 20678 diff changeset	139 EIG eigH;
7578 91f8446ce4ae handle possible error from EIG John W. Eaton <jwe@octave.org> parents: 7361 diff changeset	140
20741 a5ab31b52ae8 eliminate more uses of error_state John W. Eaton <jwe@octave.org> parents: 20678 diff changeset	141 try
a5ab31b52ae8 eliminate more uses of error_state John W. Eaton <jwe@octave.org> parents: 20678 diff changeset	142 {
a5ab31b52ae8 eliminate more uses of error_state John W. Eaton <jwe@octave.org> parents: 20678 diff changeset	143 eigH = EIG (H);
a5ab31b52ae8 eliminate more uses of error_state John W. Eaton <jwe@octave.org> parents: 20678 diff changeset	144 }
20785 b6f2909e7f94 always throw exception after debugging with debug_on_error John W. Eaton <jwe@octave.org> parents: 20756 diff changeset	145 catch (octave_execution_exception& e)
7578 91f8446ce4ae handle possible error from EIG John W. Eaton <jwe@octave.org> parents: 7361 diff changeset	146 {
20756 200ae1d650b7 propagate octave_execution_exception objects through try/catch blocks John W. Eaton <jwe@octave.org> parents: 20741 diff changeset	147 error (e, "qp: failed to compute eigenvalues of H");
7578 91f8446ce4ae handle possible error from EIG John W. Eaton <jwe@octave.org> parents: 7361 diff changeset	148 }
91f8446ce4ae handle possible error from EIG John W. Eaton <jwe@octave.org> parents: 7361 diff changeset	149
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	150 ColumnVector eigenvalH = real (eigH.eigenvalues ());
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	151 Matrix eigenvecH = real (eigH.eigenvectors ());
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	152 double minReal = eigenvalH.min ();
5295 015b15716cbe [project @ 2005-04-21 16:38:38 by jwe] jwe parents: 5293 diff changeset	153 octave_idx_type indminR = 0;
015b15716cbe [project @ 2005-04-21 16:38:38 by jwe] jwe parents: 5293 diff changeset	154 for (octave_idx_type i = 0; i < n; i++)
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	155 {
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	156 if (minReal == eigenvalH(i))
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	157 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	158 indminR = i;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	159 break;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	160 }
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	161 }
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	162
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	163 bool done = false;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	164
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	165 double alpha = 0.0;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	166
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	167 Matrix R;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	168 Matrix Y (n, 0, 0.0);
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	169
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	170 ColumnVector g (n, 0.0);
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	171 ColumnVector p (n, 0.0);
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 ColumnVector lambda_tmp (n_in, 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 while (! done)
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	176 {
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	177 iter++;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	178
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	179 // Current Gradient
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	180 // g = q + H * x;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	181
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	182 g = q + H * x;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	183
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	184 if (n_act == 0)
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	185 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	186 // There are no active constraints.
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	187
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	188 if (minReal > 0.0)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	189 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	190 // Inverting the Hessian. Using the Cholesky
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	191 // factorization since the Hessian is positive
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	192 // definite.
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	193
21269 3c8a3d35661a better use of templates for Cholesky factorization John W. Eaton <jwe@octave.org> parents: 21200 diff changeset	194 chol<Matrix> cholH (H);
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	195
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	196 R = cholH.chol_matrix ();
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	197
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	198 Matrix Hinv = chol2inv (R);
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	199
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	200 // Computing the unconstrained step.
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	201 // p = -Hinv * g;
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	202
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	203 p = -Hinv * g;
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	204
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	205 info = 0;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	206 }
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	207 else
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	208 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	209 // Finding the negative curvature of H.
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	210
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	211 p = eigenvecH.column (indminR);
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 // Following the negative curvature of H.
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	214
15220 61822c866ba1 use std::numeric_limits<T>::epsilon in C++ code John W. Eaton <jwe@octave.org> parents: 15195 diff changeset	215 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	216 p = -p;
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	217
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	218 info = 1;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	219 }
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 // Multipliers are zero.
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	222 lambda_tmp.fill (0.0);
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	223 }
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	224 else
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	225 {
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	226 // There are active constraints.
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 // Computing the null space.
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	229
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	230 octave_idx_type rank;
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	231
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	232 Matrix Z = null (Aact, rank);
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	233
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	234 octave_idx_type dimZ = n - rank;
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	235
17787 175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	236 // FIXME: still remain to handle the case of
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	237 // non-full rank active set matrix.
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	238
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	239 // Computing the Y matrix (orthogonal to Z)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	240 Y = Aact.pseudo_inverse ();
7361 a2870fd8ac58 [project @ 2008-01-12 07:50:54 by jwe] jwe parents: 7035 diff changeset	241
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	242 // Reduced Hessian
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	243 Matrix Zt = Z.transpose ();
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	244 Matrix rH = Zt * H * Z;
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	245
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	246 octave_idx_type pR = 0;
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	247
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	248 if (dimZ > 0)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	249 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	250 // Computing the Cholesky factorization (pR = 0 means
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	251 // that the reduced Hessian was positive definite).
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	252
21269 3c8a3d35661a better use of templates for Cholesky factorization John W. Eaton <jwe@octave.org> parents: 21200 diff changeset	253 chol<Matrix> cholrH (rH, pR);
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	254 Matrix tR = cholrH.chol_matrix ();
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	255 if (pR == 0)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	256 R = tR;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	257 }
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	258
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	259 if (pR == 0)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	260 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	261 info = 0;
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	262
11586 12df7854fa7c strip trailing whitespace from source files John W. Eaton <jwe@octave.org> parents: 11523 diff changeset	263 // Computing the step pz.
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	264 if (dimZ > 0)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	265 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	266 // Using the Cholesky factorization to invert rH
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	267
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	268 Matrix rHinv = chol2inv (R);
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	269
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	270 ColumnVector pz = -rHinv * Zt * g;
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	271
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	272 // Global step.
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	273 p = Z * pz;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	274 }
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	275 else
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	276 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	277 // Global step.
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	278 p.fill (0.0);
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	279 }
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	280 }
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	281 else
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	282 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	283 info = 1;
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	284
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	285 // Searching for the most negative curvature.
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	286
20741 a5ab31b52ae8 eliminate more uses of error_state John W. Eaton <jwe@octave.org> parents: 20678 diff changeset	287 EIG eigrH;
7578 91f8446ce4ae handle possible error from EIG John W. Eaton <jwe@octave.org> parents: 7361 diff changeset	288
20741 a5ab31b52ae8 eliminate more uses of error_state John W. Eaton <jwe@octave.org> parents: 20678 diff changeset	289 try
a5ab31b52ae8 eliminate more uses of error_state John W. Eaton <jwe@octave.org> parents: 20678 diff changeset	290 {
a5ab31b52ae8 eliminate more uses of error_state John W. Eaton <jwe@octave.org> parents: 20678 diff changeset	291 eigrH = EIG (rH);
a5ab31b52ae8 eliminate more uses of error_state John W. Eaton <jwe@octave.org> parents: 20678 diff changeset	292 }
20785 b6f2909e7f94 always throw exception after debugging with debug_on_error John W. Eaton <jwe@octave.org> parents: 20756 diff changeset	293 catch (octave_execution_exception& e)
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	294 {
20756 200ae1d650b7 propagate octave_execution_exception objects through try/catch blocks John W. Eaton <jwe@octave.org> parents: 20741 diff changeset	295 error (e, "qp: failed to compute eigenvalues of rH");
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	296 }
7578 91f8446ce4ae handle possible error from EIG John W. Eaton <jwe@octave.org> parents: 7361 diff changeset	297
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	298 ColumnVector eigenvalrH = real (eigrH.eigenvalues ());
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	299 Matrix eigenvecrH = real (eigrH.eigenvectors ());
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	300 double mRrH = eigenvalrH.min ();
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	301 indminR = 0;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	302 for (octave_idx_type i = 0; i < n; i++)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	303 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	304 if (mRrH == eigenvalH(i))
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	305 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	306 indminR = i;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	307 break;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	308 }
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	309 }
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	310
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	311 ColumnVector eVrH = eigenvecrH.column (indminR);
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	312
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	313 // Computing the step pz.
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	314 p = Z * eVrH;
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	315
15220 61822c866ba1 use std::numeric_limits<T>::epsilon in C++ code John W. Eaton <jwe@octave.org> parents: 15195 diff changeset	316 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	317 p = -p;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	318 }
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	319 }
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	320
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	321 // Checking the step-size.
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	322 ColumnVector abs_p (n);
5295 015b15716cbe [project @ 2005-04-21 16:38:38 by jwe] jwe parents: 5293 diff changeset	323 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	324 abs_p(i) = std::abs (p(i));
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	325 double max_p = abs_p.max ();
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	326
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	327 if (max_p < rtol)
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	328 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	329 // The step is null. Checking constraints.
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	330 if (n_act - n_eq == 0)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	331 // Solution is found because no inequality
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	332 // constraints are active.
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	333 done = true;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	334 else
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	335 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	336 // Computing the multipliers only for the inequality
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	337 // constraints that are active. We do NOT compute
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	338 // multipliers for the equality constraints.
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	339 Matrix Yt = Y.transpose ();
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	340 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	341 lambda_tmp = Yt * (g + H * p);
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	342
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	343 // Checking the multipliers. We remove the most
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	344 // negative from the set (if any).
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	345 double min_lambda = lambda_tmp.min ();
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	346 if (min_lambda >= 0)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	347 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	348 // Solution is found.
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	349 done = true;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	350 }
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	351 else
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	352 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	353 octave_idx_type which_eig = 0;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	354 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	355 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	356 if (lambda_tmp(i) == min_lambda)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	357 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	358 which_eig = i;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	359 break;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	360 }
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	361 }
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	362
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	363 // At least one multiplier is negative, we
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	364 // remove it from the set.
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	365
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	366 n_act--;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	367 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	368 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	369 Wact(i) = Wact(i+1);
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	370 for (octave_idx_type j = 0; j < n; j++)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	371 Aact(n_eq+i,j) = Aact(n_eq+i+1,j);
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	372 bact(n_eq+i) = bact(n_eq+i+1);
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	373 }
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	374
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	375 // Resizing the active set.
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	376 Wact.resize (n_act-n_eq);
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	377 bact.resize (n_act);
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	378 Aact.resize (n_act, n);
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	379 }
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 }
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	382 else
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	383 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	384 // The step is not null.
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	385 if (n_act - n_eq == n_in)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	386 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	387 // All inequality constraints were active. We can
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	388 // add the whole step.
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	389 x += p;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	390 }
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	391 else
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	392 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	393 // Some constraints were not active. Checking if
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	394 // there is a blocking constraint.
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	395 alpha = 1.0;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	396 octave_idx_type is_block = -1;
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	397
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	398 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	399 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	400 bool found = false;
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	401
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	402 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	403 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	404 if (Wact(j) == i)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	405 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	406 found = true;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	407 break;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	408 }
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	409 }
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	410
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	411 if (! found)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	412 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	413 // 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	414 // blocking constraint?
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	415
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	416 RowVector tmp_row = Ain.row (i);
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	417 double tmp = tmp_row * p;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	418 double res = tmp_row * x;
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	419
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	420 if (tmp < 0.0)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	421 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	422 double alpha_tmp = (bin(i) - res) / tmp;
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	423
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	424 if (alpha_tmp < alpha)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	425 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	426 alpha = alpha_tmp;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	427 is_block = i;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	428 }
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 }
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	432
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	433 // 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	434 // constraint (if any).
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	435 if (is_block >= 0)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	436 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	437 // There is a blocking constraint (index in
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	438 // is_block) which is added to the active set.
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	439 n_act++;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	440 Aact = Aact.stack (Ain.row (is_block));
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	441 bact.resize (n_act, bin(is_block));
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	442 Wact.resize (n_act-n_eq, is_block);
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	443
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	444 // Adding the reduced step
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	445 x += alpha * p;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	446 }
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	447 else
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	448 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	449 // There are no blocking constraints. Adding the
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	450 // whole step.
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	451 x += alpha * p;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	452 }
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 }
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	455
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	456 if (iter == maxit)
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	457 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	458 done = true;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	459 // warning ("qp_main: maximum number of iteration reached");
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	460 info = 3;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	461 }
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	462 }
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 lambda_tmp = Y.transpose () * (g + H * p);
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	465
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	466 // Reordering the Lagrange multipliers.
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	467
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	468 lambda.resize (n_tot);
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	469 lambda.fill (0.0);
5295 015b15716cbe [project @ 2005-04-21 16:38:38 by jwe] jwe parents: 5293 diff changeset	470 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	471 lambda(i) = lambda_tmp(i);
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	472
5295 015b15716cbe [project @ 2005-04-21 16:38:38 by jwe] jwe parents: 5293 diff changeset	473 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	474 {
5295 015b15716cbe [project @ 2005-04-21 16:38:38 by jwe] jwe parents: 5293 diff changeset	475 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	476 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	477 if (Wact(j) == i - n_eq)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	478 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	479 lambda(i) = lambda_tmp(n_eq+j);
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	480 break;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	481 }
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	482 }
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	483 }
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 return info;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	486 }
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	487
15039 e753177cde93 maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory Rik <rik@octave.org> parents: 14501 diff changeset	488 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	489 doc: /* -- texinfo --
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	490 @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})
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21724 diff changeset	491 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	492 @end deftypefn */)
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	493 {
20802 8bb38ba1bad6 eliminate return statements after calls to print_usage John W. Eaton <jwe@octave.org> parents: 20785 diff changeset	494 if (args.length () != 8)
8bb38ba1bad6 eliminate return statements after calls to print_usage John W. Eaton <jwe@octave.org> parents: 20785 diff changeset	495 print_usage ();
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	496
20892 c07bee629973 2015 Code Sprint: use ovl (). Rik <rik@octave.org> parents: 20853 diff changeset	497 const ColumnVector x0 (args(0).vector_value ());
c07bee629973 2015 Code Sprint: use ovl (). Rik <rik@octave.org> parents: 20853 diff changeset	498 const Matrix H (args(1).matrix_value ());
c07bee629973 2015 Code Sprint: use ovl (). Rik <rik@octave.org> parents: 20853 diff changeset	499 const ColumnVector q (args(2).vector_value ());
c07bee629973 2015 Code Sprint: use ovl (). Rik <rik@octave.org> parents: 20853 diff changeset	500 const Matrix Aeq (args(3).matrix_value ());
c07bee629973 2015 Code Sprint: use ovl (). Rik <rik@octave.org> parents: 20853 diff changeset	501 const ColumnVector beq (args(4).vector_value ());
c07bee629973 2015 Code Sprint: use ovl (). Rik <rik@octave.org> parents: 20853 diff changeset	502 const Matrix Ain (args(5).matrix_value ());
c07bee629973 2015 Code Sprint: use ovl (). Rik <rik@octave.org> parents: 20853 diff changeset	503 const ColumnVector bin (args(6).vector_value ());
c07bee629973 2015 Code Sprint: use ovl (). Rik <rik@octave.org> parents: 20853 diff changeset	504 const int maxit (args(7).int_value ());
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	505
20802 8bb38ba1bad6 eliminate return statements after calls to print_usage John W. Eaton <jwe@octave.org> parents: 20785 diff changeset	506 int iter = 0;
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	507
20892 c07bee629973 2015 Code Sprint: use ovl (). Rik <rik@octave.org> parents: 20853 diff changeset	508 // 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	509 ColumnVector x = x0;
20678 4b00afb5e9c3 eliminate more uses of error_state John W. Eaton <jwe@octave.org> parents: 20232 diff changeset	510
20802 8bb38ba1bad6 eliminate return statements after calls to print_usage John W. Eaton <jwe@octave.org> parents: 20785 diff changeset	511 // Reordering the Lagrange multipliers
8bb38ba1bad6 eliminate return statements after calls to print_usage John W. Eaton <jwe@octave.org> parents: 20785 diff changeset	512 ColumnVector lambda;
8bb38ba1bad6 eliminate return statements after calls to print_usage John W. Eaton <jwe@octave.org> parents: 20785 diff changeset	513
8bb38ba1bad6 eliminate return statements after calls to print_usage John W. Eaton <jwe@octave.org> parents: 20785 diff changeset	514 int info = qp (H, q, Aeq, beq, Ain, bin, maxit, x, lambda, iter);
8bb38ba1bad6 eliminate return statements after calls to print_usage John W. Eaton <jwe@octave.org> parents: 20785 diff changeset	515
20892 c07bee629973 2015 Code Sprint: use ovl (). Rik <rik@octave.org> parents: 20853 diff changeset	516 return ovl (x, lambda, info, iter);
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	517 }
12805 3641167e5b75 codesprint: .cc helper functions do not need tests Rik <octave@nomad.inbox5.com>* parents: 11586 diff changeset	518
3641167e5b75 codesprint: .cc helper functions do not need tests Rik <octave@nomad.inbox5.com>* parents: 11586 diff changeset	519 /*
3641167e5b75 codesprint: .cc helper functions do not need tests Rik <octave@nomad.inbox5.com>* parents: 11586 diff changeset	520 ## 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	521 %!assert (1)
3641167e5b75 codesprint: .cc helper functions do not need tests Rik <octave@nomad.inbox5.com>* parents: 11586 diff changeset	522 */

Mercurial > octave

annotate libinterp/corefcn/__qp__.cc @ 21966:112b20240c87

annotate libinterp/corefcn/qp.cc @ 21966:112b20240c87