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