octave: src/DLD-FUNCTIONS/sqrtm.cc annotate

annotate src/DLD-FUNCTIONS/sqrtm.cc @ 7795:df9519e9990c

Handle single precision eps values

author	David Bateman <dbateman@free.fr>
date	Mon, 12 May 2008 22:57:11 +0200
parents	82be108cc558
children	6f2d95255911

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
7017 a1dbe9d80eee [project @ 2007-10-12 21:27:11 by jwe] jwe parents: 7016 diff changeset	3 Copyright (C) 2001, 2003, 2005, 2006, 2007 Ross Lippert and Paul Kienzle
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	4
7016 93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5823 diff changeset	5 This file is part of Octave.
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	6
7016 93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5823 diff changeset	7 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	8 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	9 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	10 option) any later version.
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5823 diff changeset	11
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5823 diff changeset	12 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	13 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	14 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	15 for more details.
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	16
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	17 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	18 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	19 <http://www.gnu.org/licenses/>.
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	20
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 #ifdef HAVE_CONFIG_H
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	24 #include <config.h>
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	25 #endif
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	26
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	27 #include <float.h>
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	28
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	29 #include "CmplxSCHUR.h"
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	30 #include "fCmplxSCHUR.h"
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	31 #include "lo-ieee.h"
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	32 #include "lo-mappers.h"
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	33
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	34 #include "defun-dld.h"
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	35 #include "error.h"
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	36 #include "gripes.h"
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	37 #include "utils.h"
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	38
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	39 template <class T>
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	40 static inline T
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	41 getmin (T x, T y)
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	42 {
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	43 return x < y ? x : y;
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	44 }
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	45
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	46 template <class T>
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	47 static inline T
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	48 getmax (T x, T y)
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	49 {
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	50 return x > y ? x : y;
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	51 }
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	52
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	53 static double
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	54 frobnorm (const ComplexMatrix& A)
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	55 {
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	56 double sum = 0;
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	57
5275 23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe] jwe parents: 4331 diff changeset	58 for (octave_idx_type i = 0; i < A.rows (); i++)
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe] jwe parents: 4331 diff changeset	59 for (octave_idx_type j = 0; j < A.columns (); j++)
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	60 sum += real (A(i,j) * conj (A(i,j)));
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	61
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	62 return sqrt (sum);
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	63 }
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	64
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	65 static double
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	66 frobnorm (const Matrix& A)
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	67 {
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	68 double sum = 0;
5275 23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe] jwe parents: 4331 diff changeset	69 for (octave_idx_type i = 0; i < A.rows (); i++)
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe] jwe parents: 4331 diff changeset	70 for (octave_idx_type j = 0; j < A.columns (); j++)
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	71 sum += A(i,j) * A(i,j);
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	72
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	73 return sqrt (sum);
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	74 }
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	75
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	76 static float
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	77 frobnorm (const FloatComplexMatrix& A)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	78 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	79 float sum = 0;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	80
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	81 for (octave_idx_type i = 0; i < A.rows (); i++)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	82 for (octave_idx_type j = 0; j < A.columns (); j++)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	83 sum += real (A(i,j) * conj (A(i,j)));
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	84
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	85 return sqrt (sum);
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	86 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	87
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	88 static float
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	89 frobnorm (const FloatMatrix& A)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	90 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	91 float sum = 0;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	92 for (octave_idx_type i = 0; i < A.rows (); i++)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	93 for (octave_idx_type j = 0; j < A.columns (); j++)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	94 sum += A(i,j) * A(i,j);
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	95
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	96 return sqrt (sum);
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	97 }
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	98
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	99 static ComplexMatrix
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	100 sqrtm_from_schur (const ComplexMatrix& U, const ComplexMatrix& T)
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	101 {
5275 23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe] jwe parents: 4331 diff changeset	102 const octave_idx_type n = U.rows ();
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	103
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	104 ComplexMatrix R (n, n, 0.0);
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	105
5275 23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe] jwe parents: 4331 diff changeset	106 for (octave_idx_type j = 0; j < n; j++)
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	107 R(j,j) = sqrt (T(j,j));
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	108
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	109 const double fudge = sqrt (DBL_MIN);
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	110
5275 23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe] jwe parents: 4331 diff changeset	111 for (octave_idx_type p = 0; p < n-1; p++)
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	112 {
5275 23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe] jwe parents: 4331 diff changeset	113 for (octave_idx_type i = 0; i < n-(p+1); i++)
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	114 {
5275 23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe] jwe parents: 4331 diff changeset	115 const octave_idx_type j = i + p + 1;
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	116
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	117 Complex s = T(i,j);
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	118
5275 23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe] jwe parents: 4331 diff changeset	119 for (octave_idx_type k = i+1; k < j; k++)
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	120 s -= R(i,k) * R(k,j);
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	121
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	122 // dividing
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	123 // R(i,j) = s/(R(i,i)+R(j,j));
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	124 // screwing around to not / 0
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	125
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	126 const Complex d = R(i,i) + R(j,j) + fudge;
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	127 const Complex conjd = conj (d);
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	128
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	129 R(i,j) = (sconjd)/(dconjd);
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	130 }
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	131 }
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	132
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	133 return U * R * U.hermitian ();
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	134 }
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	135
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	136 static FloatComplexMatrix
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	137 sqrtm_from_schur (const FloatComplexMatrix& U, const FloatComplexMatrix& T)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	138 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	139 const octave_idx_type n = U.rows ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	140
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	141 FloatComplexMatrix R (n, n, 0.0);
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	142
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	143 for (octave_idx_type j = 0; j < n; j++)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	144 R(j,j) = sqrt (T(j,j));
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	145
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	146 const float fudge = sqrt (FLT_MIN);
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	147
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	148 for (octave_idx_type p = 0; p < n-1; p++)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	149 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	150 for (octave_idx_type i = 0; i < n-(p+1); i++)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	151 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	152 const octave_idx_type j = i + p + 1;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	153
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	154 FloatComplex s = T(i,j);
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	155
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	156 for (octave_idx_type k = i+1; k < j; k++)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	157 s -= R(i,k) * R(k,j);
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	158
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	159 // dividing
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	160 // R(i,j) = s/(R(i,i)+R(j,j));
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	161 // screwing around to not / 0
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	162
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	163 const FloatComplex d = R(i,i) + R(j,j) + fudge;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	164 const FloatComplex conjd = conj (d);
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	165
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	166 R(i,j) = (sconjd)/(dconjd);
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	167 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	168 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	169
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	170 return U * R * U.hermitian ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	171 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	172
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	173 DEFUN_DLD (sqrtm, args, nargout,
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	174 "-- texinfo --\n\
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	175 @deftypefn {Loadable Function} {[@var{result}, @var{error_estimate}] =} sqrtm (@var{a})\n\
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	176 Compute the matrix square root of the square matrix @var{a}.\n\
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	177 \n\
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	178 Ref: Nicholas J. Higham. A new sqrtm for MATLAB. Numerical Analysis\n\
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	179 Report No. 336, Manchester Centre for Computational Mathematics,\n\
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	180 Manchester, England, January 1999.\n\
5642 2618a0750ae6 [project @ 2006-03-06 21:26:48 by jwe] jwe parents: 5307 diff changeset	181 @seealso{expm, logm, funm}\n\
2618a0750ae6 [project @ 2006-03-06 21:26:48 by jwe] jwe parents: 5307 diff changeset	182 @end deftypefn")
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	183 {
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	184 octave_value_list retval;
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	185
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	186 int nargin = args.length ();
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	187
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	188 if (nargin != 1)
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	189 {
5823 080c08b192d8 [project @ 2006-05-19 05:32:17 by jwe] jwe parents: 5775 diff changeset	190 print_usage ();
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	191 return retval;
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	192 }
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	193
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	194 octave_value arg = args(0);
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	195
5275 23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe] jwe parents: 4331 diff changeset	196 octave_idx_type n = arg.rows ();
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe] jwe parents: 4331 diff changeset	197 octave_idx_type nc = arg.columns ();
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	198
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	199 int arg_is_empty = empty_arg ("sqrtm", n, nc);
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	200
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	201 if (arg_is_empty < 0)
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	202 return retval;
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	203 else if (arg_is_empty > 0)
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	204 return octave_value (Matrix ());
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	205
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	206 if (n != nc)
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	207 {
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	208 gripe_square_matrix_required ("sqrtm");
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	209 return retval;
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	210 }
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	211
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	212 retval(1) = lo_ieee_inf_value ();
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	213 retval(0) = lo_ieee_nan_value ();
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	214
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	215
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	216 if (arg.is_single_type ())
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	217 {
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	218 if (arg.is_real_scalar ())
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	219 {
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	220 float d = arg.float_value ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	221 if (d > 0.0)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	222 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	223 retval(0) = sqrt (d);
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	224 retval(1) = 0.0;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	225 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	226 else
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	227 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	228 retval(0) = FloatComplex (0.0, sqrt (d));
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	229 retval(1) = 0.0;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	230 }
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	231 }
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	232 else if (arg.is_complex_scalar ())
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	233 {
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	234 FloatComplex c = arg.float_complex_value ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	235 retval(0) = sqrt (c);
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	236 retval(1) = 0.0;
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	237 }
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	238 else if (arg.is_matrix_type ())
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	239 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	240 float err, minT;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	241
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	242 if (arg.is_real_matrix ())
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	243 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	244 FloatMatrix A = arg.float_matrix_value();
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	245
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	246 if (error_state)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	247 return retval;
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	248
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	249 // FIXME -- eventually, FloatComplexSCHUR will accept a
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	250 // real matrix arg.
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	251
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	252 FloatComplexMatrix Ac (A);
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	253
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	254 const FloatComplexSCHUR schur (Ac, std::string ());
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	255
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	256 if (error_state)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	257 return retval;
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	258
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	259 const FloatComplexMatrix U (schur.unitary_matrix ());
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	260 const FloatComplexMatrix T (schur.schur_matrix ());
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	261 const FloatComplexMatrix X (sqrtm_from_schur (U, T));
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	262
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	263 // Check for minimal imaginary part
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	264 float normX = 0.0;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	265 float imagX = 0.0;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	266 for (octave_idx_type i = 0; i < n; i++)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	267 for (octave_idx_type j = 0; j < n; j++)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	268 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	269 imagX = getmax (imagX, imag (X(i,j)));
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	270 normX = getmax (normX, abs (X(i,j)));
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	271 }
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	272
7795 df9519e9990c Handle single precision eps values David Bateman <dbateman@free.fr> parents: 7789 diff changeset	273 if (imagX < normX * 100 * FLT_EPSILON)
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	274 retval(0) = real (X);
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	275 else
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	276 retval(0) = X;
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	277
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	278 // Compute error
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	279 // FIXME can we estimate the error without doing the
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	280 // matrix multiply?
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	281
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	282 err = frobnorm (X*X - FloatComplexMatrix (A)) / frobnorm (A);
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	283
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	284 if (xisnan (err))
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	285 err = lo_ieee_float_inf_value ();
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	286
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	287 // Find min diagonal
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	288 minT = lo_ieee_float_inf_value ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	289 for (octave_idx_type i=0; i < n; i++)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	290 minT = getmin(minT, abs(T(i,i)));
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	291 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	292 else
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	293 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	294 FloatComplexMatrix A = arg.float_complex_matrix_value ();
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	295
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	296 if (error_state)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	297 return retval;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	298
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	299 const FloatComplexSCHUR schur (A, std::string ());
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	300
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	301 if (error_state)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	302 return retval;
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	303
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	304 const FloatComplexMatrix U (schur.unitary_matrix ());
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	305 const FloatComplexMatrix T (schur.schur_matrix ());
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	306 const FloatComplexMatrix X (sqrtm_from_schur (U, T));
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	307
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	308 retval(0) = X;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	309
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	310 err = frobnorm (X*X - A) / frobnorm (A);
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	311
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	312 if (xisnan (err))
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	313 err = lo_ieee_float_inf_value ();
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	314
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	315 minT = lo_ieee_float_inf_value ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	316 for (octave_idx_type i = 0; i < n; i++)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	317 minT = getmin (minT, abs (T(i,i)));
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	318 }
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	319
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	320 retval(1) = err;
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	321
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	322 if (nargout < 2)
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	323 {
7795 df9519e9990c Handle single precision eps values David Bateman <dbateman@free.fr> parents: 7789 diff changeset	324 if (err > 100(minT+FLT_EPSILON)n)
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	325 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	326 if (minT == 0.0)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	327 error ("sqrtm: A is singular, sqrt may not exist");
7795 df9519e9990c Handle single precision eps values David Bateman <dbateman@free.fr> parents: 7789 diff changeset	328 else if (minT <= sqrt (FLT_MIN))
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	329 error ("sqrtm: A is nearly singular, failed to find sqrt");
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	330 else
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	331 error ("sqrtm: failed to find sqrt");
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	332 }
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	333 }
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	334 }
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	335 }
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	336 else
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	337 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	338 if (arg.is_real_scalar ())
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	339 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	340 double d = arg.double_value ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	341 if (d > 0.0)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	342 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	343 retval(0) = sqrt (d);
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	344 retval(1) = 0.0;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	345 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	346 else
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	347 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	348 retval(0) = Complex (0.0, sqrt (d));
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	349 retval(1) = 0.0;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	350 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	351 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	352 else if (arg.is_complex_scalar ())
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	353 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	354 Complex c = arg.complex_value ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	355 retval(0) = sqrt (c);
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	356 retval(1) = 0.0;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	357 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	358 else if (arg.is_matrix_type ())
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	359 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	360 double err, minT;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	361
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	362 if (arg.is_real_matrix ())
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	363 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	364 Matrix A = arg.matrix_value();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	365
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	366 if (error_state)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	367 return retval;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	368
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	369 // FIXME -- eventually, ComplexSCHUR will accept a
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	370 // real matrix arg.
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	371
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	372 ComplexMatrix Ac (A);
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	373
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	374 const ComplexSCHUR schur (Ac, std::string ());
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	375
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	376 if (error_state)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	377 return retval;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	378
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	379 const ComplexMatrix U (schur.unitary_matrix ());
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	380 const ComplexMatrix T (schur.schur_matrix ());
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	381 const ComplexMatrix X (sqrtm_from_schur (U, T));
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	382
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	383 // Check for minimal imaginary part
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	384 double normX = 0.0;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	385 double imagX = 0.0;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	386 for (octave_idx_type i = 0; i < n; i++)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	387 for (octave_idx_type j = 0; j < n; j++)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	388 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	389 imagX = getmax (imagX, imag (X(i,j)));
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	390 normX = getmax (normX, abs (X(i,j)));
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	391 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	392
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	393 if (imagX < normX * 100 * DBL_EPSILON)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	394 retval(0) = real (X);
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	395 else
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	396 retval(0) = X;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	397
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	398 // Compute error
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	399 // FIXME can we estimate the error without doing the
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	400 // matrix multiply?
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	401
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	402 err = frobnorm (X*X - ComplexMatrix (A)) / frobnorm (A);
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	403
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	404 if (xisnan (err))
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	405 err = lo_ieee_inf_value ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	406
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	407 // Find min diagonal
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	408 minT = lo_ieee_inf_value ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	409 for (octave_idx_type i=0; i < n; i++)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	410 minT = getmin(minT, abs(T(i,i)));
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	411 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	412 else
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	413 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	414 ComplexMatrix A = arg.complex_matrix_value ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	415
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	416 if (error_state)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	417 return retval;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	418
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	419 const ComplexSCHUR schur (A, std::string ());
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	420
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	421 if (error_state)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	422 return retval;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	423
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	424 const ComplexMatrix U (schur.unitary_matrix ());
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	425 const ComplexMatrix T (schur.schur_matrix ());
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	426 const ComplexMatrix X (sqrtm_from_schur (U, T));
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	427
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	428 retval(0) = X;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	429
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	430 err = frobnorm (X*X - A) / frobnorm (A);
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	431
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	432 if (xisnan (err))
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	433 err = lo_ieee_inf_value ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	434
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	435 minT = lo_ieee_inf_value ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	436 for (octave_idx_type i = 0; i < n; i++)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	437 minT = getmin (minT, abs (T(i,i)));
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	438 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	439
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	440 retval(1) = err;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	441
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	442 if (nargout < 2)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	443 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	444 if (err > 100(minT+DBL_EPSILON)n)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	445 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	446 if (minT == 0.0)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	447 error ("sqrtm: A is singular, sqrt may not exist");
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	448 else if (minT <= sqrt (DBL_MIN))
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	449 error ("sqrtm: A is nearly singular, failed to find sqrt");
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	450 else
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	451 error ("sqrtm: failed to find sqrt");
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	452 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	453 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	454 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	455 else
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	456 gripe_wrong_type_arg ("sqrtm", arg);
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7017 diff changeset	457 }
4331 d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	458
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	459 return retval;
d3278845f764 [project @ 2003-02-19 00:12:31 by jwe] jwe parents: diff changeset	460 }
5307 4c8a2e4e0717 [project @ 2005-04-26 19:24:27 by jwe] jwe parents: 5275 diff changeset	461
4c8a2e4e0717 [project @ 2005-04-26 19:24:27 by jwe] jwe parents: 5275 diff changeset	462 /*
4c8a2e4e0717 [project @ 2005-04-26 19:24:27 by jwe] jwe parents: 5275 diff changeset	463 ;;; Local Variables: ***
4c8a2e4e0717 [project @ 2005-04-26 19:24:27 by jwe] jwe parents: 5275 diff changeset	464 ;;; mode: C++ ***
4c8a2e4e0717 [project @ 2005-04-26 19:24:27 by jwe] jwe parents: 5275 diff changeset	465 ;;; End: ***
4c8a2e4e0717 [project @ 2005-04-26 19:24:27 by jwe] jwe parents: 5275 diff changeset	466 */

Mercurial > octave

annotate src/DLD-FUNCTIONS/sqrtm.cc @ 7795:df9519e9990c