Mercurial > octave-nkf
annotate src/DLD-FUNCTIONS/sqrtm.cc @ 11553:01f703952eff
Improve docstrings for functions in DLD-FUNCTIONS directory.
Use same variable names in error() strings and in documentation.
| author | Rik <octave@nomad.inbox5.com> |
|---|---|
| date | Sun, 16 Jan 2011 22:13:23 -0800 |
| parents | fd0a3ac60b0e |
| children | 12df7854fa7c |
| rev | line source |
|---|---|
| 4331 | 1 /* |
| 2 | |
| 11523 | 3 Copyright (C) 2001-2011 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\ | |
|
11553
01f703952eff
Improve docstrings for functions in DLD-FUNCTIONS directory.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
217 @deftypefn {Loadable Function} {[@var{result}, @var{error_estimate}] =} sqrtm (@var{A})\n\ |
|
01f703952eff
Improve docstrings for functions in DLD-FUNCTIONS directory.
Rik <octave@nomad.inbox5.com>
parents:
11523
diff
changeset
|
218 Compute the matrix square root of the square matrix @var{A}.\n\ |
| 4331 | 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 } |