GNU Octave  4.0.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
sylvester.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1996-2015 John W. Eaton
4 
5 This file is part of Octave.
6 
7 Octave is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11 
12 Octave is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with Octave; see the file COPYING. If not, see
19 <http://www.gnu.org/licenses/>.
20 
21 */
22 
23 // Author: A. S. Hodel <scotte@eng.auburn.edu>
24 
25 #ifdef HAVE_CONFIG_H
26 #include <config.h>
27 #endif
28 
29 #include "defun.h"
30 #include "error.h"
31 #include "gripes.h"
32 #include "oct-obj.h"
33 #include "utils.h"
34 
35 DEFUN (sylvester, args, nargout,
36  "-*- texinfo -*-\n\
37 @deftypefn {Built-in Function} {@var{X} =} syl (@var{A}, @var{B}, @var{C})\n\
38 Solve the Sylvester equation\n\
39 @tex\n\
40 $$\n\
41  A X + X B = C\n\
42 $$\n\
43 @end tex\n\
44 @ifnottex\n\
45 \n\
46 @example\n\
47 A X + X B = C\n\
48 @end example\n\
49 \n\
50 @end ifnottex\n\
51 using standard @sc{lapack} subroutines.\n\
52 \n\
53 For example:\n\
54 \n\
55 @example\n\
56 @group\n\
57 sylvester ([1, 2; 3, 4], [5, 6; 7, 8], [9, 10; 11, 12])\n\
58  @result{} [ 0.50000, 0.66667; 0.66667, 0.50000 ]\n\
59 @end group\n\
60 @end example\n\
61 @end deftypefn")
62 {
63  octave_value retval;
64 
65  int nargin = args.length ();
66 
67  if (nargin != 3 || nargout > 1)
68  {
69  print_usage ();
70  return retval;
71  }
72 
73  octave_value arg_a = args(0);
74  octave_value arg_b = args(1);
75  octave_value arg_c = args(2);
76 
77  octave_idx_type a_nr = arg_a.rows ();
78  octave_idx_type a_nc = arg_a.columns ();
79 
80  octave_idx_type b_nr = arg_b.rows ();
81  octave_idx_type b_nc = arg_b.columns ();
82 
83  octave_idx_type c_nr = arg_c.rows ();
84  octave_idx_type c_nc = arg_c.columns ();
85 
86  int arg_a_is_empty = empty_arg ("sylvester", a_nr, a_nc);
87  int arg_b_is_empty = empty_arg ("sylvester", b_nr, b_nc);
88  int arg_c_is_empty = empty_arg ("sylvester", c_nr, c_nc);
89 
90  bool isfloat = arg_a.is_single_type ()
91  || arg_b.is_single_type ()
92  || arg_c.is_single_type ();
93 
94  if (arg_a_is_empty > 0 && arg_b_is_empty > 0 && arg_c_is_empty > 0)
95  if (isfloat)
96  return octave_value (FloatMatrix ());
97  else
98  return octave_value (Matrix ());
99  else if (arg_a_is_empty || arg_b_is_empty || arg_c_is_empty)
100  return retval;
101 
102  // Arguments are not empty, so check for correct dimensions.
103 
104  if (a_nr != a_nc)
105  {
106  gripe_square_matrix_required ("sylvester: input A");
107  return retval;
108  }
109  else if (b_nr != b_nc)
110  {
111  gripe_square_matrix_required ("sylvester: input B");
112  return retval;
113  }
114  else if (a_nr != c_nr || b_nr != c_nc)
115  {
117  return retval;
118  }
119 
120  if (isfloat)
121  {
122  if (arg_a.is_complex_type ()
123  || arg_b.is_complex_type ()
124  || arg_c.is_complex_type ())
125  {
126  // Do everything in complex arithmetic;
127 
129 
130  if (error_state)
131  return retval;
132 
134 
135  if (error_state)
136  return retval;
137 
139 
140  if (error_state)
141  return retval;
142 
143  retval = Sylvester (ca, cb, cc);
144  }
145  else
146  {
147  // Do everything in real arithmetic.
148 
149  FloatMatrix ca = arg_a.float_matrix_value ();
150 
151  if (error_state)
152  return retval;
153 
154  FloatMatrix cb = arg_b.float_matrix_value ();
155 
156  if (error_state)
157  return retval;
158 
159  FloatMatrix cc = arg_c.float_matrix_value ();
160 
161  if (error_state)
162  return retval;
163 
164  retval = Sylvester (ca, cb, cc);
165  }
166  }
167  else
168  {
169  if (arg_a.is_complex_type ()
170  || arg_b.is_complex_type ()
171  || arg_c.is_complex_type ())
172  {
173  // Do everything in complex arithmetic;
174 
175  ComplexMatrix ca = arg_a.complex_matrix_value ();
176 
177  if (error_state)
178  return retval;
179 
180  ComplexMatrix cb = arg_b.complex_matrix_value ();
181 
182  if (error_state)
183  return retval;
184 
185  ComplexMatrix cc = arg_c.complex_matrix_value ();
186 
187  if (error_state)
188  return retval;
189 
190  retval = Sylvester (ca, cb, cc);
191  }
192  else
193  {
194  // Do everything in real arithmetic.
195 
196  Matrix ca = arg_a.matrix_value ();
197 
198  if (error_state)
199  return retval;
200 
201  Matrix cb = arg_b.matrix_value ();
202 
203  if (error_state)
204  return retval;
205 
206  Matrix cc = arg_c.matrix_value ();
207 
208  if (error_state)
209  return retval;
210 
211  retval = Sylvester (ca, cb, cc);
212  }
213  }
214 
215  return retval;
216 }
217 
218 /*
219 %!assert (sylvester ([1, 2; 3, 4], [5, 6; 7, 8], [9, 10; 11, 12]), [1/2, 2/3; 2/3, 1/2], sqrt (eps))
220 %!assert (sylvester (single ([1, 2; 3, 4]), single ([5, 6; 7, 8]), single ([9, 10; 11, 12])), single ([1/2, 2/3; 2/3, 1/2]), sqrt (eps ("single")))
221 
222 %% Test input validation
223 %!error sylvester ()
224 %!error sylvester (1)
225 %!error sylvester (1,2)
226 %!error sylvester (1, 2, 3, 4)
227 %!error <input A: .* must be a square matrix> sylvester (ones (2,3), ones (2,2), ones (2,2))
228 %!error <input B: .* must be a square matrix> sylvester (ones (2,2), ones (2,3), ones (2,2))
229 %!error <nonconformant matrices> sylvester (ones (2,2), ones (2,2), ones (3,3))
230 */
void gripe_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
octave_idx_type rows(void) const
Definition: ov.h:473
FloatComplexMatrix float_complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:795
OCTINTERP_API void print_usage(void)
Definition: defun.cc:51
ComplexMatrix Sylvester(const ComplexMatrix &a, const ComplexMatrix &b, const ComplexMatrix &c)
Definition: CMatrix.cc:3588
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:44
int empty_arg(const char *, octave_idx_type nr, octave_idx_type nc)
Definition: utils.cc:246
void gripe_square_matrix_required(const char *name)
Definition: gripes.cc:81
octave_idx_type columns(void) const
Definition: ov.h:475
int error_state
Definition: error.cc:101
bool is_complex_type(void) const
Definition: ov.h:654
octave_idx_type length(void) const
Definition: ov.cc:1525
Definition: dMatrix.h:35
Matrix matrix_value(bool frc_str_conv=false) const
Definition: ov.h:773
ComplexMatrix complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:791
FloatMatrix float_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:776
bool is_single_type(void) const
Definition: ov.h:611
return octave_value(v1.char_array_value().concat(v2.char_array_value(), ra_idx),((a1.is_sq_string()||a2.is_sq_string())? '\'': '"'))