|
515
|
1 // f-log.cc -*- C++ -*- |
|
|
2 /* |
|
|
3 |
|
|
4 Copyright (C) 1994 John W. Eaton |
|
|
5 |
|
|
6 This file is part of Octave. |
|
|
7 |
|
|
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 2, or (at your option) any |
|
|
11 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. |
|
|
17 |
|
|
18 You should have received a copy of the GNU General Public License |
|
|
19 along with Octave; see the file COPYING. If not, write to the Free |
|
|
20 Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. |
|
|
21 |
|
|
22 */ |
|
|
23 |
|
|
24 #ifdef HAVE_CONFIG_H |
|
|
25 #include "config.h" |
|
|
26 #endif |
|
|
27 |
|
|
28 #include "EIG.h" |
|
|
29 |
|
|
30 #include "tree-const.h" |
|
|
31 #include "user-prefs.h" |
|
|
32 #include "error.h" |
|
|
33 #include "gripes.h" |
|
|
34 #include "f-log.h" |
|
|
35 |
|
|
36 // XXX FIXME XXX -- the next two functions (and expm) should really be just |
|
|
37 // one... |
|
|
38 |
|
|
39 Octave_object |
|
|
40 matrix_log (const tree_constant& a) |
|
|
41 { |
|
|
42 Octave_object retval (1); |
|
|
43 |
|
|
44 tree_constant tmp = a.make_numeric ();; |
|
|
45 |
|
|
46 if (tmp.rows () == 0 || tmp.columns () == 0) |
|
|
47 { |
|
|
48 int flag = user_pref.propagate_empty_matrices; |
|
|
49 if (flag != 0) |
|
|
50 { |
|
|
51 if (flag < 0) |
|
|
52 gripe_empty_arg ("logm", 0); |
|
|
53 Matrix m; |
|
|
54 retval(0) = m; |
|
|
55 return retval; |
|
|
56 } |
|
|
57 else |
|
|
58 gripe_empty_arg ("logm", 1); |
|
|
59 } |
|
|
60 |
|
|
61 switch (tmp.const_type ()) |
|
|
62 { |
|
|
63 case tree_constant_rep::matrix_constant: |
|
|
64 { |
|
|
65 Matrix m = tmp.matrix_value (); |
|
|
66 |
|
|
67 int nr = m.rows (); |
|
|
68 int nc = m.columns (); |
|
|
69 |
|
|
70 if (nr == 0 || nc == 0 || nr != nc) |
|
|
71 gripe_square_matrix_required ("logm"); |
|
|
72 else |
|
|
73 { |
|
|
74 EIG m_eig (m); |
|
|
75 ComplexColumnVector lambda (m_eig.eigenvalues ()); |
|
|
76 ComplexMatrix Q (m_eig.eigenvectors ()); |
|
|
77 |
|
|
78 for (int i = 0; i < nr; i++) |
|
|
79 { |
|
|
80 Complex elt = lambda.elem (i); |
|
|
81 if (imag (elt) == 0.0 && real (elt) > 0.0) |
|
|
82 lambda.elem (i) = log (real (elt)); |
|
|
83 else |
|
|
84 lambda.elem (i) = log (elt); |
|
|
85 } |
|
|
86 |
|
|
87 ComplexDiagMatrix D (lambda); |
|
|
88 ComplexMatrix result = Q * D * Q.inverse (); |
|
|
89 |
|
|
90 retval(0) = result; |
|
|
91 } |
|
|
92 } |
|
|
93 break; |
|
|
94 case tree_constant_rep::complex_matrix_constant: |
|
|
95 { |
|
|
96 ComplexMatrix m = tmp.complex_matrix_value (); |
|
|
97 |
|
|
98 int nr = m.rows (); |
|
|
99 int nc = m.columns (); |
|
|
100 |
|
|
101 if (nr == 0 || nc == 0 || nr != nc) |
|
|
102 gripe_square_matrix_required ("logm"); |
|
|
103 else |
|
|
104 { |
|
|
105 EIG m_eig (m); |
|
|
106 ComplexColumnVector lambda (m_eig.eigenvalues ()); |
|
|
107 ComplexMatrix Q (m_eig.eigenvectors ()); |
|
|
108 |
|
|
109 for (int i = 0; i < nr; i++) |
|
|
110 { |
|
|
111 Complex elt = lambda.elem (i); |
|
|
112 if (imag (elt) == 0.0 && real (elt) > 0.0) |
|
|
113 lambda.elem (i) = log (real (elt)); |
|
|
114 else |
|
|
115 lambda.elem (i) = log (elt); |
|
|
116 } |
|
|
117 |
|
|
118 ComplexDiagMatrix D (lambda); |
|
|
119 ComplexMatrix result = Q * D * Q.inverse (); |
|
|
120 |
|
|
121 retval(0) = result; |
|
|
122 } |
|
|
123 } |
|
|
124 break; |
|
|
125 case tree_constant_rep::scalar_constant: |
|
|
126 { |
|
|
127 double d = tmp.double_value (); |
|
|
128 if (d > 0.0) |
|
|
129 retval(0) = log (d); |
|
|
130 else |
|
|
131 { |
|
|
132 Complex dtmp (d); |
|
|
133 retval(0) = log (dtmp); |
|
|
134 } |
|
|
135 } |
|
|
136 break; |
|
|
137 case tree_constant_rep::complex_scalar_constant: |
|
|
138 { |
|
|
139 Complex c = tmp.complex_value (); |
|
|
140 retval(0) = log (c); |
|
|
141 } |
|
|
142 break; |
|
|
143 default: |
|
|
144 break; |
|
|
145 } |
|
|
146 return retval; |
|
|
147 } |
|
|
148 |
|
|
149 Octave_object |
|
|
150 matrix_sqrt (const tree_constant& a) |
|
|
151 { |
|
|
152 Octave_object retval (1); |
|
|
153 |
|
|
154 tree_constant tmp = a.make_numeric ();; |
|
|
155 |
|
|
156 if (tmp.rows () == 0 || tmp.columns () == 0) |
|
|
157 { |
|
|
158 int flag = user_pref.propagate_empty_matrices; |
|
|
159 if (flag != 0) |
|
|
160 { |
|
|
161 if (flag < 0) |
|
|
162 gripe_empty_arg ("sqrtm", 0); |
|
|
163 Matrix m; |
|
|
164 retval(0) = m; |
|
|
165 return retval; |
|
|
166 } |
|
|
167 else |
|
|
168 gripe_empty_arg ("sqrtm", 1); |
|
|
169 } |
|
|
170 |
|
|
171 switch (tmp.const_type ()) |
|
|
172 { |
|
|
173 case tree_constant_rep::matrix_constant: |
|
|
174 { |
|
|
175 Matrix m = tmp.matrix_value (); |
|
|
176 |
|
|
177 int nr = m.rows (); |
|
|
178 int nc = m.columns (); |
|
|
179 |
|
|
180 if (nr == 0 || nc == 0 || nr != nc) |
|
|
181 gripe_square_matrix_required ("sqrtm"); |
|
|
182 else |
|
|
183 { |
|
|
184 EIG m_eig (m); |
|
|
185 ComplexColumnVector lambda (m_eig.eigenvalues ()); |
|
|
186 ComplexMatrix Q (m_eig.eigenvectors ()); |
|
|
187 |
|
|
188 for (int i = 0; i < nr; i++) |
|
|
189 { |
|
|
190 Complex elt = lambda.elem (i); |
|
|
191 if (imag (elt) == 0.0 && real (elt) > 0.0) |
|
|
192 lambda.elem (i) = sqrt (real (elt)); |
|
|
193 else |
|
|
194 lambda.elem (i) = sqrt (elt); |
|
|
195 } |
|
|
196 |
|
|
197 ComplexDiagMatrix D (lambda); |
|
|
198 ComplexMatrix result = Q * D * Q.inverse (); |
|
|
199 |
|
|
200 retval(0) = result; |
|
|
201 } |
|
|
202 } |
|
|
203 break; |
|
|
204 case tree_constant_rep::complex_matrix_constant: |
|
|
205 { |
|
|
206 ComplexMatrix m = tmp.complex_matrix_value (); |
|
|
207 |
|
|
208 int nr = m.rows (); |
|
|
209 int nc = m.columns (); |
|
|
210 |
|
|
211 if (nr == 0 || nc == 0 || nr != nc) |
|
|
212 gripe_square_matrix_required ("sqrtm"); |
|
|
213 else |
|
|
214 { |
|
|
215 EIG m_eig (m); |
|
|
216 ComplexColumnVector lambda (m_eig.eigenvalues ()); |
|
|
217 ComplexMatrix Q (m_eig.eigenvectors ()); |
|
|
218 |
|
|
219 for (int i = 0; i < nr; i++) |
|
|
220 { |
|
|
221 Complex elt = lambda.elem (i); |
|
|
222 if (imag (elt) == 0.0 && real (elt) > 0.0) |
|
|
223 lambda.elem (i) = sqrt (real (elt)); |
|
|
224 else |
|
|
225 lambda.elem (i) = sqrt (elt); |
|
|
226 } |
|
|
227 |
|
|
228 ComplexDiagMatrix D (lambda); |
|
|
229 ComplexMatrix result = Q * D * Q.inverse (); |
|
|
230 |
|
|
231 retval(0) = result; |
|
|
232 } |
|
|
233 } |
|
|
234 break; |
|
|
235 case tree_constant_rep::scalar_constant: |
|
|
236 { |
|
|
237 double d = tmp.double_value (); |
|
|
238 if (d > 0.0) |
|
|
239 retval(0) = sqrt (d); |
|
|
240 else |
|
|
241 { |
|
|
242 Complex dtmp (d); |
|
|
243 retval(0) = sqrt (dtmp); |
|
|
244 } |
|
|
245 } |
|
|
246 break; |
|
|
247 case tree_constant_rep::complex_scalar_constant: |
|
|
248 { |
|
|
249 Complex c = tmp.complex_value (); |
|
|
250 retval(0) = log (c); |
|
|
251 } |
|
|
252 break; |
|
|
253 default: |
|
|
254 break; |
|
|
255 } |
|
|
256 return retval; |
|
|
257 } |
|
|
258 |
|
|
259 /* |
|
|
260 ;;; Local Variables: *** |
|
|
261 ;;; mode: C++ *** |
|
|
262 ;;; page-delimiter: "^/\\*" *** |
|
|
263 ;;; End: *** |
|
|
264 */ |
