Mercurial > octave-nkf
annotate src/DLD-FUNCTIONS/sqrtm.cc @ 10846:a4f482e66b65
Grammarcheck more of the documentation.
Use @noindent macro appropriately.
Limit line length to 80 characters.
author | Rik <octave@nomad.inbox5.com> |
---|---|
date | Sun, 01 Aug 2010 20:22:17 -0700 |
parents | 3140cb7a05a1 |
children | fd0a3ac60b0e |
rev | line source |
---|---|
4331 | 1 /* |
2 | |
8920 | 3 Copyright (C) 2001, 2003, 2005, 2006, 2007, 2008 Ross Lippert and Paul Kienzle |
10608 | 4 Copyright (C) 2010 VZLU Prague |
4331 | 5 |
7016 | 6 This file is part of Octave. |
4331 | 7 |
7016 | 8 Octave is free software; you can redistribute it and/or modify it |
9 under the terms of the GNU General Public License as published by the | |
10 Free Software Foundation; either version 3 of the License, or (at your | |
11 option) any later version. | |
12 | |
13 Octave is distributed in the hope that it will be useful, but WITHOUT | |
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
16 for more details. | |
4331 | 17 |
18 You should have received a copy of the GNU General Public License | |
7016 | 19 along with Octave; see the file COPYING. If not, see |
20 <http://www.gnu.org/licenses/>. | |
4331 | 21 |
22 */ | |
23 | |
24 #ifdef HAVE_CONFIG_H | |
25 #include <config.h> | |
26 #endif | |
27 | |
28 #include <float.h> | |
29 | |
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 | 32 #include "lo-ieee.h" |
33 #include "lo-mappers.h" | |
10608 | 34 #include "oct-norm.h" |
4331 | 35 |
36 #include "defun-dld.h" | |
37 #include "error.h" | |
38 #include "gripes.h" | |
39 #include "utils.h" | |
10608 | 40 #include "xnorm.h" |
4331 | 41 |
10608 | 42 template <class Matrix> |
43 static void | |
44 sqrtm_utri_inplace (Matrix& T) | |
4331 | 45 { |
10608 | 46 typedef typename Matrix::element_type element_type; |
4331 | 47 |
10608 | 48 const element_type zero = element_type (); |
4331 | 49 |
10608 | 50 bool singular = false; |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
51 |
10608 | 52 /* |
53 * the following code is equivalent to this triple loop: | |
54 * | |
55 * n = rows (T); | |
56 * for j = 1:n | |
57 * T(j,j) = sqrt (T(j,j)); | |
58 * for i = j-1:-1:1 | |
59 * T(i,j) /= (T(i,i) + T(j,j)); | |
60 * k = 1:i-1; | |
61 * T(k,j) -= T(k,i) * T(i,j); | |
62 * endfor | |
63 * endfor | |
64 * | |
65 * this is an in-place, cache-aligned variant of the code | |
66 * given in Higham's paper. | |
67 */ | |
4331 | 68 |
10608 | 69 const octave_idx_type n = T.rows (); |
70 element_type *Tp = T.fortran_vec (); | |
5275 | 71 for (octave_idx_type j = 0; j < n; j++) |
4331 | 72 { |
10608 | 73 element_type *colj = Tp + n*j; |
74 if (colj[j] != zero) | |
75 colj[j] = sqrt (colj[j]); | |
76 else | |
77 singular = true; | |
4331 | 78 |
10608 | 79 for (octave_idx_type i = j-1; i >= 0; i--) |
80 { | |
81 const element_type *coli = Tp + n*i; | |
82 const element_type colji = colj[i] /= (coli[i] + colj[j]); | |
83 for (octave_idx_type k = 0; k < i; k++) | |
84 colj[k] -= coli[k] * colji; | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9064
diff
changeset
|
85 } |
4331 | 86 } |
87 | |
10608 | 88 if (singular) |
10629
b7c4954e1c00
add specific ID for sqrtm singularity warning
Jaroslav Hajek <highegg@gmail.com>
parents:
10619
diff
changeset
|
89 warning_with_id ("Octave:sqrtm:SingularMatrix", |
b7c4954e1c00
add specific ID for sqrtm singularity warning
Jaroslav Hajek <highegg@gmail.com>
parents:
10619
diff
changeset
|
90 "sqrtm: matrix is singular, may not have a square root"); |
4331 | 91 } |
92 | |
10608 | 93 template <class Matrix, class ComplexMatrix, class ComplexSCHUR> |
94 static octave_value | |
95 do_sqrtm (const octave_value& arg) | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
96 { |
10608 | 97 |
98 octave_value retval; | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
99 |
10608 | 100 MatrixType mt = arg.matrix_type (); |
101 | |
102 bool iscomplex = arg.is_complex_type (); | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
103 |
10608 | 104 typedef typename Matrix::element_type real_type; |
105 | |
106 real_type cutoff = 0, one = 1; | |
107 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
|
108 |
10608 | 109 if (! iscomplex) |
110 { | |
111 Matrix x = octave_value_extract<Matrix> (arg); | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
112 |
10608 | 113 if (mt.is_unknown ()) // if type is not known, compute it now. |
114 arg.matrix_type (mt = MatrixType (x)); | |
115 | |
116 switch (mt.type ()) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9064
diff
changeset
|
117 { |
10608 | 118 case MatrixType::Upper: |
119 case MatrixType::Diagonal: | |
120 { | |
121 if (! x.diag ().any_element_is_negative ()) | |
122 { | |
123 // Do it in real arithmetic. | |
124 sqrtm_utri_inplace (x); | |
125 retval = x; | |
10619
9f0a264d2f60
mark sqrtm(triangular) as triangular
Jaroslav Hajek <highegg@gmail.com>
parents:
10608
diff
changeset
|
126 retval.matrix_type (mt); |
10608 | 127 } |
128 else | |
129 iscomplex = true; | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
130 |
10608 | 131 break; |
132 } | |
133 case MatrixType::Lower: | |
134 { | |
135 if (! x.diag ().any_element_is_negative ()) | |
136 { | |
137 x = x.transpose (); | |
138 sqrtm_utri_inplace (x); | |
139 retval = x.transpose (); | |
10619
9f0a264d2f60
mark sqrtm(triangular) as triangular
Jaroslav Hajek <highegg@gmail.com>
parents:
10608
diff
changeset
|
140 retval.matrix_type (mt); |
10608 | 141 } |
142 else | |
143 iscomplex = true; | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
144 |
10608 | 145 break; |
146 } | |
147 default: | |
148 { | |
149 iscomplex = true; | |
150 break; | |
151 } | |
152 } | |
153 | |
154 if (iscomplex) | |
155 cutoff = 10 * x.rows () * eps * xnorm (x, one); | |
156 } | |
157 | |
158 if (iscomplex) | |
159 { | |
160 ComplexMatrix x = octave_value_extract<ComplexMatrix> (arg); | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
161 |
10608 | 162 if (mt.is_unknown ()) // if type is not known, compute it now. |
163 arg.matrix_type (mt = MatrixType (x)); | |
164 | |
165 switch (mt.type ()) | |
166 { | |
167 case MatrixType::Upper: | |
168 case MatrixType::Diagonal: | |
169 { | |
170 sqrtm_utri_inplace (x); | |
171 retval = x; | |
10619
9f0a264d2f60
mark sqrtm(triangular) as triangular
Jaroslav Hajek <highegg@gmail.com>
parents:
10608
diff
changeset
|
172 retval.matrix_type (mt); |
10608 | 173 |
174 break; | |
175 } | |
176 case MatrixType::Lower: | |
177 { | |
178 x = x.transpose (); | |
179 sqrtm_utri_inplace (x); | |
180 retval = x.transpose (); | |
10619
9f0a264d2f60
mark sqrtm(triangular) as triangular
Jaroslav Hajek <highegg@gmail.com>
parents:
10608
diff
changeset
|
181 retval.matrix_type (mt); |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
182 |
10608 | 183 break; |
184 } | |
185 default: | |
186 { | |
187 ComplexMatrix u; | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
188 |
10608 | 189 do |
190 { | |
191 ComplexSCHUR schur (x, std::string (), true); | |
192 x = schur.schur_matrix (); | |
193 u = schur.unitary_matrix (); | |
194 } | |
195 while (0); // schur no longer needed. | |
196 | |
197 sqrtm_utri_inplace (x); | |
198 | |
199 x = u * x; // original x no longer needed. | |
200 ComplexMatrix res = xgemm (x, u, blas_no_trans, blas_conj_trans); | |
201 | |
202 if (cutoff > 0 && xnorm (imag (res), one) <= cutoff) | |
203 retval = real (res); | |
204 else | |
205 retval = res; | |
206 | |
207 break; | |
208 } | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9064
diff
changeset
|
209 } |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
210 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
211 |
10608 | 212 return retval; |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
213 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
214 |
4331 | 215 DEFUN_DLD (sqrtm, args, nargout, |
216 "-*- texinfo -*-\n\ | |
217 @deftypefn {Loadable Function} {[@var{result}, @var{error_estimate}] =} sqrtm (@var{a})\n\ | |
218 Compute the matrix square root of the square matrix @var{a}.\n\ | |
219 \n\ | |
10791
3140cb7a05a1
Add spellchecker scripts for Octave and run spellcheck of documentation
Rik <octave@nomad.inbox5.com>
parents:
10629
diff
changeset
|
220 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
|
221 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
|
222 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
|
223 @seealso{expm, logm}\n\ |
5642 | 224 @end deftypefn") |
4331 | 225 { |
226 octave_value_list retval; | |
227 | |
228 int nargin = args.length (); | |
229 | |
230 if (nargin != 1) | |
231 { | |
5823 | 232 print_usage (); |
4331 | 233 return retval; |
234 } | |
235 | |
236 octave_value arg = args(0); | |
237 | |
5275 | 238 octave_idx_type n = arg.rows (); |
239 octave_idx_type nc = arg.columns (); | |
4331 | 240 |
10608 | 241 if (n != nc || arg.ndims () > 2) |
4331 | 242 { |
243 gripe_square_matrix_required ("sqrtm"); | |
244 return retval; | |
245 } | |
246 | |
10608 | 247 if (arg.is_diag_matrix ()) |
248 { | |
249 // sqrtm of a diagonal matrix is just sqrt. | |
250 retval(0) = arg.sqrt (); | |
251 } | |
252 else if (arg.is_single_type ()) | |
4331 | 253 { |
10608 | 254 retval(0) = do_sqrtm<FloatMatrix, FloatComplexMatrix, FloatComplexSCHUR> (arg); |
255 } | |
256 else if (arg.is_numeric_type ()) | |
257 { | |
258 retval(0) = do_sqrtm<Matrix, ComplexMatrix, ComplexSCHUR> (arg); | |
259 } | |
4331 | 260 |
10608 | 261 if (nargout > 1 && ! error_state) |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
262 { |
10608 | 263 // This corresponds to generic code |
264 // 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
|
265 |
10608 | 266 octave_value s = retval(0); |
267 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
|
268 } |
4331 | 269 |
270 return retval; | |
271 } |