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
tril.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2004-2017 David Bateman
4 Copyright (C) 2009 VZLU Prague
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 3 of the License, or (at your
11 option) any later version.
12 
13 Octave is distributed in the hope that it will be useful, but WITHOUT
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 for more details.
17 
18 You should have received a copy of the GNU General Public License
19 along with Octave; see the file COPYING. If not, see
20 <http://www.gnu.org/licenses/>.
21 
22 */
23 
24 #if defined (HAVE_CONFIG_H)
25 # include "config.h"
26 #endif
27 
28 #include <algorithm>
29 #include "Array.h"
30 #include "Sparse.h"
31 #include "mx-base.h"
32 
33 #include "ov.h"
34 #include "Cell.h"
35 
36 #include "defun.h"
37 #include "error.h"
38 #include "ovl.h"
39 
40 // The bulk of the work.
41 template <typename T>
42 static Array<T>
43 do_tril (const Array<T>& a, octave_idx_type k, bool pack)
44 {
45  octave_idx_type nr = a.rows ();
46  octave_idx_type nc = a.columns ();
47  const T *avec = a.fortran_vec ();
49 
50  if (pack)
51  {
52  octave_idx_type j1 = std::min (std::max (zero, k), nc);
53  octave_idx_type j2 = std::min (std::max (zero, nr + k), nc);
54  octave_idx_type n = j1 * nr + ((j2 - j1) * (nr-(j1-k) + nr-(j2-1-k))) / 2;
55  Array<T> r (dim_vector (n, 1));
56  T *rvec = r.fortran_vec ();
57  for (octave_idx_type j = 0; j < nc; j++)
58  {
59  octave_idx_type ii = std::min (std::max (zero, j - k), nr);
60  rvec = std::copy (avec + ii, avec + nr, rvec);
61  avec += nr;
62  }
63 
64  return r;
65  }
66  else
67  {
68  Array<T> r (a.dims ());
69  T *rvec = r.fortran_vec ();
70  for (octave_idx_type j = 0; j < nc; j++)
71  {
72  octave_idx_type ii = std::min (std::max (zero, j - k), nr);
73  std::fill (rvec, rvec + ii, T ());
74  std::copy (avec + ii, avec + nr, rvec + ii);
75  avec += nr;
76  rvec += nr;
77  }
78 
79  return r;
80  }
81 }
82 
83 template <typename T>
84 static Array<T>
85 do_triu (const Array<T>& a, octave_idx_type k, bool pack)
86 {
87  octave_idx_type nr = a.rows ();
88  octave_idx_type nc = a.columns ();
89  const T *avec = a.fortran_vec ();
91 
92  if (pack)
93  {
94  octave_idx_type j1 = std::min (std::max (zero, k), nc);
95  octave_idx_type j2 = std::min (std::max (zero, nr + k), nc);
97  = ((j2 - j1) * ((j1+1-k) + (j2-k))) / 2 + (nc - j2) * nr;
98  Array<T> r (dim_vector (n, 1));
99  T *rvec = r.fortran_vec ();
100  for (octave_idx_type j = 0; j < nc; j++)
101  {
102  octave_idx_type ii = std::min (std::max (zero, j + 1 - k), nr);
103  rvec = std::copy (avec, avec + ii, rvec);
104  avec += nr;
105  }
106 
107  return r;
108  }
109  else
110  {
111  NoAlias<Array<T> > r (a.dims ());
112  T *rvec = r.fortran_vec ();
113  for (octave_idx_type j = 0; j < nc; j++)
114  {
115  octave_idx_type ii = std::min (std::max (zero, j + 1 - k), nr);
116  std::copy (avec, avec + ii, rvec);
117  std::fill (rvec + ii, rvec + nr, T ());
118  avec += nr;
119  rvec += nr;
120  }
121 
122  return r;
123  }
124 }
125 
126 // These two are by David Bateman.
127 // FIXME: optimizations possible. "pack" support missing.
128 
129 template <typename T>
130 static Sparse<T>
131 do_tril (const Sparse<T>& a, octave_idx_type k, bool pack)
132 {
133  if (pack) // FIXME
134  error ("tril: \"pack\" not implemented for sparse matrices");
135 
136  Sparse<T> m = a;
137  octave_idx_type nc = m.cols ();
138 
139  for (octave_idx_type j = 0; j < nc; j++)
140  for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
141  if (m.ridx (i) < j-k)
142  m.data(i) = 0.;
143 
144  m.maybe_compress (true);
145 
146  return m;
147 }
148 
149 template <typename T>
150 static Sparse<T>
151 do_triu (const Sparse<T>& a, octave_idx_type k, bool pack)
152 {
153  if (pack) // FIXME
154  error ("triu: \"pack\" not implemented for sparse matrices");
155 
156  Sparse<T> m = a;
157  octave_idx_type nc = m.cols ();
158 
159  for (octave_idx_type j = 0; j < nc; j++)
160  for (octave_idx_type i = m.cidx (j); i < m.cidx (j+1); i++)
161  if (m.ridx (i) > j-k)
162  m.data(i) = 0.;
163 
164  m.maybe_compress (true);
165  return m;
166 }
167 
168 // Convenience dispatchers.
169 template <typename T>
170 static Array<T>
171 do_trilu (const Array<T>& a, octave_idx_type k, bool lower, bool pack)
172 {
173  return lower ? do_tril (a, k, pack) : do_triu (a, k, pack);
174 }
175 
176 template <typename T>
177 static Sparse<T>
178 do_trilu (const Sparse<T>& a, octave_idx_type k, bool lower, bool pack)
179 {
180  return lower ? do_tril (a, k, pack) : do_triu (a, k, pack);
181 }
182 
183 static octave_value
185  const octave_value_list& args)
186 {
187  bool lower = (name == "tril");
188 
189  int nargin = args.length ();
190  bool pack = false;
191 
192  if (nargin >= 2 && args(nargin-1).is_string ())
193  {
194  pack = (args(nargin-1).string_value () == "pack");
195  nargin--;
196  }
197 
198  if (nargin < 1 || nargin > 2)
199  print_usage ();
200 
201  octave_idx_type k = 0;
202  if (nargin == 2)
203  k = args(1).idx_type_value (true);
204 
205  octave_value arg = args(0);
206 
207  dim_vector dims = arg.dims ();
208  if (dims.ndims () != 2)
209  error ("%s: need a 2-D matrix", name.c_str ());
210  else if (k < -dims(0) || k > dims(1))
211  error ("%s: requested diagonal out of range", name.c_str ());
212 
214 
215  switch (arg.builtin_type ())
216  {
217  case btyp_double:
218  if (arg.is_sparse_type ())
219  retval = do_trilu (arg.sparse_matrix_value (), k, lower, pack);
220  else
221  retval = do_trilu (arg.array_value (), k, lower, pack);
222  break;
223 
224  case btyp_complex:
225  if (arg.is_sparse_type ())
226  retval = do_trilu (arg.sparse_complex_matrix_value (), k, lower,
227  pack);
228  else
229  retval = do_trilu (arg.complex_array_value (), k, lower, pack);
230  break;
231 
232  case btyp_bool:
233  if (arg.is_sparse_type ())
234  retval = do_trilu (arg.sparse_bool_matrix_value (), k, lower,
235  pack);
236  else
237  retval = do_trilu (arg.bool_array_value (), k, lower, pack);
238  break;
239 
240 #define ARRAYCASE(TYP) \
241  case btyp_ ## TYP: \
242  retval = do_trilu (arg.TYP ## _array_value (), k, lower, pack); \
243  break
244 
245  ARRAYCASE (float);
246  ARRAYCASE (float_complex);
247  ARRAYCASE (int8);
248  ARRAYCASE (int16);
249  ARRAYCASE (int32);
250  ARRAYCASE (int64);
251  ARRAYCASE (uint8);
252  ARRAYCASE (uint16);
253  ARRAYCASE (uint32);
254  ARRAYCASE (uint64);
255  ARRAYCASE (char);
256 
257 #undef ARRAYCASE
258 
259  default:
260  {
261  // Generic code that works on octave-values, that is slow
262  // but will also work on arbitrary user types
263  if (pack) // FIXME
264  error ("%s: \"pack\" not implemented for class %s",
265  name.c_str (), arg.class_name ().c_str ());
266 
267  octave_value tmp = arg;
268  if (arg.is_empty ())
269  return arg;
270 
271  octave_idx_type nr = dims(0);
272  octave_idx_type nc = dims(1);
273 
274  // The sole purpose of this code is to force the correct matrix size.
275  // This would not be necessary if the octave_value resize function
276  // allowed a fill_value. It also allows odd attributes in some user
277  // types to be handled. With a fill_value, it should be replaced with
278  //
279  // octave_value_list ov_idx;
280  // tmp = tmp.resize(dim_vector (0,0)).resize (dims, fill_value);
281 
282  octave_value_list ov_idx;
283  std::list<octave_value_list> idx_tmp;
284  ov_idx(1) = static_cast<double> (nc+1);
285  ov_idx(0) = Range (1, nr);
286  idx_tmp.push_back (ov_idx);
287  ov_idx(1) = static_cast<double> (nc);
288  tmp = tmp.resize (dim_vector (0,0));
289  tmp = tmp.subsasgn ("(",idx_tmp, arg.do_index_op (ov_idx));
290  tmp = tmp.resize (dims);
291 
292  if (lower)
293  {
294  octave_idx_type st = nc < nr + k ? nc : nr + k;
295 
296  for (octave_idx_type j = 1; j <= st; j++)
297  {
298  octave_idx_type nr_limit = 1 > j - k ? 1 : j - k;
299  ov_idx(1) = static_cast<double> (j);
300  ov_idx(0) = Range (nr_limit, nr);
301  std::list<octave_value_list> idx;
302  idx.push_back (ov_idx);
303 
304  tmp = tmp.subsasgn ("(", idx, arg.do_index_op (ov_idx));
305  }
306  }
307  else
308  {
309  octave_idx_type st = k + 1 > 1 ? k + 1 : 1;
310 
311  for (octave_idx_type j = st; j <= nc; j++)
312  {
313  octave_idx_type nr_limit = nr < j - k ? nr : j - k;
314  ov_idx(1) = static_cast<double> (j);
315  ov_idx(0) = Range (1, nr_limit);
316  std::list<octave_value_list> idx;
317  idx.push_back (ov_idx);
318 
319  tmp = tmp.subsasgn ("(", idx, arg.do_index_op (ov_idx));
320  }
321  }
322 
323  retval = tmp;
324  }
325  }
326 
327  return retval;
328 }
329 
330 DEFUN (tril, args, ,
331  doc: /* -*- texinfo -*-
332 @deftypefn {} {} tril (@var{A})
333 @deftypefnx {} {} tril (@var{A}, @var{k})
334 @deftypefnx {} {} tril (@var{A}, @var{k}, @var{pack})
335 @deftypefnx {} {} triu (@var{A})
336 @deftypefnx {} {} triu (@var{A}, @var{k})
337 @deftypefnx {} {} triu (@var{A}, @var{k}, @var{pack})
338 Return a new matrix formed by extracting the lower (@code{tril})
339 or upper (@code{triu}) triangular part of the matrix @var{A}, and
340 setting all other elements to zero.
341 
342 The second argument is optional, and specifies how many diagonals above or
343 below the main diagonal should also be set to zero.
344 
345 The default value of @var{k} is zero, so that @code{triu} and @code{tril}
346 normally include the main diagonal as part of the result.
347 
348 If the value of @var{k} is nonzero integer, the selection of elements starts
349 at an offset of @var{k} diagonals above or below the main diagonal; above
350 for positive @var{k} and below for negative @var{k}.
351 
352 The absolute value of @var{k} must not be greater than the number of
353 subdiagonals or superdiagonals.
354 
355 For example:
356 
357 @example
358 @group
359 tril (ones (3), -1)
360  @result{} 0 0 0
361  1 0 0
362  1 1 0
363 @end group
364 @end example
365 
366 @noindent
367 and
368 
369 @example
370 @group
371 tril (ones (3), 1)
372  @result{} 1 1 0
373  1 1 1
374  1 1 1
375 @end group
376 @end example
377 
378 If the option @qcode{"pack"} is given as third argument, the extracted
379 elements are not inserted into a matrix, but rather stacked column-wise one
380 above other.
381 @seealso{diag}
382 @end deftypefn */)
383 {
384  return do_trilu ("tril", args);
385 }
386 
387 DEFUN (triu, args, ,
388  doc: /* -*- texinfo -*-
389 @deftypefn {} {} triu (@var{A})
390 @deftypefnx {} {} triu (@var{A}, @var{k})
391 @deftypefnx {} {} triu (@var{A}, @var{k}, @var{pack})
392 See the documentation for the @code{tril} function (@pxref{tril}).
393 @seealso{tril}
394 @end deftypefn */)
395 {
396  return do_trilu ("triu", args);
397 }
398 
399 /*
400 %!test
401 %! a = [1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12];
402 %!
403 %! l0 = [1, 0, 0; 4, 5, 0; 7, 8, 9; 10, 11, 12];
404 %! l1 = [1, 2, 0; 4, 5, 6; 7, 8, 9; 10, 11, 12];
405 %! l2 = [1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12];
406 %! lm1 = [0, 0, 0; 4, 0, 0; 7, 8, 0; 10, 11, 12];
407 %! lm2 = [0, 0, 0; 0, 0, 0; 7, 0, 0; 10, 11, 0];
408 %! lm3 = [0, 0, 0; 0, 0, 0; 0, 0, 0; 10, 0, 0];
409 %! lm4 = [0, 0, 0; 0, 0, 0; 0, 0, 0; 0, 0, 0];
410 %!
411 %! assert (tril (a, -4), lm4);
412 %! assert (tril (a, -3), lm3);
413 %! assert (tril (a, -2), lm2);
414 %! assert (tril (a, -1), lm1);
415 %! assert (tril (a), l0);
416 %! assert (tril (a, 1), l1);
417 %! assert (tril (a, 2), l2);
418 
419 %!error tril ()
420 */
ComplexNDArray complex_array_value(bool frc_str_conv=false) const
Definition: ov.h:812
octave_idx_type cols(void) const
Definition: Sparse.h:272
T * data(void)
Definition: Sparse.h:521
OCTINTERP_API void print_usage(void)
Definition: defun.cc:52
OCTAVE_EXPORT octave_value_list uint16
Definition: ov.cc:1258
octave_idx_type length(void) const
Definition: ovl.h:96
for large enough k
Definition: lu.cc:606
Definition: Range.h:33
#define DEFUN(name, args_name, nargout_name, doc)
Definition: defun.h:46
void error(const char *fmt,...)
Definition: error.cc:570
is greater than zero
Definition: load-path.cc:2339
octave_idx_type * cidx(void)
Definition: Sparse.h:543
static Array< T > do_tril(const Array< T > &a, octave_idx_type k, bool pack)
Definition: tril.cc:43
boolNDArray bool_array_value(bool warn=false) const
Definition: ov.h:825
octave_idx_type rows(void) const
Definition: Array.h:401
octave_value arg
Definition: pr-output.cc:3440
static Array< T > do_triu(const Array< T > &a, octave_idx_type k, bool pack)
Definition: tril.cc:85
calling an anonymous function involves an overhead quite comparable to the overhead of an m file function Passing a handle to a built in function is because the interpreter is not involved in the internal loop For a
Definition: cellfun.cc:398
JNIEnv void * args
Definition: ov-java.cc:67
const dim_vector & dims(void) const
Return a const-reference so that dims ()(i) works efficiently.
Definition: Array.h:439
#define ARRAYCASE(TYP)
OCTAVE_EXPORT octave_value_list any number nd example oindent prints the prompt xample Pick a any number!nd example oindent and waits for the user to enter a value The string entered by the user is evaluated as an so it may be a literal a variable name
Definition: input.cc:871
bool is_sparse_type(void) const
Definition: ov.h:682
OCTAVE_EXPORT octave_value_list uint32
Definition: ov.cc:1258
nd deftypefn *octave_map m
Definition: ov-struct.cc:2058
Sparse< T > maybe_compress(bool remove_zeros=false)
Definition: Sparse.h:481
int nargin
Definition: graphics.cc:10115
OCTAVE_EXPORT octave_value_list int16
Definition: ov.cc:1302
SparseBoolMatrix sparse_bool_matrix_value(bool warn=false) const
Definition: ov.h:841
double tmp
Definition: data.cc:6300
octave_value retval
Definition: data.cc:6294
OCTAVE_EXPORT octave_value_list int32
Definition: ov.cc:1258
SparseComplexMatrix sparse_complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:838
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
dim_vector dims(void) const
Definition: ov.h:486
OCTAVE_EXPORT octave_value_list int64
Definition: ov.cc:1258
N Dimensional Array with copy-on-write semantics.
Definition: Array.h:126
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:228
bool is_empty(void) const
Definition: ov.h:542
This is a simple wrapper template that will subclass an Array type or any later type derived from ...
Definition: Array.h:888
NDArray array_value(bool frc_str_conv=false) const
Definition: ov.h:793
octave_idx_type * ridx(void)
Definition: Sparse.h:530
=val(i)}if ode{val(i)}occurs in table i
Definition: lookup.cc:239
octave_idx_type ndims(void) const
Number of dimensions.
Definition: dim-vector.h:301
std::string class_name(void) const
Definition: ov.h:1234
void resize(octave_idx_type n, const octave_value &rfv=octave_value())
Definition: ovl.h:100
SparseMatrix sparse_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:834
const T * fortran_vec(void) const
Definition: Array.h:584
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
builtin_type_t builtin_type(void) const
Definition: ov.h:619
octave_idx_type columns(void) const
Definition: Array.h:410
static Array< T > do_trilu(const Array< T > &a, octave_idx_type k, bool lower, bool pack)
Definition: tril.cc:171
charNDArray min(char d, const charNDArray &m)
Definition: chNDArray.cc:205
octave_value do_index_op(const octave_value_list &idx, bool resize_ok=false)
Definition: ov.h:454