octave-jordi: src/DLD-FUNCTIONS/__qp_

annotate src/DLD-FUNCTIONS/qp.cc @ 10315:57a59eae83cc

untabify src C++ source files

author	John W. Eaton <jwe@octave.org>
date	Thu, 11 Feb 2010 12:41:46 -0500
parents	40dfc0c99116
children	fd0a3ac60b0e

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
7361 a2870fd8ac58 [project @ 2008-01-12 07:50:54 by jwe] jwe parents: 7035 diff changeset	3 Copyright (C) 2000, 2001, 2004, 2005, 2006, 2007, 2008 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
5293 7a51b27574bc [project @ 2005-04-21 16:03:29 by jwe] jwe parents: 5290 diff changeset	23 #ifdef HAVE_CONFIG_H
7a51b27574bc [project @ 2005-04-21 16:03:29 by jwe] jwe parents: 5290 diff changeset	24 #include <config.h>
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
5293 7a51b27574bc [project @ 2005-04-21 16:03:29 by jwe] jwe parents: 5290 diff changeset	29 #include "dbleCHOL.h"
7a51b27574bc [project @ 2005-04-21 16:03:29 by jwe] jwe parents: 5290 diff changeset	30 #include "dbleSVD.h"
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
5293 7a51b27574bc [project @ 2005-04-21 16:03:29 by jwe] jwe parents: 5290 diff changeset	34 #include "defun-dld.h"
7a51b27574bc [project @ 2005-04-21 16:03:29 by jwe] jwe parents: 5290 diff changeset	35 #include "error.h"
7a51b27574bc [project @ 2005-04-21 16:03:29 by jwe] jwe parents: 5290 diff changeset	36 #include "gripes.h"
7a51b27574bc [project @ 2005-04-21 16:03:29 by jwe] jwe parents: 5290 diff changeset	37 #include "oct-obj.h"
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 {
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	50 SVD A_svd (A);
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	54 ColumnVector s = S.diag ();
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	63 double tol = tmp * s(0) * DBL_EPSILON;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	64
5295 015b15716cbe [project @ 2005-04-21 16:38:38 by jwe] jwe parents: 5293 diff changeset	65 octave_idx_type n = s.length ();
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++)
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	79 if (std::abs (retval(i)) < DBL_EPSILON)
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	97 double rtol = sqrt (DBL_EPSILON);
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.
5295 015b15716cbe [project @ 2005-04-21 16:38:38 by jwe] jwe parents: 5293 diff changeset	100 octave_idx_type n = x.length ();
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.
5295 015b15716cbe [project @ 2005-04-21 16:38:38 by jwe] jwe parents: 5293 diff changeset	103 octave_idx_type n_eq = beq.length ();
015b15716cbe [project @ 2005-04-21 16:38:38 by jwe] jwe parents: 5293 diff changeset	104 octave_idx_type n_in = bin.length ();
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	139 EIG eigH (H);
7578 91f8446ce4ae handle possible error from EIG John W. Eaton <jwe@octave.org> parents: 7361 diff changeset	140
91f8446ce4ae handle possible error from EIG John W. Eaton <jwe@octave.org> parents: 7361 diff changeset	141 if (error_state)
91f8446ce4ae handle possible error from EIG John W. Eaton <jwe@octave.org> parents: 7361 diff changeset	142 {
91f8446ce4ae handle possible error from EIG John W. Eaton <jwe@octave.org> parents: 7361 diff changeset	143 error ("qp: failed to compute eigenvalues of H");
91f8446ce4ae handle possible error from EIG John W. Eaton <jwe@octave.org> parents: 7361 diff changeset	144 return -1;
91f8446ce4ae handle possible error from EIG John W. Eaton <jwe@octave.org> parents: 7361 diff changeset	145 }
91f8446ce4ae handle possible error from EIG John W. Eaton <jwe@octave.org> parents: 7361 diff changeset	146
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	147 ColumnVector eigenvalH = real (eigH.eigenvalues ());
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	148 Matrix eigenvecH = real (eigH.eigenvectors ());
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	149 double minReal = eigenvalH.min ();
5295 015b15716cbe [project @ 2005-04-21 16:38:38 by jwe] jwe parents: 5293 diff changeset	150 octave_idx_type indminR = 0;
015b15716cbe [project @ 2005-04-21 16:38:38 by jwe] jwe parents: 5293 diff changeset	151 for (octave_idx_type i = 0; i < n; i++)
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 if (minReal == eigenvalH(i))
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	154 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	155 indminR = i;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	156 break;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	157 }
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	158 }
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	159
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	160 bool done = false;
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 double alpha = 0.0;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	163
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	164 Matrix R;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	165 Matrix Y (n, 0, 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 ColumnVector g (n, 0.0);
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	168 ColumnVector p (n, 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 lambda_tmp (n_in, 0.0);
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	171
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	172 while (! done)
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	173 {
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	174 iter++;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	175
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	176 // Current Gradient
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	177 // g = q + H * x;
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 g = q + H * x;
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 if (n_act == 0)
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	182 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	183 // There are no active constraints.
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	184
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	185 if (minReal > 0.0)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	186 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	187 // Inverting the Hessian. Using the Cholesky
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	188 // factorization since the Hessian is positive
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	189 // definite.
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	190
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	191 CHOL cholH (H);
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	192
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	193 R = cholH.chol_matrix ();
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	194
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	195 Matrix Hinv = chol2inv (R);
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	196
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	197 // Computing the unconstrained step.
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	198 // p = -Hinv * g;
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 p = -Hinv * g;
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	201
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	202 info = 0;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	203 }
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	204 else
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	205 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	206 // Finding the negative curvature of H.
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	207
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	208 p = eigenvecH.column (indminR);
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 // Following the negative curvature of H.
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	211
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	212 if (p.transpose () * g > DBL_EPSILON)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	213 p = -p;
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 = 1;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	216 }
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 // Multipliers are zero.
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	219 lambda_tmp.fill (0.0);
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	220 }
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	221 else
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 // There are active constraints.
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	224
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	225 // Computing the null space.
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	226
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	227 octave_idx_type rank;
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	228
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	229 Matrix Z = null (Aact, rank);
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 octave_idx_type dimZ = n - rank;
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	232
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	233 // FIXME -- still remain to handle the case of
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	234 // non-full rank active set matrix.
5290 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 // Computing the Y matrix (orthogonal to Z)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	237 Y = Aact.pseudo_inverse ();
7361 a2870fd8ac58 [project @ 2008-01-12 07:50:54 by jwe] jwe parents: 7035 diff changeset	238
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	239 // Reduced Hessian
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	240 Matrix Zt = Z.transpose ();
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	241 Matrix rH = Zt * H * Z;
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	242
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	243 octave_idx_type pR = 0;
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	244
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	245 if (dimZ > 0)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	246 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	247 // Computing the Cholesky factorization (pR = 0 means
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	248 // that the reduced Hessian was positive definite).
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	249
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	250 CHOL cholrH (rH, pR);
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	251 Matrix tR = cholrH.chol_matrix ();
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	252 if (pR == 0)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	253 R = tR;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	254 }
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 if (pR == 0)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	257 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	258 info = 0;
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	259
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	260 // Computing the step pz.
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	261 if (dimZ > 0)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	262 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	263 // Using the Cholesky factorization to invert rH
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	264
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	265 Matrix rHinv = chol2inv (R);
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	266
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	267 ColumnVector pz = -rHinv * Zt * g;
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 // Global step.
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	270 p = Z * pz;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	271 }
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	272 else
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	273 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	274 // Global step.
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	275 p.fill (0.0);
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 }
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	278 else
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 info = 1;
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 // Searching for the most negative curvature.
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	283
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	284 EIG eigrH (rH);
7578 91f8446ce4ae handle possible error from EIG John W. Eaton <jwe@octave.org> parents: 7361 diff changeset	285
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	286 if (error_state)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	287 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	288 error ("qp: failed to compute eigenvalues of rH");
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	289 return -1;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	290 }
7578 91f8446ce4ae handle possible error from EIG John W. Eaton <jwe@octave.org> parents: 7361 diff changeset	291
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	292 ColumnVector eigenvalrH = real (eigrH.eigenvalues ());
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	293 Matrix eigenvecrH = real (eigrH.eigenvectors ());
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	294 double mRrH = eigenvalrH.min ();
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	295 indminR = 0;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	296 for (octave_idx_type i = 0; i < n; i++)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	297 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	298 if (mRrH == eigenvalH(i))
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	299 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	300 indminR = i;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	301 break;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	302 }
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	303 }
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	304
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	305 ColumnVector eVrH = eigenvecrH.column (indminR);
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	306
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	307 // Computing the step pz.
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	308 p = Z * eVrH;
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	309
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	310 if (p.transpose () * g > DBL_EPSILON)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	311 p = -p;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	312 }
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	313 }
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	314
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	315 // Checking the step-size.
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	316 ColumnVector abs_p (n);
5295 015b15716cbe [project @ 2005-04-21 16:38:38 by jwe] jwe parents: 5293 diff changeset	317 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	318 abs_p(i) = std::abs (p(i));
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	319 double max_p = abs_p.max ();
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 if (max_p < rtol)
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	322 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	323 // The step is null. Checking constraints.
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	324 if (n_act - n_eq == 0)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	325 // Solution is found because no inequality
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	326 // constraints are active.
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	327 done = true;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	328 else
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 // Computing the multipliers only for the inequality
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	331 // constraints that are active. We do NOT compute
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	332 // multipliers for the equality constraints.
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	333 Matrix Yt = Y.transpose ();
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	334 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	335 lambda_tmp = Yt * (g + H * p);
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	336
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	337 // Checking the multipliers. We remove the most
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	338 // negative from the set (if any).
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	339 double min_lambda = lambda_tmp.min ();
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	340 if (min_lambda >= 0)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	341 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	342 // Solution is found.
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	343 done = true;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	344 }
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	345 else
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	346 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	347 octave_idx_type which_eig = 0;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	348 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	349 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	350 if (lambda_tmp(i) == min_lambda)
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 which_eig = i;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	353 break;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	354 }
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	355 }
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	356
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	357 // At least one multiplier is negative, we
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	358 // remove it from the set.
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	359
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	360 n_act--;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	361 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	362 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	363 Wact(i) = Wact(i+1);
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	364 for (octave_idx_type j = 0; j < n; j++)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	365 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	366 bact(n_eq+i) = bact(n_eq+i+1);
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	367 }
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	368
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	369 // Resizing the active set.
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	370 Wact.resize (n_act-n_eq);
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	371 bact.resize (n_act);
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	372 Aact.resize (n_act, n);
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	373 }
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	374 }
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	375 }
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	376 else
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	377 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	378 // The step is not null.
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	379 if (n_act - n_eq == n_in)
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 // All inequality constraints were active. We can
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	382 // add the whole step.
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	383 x += p;
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 else
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 // Some constraints were not active. Checking if
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	388 // there is a blocking constraint.
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	389 alpha = 1.0;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	390 octave_idx_type is_block = -1;
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	391
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	392 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	393 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	394 bool found = false;
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	395
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	396 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	397 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	398 if (Wact(j) == 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 found = true;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	401 break;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	402 }
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	403 }
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	404
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	405 if (! found)
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 // 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	408 // blocking constraint?
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	409
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	410 RowVector tmp_row = Ain.row (i);
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	411 double tmp = tmp_row * p;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	412 double res = tmp_row * x;
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	413
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	414 if (tmp < 0.0)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	415 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	416 double alpha_tmp = (bin(i) - res) / tmp;
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	417
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	418 if (alpha_tmp < alpha)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	419 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	420 alpha = alpha_tmp;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	421 is_block = i;
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 }
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	424 }
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	425 }
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	426
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	427 // 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	428 // constraint (if any).
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	429 if (is_block >= 0)
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 // There is a blocking constraint (index in
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	432 // is_block) which is added to the active set.
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	433 n_act++;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	434 Aact = Aact.stack (Ain.row (is_block));
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	435 bact.resize (n_act, bin(is_block));
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	436 Wact.resize (n_act-n_eq, is_block);
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	437
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	438 // Adding the reduced step
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	439 x += alpha * p;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	440 }
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	441 else
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	442 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	443 // There are no blocking constraints. Adding the
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	444 // whole 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 }
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	448 }
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	449
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	450 if (iter == maxit)
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	451 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	452 done = true;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	453 // warning ("qp_main: maximum number of iteration reached");
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	454 info = 3;
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
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	458 lambda_tmp = Y.transpose () * (g + H * p);
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	459
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	460 // Reordering the Lagrange multipliers.
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	461
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	462 lambda.resize (n_tot);
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	463 lambda.fill (0.0);
5295 015b15716cbe [project @ 2005-04-21 16:38:38 by jwe] jwe parents: 5293 diff changeset	464 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	465 lambda(i) = lambda_tmp(i);
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	466
5295 015b15716cbe [project @ 2005-04-21 16:38:38 by jwe] jwe parents: 5293 diff changeset	467 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	468 {
5295 015b15716cbe [project @ 2005-04-21 16:38:38 by jwe] jwe parents: 5293 diff changeset	469 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	470 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	471 if (Wact(j) == i - n_eq)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	472 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	473 lambda(i) = lambda_tmp(n_eq+j);
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	474 break;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	475 }
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	476 }
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	477 }
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	478
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	479 return info;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	480 }
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	481
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	482 DEFUN_DLD (__qp__, args, ,
6945 6bbf56a9718a [project @ 2007-10-02 20:47:22 by jwe] jwe parents: 6860 diff changeset	483 "-- texinfo --\n\
6bbf56a9718a [project @ 2007-10-02 20:47:22 by jwe] jwe parents: 6860 diff changeset	484 @deftypefn {Loadable Function} {[@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})\n\
6bbf56a9718a [project @ 2007-10-02 20:47:22 by jwe] jwe parents: 6860 diff changeset	485 Undocumented internal function.\n\
6bbf56a9718a [project @ 2007-10-02 20:47:22 by jwe] jwe parents: 6860 diff changeset	486 @end deftypefn")
5290 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 octave_value_list retval;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	489
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	490 if (args.length () == 8)
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	491 {
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	492 const ColumnVector x0 (args(0) . vector_value ());
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	493 const Matrix H (args(1) . matrix_value ());
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	494 const ColumnVector q (args(2) . vector_value ());
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	495 const Matrix Aeq (args(3) . matrix_value ());
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	496 const ColumnVector beq (args(4) . vector_value ());
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	497 const Matrix Ain (args(5) . matrix_value ());
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	498 const ColumnVector bin (args(6) . vector_value ());
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	499 const int maxit (args(7) . int_value ());
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	500
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	501 if (! error_state)
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	502 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	503 int iter = 0;
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	504
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	505 // Copying the initial guess in the working variable
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	506 ColumnVector x = x0;
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	507
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	508 // Reordering the Lagrange multipliers
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	509 ColumnVector lambda;
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	510
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	511 int info = qp (H, q, Aeq, beq, Ain, bin, maxit, x, lambda, iter);
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	512
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	513 if (! error_state)
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	514 {
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	515 retval(3) = iter;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	516 retval(2) = info;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	517 retval(1) = lambda;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	518 retval(0) = x;
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	519 }
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	520 else
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	521 error ("qp: internal error");
40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	522 }
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	523 else
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 7729 diff changeset	524 error ("__qp__: invalid arguments");
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	525 }
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	526 else
5823 080c08b192d8 [project @ 2006-05-19 05:32:17 by jwe] jwe parents: 5775 diff changeset	527 print_usage ();
5290 41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	528
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	529 return retval;
41273fff034d [project @ 2005-04-21 14:46:50 by jwe] jwe parents: diff changeset	530 }

Mercurial > hg > octave-jordi

annotate src/DLD-FUNCTIONS/__qp__.cc @ 10315:57a59eae83cc

annotate src/DLD-FUNCTIONS/qp.cc @ 10315:57a59eae83cc