GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
PermMatrix.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (C) 2008-2021 The Octave Project Developers
4 //
5 // See the file COPYRIGHT.md in the top-level directory of this
6 // distribution or <https://octave.org/copyright/>.
7 //
8 // This file is part of Octave.
9 //
10 // Octave is free software: you can redistribute it and/or modify it
11 // under the terms of the GNU General Public License as published by
12 // the Free Software Foundation, either version 3 of the License, or
13 // (at your option) any later version.
14 //
15 // Octave is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 // GNU General Public License for more details.
19 //
20 // You should have received a copy of the GNU General Public License
21 // along with Octave; see the file COPYING. If not, see
22 // <https://www.gnu.org/licenses/>.
23 //
24 ////////////////////////////////////////////////////////////////////////
25 
26 #if defined (HAVE_CONFIG_H)
27 # include "config.h"
28 #endif
29 
30 #include "PermMatrix.h"
31 #include "idx-vector.h"
32 #include "Array-util.h"
33 #include "oct-locbuf.h"
34 
35 OCTAVE_NORETURN static
36 void
38 {
39  (*current_liboctave_error_handler) ("PermMatrix: invalid permutation vector");
40 }
41 
42 void
43 PermMatrix::setup (const Array<octave_idx_type>& p, bool colp, bool check)
44 {
45  if (check)
46  {
47  if (! idx_vector (p).is_permutation (p.numel ()))
49  }
50 
51  if (! colp)
52  *this = this->transpose ();
53 }
54 
55 PermMatrix::PermMatrix (const Array<octave_idx_type>& p, bool colp, bool check)
56  : Array<octave_idx_type> (p)
57 {
58  setup (p, colp, check);
59 }
60 
61 void
63 {
64  octave_idx_type len = idx.length (n);
65 
66  if (! idx.is_permutation (len))
68 
70  for (octave_idx_type i = 0; i < len; i++) idxa(i) = idx(i);
72 
73  if (! colp)
74  *this = this->transpose ();
75 }
76 
79 {
80  setup (idx, colp, n);
81 }
82 
85 {
86  for (octave_idx_type i = 0; i < n; i++)
87  xelem (i) = i;
88 }
89 
92 {
94  if (i < 0 || j < 0 || i > len || j > len)
95  (*current_liboctave_error_handler) ("index out of range");
96 
97  return elem (i, j);
98 }
99 
102 {
104 
106 
107  for (octave_idx_type i = 0; i < len; ++i)
108  retval.xelem (xelem (i)) = i;
109 
110  return retval;
111 }
112 
115 {
116  return transpose ();
117 }
118 
121 {
122  // Determine the sign of a permutation in linear time.
123  // Is this widely known?
124 
126  const octave_idx_type *pa = data ();
127 
130 
131  for (octave_idx_type i = 0; i < len; i++)
132  {
133  p[i] = pa[i];
134  q[p[i]] = i;
135  }
136 
137  bool neg = false;
138 
139  for (octave_idx_type i = 0; i < len; i++)
140  {
141  octave_idx_type j = p[i];
142  octave_idx_type k = q[i];
143  if (j != i)
144  {
145  p[k] = p[i];
146  q[j] = q[i];
147  neg = ! neg;
148  }
149  }
150 
151  return neg ? -1 : 1;
152 }
153 
156 {
157  if (m < 0)
158  return inverse ().pos_power (-m);
159  else if (m > 0)
160  return pos_power (m);
161  else
162  return PermMatrix (rows ());
163 }
164 
167 {
168  octave_idx_type n = rows ();
169 
170  const octave_idx_type *p = data ();
171  Array<octave_idx_type> res_pvec (dim_vector (n, 1), -1);
172  octave_idx_type *q = res_pvec.fortran_vec ();
173 
174  for (octave_idx_type ics = 0; ics < n; ics++)
175  {
176  if (q[ics] > 0)
177  continue;
178 
179  // go forward m steps or until a cycle is found.
180  octave_idx_type ic, j;
181  for (j = 1, ic = p[ics]; j != m && ic != ics; j++, ic = p[ic]) ;
182  if (ic == ics)
183  {
184  // reduce power.
185  octave_idx_type mm = m % j;
186  // go forward mm steps.
187  for (j = 0, ic = ics; j != mm; j++, ic = p[ic]) ;
188  }
189 
190  // now ic = p^m[ics]. Loop through the whole cycle.
191  octave_idx_type jcs = ics;
192  do
193  {
194  q[jcs] = ic;
195  jcs = p[jcs]; ic = p[ic];
196  }
197  while (jcs != ics);
198 
199  }
200 
201  return PermMatrix (res_pvec, true, false);
202 }
203 
206 {
207  return PermMatrix (n);
208 }
209 
211 operator *(const PermMatrix& a, const PermMatrix& b)
212 {
213  PermMatrix r;
214 
215  const Array<octave_idx_type> ia = a.col_perm_vec ();
216  const Array<octave_idx_type> ib = b.col_perm_vec ();
217 
218  octave_idx_type n = a.columns ();
219 
220  if (n != b.rows ())
221  octave::err_nonconformant ("operator *", n, n, b.rows (), b.rows ());
222 
223  r = PermMatrix (ia.index (idx_vector (ib)), true, false);
224 
225  return r;
226 }
static OCTAVE_NORETURN void err_invalid_permutation(void)
Definition: PermMatrix.cc:37
PermMatrix operator*(const PermMatrix &a, const PermMatrix &b)
Definition: PermMatrix.cc:211
octave_idx_type & xelem(octave_idx_type n)
Definition: Array.h:469
octave_idx_type numel(void) const
Number of elements in the array.
Definition: Array.h:377
Array< T > & operator=(const Array< T > &a)
Definition: Array.h:325
const octave_idx_type * data(void) const
Definition: Array.h:581
Array< T > index(const idx_vector &i) const
Indexing without resizing.
Definition: Array.cc:698
const T * fortran_vec(void) const
Size of the specified dimension.
Definition: Array.h:583
octave_idx_type perm_length(void) const
Definition: PermMatrix.h:64
const Array< octave_idx_type > & col_perm_vec(void) const
Definition: PermMatrix.h:81
octave_idx_type columns(void) const
Definition: PermMatrix.h:62
PermMatrix inverse(void) const
Definition: PermMatrix.cc:114
octave_idx_type elem(octave_idx_type i, octave_idx_type j) const
Definition: PermMatrix.h:85
PermMatrix pos_power(octave_idx_type m) const
Definition: PermMatrix.cc:166
octave_idx_type rows(void) const
Definition: PermMatrix.h:60
octave_idx_type determinant(void) const
Definition: PermMatrix.cc:120
void setup(const Array< octave_idx_type > &p, bool colp, bool check)
Definition: PermMatrix.cc:43
octave_idx_type checkelem(octave_idx_type i, octave_idx_type j) const
Definition: PermMatrix.cc:91
PermMatrix(void)=default
static PermMatrix eye(octave_idx_type n)
Definition: PermMatrix.cc:205
PermMatrix power(octave_idx_type n) const
Definition: PermMatrix.cc:155
PermMatrix transpose(void) const
Definition: PermMatrix.cc:101
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:95
octave_idx_type length(octave_idx_type n=0) const
Definition: idx-vector.h:558
bool is_permutation(octave_idx_type n) const
Definition: idx-vector.cc:1141
T octave_idx_type m
Definition: mx-inlines.cc:773
octave_idx_type n
Definition: mx-inlines.cc:753
T * r
Definition: mx-inlines.cc:773
void err_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:44
octave_value::octave_value(const Array< char > &chm, char type) return retval
Definition: ov.cc:811
F77_RET_T len
Definition: xerbla.cc:61