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
PermMatrix.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2008-2015 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 #ifdef HAVE_CONFIG_H
24 #include <config.h>
25 #endif
26 
27 #include "PermMatrix.h"
28 #include "idx-vector.h"
29 #include "Array-util.h"
30 #include "oct-locbuf.h"
31 
32 static void
34 {
35  (*current_liboctave_error_handler)
36  ("PermMatrix: invalid permutation vector");
37 }
38 
39 void
40 PermMatrix::setup (const Array<octave_idx_type>& p, bool colp, bool check)
41 {
42  if (check)
43  {
44  if (! idx_vector (p).is_permutation (p.length ()))
45  {
48  }
49  }
50 
51  if (!colp)
52  *this = this->transpose ();
53 }
54 
56  : Array<octave_idx_type> (p)
57 {
58  setup (p, false, true);
59 }
60 
61 PermMatrix::PermMatrix (const Array<octave_idx_type>& p, bool colp, bool check)
62  : Array<octave_idx_type> (p)
63 {
64  setup (p, colp, check);
65 }
66 
67 void
68 PermMatrix::setup (const idx_vector& idx, bool colp, octave_idx_type n)
69 {
70  octave_idx_type len = idx.length (n);
71 
72  if (! idx.is_permutation (len))
74  else
75  {
76  Array<octave_idx_type> idxa (dim_vector (len, 1));
77  for (octave_idx_type i = 0; i < len; i++) idxa(i) = idx(i);
79  }
80 
81  if (!colp)
82  *this = this->transpose ();
83 }
84 
87 {
88  setup (idx, false, 0);
89 }
90 
93 {
94  setup (idx, colp, n);
95 }
96 
98  : Array<octave_idx_type> (dim_vector (n, 1))
99 {
100  for (octave_idx_type i = 0; i < n; i++)
101  xelem (i) = i;
102 }
103 
106 {
108  if (i < 0 || j < 0 || i > len || j > len)
109  {
110  (*current_liboctave_error_handler) ("index out of range");
111  return 0;
112  }
113  else
114  return elem (i, j);
115 }
116 
117 
120 {
122 
123  PermMatrix retval (len);
124 
125  for (octave_idx_type i = 0; i < len; ++i)
126  retval.xelem (xelem (i)) = i;
127 
128  return retval;
129 }
130 
133 {
134  return transpose ();
135 }
136 
139 {
140  // Determine the sign of a permutation in linear time.
141  // Is this widely known?
142 
143  octave_idx_type len = perm_length ();
144  const octave_idx_type *pa = data ();
145 
148 
149  for (octave_idx_type i = 0; i < len; i++)
150  {
151  p[i] = pa[i];
152  q[p[i]] = i;
153  }
154 
155  bool neg = false;
156 
157  for (octave_idx_type i = 0; i < len; i++)
158  {
159  octave_idx_type j = p[i];
160  octave_idx_type k = q[i];
161  if (j != i)
162  {
163  p[k] = p[i];
164  q[j] = q[i];
165  neg = ! neg;
166  }
167  }
168 
169  return neg ? -1 : 1;
170 }
171 
174 {
175  if (m < 0)
176  return inverse ().pos_power (-m);
177  else if (m > 0)
178  return pos_power (m);
179  else
180  return PermMatrix (rows ());
181 }
182 
185 {
186  octave_idx_type n = rows ();
187 
188  const octave_idx_type *p = data ();
189  Array<octave_idx_type> res_pvec (dim_vector (n, 1), -1);
190  octave_idx_type *q = res_pvec.fortran_vec ();
191 
192  for (octave_idx_type ics = 0; ics < n; ics++)
193  {
194  if (q[ics] > 0)
195  continue;
196 
197  // go forward m steps or until a cycle is found.
198  octave_idx_type ic, j;
199  for (j = 1, ic = p[ics]; j != m && ic != ics; j++, ic = p[ic]) ;
200  if (ic == ics)
201  {
202  // reduce power.
203  octave_idx_type mm = m % j;
204  // go forward mm steps.
205  for (j = 0, ic = ics; j != mm; j++, ic = p[ic]) ;
206  }
207 
208  // now ic = p^m[ics]. Loop through the whole cycle.
209  octave_idx_type jcs = ics;
210  do
211  {
212  q[jcs] = ic;
213  jcs = p[jcs]; ic = p[ic];
214  }
215  while (jcs != ics);
216 
217  }
218 
219  return PermMatrix (res_pvec, true, false);
220 }
221 
224 {
225  return PermMatrix (n);
226 }
227 
229 operator *(const PermMatrix& a, const PermMatrix& b)
230 {
231  PermMatrix r;
232 
233  const Array<octave_idx_type> ia = a.col_perm_vec ();
234  const Array<octave_idx_type> ib = b.col_perm_vec ();
235 
236  octave_idx_type n = a.columns ();
237 
238  if (n != b.rows ())
239  gripe_nonconformant ("operator *", n, n, b.rows (), b.rows ());
240  else
241  r = PermMatrix (ia.index (idx_vector (ib)), true, false);
242 
243  return r;
244 }
void gripe_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
octave_idx_type length(octave_idx_type n=0) const
Definition: idx-vector.h:551
octave_idx_type rows(void) const
Definition: PermMatrix.h:55
static void gripe_invalid_permutation(void)
Definition: PermMatrix.cc:33
octave_idx_type perm_length(void) const
Definition: PermMatrix.h:59
PermMatrix(void)
Definition: PermMatrix.h:36
PermMatrix transpose(void) const
Definition: PermMatrix.cc:119
octave_idx_type checkelem(octave_idx_type i, octave_idx_type j) const
Definition: PermMatrix.cc:105
octave_idx_type columns(void) const
Definition: PermMatrix.h:57
octave_idx_type determinant(void) const
Definition: PermMatrix.cc:138
PermMatrix operator*(const PermMatrix &a, const PermMatrix &b)
Definition: PermMatrix.cc:229
PermMatrix inverse(void) const
Definition: PermMatrix.cc:132
const Array< octave_idx_type > & col_perm_vec(void) const
Definition: PermMatrix.h:72
const octave_idx_type * data(void) const
Definition: Array.h:479
octave_idx_type elem(octave_idx_type i, octave_idx_type j) const
Definition: PermMatrix.h:76
bool is_permutation(octave_idx_type n) const
Definition: idx-vector.cc:1145
PermMatrix pos_power(octave_idx_type m) const
Definition: PermMatrix.cc:184
octave_idx_type & xelem(octave_idx_type n)
Definition: Array.h:353
octave_idx_type length(void) const
Number of elements in the array.
Definition: Array.h:267
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:197
const T * fortran_vec(void) const
Definition: Array.h:481
static PermMatrix eye(octave_idx_type n)
Definition: PermMatrix.cc:223
PermMatrix power(octave_idx_type n) const
Definition: PermMatrix.cc:173
Array< T > & operator=(const Array< T > &a)
Definition: Array.h:226
void setup(const Array< octave_idx_type > &p, bool colp, bool check)
Definition: PermMatrix.cc:40
Array< T > index(const idx_vector &i) const
Indexing without resizing.
Definition: Array.cc:716