Mercurial > octave
view scripts/sparse/spaugment.m @ 20711:7b608fadc663
Make error messages more specific about the variable and problem encountered.
* besselj.cc, bitfcns.cc, colloc.cc, daspk.cc, dasrt.cc, dassl.cc, data.cc,
dirfns.cc, ellipj.cc, error.cc, gl-render.cc, graphics.cc, graphics.in.h,
load-path.cc, lsode.cc, lu.cc, luinc.cc, oct-hist.cc, oct-obj.cc,
octave-link.cc, quad.cc, rand.cc, symtab.cc, sysdep.cc, toplev.cc, utils.cc,
variables.cc, __init_fltk__.cc, chol.cc, fftw.cc, ov-cell.cc, ov-ch-mat.cc,
ov-class.cc, ov-classdef.cc, ov-cx-mat.cc, ov-fcn-inline.cc, ov-struct.cc,
ov-usr-fcn.cc, CMatrix.cc, dMatrix.cc, fCMatrix.cc, fMatrix.cc, lo-specfun.cc,
curl.m, divergence.m, __fltk_file_filter__.m, __uiobject_split_args__.m,
uigetfile.m, doc.m, imshow.m, rref.m, subspace.m, edit.m, fileattrib.m, open.m,
substruct.m, annotation.m, axis.m, caxis.m, datetick.m, hidden.m, legend.m,
whitebg.m, colorbar.m, __add_datasource__.m, __ezplot__.m, __pie__.m,
__plt_get_axis_arg__.m, pan.m, __print_parse_opts__.m, rotate3d.m, subplot.m,
zoom.m, compan.m, addpref.m, getpref.m, setpref.m, powerset.m, bicg.m,
bicgstab.m, cgs.m, qmr.m, spaugment.m, pascal.m, moment.m, cstrcat.m,
system.tst:
Make error messages more specific about the variable and problem encountered.
author | Rik <rik@octave.org> |
---|---|
date | Wed, 18 Nov 2015 10:40:26 -0800 |
parents | df437a52bcaf |
children | 516bb87ea72e |
line wrap: on
line source
## Copyright (C) 2008-2015 David Bateman ## ## This file is part of Octave. ## ## Octave is free software; you can redistribute it and/or modify it ## under the terms of the GNU General Public License as published by ## the Free Software Foundation; either version 3 of the License, or (at ## your option) any later version. ## ## Octave is distributed in the hope that it will be useful, but ## WITHOUT ANY WARRANTY; without even the implied warranty of ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ## General Public License for more details. ## ## You should have received a copy of the GNU General Public License ## along with Octave; see the file COPYING. If not, see ## <http://www.gnu.org/licenses/>. ## -*- texinfo -*- ## @deftypefn {Function File} {@var{s} =} spaugment (@var{A}, @var{c}) ## Create the augmented matrix of @var{A}. ## ## This is given by ## ## @example ## @group ## [@var{c} * eye(@var{m}, @var{m}), @var{A}; ## @var{A}', zeros(@var{n}, @var{n})] ## @end group ## @end example ## ## @noindent ## This is related to the least squares solution of ## @code{@var{A} \ @var{b}}, by ## ## @example ## @group ## @var{s} * [ @var{r} / @var{c}; x] = [ @var{b}, zeros(@var{n}, columns(@var{b})) ] ## @end group ## @end example ## ## @noindent ## where @var{r} is the residual error ## ## @example ## @var{r} = @var{b} - @var{A} * @var{x} ## @end example ## ## As the matrix @var{s} is symmetric indefinite it can be factorized with ## @code{lu}, and the minimum norm solution can therefore be found without the ## need for a @code{qr} factorization. As the residual error will be ## @code{zeros (@var{m}, @var{m})} for underdetermined problems, and example ## can be ## ## @example ## @group ## m = 11; n = 10; mn = max (m, n); ## A = spdiags ([ones(mn,1), 10*ones(mn,1), -ones(mn,1)], ## [-1, 0, 1], m, n); ## x0 = A \ ones (m,1); ## s = spaugment (A); ## [L, U, P, Q] = lu (s); ## x1 = Q * (U \ (L \ (P * [ones(m,1); zeros(n,1)]))); ## x1 = x1(end - n + 1 : end); ## @end group ## @end example ## ## To find the solution of an overdetermined problem needs an estimate of the ## residual error @var{r} and so it is more complex to formulate a minimum norm ## solution using the @code{spaugment} function. ## ## In general the left division operator is more stable and faster than using ## the @code{spaugment} function. ## @seealso{mldivide} ## @end deftypefn function s = spaugment (A, c) if (nargin < 2) if (issparse (A)) c = max (max (abs (A))) / 1000; else if (ndims (A) != 2) error ("spaugment: A must be a 2-D matrix"); else c = max (abs (A(:))) / 1000; endif endif elseif (! isscalar (c)) error ("spaugment: C must be a scalar"); endif [m, n] = size (A); s = [ c * speye(m, m), A; A', sparse(n, n)]; endfunction %!testif HAVE_UMFPACK %! m = 11; n = 10; mn = max (m ,n); %! A = spdiags ([ones(mn,1), 10*ones(mn,1), -ones(mn,1)],[-1,0,1], m, n); %! x0 = A \ ones (m,1); %! s = spaugment (A); %! [L, U, P, Q] = lu (s); %! x1 = Q * (U \ (L \ (P * [ones(m,1); zeros(n,1)]))); %! x1 = x1(end - n + 1 : end); %! assert (x1, x0, 1e-6);