octave-nkf: libinterp/corefcn/sqrtm.cc annotate

annotate libinterp/corefcn/sqrtm.cc @ 20587:f90c8372b7ba

eliminate many more simple uses of error_state * Cell.cc, __ichol__.cc, __ilu__.cc, balance.cc, bsxfun.cc, colloc.cc, det.cc, dlmread.cc, dynamic-ld.cc, eig.cc, fft.cc, fft2.cc, fftn.cc, gcd.cc, getgrent.cc, getpwent.cc, givens.cc, hess.cc, input.cc, levenshtein.cc, load-path.cc, lookup.cc, ls-mat-ascii.cc, ls-mat4.cc, lsode.cc, lu.cc, max.cc, md5sum.cc, mex.cc, pager.cc, pinv.cc, pr-output.cc, qz.cc, schur.cc, sparse.cc, sqrtm.cc, str2double.cc, strfns.cc, sub2ind.cc, sysdep.cc, time.cc, toplev.cc, tril.cc, tsearch.cc, typecast.cc, __init_gnuplot__.cc, __magick_read__.cc, __osmesa_print__.cc, amd.cc, audiodevinfo.cc, dmperm.cc, fftw.cc, symrcm.cc, ov-base-diag.cc, ov-base-sparse.cc, ov-base.cc, ov-bool-sparse.cc, ov-builtin.cc, ov-complex.cc, ov-cx-diag.cc, ov-cx-mat.cc, ov-cx-sparse.cc, ov-fcn-handle.cc, ov-fcn-inline.cc, ov-float.cc, ov-flt-complex.cc, ov-flt-cx-diag.cc, ov-flt-cx-mat.cc, ov-flt-re-diag.cc, ov-flt-re-mat.cc, ov-lazy-idx.cc, ov-mex-fcn.cc, ov-perm.cc, ov-range.cc, ov-re-diag.cc, ov-re-mat.cc, ov-re-sparse.cc, ov-scalar.cc, ov-str-mat.cc, op-bm-b.cc, op-bm-bm.cc, op-sbm-b.cc, op-sbm-bm.cc, op-str-m.cc, op-str-s.cc, oct-parse.in.yy, pt-cbinop.cc, pt-colon.cc, pt-decl.cc, pt-exp.cc, pt-id.cc, pt-misc.cc, pt-select.cc, pt-unop.cc: Eliminate simple uses of error_state.

author	John W. Eaton <jwe@octave.org>
date	Mon, 05 Oct 2015 19:29:36 -0400
parents	4197fc428c7d
children

rev	line source
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	1 /*
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	2
19731 4197fc428c7d maint: Update copyright notices for 2015. John W. Eaton <jwe@octave.org> parents: 19074 diff changeset	3 Copyright (C) 2001-2015 Ross Lippert and Paul Kienzle
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	4 Copyright (C) 2010 VZLU Prague
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	5
7016 93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5823 diff changeset	6 This file is part of Octave.
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	7
7016 93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5823 diff changeset	8 Octave is free software; you can redistribute it and/or modify it
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5823 diff changeset	9 under the terms of the GNU General Public License as published by the
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5823 diff changeset	10 Free Software Foundation; either version 3 of the License, or (at your
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5823 diff changeset	11 option) any later version.
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5823 diff changeset	12
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5823 diff changeset	13 Octave is distributed in the hope that it will be useful, but WITHOUT
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5823 diff changeset	14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5823 diff changeset	15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5823 diff changeset	16 for more details.
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	17
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	18 You should have received a copy of the GNU General Public License
7016 93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5823 diff changeset	19 along with Octave; see the file COPYING. If not, see
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5823 diff changeset	20 <http://www.gnu.org/licenses/>.
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	21
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	22 */
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	23
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	24 #ifdef HAVE_CONFIG_H
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	25 #include <config.h>
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	26 #endif
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	27
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	28 #include <float.h>
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	29
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	30 #include "CmplxSCHUR.h"
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	31 #include "fCmplxSCHUR.h"
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	32 #include "lo-ieee.h"
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	33 #include "lo-mappers.h"
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	34 #include "oct-norm.h"
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	35
15039 e753177cde93 maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory Rik <rik@octave.org> parents: 14501 diff changeset	36 #include "defun.h"
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	37 #include "error.h"
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	38 #include "gripes.h"
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	39 #include "utils.h"
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	40 #include "xnorm.h"
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	41
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	42 template <class Matrix>
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	43 static void
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	44 sqrtm_utri_inplace (Matrix& T)
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	45 {
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	46 typedef typename Matrix::element_type element_type;
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	47
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	48 const element_type zero = element_type ();
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	49
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	50 bool singular = false;
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	51
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	52 // The following code is equivalent to this triple loop:
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	53 //
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	54 // n = rows (T);
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	55 // for j = 1:n
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	56 // T(j,j) = sqrt (T(j,j));
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	57 // for i = j-1:-1:1
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	58 // T(i,j) /= (T(i,i) + T(j,j));
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	59 // k = 1:i-1;
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	60 // T(k,j) -= T(k,i) * T(i,j);
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	61 // endfor
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	62 // endfor
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	63 //
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	64 // this is an in-place, cache-aligned variant of the code
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	65 // given in Higham's paper.
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	66
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	67 const octave_idx_type n = T.rows ();
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	68 element_type *Tp = T.fortran_vec ();
5275 23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe] jwe parents: 4331 diff changeset	69 for (octave_idx_type j = 0; j < n; j++)
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	70 {
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	71 element_type colj = Tp + nj;
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	72 if (colj[j] != zero)
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	73 colj[j] = sqrt (colj[j]);
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	74 else
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	75 singular = true;
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	76
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	77 for (octave_idx_type i = j-1; i >= 0; i--)
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	78 {
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	79 const element_type coli = Tp + ni;
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	80 const element_type colji = colj[i] /= (coli[i] + colj[j]);
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	81 for (octave_idx_type k = 0; k < i; k++)
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	82 colj[k] -= coli[k] * colji;
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9064 diff changeset	83 }
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	84 }
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	85
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	86 if (singular)
10629 b7c4954e1c00 add specific ID for sqrtm singularity warning Jaroslav Hajek <highegg@gmail.com> parents: 10619 diff changeset	87 warning_with_id ("Octave:sqrtm:SingularMatrix",
b7c4954e1c00 add specific ID for sqrtm singularity warning Jaroslav Hajek <highegg@gmail.com> parents: 10619 diff changeset	88 "sqrtm: matrix is singular, may not have a square root");
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	89 }
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	90
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	91 template <class Matrix, class ComplexMatrix, class ComplexSCHUR>
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	92 static octave_value
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	93 do_sqrtm (const octave_value& arg)
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	94 {
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	95
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	96 octave_value retval;
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	97
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	98 MatrixType mt = arg.matrix_type ();
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	99
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	100 bool iscomplex = arg.is_complex_type ();
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	101
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	102 typedef typename Matrix::element_type real_type;
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	103
18099 6a71e5030df5 Follow coding convention of defining and initializing only 1 variable per line in liboctinterp. Rik <rik@octave.org> parents: 17787 diff changeset	104 real_type cutoff = 0;
6a71e5030df5 Follow coding convention of defining and initializing only 1 variable per line in liboctinterp. Rik <rik@octave.org> parents: 17787 diff changeset	105 real_type one = 1;
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	106 real_type eps = std::numeric_limits<real_type>::epsilon ();
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	107
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	108 if (! iscomplex)
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	109 {
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	110 Matrix x = octave_value_extract<Matrix> (arg);
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	111
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	112 if (mt.is_unknown ()) // if type is not known, compute it now.
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	113 arg.matrix_type (mt = MatrixType (x));
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	114
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	115 switch (mt.type ())
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9064 diff changeset	116 {
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	117 case MatrixType::Upper:
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	118 case MatrixType::Diagonal:
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	119 if (! x.diag ().any_element_is_negative ())
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	120 {
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	121 // Do it in real arithmetic.
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	122 sqrtm_utri_inplace (x);
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	123 retval = x;
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	124 retval.matrix_type (mt);
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	125 }
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	126 else
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	127 iscomplex = true;
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	128 break;
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	129
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	130 case MatrixType::Lower:
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	131 if (! x.diag ().any_element_is_negative ())
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	132 {
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	133 x = x.transpose ();
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	134 sqrtm_utri_inplace (x);
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	135 retval = x.transpose ();
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	136 retval.matrix_type (mt);
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	137 }
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	138 else
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	139 iscomplex = true;
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	140 break;
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	141
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	142 default:
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	143 iscomplex = true;
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	144 break;
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	145 }
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	146
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	147 if (iscomplex)
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	148 cutoff = 10 * x.rows () * eps * xnorm (x, one);
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	149 }
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	150
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	151 if (iscomplex)
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	152 {
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	153 ComplexMatrix x = octave_value_extract<ComplexMatrix> (arg);
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	154
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	155 if (mt.is_unknown ()) // if type is not known, compute it now.
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	156 arg.matrix_type (mt = MatrixType (x));
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	157
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	158 switch (mt.type ())
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	159 {
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	160 case MatrixType::Upper:
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	161 case MatrixType::Diagonal:
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	162 sqrtm_utri_inplace (x);
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	163 retval = x;
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	164 retval.matrix_type (mt);
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	165 break;
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	166
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	167 case MatrixType::Lower:
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	168 x = x.transpose ();
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	169 sqrtm_utri_inplace (x);
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	170 retval = x.transpose ();
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	171 retval.matrix_type (mt);
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	172 break;
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	173
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	174 default:
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	175 {
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	176 ComplexMatrix u;
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	177
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	178 do
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	179 {
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	180 ComplexSCHUR schur (x, std::string (), true);
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	181 x = schur.schur_matrix ();
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	182 u = schur.unitary_matrix ();
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	183 }
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	184 while (0); // schur no longer needed.
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	185
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	186 sqrtm_utri_inplace (x);
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	187
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	188 x = u * x; // original x no longer needed.
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	189 ComplexMatrix res = xgemm (x, u, blas_no_trans, blas_conj_trans);
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	190
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	191 if (cutoff > 0 && xnorm (imag (res), one) <= cutoff)
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	192 retval = real (res);
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	193 else
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	194 retval = res;
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	195 }
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	196 break;
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9064 diff changeset	197 }
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	198 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	199
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	200 return retval;
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	201 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	202
15039 e753177cde93 maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory Rik <rik@octave.org> parents: 14501 diff changeset	203 DEFUN (sqrtm, args, nargout,
17787 175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	204 "-- texinfo --\n\
15039 e753177cde93 maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory Rik <rik@octave.org> parents: 14501 diff changeset	205 @deftypefn {Built-in Function} {@var{s} =} sqrtm (@var{A})\n\
e753177cde93 maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory Rik <rik@octave.org> parents: 14501 diff changeset	206 @deftypefnx {Built-in Function} {[@var{s}, @var{error_estimate}] =} sqrtm (@var{A})\n\
11553 01f703952eff Improve docstrings for functions in DLD-FUNCTIONS directory. Rik <octave@nomad.inbox5.com> parents: 11523 diff changeset	207 Compute the matrix square root of the square matrix @var{A}.\n\
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	208 \n\
19074 0850b5212619 doc: Add @nospell macro around proper names in documentation. Rik <rik@octave.org> parents: 18099 diff changeset	209 Ref: @nospell{N.J. Higham}. @cite{A New sqrtm for @sc{matlab}}. Numerical\n\
10791 3140cb7a05a1 Add spellchecker scripts for Octave and run spellcheck of documentation Rik <octave@nomad.inbox5.com> parents: 10629 diff changeset	210 Analysis Report No. 336, Manchester @nospell{Centre} for Computational\n\
3140cb7a05a1 Add spellchecker scripts for Octave and run spellcheck of documentation Rik <octave@nomad.inbox5.com> parents: 10629 diff changeset	211 Mathematics, Manchester, England, January 1999.\n\
8286 6f2d95255911 fix @seealso references to point to existing anchors Thorsten Meyer <thorsten.meyier@gmx.de> parents: 7795 diff changeset	212 @seealso{expm, logm}\n\
5642 2618a0750ae6 [project @ 2006-03-06 21:26:48 by jwe] jwe parents: 5307 diff changeset	213 @end deftypefn")
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	214 {
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	215 octave_value_list retval;
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	216
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	217 int nargin = args.length ();
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	218
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	219 if (nargin != 1)
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	220 {
5823 080c08b192d8 [project @ 2006-05-19 05:32:17 by jwe] jwe parents: 5775 diff changeset	221 print_usage ();
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	222 return retval;
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	223 }
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	224
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	225 octave_value arg = args(0);
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	226
5275 23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe] jwe parents: 4331 diff changeset	227 octave_idx_type n = arg.rows ();
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe] jwe parents: 4331 diff changeset	228 octave_idx_type nc = arg.columns ();
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	229
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	230 if (n != nc \|\| arg.ndims () > 2)
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	231 {
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	232 gripe_square_matrix_required ("sqrtm");
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	233 return retval;
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	234 }
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	235
14022 de90542b7afc Return retval(1) before retval(0) to avoid re-sizing call on retval. Rik <octave@nomad.inbox5.com> parents: 12926 diff changeset	236 if (nargout > 1)
de90542b7afc Return retval(1) before retval(0) to avoid re-sizing call on retval. Rik <octave@nomad.inbox5.com> parents: 12926 diff changeset	237 {
de90542b7afc Return retval(1) before retval(0) to avoid re-sizing call on retval. Rik <octave@nomad.inbox5.com> parents: 12926 diff changeset	238 retval.resize (1, 2);
de90542b7afc Return retval(1) before retval(0) to avoid re-sizing call on retval. Rik <octave@nomad.inbox5.com> parents: 12926 diff changeset	239 retval(2) = -1.0;
de90542b7afc Return retval(1) before retval(0) to avoid re-sizing call on retval. Rik <octave@nomad.inbox5.com> parents: 12926 diff changeset	240 }
de90542b7afc Return retval(1) before retval(0) to avoid re-sizing call on retval. Rik <octave@nomad.inbox5.com> parents: 12926 diff changeset	241
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	242 if (arg.is_diag_matrix ())
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	243 // sqrtm of a diagonal matrix is just sqrt.
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	244 retval(0) = arg.sqrt ();
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	245 else if (arg.is_single_type ())
17787 175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	246 retval(0) = do_sqrtm<FloatMatrix, FloatComplexMatrix, FloatComplexSCHUR>
175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	247 (arg);
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	248 else if (arg.is_numeric_type ())
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	249 retval(0) = do_sqrtm<Matrix, ComplexMatrix, ComplexSCHUR> (arg);
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	250
20587 f90c8372b7ba eliminate many more simple uses of error_state John W. Eaton <jwe@octave.org> parents: 19731 diff changeset	251 if (nargout > 1)
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	252 {
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	253 // This corresponds to generic code
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	254 //
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	255 // norm (s*s - x, "fro") / norm (x, "fro");
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	256
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	257 octave_value s = retval(0);
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	258 retval(1) = xfrobnorm (s*s - arg) / xfrobnorm (arg);
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	259 }
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	260
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	261 return retval;
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	262 }
12829 ab529124f516 codesprint: Wrote 3 tests for sqrtm.cc Jordi Gutiérrez Hermoso <jordigh@gmail.com> parents: 12584 diff changeset	263
ab529124f516 codesprint: Wrote 3 tests for sqrtm.cc Jordi Gutiérrez Hermoso <jordigh@gmail.com> parents: 12584 diff changeset	264 /*
ab529124f516 codesprint: Wrote 3 tests for sqrtm.cc Jordi Gutiérrez Hermoso <jordigh@gmail.com> parents: 12584 diff changeset	265 %!assert (sqrtm (2ones (2)), ones (2), 3eps)
ab529124f516 codesprint: Wrote 3 tests for sqrtm.cc Jordi Gutiérrez Hermoso <jordigh@gmail.com> parents: 12584 diff changeset	266
ab529124f516 codesprint: Wrote 3 tests for sqrtm.cc Jordi Gutiérrez Hermoso <jordigh@gmail.com> parents: 12584 diff changeset	267 ## The following two tests are from the reference in the docstring above.
12926 838ada8779a5 Fix brain-o in sqrtm test Jordi Gutiérrez Hermoso <jordigh@gmail.com> parents: 12917 diff changeset	268 %!test
838ada8779a5 Fix brain-o in sqrtm test Jordi Gutiérrez Hermoso <jordigh@gmail.com> parents: 12917 diff changeset	269 %! x = [0 1; 0 0];
14501 60e5cf354d80 Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	270 %! assert (any (isnan (sqrtm (x))(:)));
12829 ab529124f516 codesprint: Wrote 3 tests for sqrtm.cc Jordi Gutiérrez Hermoso <jordigh@gmail.com> parents: 12584 diff changeset	271
ab529124f516 codesprint: Wrote 3 tests for sqrtm.cc Jordi Gutiérrez Hermoso <jordigh@gmail.com> parents: 12584 diff changeset	272 %!test
14501 60e5cf354d80 Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	273 %! x = eye (4); x(2,2) = x(3,3) = 2^-26; x(1,4) = 1;
60e5cf354d80 Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	274 %! z = eye (4); z(2,2) = z(3,3) = 2^-13; z(1,4) = 0.5;
60e5cf354d80 Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	275 %! [y, err] = sqrtm (x);
60e5cf354d80 Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	276 %! assert (y, z);
17336 b81b9d079515 Use '##' for comments which stand alone on a line. Rik <rik@octave.org> parents: 15195 diff changeset	277 %! assert (err, 0); # Yes, this one has to hold exactly
12829 ab529124f516 codesprint: Wrote 3 tests for sqrtm.cc Jordi Gutiérrez Hermoso <jordigh@gmail.com> parents: 12584 diff changeset	278 */

Mercurial > octave-nkf

annotate libinterp/corefcn/sqrtm.cc @ 20587:f90c8372b7ba