Mercurial > octave-nkf
diff liboctave/UMFPACK/UMFPACK/OCTAVE/umfpack_solve.m @ 5164:57077d0ddc8e
[project @ 2005-02-25 19:55:24 by jwe]
author | jwe |
---|---|
date | Fri, 25 Feb 2005 19:55:28 +0000 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/UMFPACK/UMFPACK/OCTAVE/umfpack_solve.m Fri Feb 25 19:55:28 2005 +0000 @@ -0,0 +1,97 @@ +function x = umfpack_solve (arg1, op, arg2, Control) +% UMFPACK_SOLVE +% +% x = umfpack_solve (A, '\', b, Control) +% x = umfpack_solve (b, '/', A, Control) +% +% Computes x = A\b, or b/A, where A is square. Uses UMFPACK if A is sparse. +% The Control argument is optional. +% +% See also umfpack, umfpack_make, umfpack_details, umfpack_report, +% and umfpack_simple. + +% UMFPACK Version 4.3 (Jan. 16, 2004), Copyright (c) 2004 by Timothy A. +% Davis. All Rights Reserved. Type umfpack_details for License. + +%------------------------------------------------------------------------------- +% check inputs and get default control parameters +%------------------------------------------------------------------------------- + +if (op == '\\') + A = arg1 ; + b = arg2 ; +elseif (op == '/') + A = arg2 ; + b = arg1 ; +else + help umfack_solve + error ('umfpack_solve: unrecognized operator') ; +end + +[m n] = size (A) ; +if (m ~= n) + help umfpack_solve + error ('umfpack_solve: A must be square') ; +end + +[m1 n1] = size (b) ; +if ((op == '\\' & n ~= m1) | (op == '/' & n1 ~= m)) + help umfpack_solve + error ('umfpack_solve: b has the wrong dimensions') ; +end + +if (nargin < 4) + Control = umfpack ; +end + +%------------------------------------------------------------------------------- +% solve the system +%------------------------------------------------------------------------------- + +if (op == '\\') + + if (~issparse (A)) + + % A is not sparse, so just use MATLAB + x = A\b ; + + elseif (n1 == 1 & ~issparse (b)) + + % the UMFPACK '\' requires b to be a dense column vector + x = umfpack (A, '\\', b, Control) ; + + else + + % factorize with UMFPACK and do the forward/back solves in MATLAB + [L, U, P, Q, R] = umfpack (A, Control) ; + keyboard + x = Q * (U \ (L \ (P * (R \ b)))) ; + + end + +else + + if (~issparse (A)) + + % A is not sparse, so just use MATLAB + x = b/A ; + + elseif (m1 == 1 & ~issparse (b)) + + % the UMFPACK '/' requires b to be a dense column vector + x = umfpack (b, '/', A, Control) ; + + else + + % factorize with UMFPACK and do the forward/back solves in MATLAB + % this mimics the behavior of x = b/A, except for the row scaling + [L, U, P, Q, R] = umfpack (A.', Control) ; + x = (Q * (U \ (L \ (P * (R \ (b.')))))).' ; + + % an alternative method: + % [L, U, P, Q, r] = umfpack (A, Control) ; + % x = (R \ (P' * (L.' \ (U.' \ (Q' * b.'))))).' ; + + end + +end