octave: libinterp/corefcn/sqrtm.cc annotate

annotate libinterp/corefcn/sqrtm.cc @ 21266:e69eaee28737

make better use of templates for Schur decomposition * liboctave/numeric/schur.h, liboctave/numeric/schur.cc: New files generated from SCHUR.h, SCHUR.cc, CmplxSCHUR.h, CmplxSCHUR.cc, dbleSCHUR.h, dbleSCHUR.cc, fCmplxSCHUR.h, fCmplxSCHUR.cc, floatSCHUR.h, and floatSCHUR.cc and making them templates. * liboctave/numeric/module.mk: Update. * libinterp/corefcn/schur.cc, sqrtm.cc, CMatrix.cc, dMatrix.cc, fCMatrix.cc, fMatrix.cc, mx-defs.h, mx-ext.h: Use new template classes and header file.

author	John W. Eaton <jwe@octave.org>
date	Mon, 15 Feb 2016 20:06:12 -0500
parents	fcac5dbbf9ed
children	40de9f8f23a6

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
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	24 #ifdef HAVE_CONFIG_H
21200 fcac5dbbf9ed maint: Indent #ifdef blocks in libinterp. Rik <rik@octave.org> parents: 21139 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
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
21266 e69eaee28737 make better use of templates for Schur decomposition John W. Eaton <jwe@octave.org> parents: 21200 diff changeset	30 #include "schur.h"
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	31 #include "lo-ieee.h"
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	32 #include "lo-mappers.h"
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	33 #include "oct-norm.h"
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	34
15039 e753177cde93 maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory Rik <rik@octave.org> parents: 14501 diff changeset	35 #include "defun.h"
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	36 #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	37 #include "errwarn.h"
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	38 #include "utils.h"
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	39 #include "xnorm.h"
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	40
21139 538b57866b90 consistently use "typename" intead of "class" in template declarations John W. Eaton <jwe@octave.org> parents: 21110 diff changeset	41 template <typename Matrix>
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	42 static void
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	43 sqrtm_utri_inplace (Matrix& T)
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 typedef typename Matrix::element_type 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 const element_type zero = element_type ();
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	48
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	49 bool singular = false;
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	50
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	51 // 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	52 //
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	53 // n = rows (T);
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	54 // for j = 1:n
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	55 // T(j,j) = sqrt (T(j,j));
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	56 // for i = j-1:-1:1
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	57 // 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	58 // k = 1:i-1;
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	59 // 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	60 // endfor
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 //
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	63 // 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	64 // given in Higham's paper.
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	65
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	66 const octave_idx_type n = T.rows ();
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	67 element_type *Tp = T.fortran_vec ();
5275 23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe] jwe parents: 4331 diff changeset	68 for (octave_idx_type j = 0; j < n; j++)
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	69 {
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	70 element_type colj = Tp + nj;
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	71 if (colj[j] != zero)
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	72 colj[j] = sqrt (colj[j]);
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	73 else
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	74 singular = true;
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	75
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	76 for (octave_idx_type i = j-1; i >= 0; i--)
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	77 {
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	78 const element_type coli = Tp + ni;
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	79 const element_type colji = colj[i] /= (coli[i] + colj[j]);
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	80 for (octave_idx_type k = 0; k < i; k++)
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	81 colj[k] -= coli[k] * colji;
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9064 diff changeset	82 }
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	83 }
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	84
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	85 if (singular)
10629 b7c4954e1c00 add specific ID for sqrtm singularity warning Jaroslav Hajek <highegg@gmail.com> parents: 10619 diff changeset	86 warning_with_id ("Octave:sqrtm:SingularMatrix",
b7c4954e1c00 add specific ID for sqrtm singularity warning Jaroslav Hajek <highegg@gmail.com> parents: 10619 diff changeset	87 "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	88 }
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	89
21139 538b57866b90 consistently use "typename" intead of "class" in template declarations John W. Eaton <jwe@octave.org> parents: 21110 diff changeset	90 template <typename Matrix, typename ComplexMatrix, typename ComplexSCHUR>
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	91 static octave_value
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	92 do_sqrtm (const octave_value& arg)
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	93 {
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	94
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	95 octave_value retval;
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	96
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	97 MatrixType mt = arg.matrix_type ();
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	98
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	99 bool iscomplex = arg.is_complex_type ();
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	100
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	101 typedef typename Matrix::element_type real_type;
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	102
18100 6a71e5030df5 Follow coding convention of defining and initializing only 1 variable per line in liboctinterp. Rik <rik@octave.org> parents: 17787 diff changeset	103 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	104 real_type one = 1;
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	105 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	106
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	107 if (! iscomplex)
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	108 {
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	109 Matrix x = octave_value_extract<Matrix> (arg);
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	110
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	111 if (mt.is_unknown ()) // if type is not known, compute it now.
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	112 arg.matrix_type (mt = MatrixType (x));
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	113
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	114 switch (mt.type ())
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9064 diff changeset	115 {
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	116 case MatrixType::Upper:
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	117 case MatrixType::Diagonal:
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	118 if (! x.diag ().any_element_is_negative ())
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	119 {
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	120 // Do it in real arithmetic.
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	121 sqrtm_utri_inplace (x);
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	122 retval = x;
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	123 retval.matrix_type (mt);
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	124 }
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	125 else
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	126 iscomplex = true;
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	127 break;
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	128
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	129 case MatrixType::Lower:
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	130 if (! x.diag ().any_element_is_negative ())
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	131 {
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	132 x = x.transpose ();
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	133 sqrtm_utri_inplace (x);
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	134 retval = x.transpose ();
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	135 retval.matrix_type (mt);
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	136 }
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	137 else
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	138 iscomplex = true;
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	139 break;
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	140
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	141 default:
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	142 iscomplex = true;
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	143 break;
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	144 }
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 if (iscomplex)
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	147 cutoff = 10 * x.rows () * eps * xnorm (x, one);
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	148 }
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 if (iscomplex)
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	151 {
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	152 ComplexMatrix x = octave_value_extract<ComplexMatrix> (arg);
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	153
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	154 if (mt.is_unknown ()) // if type is not known, compute it now.
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	155 arg.matrix_type (mt = MatrixType (x));
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 switch (mt.type ())
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	158 {
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	159 case MatrixType::Upper:
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	160 case MatrixType::Diagonal:
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	161 sqrtm_utri_inplace (x);
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	162 retval = x;
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	163 retval.matrix_type (mt);
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	164 break;
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	165
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	166 case MatrixType::Lower:
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	167 x = x.transpose ();
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	168 sqrtm_utri_inplace (x);
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	169 retval = x.transpose ();
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	170 retval.matrix_type (mt);
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	171 break;
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	172
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	173 default:
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	174 {
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	175 ComplexMatrix u;
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	176
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	177 do
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	178 {
21266 e69eaee28737 make better use of templates for Schur decomposition John W. Eaton <jwe@octave.org> parents: 21200 diff changeset	179 ComplexSCHUR schur_fact (x, "", true);
e69eaee28737 make better use of templates for Schur decomposition John W. Eaton <jwe@octave.org> parents: 21200 diff changeset	180 x = schur_fact.schur_matrix ();
e69eaee28737 make better use of templates for Schur decomposition John W. Eaton <jwe@octave.org> parents: 21200 diff changeset	181 u = schur_fact.unitary_matrix ();
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	182 }
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	183 while (0); // schur no longer needed.
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 sqrtm_utri_inplace (x);
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	186
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	187 x = u * x; // original x no longer needed.
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	188 ComplexMatrix res = xgemm (x, u, blas_no_trans, blas_conj_trans);
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	189
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	190 if (cutoff > 0 && xnorm (imag (res), one) <= cutoff)
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	191 retval = real (res);
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	192 else
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	193 retval = res;
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	194 }
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	195 break;
10154 40dfc0c99116 DLD-FUNCTIONS/.cc: untabify John W. Eaton <jwe@octave.org>* parents: 9064 diff changeset	196 }
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	197 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	198
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	199 return retval;
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	200 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	201
15039 e753177cde93 maint: Move non-dynamically linked functions from DLD-FUNCTIONS/ to corefcn/ directory Rik <rik@octave.org> parents: 14501 diff changeset	202 DEFUN (sqrtm, args, nargout,
17787 175b392e91fe Use GNU style coding conventions for code in libinterp/ Rik <rik@octave.org> parents: 17744 diff changeset	203 "-- texinfo --\n\
20853 1142cf6abc0d 2015 Code Sprint: remove class of function from docstring for all C++ files. Rik <rik@octave.org> parents: 20819 diff changeset	204 @deftypefn {} {@var{s} =} sqrtm (@var{A})\n\
1142cf6abc0d 2015 Code Sprint: remove class of function from docstring for all C++ files. Rik <rik@octave.org> parents: 20819 diff changeset	205 @deftypefnx {} {[@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	206 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	207 \n\
19040 0850b5212619 doc: Add @nospell macro around proper names in documentation. Rik <rik@octave.org> parents: 18100 diff changeset	208 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	209 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	210 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	211 @seealso{expm, logm}\n\
5642 2618a0750ae6 [project @ 2006-03-06 21:26:48 by jwe] jwe parents: 5307 diff changeset	212 @end deftypefn")
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	213 {
20819 f428cbe7576f eliminate unnecessary uses of nargin John W. Eaton <jwe@octave.org> parents: 20801 diff changeset	214 if (args.length () != 1)
20801 a542a9bf177e eliminate return statements after calls to print_usage John W. Eaton <jwe@octave.org> parents: 20555 diff changeset	215 print_usage ();
4331 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 octave_value arg = args(0);
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	218
5275 23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe] jwe parents: 4331 diff changeset	219 octave_idx_type n = arg.rows ();
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe] jwe parents: 4331 diff changeset	220 octave_idx_type nc = arg.columns ();
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	221
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	222 if (n != nc \|\| arg.ndims () > 2)
21110 3d0d84305600 Use err_square_matrix_required more widely. Rik <rik@octave.org> parents: 21100 diff changeset	223 err_square_matrix_required ("sqrtm", "A");
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	224
20892 c07bee629973 2015 Code Sprint: use ovl (). Rik <rik@octave.org> parents: 20853 diff changeset	225 octave_value_list retval (nargout > 1 ? 3 : 1);
c07bee629973 2015 Code Sprint: use ovl (). Rik <rik@octave.org> parents: 20853 diff changeset	226
14022 de90542b7afc Return retval(1) before retval(0) to avoid re-sizing call on retval. Rik <octave@nomad.inbox5.com> parents: 12926 diff changeset	227 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	228 {
20892 c07bee629973 2015 Code Sprint: use ovl (). Rik <rik@octave.org> parents: 20853 diff changeset	229 // 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	230 // 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	231 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	232 }
de90542b7afc Return retval(1) before retval(0) to avoid re-sizing call on retval. Rik <octave@nomad.inbox5.com> parents: 12926 diff changeset	233
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	234 if (arg.is_diag_matrix ())
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	235 // sqrtm of a diagonal matrix is just sqrt.
b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	236 retval(0) = arg.sqrt ();
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	237 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	238 retval(0) = do_sqrtm<FloatMatrix, FloatComplexMatrix,
e69eaee28737 make better use of templates for Schur decomposition John W. Eaton <jwe@octave.org> parents: 21200 diff changeset	239 schur<FloatComplexMatrix> > (arg);
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	240 else if (arg.is_numeric_type ())
21266 e69eaee28737 make better use of templates for Schur decomposition John W. Eaton <jwe@octave.org> parents: 21200 diff changeset	241 retval(0) = do_sqrtm<Matrix, ComplexMatrix, schur<ComplexMatrix> > (arg);
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	242
20555 f90c8372b7ba eliminate many more simple uses of error_state John W. Eaton <jwe@octave.org> parents: 19697 diff changeset	243 if (nargout > 1)
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	244 {
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	245 // This corresponds to generic code
12902 b9c66f891281 maint: style fixes for sqrtm John W. Eaton <jwe@octave.org> parents: 12829 diff changeset	246 //
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	247 // 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	248
10608 f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	249 octave_value s = retval(0);
f9860b622680 improve sqrtm Jaroslav Hajek <highegg@gmail.com> parents: 10155 diff changeset	250 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	251 }
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	252
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	253 return retval;
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	254 }
12829 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 /*
ab529124f516 codesprint: Wrote 3 tests for sqrtm.cc Jordi Gutiérrez Hermoso <jordigh@gmail.com> parents: 12584 diff changeset	257 %!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	258
ab529124f516 codesprint: Wrote 3 tests for sqrtm.cc Jordi Gutiérrez Hermoso <jordigh@gmail.com> parents: 12584 diff changeset	259 ## 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	260 %!test
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 @ 21266:e69eaee28737