Mercurial > octave
annotate src/DLD-FUNCTIONS/sqrtm.cc @ 10619:9f0a264d2f60
mark sqrtm(triangular) as triangular
author | Jaroslav Hajek <highegg@gmail.com> |
---|---|
date | Fri, 07 May 2010 12:40:06 +0200 |
parents | f9860b622680 |
children | b7c4954e1c00 |
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) |
89 warning ("sqrtm: matrix is singular, may not have a square root"); | |
4331 | 90 } |
91 | |
10608 | 92 template <class Matrix, class ComplexMatrix, class ComplexSCHUR> |
93 static octave_value | |
94 do_sqrtm (const octave_value& arg) | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
95 { |
10608 | 96 |
97 octave_value retval; | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
98 |
10608 | 99 MatrixType mt = arg.matrix_type (); |
100 | |
101 bool iscomplex = arg.is_complex_type (); | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
102 |
10608 | 103 typedef typename Matrix::element_type real_type; |
104 | |
105 real_type cutoff = 0, one = 1; | |
106 real_type eps = std::numeric_limits<real_type>::epsilon (); | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
107 |
10608 | 108 if (! iscomplex) |
109 { | |
110 Matrix x = octave_value_extract<Matrix> (arg); | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
111 |
10608 | 112 if (mt.is_unknown ()) // if type is not known, compute it now. |
113 arg.matrix_type (mt = MatrixType (x)); | |
114 | |
115 switch (mt.type ()) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9064
diff
changeset
|
116 { |
10608 | 117 case MatrixType::Upper: |
118 case MatrixType::Diagonal: | |
119 { | |
120 if (! x.diag ().any_element_is_negative ()) | |
121 { | |
122 // Do it in real arithmetic. | |
123 sqrtm_utri_inplace (x); | |
124 retval = x; | |
10619
9f0a264d2f60
mark sqrtm(triangular) as triangular
Jaroslav Hajek <highegg@gmail.com>
parents:
10608
diff
changeset
|
125 retval.matrix_type (mt); |
10608 | 126 } |
127 else | |
128 iscomplex = true; | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
129 |
10608 | 130 break; |
131 } | |
132 case MatrixType::Lower: | |
133 { | |
134 if (! x.diag ().any_element_is_negative ()) | |
135 { | |
136 x = x.transpose (); | |
137 sqrtm_utri_inplace (x); | |
138 retval = x.transpose (); | |
10619
9f0a264d2f60
mark sqrtm(triangular) as triangular
Jaroslav Hajek <highegg@gmail.com>
parents:
10608
diff
changeset
|
139 retval.matrix_type (mt); |
10608 | 140 } |
141 else | |
142 iscomplex = true; | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
143 |
10608 | 144 break; |
145 } | |
146 default: | |
147 { | |
148 iscomplex = true; | |
149 break; | |
150 } | |
151 } | |
152 | |
153 if (iscomplex) | |
154 cutoff = 10 * x.rows () * eps * xnorm (x, one); | |
155 } | |
156 | |
157 if (iscomplex) | |
158 { | |
159 ComplexMatrix x = octave_value_extract<ComplexMatrix> (arg); | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
160 |
10608 | 161 if (mt.is_unknown ()) // if type is not known, compute it now. |
162 arg.matrix_type (mt = MatrixType (x)); | |
163 | |
164 switch (mt.type ()) | |
165 { | |
166 case MatrixType::Upper: | |
167 case MatrixType::Diagonal: | |
168 { | |
169 sqrtm_utri_inplace (x); | |
170 retval = x; | |
10619
9f0a264d2f60
mark sqrtm(triangular) as triangular
Jaroslav Hajek <highegg@gmail.com>
parents:
10608
diff
changeset
|
171 retval.matrix_type (mt); |
10608 | 172 |
173 break; | |
174 } | |
175 case MatrixType::Lower: | |
176 { | |
177 x = x.transpose (); | |
178 sqrtm_utri_inplace (x); | |
179 retval = x.transpose (); | |
10619
9f0a264d2f60
mark sqrtm(triangular) as triangular
Jaroslav Hajek <highegg@gmail.com>
parents:
10608
diff
changeset
|
180 retval.matrix_type (mt); |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
181 |
10608 | 182 break; |
183 } | |
184 default: | |
185 { | |
186 ComplexMatrix u; | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
187 |
10608 | 188 do |
189 { | |
190 ComplexSCHUR schur (x, std::string (), true); | |
191 x = schur.schur_matrix (); | |
192 u = schur.unitary_matrix (); | |
193 } | |
194 while (0); // schur no longer needed. | |
195 | |
196 sqrtm_utri_inplace (x); | |
197 | |
198 x = u * x; // original x no longer needed. | |
199 ComplexMatrix res = xgemm (x, u, blas_no_trans, blas_conj_trans); | |
200 | |
201 if (cutoff > 0 && xnorm (imag (res), one) <= cutoff) | |
202 retval = real (res); | |
203 else | |
204 retval = res; | |
205 | |
206 break; | |
207 } | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9064
diff
changeset
|
208 } |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
209 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
210 |
10608 | 211 return retval; |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
212 } |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
213 |
4331 | 214 DEFUN_DLD (sqrtm, args, nargout, |
215 "-*- texinfo -*-\n\ | |
216 @deftypefn {Loadable Function} {[@var{result}, @var{error_estimate}] =} sqrtm (@var{a})\n\ | |
217 Compute the matrix square root of the square matrix @var{a}.\n\ | |
218 \n\ | |
9064
7c02ec148a3c
Check grammar on all .cc files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
219 Ref: Nicholas J. Higham. A new sqrtm for @sc{matlab}. Numerical Analysis\n\ |
4331 | 220 Report No. 336, Manchester Centre for Computational Mathematics,\n\ |
221 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
|
222 @seealso{expm, logm}\n\ |
5642 | 223 @end deftypefn") |
4331 | 224 { |
225 octave_value_list retval; | |
226 | |
227 int nargin = args.length (); | |
228 | |
229 if (nargin != 1) | |
230 { | |
5823 | 231 print_usage (); |
4331 | 232 return retval; |
233 } | |
234 | |
235 octave_value arg = args(0); | |
236 | |
5275 | 237 octave_idx_type n = arg.rows (); |
238 octave_idx_type nc = arg.columns (); | |
4331 | 239 |
10608 | 240 if (n != nc || arg.ndims () > 2) |
4331 | 241 { |
242 gripe_square_matrix_required ("sqrtm"); | |
243 return retval; | |
244 } | |
245 | |
10608 | 246 if (arg.is_diag_matrix ()) |
247 { | |
248 // sqrtm of a diagonal matrix is just sqrt. | |
249 retval(0) = arg.sqrt (); | |
250 } | |
251 else if (arg.is_single_type ()) | |
4331 | 252 { |
10608 | 253 retval(0) = do_sqrtm<FloatMatrix, FloatComplexMatrix, FloatComplexSCHUR> (arg); |
254 } | |
255 else if (arg.is_numeric_type ()) | |
256 { | |
257 retval(0) = do_sqrtm<Matrix, ComplexMatrix, ComplexSCHUR> (arg); | |
258 } | |
4331 | 259 |
10608 | 260 if (nargout > 1 && ! error_state) |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
261 { |
10608 | 262 // This corresponds to generic code |
263 // 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
|
264 |
10608 | 265 octave_value s = retval(0); |
266 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
|
267 } |
4331 | 268 |
269 return retval; | |
270 } |