Mercurial > hg > octave-lyh
changeset 17521:c3a3532e3d98
linsolve.m: Add new function for Matlab compatibility.
* scripts/linear-algebra/linsolve.m: New function.
* scripts/linear-algebra/module.mk: Add linsolve.m to build system.
* NEWS: Announce new function.
author | Nir Krakauer < nkrakauer@ccny.cuny.edu> |
---|---|
date | Thu, 26 Sep 2013 08:30:26 -0700 |
parents | cdeadf62663f |
children | b66f068e4468 |
files | NEWS scripts/linear-algebra/linsolve.m scripts/linear-algebra/module.mk |
diffstat | 3 files changed, 110 insertions(+), 18 deletions(-) [+] |
line wrap: on
line diff
--- a/NEWS +++ b/NEWS @@ -265,25 +265,24 @@ inputdlg msgbox warndlg ** Other new functions added in 3.8.0: - + atan2d erfi jit_startcnt - base64_decode expint lines - base64_encode findfigs polyeig - betaincinv flintmax prefdir - built_in_docstrings_file fminsearch preferences - cmpermute gallery readline_re_read_init_file - cmunique gco readline_read_init_file - colorcube hdl2struct rgbplot - copyobj history_save save_default_options - dawson imformats shrinkfaces - dblist importdata splinefit - debug_jit isaxes stemleaf - desktop iscolormap strjoin - doc_cache_create isequaln struct2hdl - ellipj jit_debug tetramesh - ellipke jit_enable waterfall - erfcinv - + base64_decode expint lines + base64_encode findfigs linsolve + betaincinv flintmax polyeig + built_in_docstrings_file fminsearch prefdir + cmpermute gallery preferences + cmunique gco readline_re_read_init_file + colorcube hdl2struct readline_read_init_file + copyobj history_save rgbplot + dawson imformats save_default_options + dblist importdata shrinkfaces + debug_jit isaxes splinefit + desktop iscolormap stemleaf + doc_cache_create isequaln strjoin + ellipj jit_debug struct2hdl + ellipke jit_enable tetramesh + erfcinv waterfall ** Deprecated functions.
new file mode 100644 --- /dev/null +++ b/scripts/linear-algebra/linsolve.m @@ -0,0 +1,92 @@ +## Copyright (C) 2013 Nir Krakauer +## +## This program 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. +## +## This program 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 this program; If not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn{Function File}{[@var{X}, @var{R}] =} csaps(@var{A}, @var{B}, @var{options}=[]) +## +## Solve a linear system@* +## +## With no options, this is the same as @code{A \ B} +## +## Possible option fields (set to true/false): +## @table @asis +## @item @var{LT} +## A is lower triangular +## @item @var{UT} +## A is upper triangular +## @item @var{UHESS} +## A is upper Hessenberg (currently makes no difference) +## @item @var{SYM} +## A is symmetric (currently makes no difference) +## @item @var{POSDEF} +## A is positive definite +## @item @var{RECT} +## A is general rectangular (currently makes no difference) +## @item @var{TRANSA} +## Compute @code{transpose(A) \ B} +## @end table +## +## The optional second output @var{R} is the inverse condition number of @var{A} (zero if matrix is singular) +## @end deftypefn + +## Author: Nir Krakauer <nkrakauer@ccny.cuny.edu> + +function [X, R] = linsolve(A, B, options) + +trans_A = false; + +#process any options +if nargin > 2 + if ~isstruct(options) + error('Third input must be a structure') + endif + if isfield(options, 'TRANSA') && options.TRANSA + trans_A = true; + A = A'; + endif + if isfield(options, 'POSDEF') && options.POSDEF + A = matrix_type (A, 'positive definite'); + endif + if isfield(options, 'LT') && options.LT + if trans_A + A = matrix_type (A, 'upper'); + else + A = matrix_type (A, 'lower'); + endif + endif + if isfield(options, 'UT') && options.UT + if trans_A + A = matrix_type (A, 'lower'); + else + A = matrix_type (A, 'upper'); + endif + endif +endif + +X = A \ B; + +if nargout > 1 + if issquare(A) + R = 1 ./ cond(A); + else + R = 0; + endif +endif + +%!shared n, A, B, x, opts +%! n = 4; A = triu(rand(n)); x = rand(n, 1); B = A' * x; +%! opts.UT = true; opts.TRANSA = true; +%!assert (linsolve(A, B, opts), A' \ B); +