GNU Octave  4.2.1
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
ov-base-diag.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2008-2017 Jaroslav Hajek
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 // This file should not include config.h. It is only included in other
24 // C++ source files that should have included config.h before including
25 // this file.
26 
27 #include <iostream>
28 
29 #include "mach-info.h"
30 #include "lo-ieee.h"
31 
32 #include "ov-base-diag.h"
33 #include "mxarray.h"
34 #include "ov-base.h"
35 #include "ov-base-mat.h"
36 #include "pr-output.h"
37 #include "error.h"
38 #include "errwarn.h"
39 #include "oct-stream.h"
40 #include "ops.h"
41 
42 #include "ls-oct-text.h"
43 
44 template <typename DMT, typename MT>
47  const std::list<octave_value_list>& idx)
48 {
50 
51  switch (type[0])
52  {
53  case '(':
54  retval = do_index_op (idx.front ());
55  break;
56 
57  case '{':
58  case '.':
59  {
60  std::string nm = type_name ();
61  error ("%s cannot be indexed with %c", nm.c_str (), type[0]);
62  }
63  break;
64 
65  default:
67  }
68 
69  return retval.next_subsref (type, idx);
70 }
71 
72 template <typename DMT, typename MT>
75 {
77  if (matrix.rows () == 1 || matrix.cols () == 1)
78  {
79  // Rather odd special case. This is a row or column vector
80  // represented as a diagonal matrix with a single nonzero entry, but
81  // Fdiag semantics are to product a diagonal matrix for vector
82  // inputs.
83  if (k == 0)
84  // Returns Diag2Array<T> with nnz <= 1.
85  retval = matrix.build_diag_matrix ();
86  else
87  // Returns Array<T> matrix
88  retval = matrix.array_value ().diag (k);
89  }
90  else
91  // Returns Array<T> vector
92  retval = matrix.extract_diag (k);
93  return retval;
94 }
95 
96 template <typename DMT, typename MT>
99  bool resize_ok)
100 {
102 
103  if (idx.length () == 2 && ! resize_ok)
104  {
105  int k = 0; // index we're accesing when index_vector throws
106  try
107  {
108  idx_vector idx0 = idx(0).index_vector ();
109  k = 1;
110  idx_vector idx1 = idx(1).index_vector ();
111 
112  if (idx0.is_scalar () && idx1.is_scalar ())
113  {
114  retval = matrix.checkelem (idx0(0), idx1(0));
115  }
116  else
117  {
118  octave_idx_type m = idx0.length (matrix.rows ());
119  octave_idx_type n = idx1.length (matrix.columns ());
120  if (idx0.is_colon_equiv (m) && idx1.is_colon_equiv (n)
121  && m <= matrix.rows () && n <= matrix.rows ())
122  {
123  DMT rm (matrix);
124  rm.resize (m, n);
125  retval = rm;
126  }
127  else
128  retval = to_dense ().do_index_op (idx, resize_ok);
129  }
130  }
131  catch (octave::index_exception& e)
132  {
133  // Rethrow to allow more info to be reported later.
134  e.set_pos_if_unset (2, k+1);
135  throw;
136  }
137  }
138  else
139  retval = to_dense ().do_index_op (idx, resize_ok);
140 
141  return retval;
142 }
143 
144 template <typename DMT, typename MT>
147  const std::list<octave_value_list>& idx,
148  const octave_value& rhs)
149 {
151 
152  switch (type[0])
153  {
154  case '(':
155  {
156  if (type.length () != 1)
157  {
158  std::string nm = type_name ();
159  error ("in indexed assignment of %s, last lhs index must be ()",
160  nm.c_str ());
161  }
162 
163  octave_value_list jdx = idx.front ();
164 
165  // FIXME: Mostly repeated code for cases 1 and 2 could be
166  // consolidated for DRY (Don't Repeat Yourself).
167  // Check for assignments to diagonal elements which should not
168  // destroy the diagonal property of the matrix.
169  // If D is a diagonal matrix then the assignment can be
170  // 1) linear, D(i) = x, where ind2sub results in case #2 below
171  // 2) subscript D(i,i) = x, where both indices are equal.
172  if (jdx.length () == 1 && jdx(0).is_scalar_type ())
173  {
174  typename DMT::element_type val;
175  int k = 0;
176  try
177  {
178  idx_vector ind = jdx(0).index_vector ();
179  k = 1;
180  dim_vector dv (matrix.rows (), matrix.cols ());
181  Array<idx_vector> ivec = ind2sub (dv, ind);
182  idx_vector i0 = ivec(0);
183  idx_vector i1 = ivec(1);
184 
185  if (i0(0) == i1(0)
186  && chk_valid_scalar (rhs, val))
187  {
188  matrix.dgelem (i0(0)) = val;
189  retval = this;
190  this->count++;
191  // invalidate cache
192  dense_cache = octave_value ();
193  }
194  }
195  catch (octave::index_exception& e)
196  {
197  // Rethrow to allow more info to be reported later.
198  e.set_pos_if_unset (2, k+1);
199  throw;
200  }
201  }
202  else if (jdx.length () == 2
203  && jdx(0).is_scalar_type () && jdx(1).is_scalar_type ())
204  {
205  typename DMT::element_type val;
206  int k = 0;
207  try
208  {
209  idx_vector i0 = jdx(0).index_vector ();
210  k = 1;
211  idx_vector i1 = jdx(1).index_vector ();
212  if (i0(0) == i1(0)
213  && i0(0) < matrix.rows () && i1(0) < matrix.cols ()
214  && chk_valid_scalar (rhs, val))
215  {
216  matrix.dgelem (i0(0)) = val;
217  retval = this;
218  this->count++;
219  // invalidate cache
220  dense_cache = octave_value ();
221  }
222  }
223  catch (octave::index_exception& e)
224  {
225  // Rethrow to allow more info to be reported later.
226  e.set_pos_if_unset (2, k+1);
227  throw;
228  }
229  }
230 
231  if (! retval.is_defined ())
232  retval = numeric_assign (type, idx, rhs);
233  }
234  break;
235 
236  case '{':
237  case '.':
238  {
239  if (! is_empty ())
240  {
241  std::string nm = type_name ();
242  error ("%s cannot be indexed with %c", nm.c_str (), type[0]);
243  }
244 
246 
247  retval = tmp.subsasgn (type, idx, rhs);
248  }
249  break;
250 
251  default:
252  panic_impossible ();
253  }
254 
255  return retval;
256 }
257 
258 template <typename DMT, typename MT>
261 {
263  if (dv.ndims () == 2)
264  {
265  DMT rm (matrix);
266  rm.resize (dv(0), dv(1));
267  retval = rm;
268  }
269  else
270  retval = to_dense ().resize (dv, fill);
271  return retval;
272 }
273 
274 // Return true if this matrix has all true elements (non-zero, not NA/NaN).
275 template <typename DMT, typename MT>
276 bool
278 {
279  if (dims ().numel () > 1)
280  {
282  // Throw error if any NaN or NA by calling is_true().
283  octave_value (matrix.extract_diag ()).is_true ();
284  return false; // > 1x1 diagonal always has zeros
285  }
286  else
287  return to_dense ().is_true (); // 0x0 or 1x1, handle NaN etc.
288 }
289 
290 // FIXME: This should be achieveable using ::real
291 template <typename T> inline T helper_getreal (T x) { return x; }
292 template <typename T> inline T helper_getreal (std::complex<T> x)
293 { return x.real (); }
294 // FIXME: We really need some traits so that ad hoc hooks like this
295 // are not necessary.
296 template <typename T> inline T helper_iscomplex (T) { return false; }
297 template <typename T> inline T helper_iscomplex (std::complex<T>) { return true; }
298 
299 template <typename DMT, typename MT>
300 double
301 octave_base_diag<DMT, MT>::double_value (bool force_conversion) const
302 {
303  typedef typename DMT::element_type el_type;
304 
305  if (helper_iscomplex (el_type ()) && ! force_conversion)
306  warn_implicit_conversion ("Octave:imag-to-real",
307  "complex matrix", "real scalar");
308 
309  if (is_empty ())
310  err_invalid_conversion (type_name (), "real scalar");
311 
312  warn_implicit_conversion ("Octave:array-to-scalar",
313  type_name (), "real scalar");
314 
315  return helper_getreal (el_type (matrix (0, 0)));
316 }
317 
318 template <typename DMT, typename MT>
319 float
320 octave_base_diag<DMT, MT>::float_value (bool force_conversion) const
321 {
322  typedef typename DMT::element_type el_type;
323 
324  if (helper_iscomplex (el_type ()) && ! force_conversion)
325  warn_implicit_conversion ("Octave:imag-to-real",
326  "complex matrix", "real scalar");
327 
328  if (! (numel () > 0))
329  err_invalid_conversion (type_name (), "real scalar");
330 
331  warn_implicit_conversion ("Octave:array-to-scalar",
332  type_name (), "real scalar");
333 
334  return helper_getreal (el_type (matrix (0, 0)));
335 }
336 
337 template <typename DMT, typename MT>
338 Complex
340 {
341  if (rows () == 0 || columns () == 0)
342  err_invalid_conversion (type_name (), "complex scalar");
343 
344  warn_implicit_conversion ("Octave:array-to-scalar",
345  type_name (), "complex scalar");
346 
347  return matrix(0, 0);
348 }
349 
350 template <typename DMT, typename MT>
353 {
354  float tmp = lo_ieee_float_nan_value ();
355 
356  FloatComplex retval (tmp, tmp);
357 
358  if (rows () == 0 || columns () == 0)
359  err_invalid_conversion (type_name (), "complex scalar");
360 
361  warn_implicit_conversion ("Octave:array-to-scalar",
362  type_name (), "complex scalar");
363 
364  retval = matrix (0, 0);
365 
366  return retval;
367 }
368 
369 template <typename DMT, typename MT>
370 Matrix
372 {
373  return Matrix (diag_matrix_value ());
374 }
375 
376 template <typename DMT, typename MT>
379 {
380  return FloatMatrix (float_diag_matrix_value ());
381 }
382 
383 template <typename DMT, typename MT>
386 {
387  return ComplexMatrix (complex_diag_matrix_value ());
388 }
389 
390 template <typename DMT, typename MT>
393 {
394  return FloatComplexMatrix (float_complex_diag_matrix_value ());
395 }
396 
397 template <typename DMT, typename MT>
398 NDArray
400 {
401  return NDArray (matrix_value ());
402 }
403 
404 template <typename DMT, typename MT>
407 {
408  return FloatNDArray (float_matrix_value ());
409 }
410 
411 template <typename DMT, typename MT>
414 {
415  return ComplexNDArray (complex_matrix_value ());
416 }
417 
418 template <typename DMT, typename MT>
421 {
422  return FloatComplexNDArray (float_complex_matrix_value ());
423 }
424 
425 template <typename DMT, typename MT>
428 {
429  return to_dense ().bool_array_value (warn);
430 }
431 
432 template <typename DMT, typename MT>
435 {
436  return to_dense ().char_array_value (warn);
437 }
438 
439 template <typename DMT, typename MT>
442 {
443  return SparseMatrix (diag_matrix_value ());
444 }
445 
446 template <typename DMT, typename MT>
449 {
450  return SparseComplexMatrix (complex_diag_matrix_value ());
451 }
452 
453 template <typename DMT, typename MT>
455 octave_base_diag<DMT, MT>::index_vector (bool require_integers) const
456 {
457  return to_dense ().index_vector (require_integers);
458 }
459 
460 template <typename DMT, typename MT>
463  char type) const
464 {
465  return to_dense ().convert_to_str_internal (pad, force, type);
466 }
467 
468 template <typename DMT, typename MT>
469 bool
471 {
472  os << "# rows: " << matrix.rows () << "\n"
473  << "# columns: " << matrix.columns () << "\n";
474 
475  os << matrix.extract_diag ();
476 
477  return true;
478 }
479 
480 template <typename DMT, typename MT>
481 bool
483 {
484  octave_idx_type r = 0;
485  octave_idx_type c = 0;
486 
487  if (! extract_keyword (is, "rows", r, true)
488  || ! extract_keyword (is, "columns", c, true))
489  error ("load: failed to extract number of rows and columns");
490 
491  octave_idx_type l = r < c ? r : c;
492  MT tmp (l, 1);
493  is >> tmp;
494 
495  if (! is)
496  error ("load: failed to load diagonal matrix constant");
497 
498  // This is a little tricky, as we have the Matrix type, but
499  // not ColumnVector type. We need to help the compiler get
500  // through the inheritance tree.
501  typedef typename DMT::element_type el_type;
503  matrix.resize (r, c);
504 
505  // Invalidate cache. Probably not necessary, but safe.
506  dense_cache = octave_value ();
507 
508  return true;
509 }
510 
511 template <typename DMT, typename MT>
512 void
514  bool pr_as_read_syntax) const
515 {
516  return octave_print_internal (os, matrix, pr_as_read_syntax,
517  current_print_indent_level ());
518 }
519 
520 template <typename DMT, typename MT>
521 mxArray *
523 {
524  return to_dense ().as_mxArray ();
525 }
526 
527 template <typename DMT, typename MT>
528 bool
530 {
531  dim_vector dv = dims ();
532 
533  return (dv.all_ones () || dv.any_zero ());
534 }
535 
536 template <typename DMT, typename MT>
537 void
538 octave_base_diag<DMT, MT>::print (std::ostream& os, bool pr_as_read_syntax)
539 {
540  print_raw (os, pr_as_read_syntax);
541  newline (os);
542 }
543 template <typename DMT, typename MT>
544 int
546  oct_data_conv::data_type output_type,
547  int skip,
549 {
550  return to_dense ().write (os, block_size, output_type, skip, flt_fmt);
551 }
552 
553 template <typename DMT, typename MT>
554 void
556  const std::string& prefix) const
557 {
558  matrix.print_info (os, prefix);
559 }
560 
561 template <typename DMT, typename MT>
564 {
565  if (n < matrix.numel ())
566  {
567  octave_idx_type nr = matrix.rows ();
568 
569  octave_idx_type r = n % nr;
570  octave_idx_type c = n / nr;
571 
572  return octave_value (matrix.elem (r, c));
573  }
574  else
575  return octave_value ();
576 }
577 
578 template <typename DMT, typename MT>
581 {
582  if (! dense_cache.is_defined ())
583  dense_cache = MT (matrix);
584 
585  return dense_cache;
586 }
void print(std::ostream &os, bool pr_as_read_syntax=false)
bool is_true(const std::string &s)
Definition: mkoctfile.cc:398
octave_value to_dense(void) const
octave_idx_type length(octave_idx_type n=0) const
Definition: idx-vector.h:541
float float_value(bool=false) const
NDArray array_value(bool=false) const
T helper_iscomplex(T)
double double_value(bool=false) const
mxArray * as_mxArray(void) const
idx_vector index_vector(bool=false) const
ind
Definition: sub2ind.cc:107
int write(octave_stream &os, int block_size, oct_data_conv::data_type output_type, int skip, octave::mach_info::float_format flt_fmt) const
identity matrix If supplied two scalar respectively For allows like xample val
Definition: data.cc:5068
bool load_ascii(std::istream &is)
octave_value resize(const dim_vector &dv, bool fill=false) const
void set_pos_if_unset(octave_idx_type nd_arg, octave_idx_type dim_arg)
octave_idx_type length(void) const
Definition: ovl.h:96
bool is_defined(void) const
Definition: ov.h:536
for large enough k
Definition: lu.cc:606
the sparsity preserving column transformation such that that defines the pivoting threshold can be given in which case it defines the then the first element defines the pivoting tolerance for the unsymmetric the values defined such that for full matrix
Definition: lu.cc:138
SparseComplexMatrix sparse_complex_matrix_value(bool=false) const
void error(const char *fmt,...)
Definition: error.cc:570
charNDArray char_array_value(bool=false) const
SparseMatrix sparse_matrix_value(bool=false) const
octave::mach_info::float_format flt_fmt
Definition: load-save.cc:723
octave_value resize(const dim_vector &dv, bool fill=false) const
Definition: ov.h:524
octave_value convert_to_str_internal(bool pad, bool force, char type) const
Definition: ov.h:1207
Template for N-dimensional array classes with like-type math operators.
Definition: MArray.h:32
bool all_ones(void) const
Definition: dim-vector.h:377
FloatComplex float_complex_value(bool=false) const
i e
Definition: data.cc:2724
octave_value do_index_op(const octave_value_list &idx, bool resize_ok=false)
Definition: ov-base-diag.cc:98
bool print_as_scalar(void) const
octave_value convert_to_str_internal(bool pad, bool force, char type) const
octave_value subsasgn(const std::string &type, const std::list< octave_value_list > &idx, const octave_value &rhs)
Definition: ov.cc:1540
OCTAVE_EXPORT octave_value_list isdir nd deftypefn *std::string nm
Definition: utils.cc:941
octave_value subsasgn(const std::string &type, const std::list< octave_value_list > &idx, const octave_value &rhs)
std::string extract_keyword(std::istream &is, const char *keyword, const bool next_only)
Definition: ls-oct-text.cc:80
nd deftypefn *octave_map m
Definition: ov-struct.cc:2058
double tmp
Definition: data.cc:6300
octave_value retval
Definition: data.cc:6294
#define panic_impossible()
Definition: error.h:40
Template for two dimensional diagonal array with math operators.
Definition: MDiagArray2.h:33
idx type
Definition: ov.cc:3129
ComplexMatrix complex_matrix_value(bool=false) const
Definition: dMatrix.h:37
bool any_zero(void) const
Definition: dim-vector.h:359
Complex complex_value(bool=false) const
the sparsity preserving column transformation such that that defines the pivoting threshold can be given in which case it defines the c
Definition: lu.cc:138
the exceeded dimensions are set to if fewer subscripts than dimensions are the exceeding dimensions are merged into the final requested dimension For consider the following dims
Definition: sub2ind.cc:255
float lo_ieee_float_nan_value(void)
Definition: lo-ieee.cc:202
void print_info(std::ostream &os, const std::string &prefix) const
bool is_colon_equiv(octave_idx_type n) const
Definition: idx-vector.h:574
N Dimensional Array with copy-on-write semantics.
Definition: Array.h:126
Array< idx_vector > ind2sub(const dim_vector &dv, const idx_vector &idx)
Definition: Array-util.cc:619
FloatComplexNDArray float_complex_array_value(bool=false) const
T::size_type numel(const T &str)
Definition: oct-string.cc:61
void err_invalid_conversion(const std::string &from, const std::string &to)
Definition: errwarn.cc:68
FloatMatrix float_matrix_value(bool=false) const
octave_value subsref(const std::string &type, const std::list< octave_value_list > &idx)
Definition: ov-base-diag.cc:46
void print_raw(std::ostream &os, bool pr_as_read_syntax=false) const
Matrix matrix_value(bool=false) const
ComplexNDArray complex_array_value(bool=false) const
FloatNDArray float_array_value(bool=false) const
void octave_print_internal(std::ostream &, char, bool)
Definition: pr-output.cc:1722
void warn_implicit_conversion(const char *id, const char *from, const char *to)
Definition: errwarn.cc:332
bool is_scalar(void) const
Definition: idx-vector.h:568
octave_value fast_elem_extract(octave_idx_type n) const
boolNDArray bool_array_value(bool warn=false) const
octave_idx_type ndims(void) const
Number of dimensions.
Definition: dim-vector.h:301
bool save_ascii(std::ostream &os)
void warn_array_as_logical(const dim_vector &dv)
Definition: errwarn.cc:276
std::complex< float > FloatComplex
Definition: oct-cmplx.h:32
FloatComplexMatrix float_complex_matrix_value(bool=false) const
std::complex< double > Complex
Definition: oct-cmplx.h:31
write the output to stdout if nargout is
Definition: load-save.cc:1576
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:87
If this string is the system will ring the terminal sometimes it is useful to be able to print the original representation of the string
Definition: utils.cc:854
dim_vector dv
Definition: sub2ind.cc:263
octave_value next_subsref(const std::string &type, const std::list< octave_value_list > &idx, size_t skip=1)
Definition: ov.cc:1462
octave_value diag(octave_idx_type k=0) const
Definition: ov-base-diag.cc:74
F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T F77_DBLE &F77_RET_T F77_REAL &F77_RET_T F77_REAL &F77_RET_T F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T F77_REAL F77_REAL &F77_RET_T F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE * x
return octave_value(v1.char_array_value().concat(v2.char_array_value(), ra_idx),((a1.is_sq_string()||a2.is_sq_string())? '\'': '"'))
static octave_value empty_conv(const std::string &type, const octave_value &rhs=octave_value())
Definition: ov.cc:2882
bool is_true(void) const
T helper_getreal(T x)
octave_value do_index_op(const octave_value_list &idx, bool resize_ok=false)
Definition: ov.h:454