58 if (U.is_square () && fact.
regular ())
64 DEFUN (lu, args, nargout,
66 @deftypefn {Built-in Function} {[@var{L}, @var{U}] =} lu (@var{A})\n\
67 @deftypefnx {Built-in Function} {[@var{L}, @var{U}, @var{P}] =} lu (@var{A})\n\
68 @deftypefnx {Built-in Function} {[@var{L}, @var{U}, @var{P}, @var{Q}] =} lu (@var{S})\n\
69 @deftypefnx {Built-in Function} {[@var{L}, @var{U}, @var{P}, @var{Q}, @var{R}] =} lu (@var{S})\n\
70 @deftypefnx {Built-in Function} {[@dots{}] =} lu (@var{S}, @var{thres})\n\
71 @deftypefnx {Built-in Function} {@var{y} =} lu (@dots{})\n\
72 @deftypefnx {Built-in Function} {[@dots{}] =} lu (@dots{}, \"vector\")\n\
73 @cindex LU decomposition\n\
74 Compute the LU@tie{}decomposition of @var{A}.\n\
76 If @var{A} is full subroutines from @sc{lapack} are used and if @var{A} is\n\
77 sparse then @sc{umfpack} is used.\n\
79 The result is returned in a permuted form, according to the optional return\n\
80 value @var{P}. For example, given the matrix @code{a = [1, 2; 3, 4]},\n\
83 [l, u, p] = lu (@var{a})\n\
108 The matrix is not required to be square.\n\
110 When called with two or three output arguments and a spare input matrix,\n\
111 @code{lu} does not attempt to perform sparsity preserving column\n\
112 permutations. Called with a fourth output argument, the sparsity\n\
113 preserving column transformation @var{Q} is returned, such that\n\
114 @code{@var{P} * @var{A} * @var{Q} = @var{L} * @var{U}}.\n\
116 Called with a fifth output argument and a sparse input matrix,\n\
117 @code{lu} attempts to use a scaling factor @var{R} on the input matrix\n\
119 @code{@var{P} * (@var{R} \\ @var{A}) * @var{Q} = @var{L} * @var{U}}.\n\
120 This typically leads to a sparser and more stable factorization.\n\
122 An additional input argument @var{thres}, that defines the pivoting\n\
123 threshold can be given. @var{thres} can be a scalar, in which case\n\
124 it defines the @sc{umfpack} pivoting tolerance for both symmetric and\n\
125 unsymmetric cases. If @var{thres} is a 2-element vector, then the first\n\
126 element defines the pivoting tolerance for the unsymmetric @sc{umfpack}\n\
127 pivoting strategy and the second for the symmetric strategy. By default,\n\
128 the values defined by @code{spparms} are used ([0.1, 0.001]).\n\
130 Given the string argument @qcode{\"vector\"}, @code{lu} returns the values\n\
131 of @var{P} and @var{Q} as vector values, such that for full matrix,\n\
132 @code{@var{A} (@var{P},:) = @var{L} * @var{U}}, and @code{@var{R}(@var{P},:)\n\
133 * @var{A} (:, @var{Q}) = @var{L} * @var{U}}.\n\
135 With two output arguments, returns the permuted forms of the upper and\n\
136 lower triangular matrices, such that @code{@var{A} = @var{L} * @var{U}}.\n\
137 With one output argument @var{y}, then the matrix returned by the @sc{lapack}\n\
138 routines is returned. If the input matrix is sparse then the matrix @var{L}\n\
139 is embedded into @var{U} to give a return value similar to the full case.\n\
140 For both full and sparse matrices, @code{lu} loses the permutation\n\
142 @seealso{luupdate, ilu, chol, hess, qr, qz, schur, svd}\n\
146 int nargin = args.
length ();
147 bool issparse = (nargin > 0 && args(0).is_sparse_type ());
148 bool scale = (nargout == 5);
150 if (nargin < 1 || (issparse && (nargin > 3 || nargout > 5))
151 || (!issparse && (nargin > 2 || nargout > 3)))
163 if (args(n).is_string ())
165 std::string tmp = args(n++).string_value ();
167 if (tmp.compare (
"vector") == 0)
170 error (
"lu: unrecognized string argument");
174 Matrix tmp = args(n++).matrix_value ();
179 error (
"lu: can not define pivoting threshold THRES for full matrices");
180 else if (tmp.
nelem () == 1)
186 else if (tmp.
nelem () == 2)
189 error (
"lu: expecting 2-element vector for THRES");
199 int arg_is_empty =
empty_arg (
"lu", nr, nc);
203 if (arg_is_empty < 0)
205 else if (arg_is_empty > 0)
220 SparseLU fact (m, Qinit, thres,
false,
true);
223 retval(0) = fact.
Y ();
244 retval(2) = fact.
Pr_vec();
258 retval(4) = fact.
R ();
262 retval(3) = fact.
Pc_vec ();
263 retval(2) = fact.
Pr_vec ();
267 retval(3) = fact.
Pc_mat ();
268 retval(2) = fact.
Pr_mat ();
292 retval(0) = fact.
Y ();
314 retval(2) = fact.
Pr_vec();
328 retval(4) = fact.
R ();
332 retval(3) = fact.
Pc_vec ();
333 retval(2) = fact.
Pr_vec ();
337 retval(3) = fact.
Pc_mat ();
338 retval(2) = fact.
Pr_mat ();
352 if (arg_is_empty < 0)
354 else if (arg_is_empty > 0)
371 retval(0) = fact.
Y ();
387 retval(2) = fact.
P_vec ();
389 retval(2) = fact.
P ();
409 retval(0) = fact.
Y ();
425 retval(2) = fact.
P_vec ();
427 retval(2) = fact.
P ();
450 retval(0) = fact.
Y ();
466 retval(2) = fact.
P_vec ();
468 retval(2) = fact.
P ();
488 retval(0) = fact.
Y ();
504 retval(2) = fact.
P_vec ();
506 retval(2) = fact.
P ();
602 DEFUN (luupdate, args, ,
604 @deftypefn {Built-in Function} {[@var{L}, @var{U}] =} luupdate (@var{L}, @var{U}, @var{x}, @var{y})\n\
605 @deftypefnx {Built-in Function} {[@var{L}, @var{U}, @var{P}] =} luupdate (@var{L}, @var{U}, @var{P}, @var{x}, @var{y})\n\
606 Given an LU@tie{}factorization of a real or complex matrix\n\
607 @w{@var{A} = @var{L}*@var{U}}, @var{L}@tie{}lower unit trapezoidal and\n\
608 @var{U}@tie{}upper trapezoidal, return the LU@tie{}factorization\n\
609 of @w{@var{A} + @var{x}*@var{y}.'}, where @var{x} and @var{y} are\n\
610 column vectors (rank-1 update) or matrices with equal number of columns\n\
613 Optionally, row-pivoted updating can be used by supplying a row permutation\n\
614 (pivoting) matrix @var{P}; in that case, an updated permutation matrix is\n\
615 returned. Note that if @var{L}, @var{U}, @var{P} is a pivoted\n\
616 LU@tie{}factorization as obtained by @code{lu}:\n\
619 [@var{L}, @var{U}, @var{P}] = lu (@var{A});\n\
623 then a factorization of @tcode{@var{A}+@var{x}*@var{y}.'} can be obtained\n\
627 [@var{L1}, @var{U1}] = lu (@var{L}, @var{U}, @var{P}*@var{x}, @var{y})\n\
634 [@var{L1}, @var{U1}, @var{P1}] = lu (@var{L}, @var{U}, @var{P}, @var{x}, @var{y})\n\
637 The first form uses the unpivoted algorithm, which is faster, but less\n\
638 stable. The second form uses a slower pivoted algorithm, which is more\n\
641 The matrix case is done as a sequence of rank-1 updates; thus, for large\n\
642 enough k, it will be both faster and more accurate to recompute the\n\
643 factorization from scratch.\n\
644 @seealso{lu, cholupdate, qrupdate}\n\
650 bool pivoted = nargin == 5;
652 if (nargin != 4 && nargin != 5)
697 retval(2) = fact.
P ();
715 retval(2) = fact.
P ();
740 retval(2) = fact.
P ();
758 retval(2) = fact.
P ();
765 error (
"luupdate: dimension mismatch");
768 error (
"luupdate: L, U, X, and Y must be numeric");
void update(const ComplexColumnVector &u, const ComplexColumnVector &v)
bool is_real_type(void) const
void gripe_wrong_type_arg(const char *name, const char *s, bool is_error)
const octave_idx_type * col_perm(void) const
octave_idx_type rows(void) const
void resize(octave_idx_type nr, octave_idx_type nc, double rfv=0)
const octave_idx_type * row_perm(void) const
FloatComplexMatrix float_complex_matrix_value(bool frc_str_conv=false) const
OCTINTERP_API void print_usage(void)
octave_idx_type length(void) const
bool is_numeric_type(void) const
#define DEFUN(name, args_name, nargout_name, doc)
void error(const char *fmt,...)
PermMatrix transpose(void) const
bool is_perm_matrix(void) const
PermMatrix Pr_mat(void) const
int empty_arg(const char *, octave_idx_type nr, octave_idx_type nc)
void update(const ColumnVector &u, const ColumnVector &v)
static bool check_lu_dims(const octave_value &l, const octave_value &u, const octave_value &p)
void update_piv(const ComplexColumnVector &u, const ComplexColumnVector &v)
SparseMatrix R(void) const
octave_idx_type nelem(void) const
Number of elements in the array.
static octave_value get_lu_l(const base_lu< MT > &fact)
octave_idx_type columns(void) const
ColumnVector Pr_vec(void) const
bool is_complex_type(void) const
void update_piv(const FloatColumnVector &u, const FloatColumnVector &v)
SparseComplexMatrix sparse_complex_matrix_value(bool frc_str_conv=false) const
Matrix matrix_value(bool frc_str_conv=false) const
void update(const FloatColumnVector &u, const FloatColumnVector &v)
static octave_value get_lu_u(const base_lu< MT > &fact)
ComplexMatrix complex_matrix_value(bool frc_str_conv=false) const
ColumnVector P_vec(void) const
FloatMatrix float_matrix_value(bool frc_str_conv=false) const
void update_piv(const ColumnVector &u, const ColumnVector &v)
PermMatrix Pc_mat(void) const
void update_piv(const FloatComplexColumnVector &u, const FloatComplexColumnVector &v)
void scale(Matrix &m, double x, double y, double z)
ColumnVector Pc_vec(void) const
bool is_undefined(void) const
PermMatrix perm_matrix_value(void) const
SparseMatrix sparse_matrix_value(bool frc_str_conv=false) const
void update(const FloatComplexColumnVector &u, const FloatComplexColumnVector &v)
bool is_single_type(void) const
static PermMatrix eye(octave_idx_type n)
return octave_value(v1.char_array_value().concat(v2.char_array_value(), ra_idx),((a1.is_sq_string()||a2.is_sq_string())? '\'': '"'))
F77_RET_T const double * x
charNDArray min(char d, const charNDArray &m)
void resize(octave_idx_type n, const double &rfv=0)