42 template <
class Matrix>
48 const element_type zero = element_type ();
50 bool singular =
false;
71 element_type *colj = Tp + n*j;
73 colj[j] =
sqrt (colj[j]);
79 const element_type *coli = Tp + n*i;
80 const element_type colji = colj[i] /= (coli[i] + colj[j]);
82 colj[k] -= coli[k] * colji;
88 "sqrtm: matrix is singular, may not have a square root");
91 template <
class Matrix,
class ComplexMatrix,
class ComplexSCHUR>
104 real_type cutoff = 0;
106 real_type
eps = std::numeric_limits<real_type>::epsilon ();
148 cutoff = 10 * x.
rows () * eps *
xnorm (x, one);
191 if (cutoff > 0 &&
xnorm (
imag (res), one) <= cutoff)
203 DEFUN (sqrtm, args, nargout,
205 @deftypefn {Built-in Function} {@var{s} =} sqrtm (@var{A})\n\
206 @deftypefnx {Built-in Function} {[@var{s}, @var{error_estimate}] =} sqrtm (@var{A})\n\
207 Compute the matrix square root of the square matrix @var{A}.\n\
209 Ref: @nospell{N.J. Higham}. @cite{A New sqrtm for @sc{matlab}}. Numerical\n\
210 Analysis Report No. 336, Manchester @nospell{Centre} for Computational\n\
211 Mathematics, Manchester, England, January 1999.\n\
212 @seealso{expm, logm}\n\
217 int nargin = args.
length ();
230 if (n != nc || arg.
ndims () > 2)
244 retval(0) = arg.
sqrt ();
246 retval(0) = do_sqrtm<FloatMatrix, FloatComplexMatrix, FloatComplexSCHUR>
249 retval(0) = do_sqrtm<Matrix, ComplexMatrix, ComplexSCHUR> (
arg);
Matrix diag(octave_idx_type k=0) const
void warning_with_id(const char *id, const char *fmt,...)
bool any_element_is_negative(bool=false) const
octave_idx_type rows(void) const
static octave_value do_sqrtm(const octave_value &arg)
bool is_unknown(void) 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)
ComplexMatrix schur_matrix(void) const
void gripe_square_matrix_required(const char *name)
octave_idx_type rows(void) const
octave_idx_type columns(void) const
int type(bool quiet=true)
bool is_complex_type(void) const
Matrix transpose(void) const
OCTAVE_API double xnorm(const ColumnVector &x, double p)
ComplexMatrix transpose(void) const
ComplexMatrix unitary_matrix(void) const
MatrixType matrix_type(void) const
ComplexMatrix xgemm(const ComplexMatrix &a, const ComplexMatrix &b, blas_trans_type transa, blas_trans_type transb)
OCTAVE_API double xfrobnorm(const Matrix &x)
void resize(octave_idx_type n, const octave_value &rfv=octave_value())
ColumnVector imag(const ComplexColumnVector &a)
static void sqrtm_utri_inplace(Matrix &T)
ComplexMatrix octave_value_extract< ComplexMatrix >(const octave_value &v)
const T * fortran_vec(void) const
bool is_single_type(void) const
ColumnVector real(const ComplexColumnVector &a)
Matrix octave_value_extract< Matrix >(const octave_value &v)
octave_value sqrt(void) const
bool is_diag_matrix(void) const
F77_RET_T const double * x