octave: libinterp/corefcn/sqrtm.cc annotate

annotate libinterp/corefcn/sqrtm.cc @ 22317:6ca3acf5fad8

move some new numeric classes to namespace octave::math * aepbalance.cc, aepbalance.h, chol.cc, chol.h, gepbalance.cc, gepbalance.h, hess.cc, hess.h, lu.cc, lu.h, qr.cc, qr.h, qrp.cc, qrp.h, schur.cc, schur.h, sparse-chol.cc, sparse-chol.h, sparse-dmsolve.cc,sparse-lu.cc, sparse-lu.h, sparse-qr.cc, sparse-qr.h, svd.cc, svd.h: Move classes to namespace octave::math. * __luinc__.cc, __qp__.cc, balance.cc, hess.cc, lu.cc, qz.cc, schur.cc, sqrtm.cc, svd.cc, chol.cc, dmperm.cc, qr.cc, lex.h, CMatrix.cc, CSparse.cc, dMatrix.cc, dSparse.cc, fCMatrix.cc, fMatrix.cc, eigs-base.cc, oct-norm.cc: Update for new namespaces.

author	John W. Eaton <jwe@octave.org>
date	Tue, 16 Aug 2016 15:48:22 -0400
parents	112b20240c87
children	bac0d6f07a3e

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
19697 4197fc428c7d maint: Update copyright notices for 2015. John W. Eaton <jwe@octave.org> parents: 19040 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
21724 aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef" John W. Eaton <jwe@octave.org> parents: 21575 diff changeset	24 #if defined (HAVE_CONFIG_H)
21301 40de9f8f23a6 Use '#include "config.h"' rather than <config.h>. Rik <rik@octave.org> parents: 21266 diff changeset	25 # include "config.h"
4331 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
21266 e69eaee28737 make better use of templates for Schur decomposition John W. Eaton <jwe@octave.org> parents: 21200 diff changeset	28 #include "schur.h"
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	29 #include "lo-ieee.h"
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	30 #include "lo-mappers.h"
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	31 #include "oct-norm.h"
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	32
15039 e753177cde93 maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory Rik <rik@octave.org> parents: 14501 diff changeset	33 #include "defun.h"
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	34 #include "error.h"
21100 e39e05d90788 Switch gripe_XXX to either err_XXX or warn_XXX naming scheme. Rik <rik@octave.org> parents: 21017 diff changeset	35 #include "errwarn.h"
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	36 #include "utils.h"
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	37 #include "xnorm.h"
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	38
21139 538b57866b90 consistently use "typename" intead of "class" in template declarations John W. Eaton <jwe@octave.org> parents: 21110 diff changeset	39 template <typename Matrix>
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	40 static void
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	41 sqrtm_utri_inplace (Matrix& T)
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	42 {
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	43 typedef typename Matrix::element_type element_type;
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	44
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	45 const element_type zero = element_type ();
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	46
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	47 bool singular = false;
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	48
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	49 // 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	50 //
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	51 // n = rows (T);
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	52 // for j = 1:n
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	53 // T(j,j) = sqrt (T(j,j));
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	54 // for i = j-1:-1:1
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	55 // 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	56 // k = 1:i-1;
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	57 // 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	58 // endfor
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	59 // endfor
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	60 //
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	61 // 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	62 // given in Higham's paper.
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	63
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	64 const octave_idx_type n = T.rows ();
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	65 element_type *Tp = T.fortran_vec ();
5275 23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe] jwe parents: 4331 diff changeset	66 for (octave_idx_type j = 0; j < n; j++)
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	67 {
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	68 element_type colj = Tp + nj;
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	69 if (colj[j] != zero)
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	70 colj[j] = sqrt (colj[j]);
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	71 else
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	72 singular = true;
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	73
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	74 for (octave_idx_type i = j-1; i >= 0; i--)
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	75 {
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	76 const element_type coli = Tp + ni;
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	77 const element_type colji = colj[i] /= (coli[i] + colj[j]);
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	78 for (octave_idx_type k = 0; k < i; k++)
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	79 colj[k] -= coli[k] * colji;
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9064 diff changeset	80 }
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	81 }
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	82
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	83 if (singular)
10629 b7c4954e1c00 add specific ID for sqrtm singularity warning Jaroslav Hajek <highegg@gmail.com> parents: 10619 diff changeset	84 warning_with_id ("Octave:sqrtm:SingularMatrix",
b7c4954e1c00 add specific ID for sqrtm singularity warning Jaroslav Hajek <highegg@gmail.com> parents: 10619 diff changeset	85 "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	86 }
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	87
21139 538b57866b90 consistently use "typename" intead of "class" in template declarations John W. Eaton <jwe@octave.org> parents: 21110 diff changeset	88 template <typename Matrix, typename ComplexMatrix, typename ComplexSCHUR>
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	89 static octave_value
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	90 do_sqrtm (const octave_value& arg)
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	91 {
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	92
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	93 octave_value retval;
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 MatrixType mt = arg.matrix_type ();
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	96
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	97 bool iscomplex = arg.is_complex_type ();
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	98
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	99 typedef typename Matrix::element_type real_type;
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	100
18100 6a71e5030df5 Follow coding convention of defining and initializing only 1 variable per line in liboctinterp. Rik <rik@octave.org> parents: 17787 diff changeset	101 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	102 real_type one = 1;
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	103 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	104
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	105 if (! iscomplex)
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	106 {
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	107 Matrix x = octave_value_extract<Matrix> (arg);
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	108
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	109 if (mt.is_unknown ()) // if type is not known, compute it now.
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	110 arg.matrix_type (mt = MatrixType (x));
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	111
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	112 switch (mt.type ())
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9064 diff changeset	113 {
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	114 case MatrixType::Upper:
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	115 case MatrixType::Diagonal:
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	116 if (! x.diag ().any_element_is_negative ())
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	117 {
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	118 // Do it in real arithmetic.
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	119 sqrtm_utri_inplace (x);
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	120 retval = x;
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	121 retval.matrix_type (mt);
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	122 }
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	123 else
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	124 iscomplex = true;
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	125 break;
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	126
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	127 case MatrixType::Lower:
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	128 if (! x.diag ().any_element_is_negative ())
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	129 {
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	130 x = x.transpose ();
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	131 sqrtm_utri_inplace (x);
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	132 retval = x.transpose ();
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	133 retval.matrix_type (mt);
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	134 }
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	135 else
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	136 iscomplex = true;
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	137 break;
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	138
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	139 default:
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	140 iscomplex = true;
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	141 break;
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	142 }
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	143
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	144 if (iscomplex)
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	145 cutoff = 10 * x.rows () * eps * xnorm (x, one);
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
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	148 if (iscomplex)
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 ComplexMatrix x = octave_value_extract<ComplexMatrix> (arg);
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	151
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	152 if (mt.is_unknown ()) // if type is not known, compute it now.
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	153 arg.matrix_type (mt = MatrixType (x));
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	154
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	155 switch (mt.type ())
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	156 {
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	157 case MatrixType::Upper:
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	158 case MatrixType::Diagonal:
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	159 sqrtm_utri_inplace (x);
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	160 retval = x;
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	161 retval.matrix_type (mt);
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	162 break;
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	163
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	164 case MatrixType::Lower:
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	165 x = x.transpose ();
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	166 sqrtm_utri_inplace (x);
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	167 retval = x.transpose ();
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	168 retval.matrix_type (mt);
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	169 break;
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	170
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	171 default:
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	172 {
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	173 ComplexMatrix u;
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	174
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	175 do
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	176 {
21266 e69eaee28737 make better use of templates for Schur decomposition John W. Eaton <jwe@octave.org> parents: 21200 diff changeset	177 ComplexSCHUR schur_fact (x, "", true);
e69eaee28737 make better use of templates for Schur decomposition John W. Eaton <jwe@octave.org> parents: 21200 diff changeset	178 x = schur_fact.schur_matrix ();
e69eaee28737 make better use of templates for Schur decomposition John W. Eaton <jwe@octave.org> parents: 21200 diff changeset	179 u = schur_fact.unitary_matrix ();
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	180 }
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	181 while (0); // schur no longer needed.
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	182
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	183 sqrtm_utri_inplace (x);
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	184
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	185 x = u * x; // original x no longer needed.
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	186 ComplexMatrix res = xgemm (x, u, blas_no_trans, blas_conj_trans);
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 if (cutoff > 0 && xnorm (imag (res), one) <= cutoff)
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	189 retval = real (res);
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	190 else
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	191 retval = res;
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	192 }
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	193 break;
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9064 diff changeset	194 }
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	195 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	196
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	197 return retval;
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
15039 e753177cde93 maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory Rik <rik@octave.org> parents: 14501 diff changeset	200 DEFUN (sqrtm, args, nargout,
21966 112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21912 diff changeset	201 doc: /* -- texinfo --
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21912 diff changeset	202 @deftypefn {} {@var{s} =} sqrtm (@var{A})
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21912 diff changeset	203 @deftypefnx {} {[@var{s}, @var{error_estimate}] =} sqrtm (@var{A})
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21912 diff changeset	204 Compute the matrix square root of the square matrix @var{A}.
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21912 diff changeset	205
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21912 diff changeset	206 Ref: @nospell{N.J. Higham}. @cite{A New sqrtm for @sc{matlab}}. Numerical
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21912 diff changeset	207 Analysis Report No. 336, Manchester @nospell{Centre} for Computational
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21912 diff changeset	208 Mathematics, Manchester, England, January 1999.
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21912 diff changeset	209 @seealso{expm, logm}
112b20240c87 move docstrings in C++ files out of C strings and into comments John W. Eaton <jwe@octave.org> parents: 21912 diff changeset	210 @end deftypefn */)
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	211 {
20819 f428cbe7576f eliminate unnecessary uses of nargin John W. Eaton <jwe@octave.org> parents: 20801 diff changeset	212 if (args.length () != 1)
20801 a542a9bf177e eliminate return statements after calls to print_usage John W. Eaton <jwe@octave.org> parents: 20555 diff changeset	213 print_usage ();
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 arg = args(0);
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	216
5275 23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe] jwe parents: 4331 diff changeset	217 octave_idx_type n = arg.rows ();
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe] jwe parents: 4331 diff changeset	218 octave_idx_type nc = arg.columns ();
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	219
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	220 if (n != nc \|\| arg.ndims () > 2)
21110 3d0d84305600 Use err_square_matrix_required more widely. Rik <rik@octave.org> parents: 21100 diff changeset	221 err_square_matrix_required ("sqrtm", "A");
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	222
20892 c07bee629973 2015 Code Sprint: use ovl (). Rik <rik@octave.org> parents: 20853 diff changeset	223 octave_value_list retval (nargout > 1 ? 3 : 1);
c07bee629973 2015 Code Sprint: use ovl (). Rik <rik@octave.org> parents: 20853 diff changeset	224
14022 de90542b7afc Return retval(1) before retval(0) to avoid re-sizing call on retval. Rik <octave@nomad.inbox5.com> parents: 12926 diff changeset	225 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	226 {
20892 c07bee629973 2015 Code Sprint: use ovl (). Rik <rik@octave.org> parents: 20853 diff changeset	227 // FIXME: Octave does not calculate a condition number with respect to
c07bee629973 2015 Code Sprint: use ovl (). Rik <rik@octave.org> parents: 20853 diff changeset	228 // sqrtm. Should this return NaN instead of -1?
14022 de90542b7afc Return retval(1) before retval(0) to avoid re-sizing call on retval. Rik <octave@nomad.inbox5.com> parents: 12926 diff changeset	229 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	230 }
de90542b7afc Return retval(1) before retval(0) to avoid re-sizing call on retval. Rik <octave@nomad.inbox5.com> parents: 12926 diff changeset	231
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	232 if (arg.is_diag_matrix ())
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	233 // sqrtm of a diagonal matrix is just sqrt.
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	234 retval(0) = arg.sqrt ();
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	235 else if (arg.is_single_type ())
21266 e69eaee28737 make better use of templates for Schur decomposition John W. Eaton <jwe@octave.org> parents: 21200 diff changeset	236 retval(0) = do_sqrtm<FloatMatrix, FloatComplexMatrix,
22317 6ca3acf5fad8 move some new numeric classes to namespace octave::math John W. Eaton <jwe@octave.org> parents: 21966 diff changeset	237 octave::math::schur<FloatComplexMatrix> > (arg);
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	238 else if (arg.is_numeric_type ())
22317 6ca3acf5fad8 move some new numeric classes to namespace octave::math John W. Eaton <jwe@octave.org> parents: 21966 diff changeset	239 retval(0) = do_sqrtm<Matrix, ComplexMatrix,
6ca3acf5fad8 move some new numeric classes to namespace octave::math John W. Eaton <jwe@octave.org> parents: 21966 diff changeset	240 octave::math::schur<ComplexMatrix> > (arg);
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	241
20555 f90c8372b7ba eliminate many more simple uses of error_state John W. Eaton <jwe@octave.org> parents: 19697 diff changeset	242 if (nargout > 1)
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	243 {
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	244 // This corresponds to generic code
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	245 //
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	246 // 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	247
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	248 octave_value s = retval(0);
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	249 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	250 }
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	251
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	252 return retval;
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	253 }
12829 ab529124f516 codesprint: Wrote 3 tests for sqrtm.cc Jordi Gutiérrez Hermoso <jordigh@gmail.com> parents: 12584 diff changeset	254
ab529124f516 codesprint: Wrote 3 tests for sqrtm.cc Jordi Gutiérrez Hermoso <jordigh@gmail.com> parents: 12584 diff changeset	255 /*
ab529124f516 codesprint: Wrote 3 tests for sqrtm.cc Jordi Gutiérrez Hermoso <jordigh@gmail.com> parents: 12584 diff changeset	256 %!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	257
ab529124f516 codesprint: Wrote 3 tests for sqrtm.cc Jordi Gutiérrez Hermoso <jordigh@gmail.com> parents: 12584 diff changeset	258 ## 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	259 %!test
21575 bc9aa534bc29 maint: Clean up BIST tests so they don't produce warnings. Rik <rik@octave.org> parents: 21301 diff changeset	260 %! warning ("off", "Octave:sqrtm:SingularMatrix", "local");
12926 838ada8779a5 Fix brain-o in sqrtm test Jordi Gutiérrez Hermoso <jordigh@gmail.com> parents: 12917 diff changeset	261 %! 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	262 %! 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	263
ab529124f516 codesprint: Wrote 3 tests for sqrtm.cc Jordi Gutiérrez Hermoso <jordigh@gmail.com> parents: 12584 diff changeset	264 %!test
14501 60e5cf354d80 Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	265 %! 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	266 %! 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	267 %! [y, err] = sqrtm (x);
60e5cf354d80 Update %!tests in DLD-FUNCTIONS/ directory with Octave coding conventions. Rik <octave@nomad.inbox5.com> parents: 14138 diff changeset	268 %! assert (y, z);
17336 b81b9d079515 Use '##' for comments which stand alone on a line. Rik <rik@octave.org> parents: 15195 diff changeset	269 %! 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	270 */

Mercurial > octave

annotate libinterp/corefcn/sqrtm.cc @ 22317:6ca3acf5fad8