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
Sparse-perm-op-defs.h
Go to the documentation of this file.
1 /* -*- C++ -*-
2 
3 Copyright (C) 2009-2017 Jason Riedy
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 #if ! defined (octave_Sparse_perm_op_defs_h)
24 #define octave_Sparse_perm_op_defs_h 1
25 
26 #include "octave-config.h"
27 
28 // Matrix multiplication
29 
30 template <typename SM>
31 SM octinternal_do_mul_colpm_sm (const octave_idx_type *pcol, const SM& a)
32 // Relabel the rows according to pcol.
33 {
34  const octave_idx_type nr = a.rows ();
35  const octave_idx_type nc = a.cols ();
36  const octave_idx_type nent = a.nnz ();
37  SM r (nr, nc, nent);
38 
40 
41  for (octave_idx_type j = 0; j <= nc; ++j)
42  r.xcidx (j) = a.cidx (j);
43 
44  for (octave_idx_type j = 0; j < nc; j++)
45  {
46  octave_quit ();
47 
48  OCTAVE_LOCAL_BUFFER (octave_idx_type, sidx, r.xcidx (j+1) - r.xcidx (j));
49  for (octave_idx_type i = r.xcidx (j), ii = 0; i < r.xcidx (j+1); i++)
50  {
51  sidx[ii++]=i;
52  r.xridx (i) = pcol[a.ridx (i)];
53  }
54  sort.sort (r.xridx () + r.xcidx (j), sidx, r.xcidx (j+1) - r.xcidx (j));
55  for (octave_idx_type i = r.xcidx (j), ii = 0; i < r.xcidx (j+1); i++)
56  r.xdata (i) = a.data (sidx[ii++]);
57  }
58 
59  return r;
60 }
61 
62 template <typename SM>
63 SM octinternal_do_mul_pm_sm (const PermMatrix& p, const SM& a)
64 {
65  const octave_idx_type nr = a.rows ();
66  if (p.cols () != nr)
67  octave::err_nonconformant ("operator *",
68  p.rows (), p.cols (), a.rows (), a.cols ());
69 
71 }
72 
73 template <typename SM>
74 SM octinternal_do_mul_sm_rowpm (const SM& a, const octave_idx_type *prow)
75 // For a row permutation, iterate across the source a and stuff the
76 // results into the correct destination column in r.
77 {
78  const octave_idx_type nr = a.rows ();
79  const octave_idx_type nc = a.cols ();
80  const octave_idx_type nent = a.nnz ();
81  SM r (nr, nc, nent);
82 
83  for (octave_idx_type j_src = 0; j_src < nc; ++j_src)
84  r.xcidx (prow[j_src]) = a.cidx (j_src+1) - a.cidx (j_src);
85  octave_idx_type k = 0;
86  for (octave_idx_type j = 0; j < nc; ++j)
87  {
88  const octave_idx_type tmp = r.xcidx (j);
89  r.xcidx (j) = k;
90  k += tmp;
91  }
92  r.xcidx (nc) = nent;
93 
94  octave_idx_type k_src = 0;
95  for (octave_idx_type j_src = 0; j_src < nc; ++j_src)
96  {
97  octave_quit ();
98  const octave_idx_type j = prow[j_src];
99  const octave_idx_type kend_src = a.cidx (j_src + 1);
100  for (k = r.xcidx (j); k_src < kend_src; ++k, ++k_src)
101  {
102  r.xridx (k) = a.ridx (k_src);
103  r.xdata (k) = a.data (k_src);
104  }
105  }
106  assert (k_src == nent);
107 
108  return r;
109 }
110 
111 template <typename SM>
112 SM octinternal_do_mul_sm_colpm (const SM& a, const octave_idx_type *pcol)
113 // For a column permutation, iterate across the destination r and pull
114 // data from the correct column of a.
115 {
116  const octave_idx_type nr = a.rows ();
117  const octave_idx_type nc = a.cols ();
118  const octave_idx_type nent = a.nnz ();
119  SM r (nr, nc, nent);
120 
121  for (octave_idx_type j = 0; j < nc; ++j)
122  {
123  const octave_idx_type j_src = pcol[j];
124  r.xcidx (j+1) = r.xcidx (j) + (a.cidx (j_src+1) - a.cidx (j_src));
125  }
126  assert (r.xcidx (nc) == nent);
127 
128  octave_idx_type k = 0;
129  for (octave_idx_type j = 0; j < nc; ++j)
130  {
131  octave_quit ();
132  const octave_idx_type j_src = pcol[j];
133  octave_idx_type k_src;
134  const octave_idx_type kend_src = a.cidx (j_src + 1);
135  for (k_src = a.cidx (j_src); k_src < kend_src; ++k_src, ++k)
136  {
137  r.xridx (k) = a.ridx (k_src);
138  r.xdata (k) = a.data (k_src);
139  }
140  }
141  assert (k == nent);
142 
143  return r;
144 }
145 
146 template <typename SM>
147 SM octinternal_do_mul_sm_pm (const SM& a, const PermMatrix& p)
148 {
149  const octave_idx_type nc = a.cols ();
150  if (p.rows () != nc)
151  octave::err_nonconformant ("operator *",
152  a.rows (), a.cols (), p.rows (), p.cols ());
153 
154  return octinternal_do_mul_sm_colpm (a, p.col_perm_vec ().data ());
155 }
156 
157 #endif
octave_idx_type rows(void) const
Definition: PermMatrix.h:59
SM octinternal_do_mul_colpm_sm(const octave_idx_type *pcol, const SM &a)
octave_idx_type cols(void) const
Definition: PermMatrix.h:60
SM octinternal_do_mul_pm_sm(const PermMatrix &p, const SM &a)
for large enough k
Definition: lu.cc:606
octave_idx_type * xcidx(void)
Definition: Sparse.h:549
SM octinternal_do_mul_sm_pm(const SM &a, const PermMatrix &p)
SM octinternal_do_mul_sm_colpm(const SM &a, const octave_idx_type *pcol)
SM octinternal_do_mul_sm_rowpm(const SM &a, const octave_idx_type *prow)
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
const Array< octave_idx_type > & col_perm_vec(void) const
Definition: PermMatrix.h:79
Sparse< bool > sort(octave_idx_type dim=0, sortmode mode=ASCENDING) const
const T * data(void) const
Definition: Array.h:582
double tmp
Definition: data.cc:6300
void err_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
void sort(T *data, octave_idx_type nel)
Definition: oct-sort.cc:1506
=val(i)}if ode{val(i)}occurs in table i
Definition: lookup.cc:239
p
Definition: lu.cc:138
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:200