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
Sparse-perm-op-defs.h
Go to the documentation of this file.
1 /* -*- C++ -*-
2 
3 Copyright (C) 2009-2015 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 // Matrix multiplication
27 
28 template <typename SM>
29 SM octinternal_do_mul_colpm_sm (const octave_idx_type *pcol, const SM& a)
30 // Relabel the rows according to pcol.
31 {
32  const octave_idx_type nr = a.rows ();
33  const octave_idx_type nc = a.cols ();
34  const octave_idx_type nent = a.nnz ();
35  SM r (nr, nc, nent);
36 
38 
39  for (octave_idx_type j = 0; j <= nc; ++j)
40  r.xcidx (j) = a.cidx (j);
41 
42  for (octave_idx_type j = 0; j < nc; j++)
43  {
44  octave_quit ();
45 
46  OCTAVE_LOCAL_BUFFER (octave_idx_type, sidx, r.xcidx (j+1) - r.xcidx (j));
47  for (octave_idx_type i = r.xcidx (j), ii = 0; i < r.xcidx (j+1); i++)
48  {
49  sidx[ii++]=i;
50  r.xridx (i) = pcol[a.ridx (i)];
51  }
52  sort.sort (r.xridx () + r.xcidx (j), sidx, r.xcidx (j+1) - r.xcidx (j));
53  for (octave_idx_type i = r.xcidx (j), ii = 0; i < r.xcidx (j+1); i++)
54  r.xdata (i) = a.data (sidx[ii++]);
55  }
56 
57  return r;
58 }
59 
60 template <typename SM>
61 SM octinternal_do_mul_pm_sm (const PermMatrix& p, const SM& a)
62 {
63  const octave_idx_type nr = a.rows ();
64  if (p.cols () != nr)
65  {
66  gripe_nonconformant ("operator *", p.rows (), p.cols (), a.rows (), a.cols ());
67  return SM ();
68  }
69 
70  return octinternal_do_mul_colpm_sm (p.col_perm_vec ().data (), a);
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  {
152  gripe_nonconformant ("operator *", a.rows (), a.cols (), p.rows (), p.cols ());
153  return SM ();
154  }
155 
156  return octinternal_do_mul_sm_colpm (a, p.col_perm_vec ().data ());
157 }
158 
159 #endif // octave_Sparse_perm_op_defs_h
void gripe_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
octave_idx_type rows(void) const
Definition: PermMatrix.h:55
SM octinternal_do_mul_colpm_sm(const octave_idx_type *pcol, const SM &a)
octave_idx_type cols(void) const
Definition: PermMatrix.h:56
SM octinternal_do_mul_pm_sm(const PermMatrix &p, const SM &a)
octave_idx_type * xcidx(void)
Definition: Sparse.h:537
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)
const Array< octave_idx_type > & col_perm_vec(void) const
Definition: PermMatrix.h:72
Sparse< bool > sort(octave_idx_type dim=0, sortmode mode=ASCENDING) const
const T * data(void) const
Definition: Array.h:479
void sort(T *data, octave_idx_type nel)
Definition: oct-sort.cc:1514
#define OCTAVE_LOCAL_BUFFER(T, buf, size)
Definition: oct-locbuf.h:197