GNU Octave  4.4.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
tril.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2004-2018 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
10 the Free Software Foundation, either version 3 of the License, or
11 (at your option) any later version.
12 
13 Octave is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License 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 <https://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 (R"(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 (R"(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 
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.issparse ())
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.issparse ())
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.issparse ())
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 (R"(%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.isempty ())
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 {} {@var{A_LO} =} tril (@var{A})
333 @deftypefnx {} {@var{A_LO} =} tril (@var{A}, @var{k})
334 @deftypefnx {} {@var{A_LO} =} tril (@var{A}, @var{k}, @var{pack})
335 Return a new matrix formed by extracting the lower triangular part of the
336 matrix @var{A}, and setting all other elements to zero.
337 
338 The optional second argument specifies how many diagonals above or below the
339 main diagonal should also be set to zero. The default value of @var{k} is
340 zero which includes the main diagonal as part of the result. If the value of
341 @var{k} is a nonzero integer then the selection of elements starts at an offset
342 of @var{k} diagonals above the main diagonal for positive @var{k} or below the
343 main diagonal for negative @var{k}. The absolute value of @var{k} may not be
344 greater than the number of subdiagonals or superdiagonals.
345 
346 Example 1 : exclude main diagonal
347 
348 @example
349 @group
350 tril (ones (3), -1)
351  @result{} 0 0 0
352  1 0 0
353  1 1 0
354 @end group
355 @end example
356 
357 @noindent
358 
359 Example 2 : include first superdiagonal
360 
361 @example
362 @group
363 tril (ones (3), 1)
364  @result{} 1 1 0
365  1 1 1
366  1 1 1
367 @end group
368 @end example
369 
370 If the optional third argument @qcode{"pack"} is given then the extracted
371 elements are not inserted into a matrix, but instead stacked column-wise one
372 above another, and returned as a column vector.
373 @seealso{triu, istril, diag}
374 @end deftypefn */)
375 {
376  return do_trilu ("tril", args);
377 }
378 
379 DEFUN (triu, args, ,
380  doc: /* -*- texinfo -*-
381 @deftypefn {} {@var{A_UP} =} triu (@var{A})
382 @deftypefnx {} {@var{A_UP} =} triu (@var{A}, @var{k})
383 @deftypefnx {} {@var{A_UP} =} triu (@var{A}, @var{k}, @var{pack})
384 Return a new matrix formed by extracting the upper triangular part of the
385 matrix @var{A}, and setting all other elements to zero.
386 
387 The optional second argument specifies how many diagonals above or below the
388 main diagonal should also be set to zero. The default value of @var{k} is
389 zero which includes the main diagonal as part of the result. If the value of
390 @var{k} is a nonzero integer then the selection of elements starts at an offset
391 of @var{k} diagonals above the main diagonal for positive @var{k} or below the
392 main diagonal for negative @var{k}. The absolute value of @var{k} may not be
393 greater than the number of subdiagonals or superdiagonals.
394 
395 Example 1 : exclude main diagonal
396 
397 @example
398 @group
399 triu (ones (3), 1)
400  @result{} 0 1 1
401  0 0 1
402  0 0 0
403 @end group
404 @end example
405 
406 @noindent
407 
408 Example 2 : include first subdiagonal
409 
410 @example
411 @group
412 triu (ones (3), -1)
413  @result{} 1 1 1
414  1 1 1
415  0 1 1
416 @end group
417 @end example
418 
419 If the optional third argument @qcode{"pack"} is given then the extracted
420 elements are not inserted into a matrix, but instead stacked column-wise one
421 above another, and returned as a column vector.
422 @seealso{tril, istriu, diag}
423 @end deftypefn */)
424 {
425  return do_trilu ("triu", args);
426 }
427 
428 /*
429 %!test
430 %! a = [1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12];
431 %!
432 %! l0 = [1, 0, 0; 4, 5, 0; 7, 8, 9; 10, 11, 12];
433 %! l1 = [1, 2, 0; 4, 5, 6; 7, 8, 9; 10, 11, 12];
434 %! l2 = [1, 2, 3; 4, 5, 6; 7, 8, 9; 10, 11, 12];
435 %! lm1 = [0, 0, 0; 4, 0, 0; 7, 8, 0; 10, 11, 12];
436 %! lm2 = [0, 0, 0; 0, 0, 0; 7, 0, 0; 10, 11, 0];
437 %! lm3 = [0, 0, 0; 0, 0, 0; 0, 0, 0; 10, 0, 0];
438 %! lm4 = [0, 0, 0; 0, 0, 0; 0, 0, 0; 0, 0, 0];
439 %!
440 %! assert (tril (a, -4), lm4);
441 %! assert (tril (a, -3), lm3);
442 %! assert (tril (a, -2), lm2);
443 %! assert (tril (a, -1), lm1);
444 %! assert (tril (a), l0);
445 %! assert (tril (a, 1), l1);
446 %! assert (tril (a, 2), l2);
447 
448 %!error tril ()
449 */
T * data(void)
Definition: Sparse.h:486
bool isempty(void) const
Definition: ov.h:529
OCTINTERP_API void print_usage(void)
Definition: defun.cc:54
OCTAVE_EXPORT octave_value_list or N dimensional array whose elements are all equal to the IEEE symbol zero divided by zero($0/0$)
octave_idx_type * cidx(void)
Definition: Sparse.h:508
for large enough k
Definition: lu.cc:617
const T * fortran_vec(void) const
Definition: Array.h:584
Definition: Range.h:33
#define DEFUN(name, args_name, nargout_name, doc)
Macro to define a builtin function.
Definition: defun.h:53
void error(const char *fmt,...)
Definition: error.cc:578
static Array< T > do_tril(const Array< T > &a, octave_idx_type k, bool pack)
Definition: tril.cc:43
octave_value arg
Definition: pr-output.cc:3244
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:400
#define ARRAYCASE(TYP)
nd deftypefn *std::string name
Definition: sysdep.cc:647
dim_vector dims(void) const
Definition: ov.h:469
double tmp
Definition: data.cc:6252
bool issparse(void) const
Definition: ov.h:730
octave_value retval
Definition: data.cc:6246
octave_idx_type * ridx(void)
Definition: Sparse.h:495
std::string class_name(void) const
Definition: ov.h:1291
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
SparseComplexMatrix sparse_complex_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:885
octave_idx_type cols(void) const
Definition: Sparse.h:259
N Dimensional Array with copy-on-write semantics.
Definition: Array.h:125
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:227
builtin_type_t builtin_type(void) const
Definition: ov.h:643
boolNDArray bool_array_value(bool warn=false) const
Definition: ov.h:872
This is a simple wrapper template that will subclass an Array<T> type or any later type derived from ...
Definition: Array.h:892
octave_idx_type length(void) const
Definition: ovl.h:96
SparseMatrix sparse_matrix_value(bool frc_str_conv=false) const
Definition: ov.h:881
Sparse< T > maybe_compress(bool remove_zeros=false)
Definition: Sparse.h:438
args.length() nargin
Definition: file-io.cc:589
for i
Definition: data.cc:5264
SparseBoolMatrix sparse_bool_matrix_value(bool warn=false) const
Definition: ov.h:888
ComplexNDArray complex_array_value(bool frc_str_conv=false) const
Definition: ov.h:859
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:888
NDArray array_value(bool frc_str_conv=false) const
Definition: ov.h:840
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:204
octave_value do_index_op(const octave_value_list &idx, bool resize_ok=false)
Definition: ov.h:444