--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/linear-algebra/linsolve.m	Thu Sep 26 08:30:26 2013 -0700
@@ -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);
+