Mercurial > octave-nkf

annotate src/DLD-FUNCTIONS/sqrtm.cc @ 14138:72c96de7a403 stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: update copyright notices for 2012
author John W. Eaton <jwe@octave.org>
date Mon, 02 Jan 2012 14:25:41 -0500
parents de90542b7afc
children 60e5cf354d80
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
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
14138
72c96de7a403 maint: update copyright notices for 2012
John W. Eaton <jwe@octave.org>
parents: 14022
diff changeset
3 Copyright (C) 2001-2012 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
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe]
jwe
parents:
diff changeset
36 #include "defun-dld.h"
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 + n*j;
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 + n*i;
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
f9860b622680 improve sqrtm
Jaroslav Hajek <highegg@gmail.com>
parents: 10155
diff changeset
104 real_type cutoff = 0, one = 1;
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 {
b9c66f891281 maint: style fixes for sqrtm
John W. Eaton <jwe@octave.org>
parents: 12829
diff changeset
179 ComplexSCHUR schur (x, std::string (), true);
b9c66f891281 maint: style fixes for sqrtm
John W. Eaton <jwe@octave.org>
parents: 12829
diff changeset
180 x = schur.schur_matrix ();
b9c66f891281 maint: style fixes for sqrtm
John W. Eaton <jwe@octave.org>
parents: 12829
diff changeset
181 u = schur.unitary_matrix ();
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
4331
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe]
jwe
parents:
diff changeset
202 DEFUN_DLD (sqrtm, args, nargout,
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe]
jwe
parents:
diff changeset
203 "-*- texinfo -*-\n\
12584
7ef7e20057fa Improve documentation strings in Linear Algebra chapter.
Rik <octave@nomad.inbox5.com>
parents: 11586
diff changeset
204 @deftypefn {Loadable Function} {@var{s} =} sqrtm (@var{A})\n\
7ef7e20057fa Improve documentation strings in Linear Algebra chapter.
Rik <octave@nomad.inbox5.com>
parents: 11586
diff changeset
205 @deftypefnx {Loadable 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
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\
10791
3140cb7a05a1 Add spellchecker scripts for Octave and run spellcheck of documentation
Rik <octave@nomad.inbox5.com>
parents: 10629
diff changeset
208 Ref: N.J. Higham. @cite{A New sqrtm for @sc{matlab}}. Numerical\n\
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 {
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe]
jwe
parents:
diff changeset
214 octave_value_list retval;
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe]
jwe
parents:
diff changeset
215
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe]
jwe
parents:
diff changeset
216 int nargin = args.length ();
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe]
jwe
parents:
diff changeset
217
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe]
jwe
parents:
diff changeset
218 if (nargin != 1)
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe]
jwe
parents:
diff changeset
219 {
5823
080c08b192d8 [project @ 2006-05-19 05:32:17 by jwe]
jwe
parents: 5775
diff changeset
220 print_usage ();
4331
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe]
jwe
parents:
diff changeset
221 return retval;
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe]
jwe
parents:
diff changeset
222 }
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 octave_value arg = args(0);
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe]
jwe
parents:
diff changeset
225
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 4331
diff changeset
226 octave_idx_type n = arg.rows ();
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 4331
diff changeset
227 octave_idx_type nc = arg.columns ();
4331
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe]
jwe
parents:
diff changeset
228
10608
f9860b622680 improve sqrtm
Jaroslav Hajek <highegg@gmail.com>
parents: 10155
diff changeset
229 if (n != nc || arg.ndims () > 2)
4331
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe]
jwe
parents:
diff changeset
230 {
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe]
jwe
parents:
diff changeset
231 gripe_square_matrix_required ("sqrtm");
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe]
jwe
parents:
diff changeset
232 return retval;
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe]
jwe
parents:
diff changeset
233 }
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe]
jwe
parents:
diff changeset
234
14022
de90542b7afc Return retval(1) before retval(0) to avoid re-sizing call on retval.
Rik <octave@nomad.inbox5.com>
parents: 12926
diff changeset
235 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
236 {
de90542b7afc Return retval(1) before retval(0) to avoid re-sizing call on retval.
Rik <octave@nomad.inbox5.com>
parents: 12926
diff changeset
237 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
238 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
239 }
de90542b7afc Return retval(1) before retval(0) to avoid re-sizing call on retval.
Rik <octave@nomad.inbox5.com>
parents: 12926
diff changeset
240
10608
f9860b622680 improve sqrtm
Jaroslav Hajek <highegg@gmail.com>
parents: 10155
diff changeset
241 if (arg.is_diag_matrix ())
12902
b9c66f891281 maint: style fixes for sqrtm
John W. Eaton <jwe@octave.org>
parents: 12829
diff changeset
242 // sqrtm of a diagonal matrix is just sqrt.
b9c66f891281 maint: style fixes for sqrtm
John W. Eaton <jwe@octave.org>
parents: 12829
diff changeset
243 retval(0) = arg.sqrt ();
10608
f9860b622680 improve sqrtm
Jaroslav Hajek <highegg@gmail.com>
parents: 10155
diff changeset
244 else if (arg.is_single_type ())
12902
b9c66f891281 maint: style fixes for sqrtm
John W. Eaton <jwe@octave.org>
parents: 12829
diff changeset
245 retval(0) = do_sqrtm<FloatMatrix, FloatComplexMatrix, FloatComplexSCHUR> (arg);
10608
f9860b622680 improve sqrtm
Jaroslav Hajek <highegg@gmail.com>
parents: 10155
diff changeset
246 else if (arg.is_numeric_type ())
12902
b9c66f891281 maint: style fixes for sqrtm
John W. Eaton <jwe@octave.org>
parents: 12829
diff changeset
247 retval(0) = do_sqrtm<Matrix, ComplexMatrix, ComplexSCHUR> (arg);
4331
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe]
jwe
parents:
diff changeset
248
10608
f9860b622680 improve sqrtm
Jaroslav Hajek <highegg@gmail.com>
parents: 10155
diff changeset
249 if (nargout > 1 && ! error_state)
7789
82be108cc558 First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
250 {
10608
f9860b622680 improve sqrtm
Jaroslav Hajek <highegg@gmail.com>
parents: 10155
diff changeset
251 // This corresponds to generic code
12902
b9c66f891281 maint: style fixes for sqrtm
John W. Eaton <jwe@octave.org>
parents: 12829
diff changeset
252 //
10608
f9860b622680 improve sqrtm
Jaroslav Hajek <highegg@gmail.com>
parents: 10155
diff changeset
253 // 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
254
10608
f9860b622680 improve sqrtm
Jaroslav Hajek <highegg@gmail.com>
parents: 10155
diff changeset
255 octave_value s = retval(0);
f9860b622680 improve sqrtm
Jaroslav Hajek <highegg@gmail.com>
parents: 10155
diff changeset
256 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
257 }
4331
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe]
jwe
parents:
diff changeset
258
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe]
jwe
parents:
diff changeset
259 return retval;
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe]
jwe
parents:
diff changeset
260 }
12829
ab529124f516 codesprint: Wrote 3 tests for sqrtm.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents: 12584
diff changeset
261
ab529124f516 codesprint: Wrote 3 tests for sqrtm.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents: 12584
diff changeset
262 /*
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 %!assert (sqrtm (2*ones (2)), ones (2), 3*eps)
ab529124f516 codesprint: Wrote 3 tests for sqrtm.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents: 12584
diff changeset
265
ab529124f516 codesprint: Wrote 3 tests for sqrtm.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents: 12584
diff changeset
266 ## The following two tests are from the reference in the docstring above.
ab529124f516 codesprint: Wrote 3 tests for sqrtm.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents: 12584
diff changeset
267
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];
838ada8779a5 Fix brain-o in sqrtm test
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents: 12917
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
ab529124f516 codesprint: Wrote 3 tests for sqrtm.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents: 12584
diff changeset
273 %! x = eye (4); x(2,2) = x(3,3) = 2^-26; x(1,4) = 1;
ab529124f516 codesprint: Wrote 3 tests for sqrtm.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents: 12584
diff changeset
274 %! z = eye (4); z(2,2) = z(3,3) = 2^-13; z(1,4) = 0.5;
ab529124f516 codesprint: Wrote 3 tests for sqrtm.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents: 12584
diff changeset
275 %! [y, err] = sqrtm(x);
ab529124f516 codesprint: Wrote 3 tests for sqrtm.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents: 12584
diff changeset
276 %! assert (y, z)
12917
d6209c354667 Change test on sqrtm.cc to use any when checking for NaN
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents: 12902
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
ab529124f516 codesprint: Wrote 3 tests for sqrtm.cc
Jordi Gutiérrez Hermoso <jordigh@gmail.com>
parents: 12584
diff changeset
279 */