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
MSparse.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2004-2015 David Bateman
4 Copyright (C) 1998-2004 Andy Adler
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 the
10 Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
12 
13 Octave is distributed in the hope that it will be useful, but WITHOUT
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 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 <http://www.gnu.org/licenses/>.
21 
22 */
23 
24 #ifdef HAVE_CONFIG_H
25 #include <config.h>
26 #endif
27 
28 #include <functional>
29 
30 #include "quit.h"
31 #include "lo-error.h"
32 #include "MArray.h"
33 #include "Array-util.h"
34 
35 #include "MSparse.h"
36 #include "MSparse-defs.h"
37 
38 // sparse array with math ops.
39 
40 // Element by element MSparse by MSparse ops.
41 
42 template <class T, class OP>
44 plus_or_minus (MSparse<T>& a, const MSparse<T>& b, OP op, const char* op_name)
45 {
46  MSparse<T> r;
47 
48  octave_idx_type a_nr = a.rows ();
49  octave_idx_type a_nc = a.cols ();
50 
51  octave_idx_type b_nr = b.rows ();
52  octave_idx_type b_nc = b.cols ();
53 
54  if (a_nr != b_nr || a_nc != b_nc)
55  gripe_nonconformant (op_name , a_nr, a_nc, b_nr, b_nc);
56  else
57  {
58  r = MSparse<T> (a_nr, a_nc, (a.nnz () + b.nnz ()));
59 
60  octave_idx_type jx = 0;
61  for (octave_idx_type i = 0 ; i < a_nc ; i++)
62  {
63  octave_idx_type ja = a.cidx (i);
64  octave_idx_type ja_max = a.cidx (i+1);
65  bool ja_lt_max= ja < ja_max;
66 
67  octave_idx_type jb = b.cidx (i);
68  octave_idx_type jb_max = b.cidx (i+1);
69  bool jb_lt_max = jb < jb_max;
70 
71  while (ja_lt_max || jb_lt_max)
72  {
73  octave_quit ();
74  if ((! jb_lt_max) || (ja_lt_max && (a.ridx (ja) < b.ridx (jb))))
75  {
76  r.ridx (jx) = a.ridx (ja);
77  r.data (jx) = op (a.data (ja), 0.);
78  jx++;
79  ja++;
80  ja_lt_max= ja < ja_max;
81  }
82  else if ((! ja_lt_max)
83  || (jb_lt_max && (b.ridx (jb) < a.ridx (ja))))
84  {
85  r.ridx (jx) = b.ridx (jb);
86  r.data (jx) = op (0., b.data (jb));
87  jx++;
88  jb++;
89  jb_lt_max= jb < jb_max;
90  }
91  else
92  {
93  if (op (a.data (ja), b.data (jb)) != 0.)
94  {
95  r.data (jx) = op (a.data (ja), b.data (jb));
96  r.ridx (jx) = a.ridx (ja);
97  jx++;
98  }
99  ja++;
100  ja_lt_max= ja < ja_max;
101  jb++;
102  jb_lt_max= jb < jb_max;
103  }
104  }
105  r.cidx (i+1) = jx;
106  }
107 
108  a = r.maybe_compress ();
109  }
110 
111  return a;
112 }
113 
114 template <typename T>
115 MSparse<T>&
117 {
118  return plus_or_minus (a, b, std::plus<T> (), "operator +=");
119 }
120 
121 template <typename T>
122 MSparse<T>&
124 {
125  return plus_or_minus (a, b, std::minus<T> (), "operator -=");
126 }
127 
128 
129 // Element by element MSparse by scalar ops.
130 
131 template <class T, class OP>
132 MArray<T>
133 plus_or_minus (const MSparse<T>& a, const T& s, OP op)
134 {
135  octave_idx_type nr = a.rows ();
136  octave_idx_type nc = a.cols ();
137 
138  MArray<T> r (dim_vector (nr, nc), op (0.0, s));
139 
140  for (octave_idx_type j = 0; j < nc; j++)
141  for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++)
142  r.elem (a.ridx (i), j) = op (a.data (i), s);
143  return r;
144 }
145 
146 template <typename T>
147 MArray<T>
148 operator + (const MSparse<T>& a, const T& s)
149 {
150  return plus_or_minus (a, s, std::plus<T> ());
151 }
152 
153 template <typename T>
154 MArray<T>
155 operator - (const MSparse<T>& a, const T& s)
156 {
157  return plus_or_minus (a, s, std::minus<T> ());
158 }
159 
160 
161 template <class T, class OP>
163 times_or_divide (const MSparse<T>& a, const T& s, OP op)
164 {
165  octave_idx_type nr = a.rows ();
166  octave_idx_type nc = a.cols ();
167  octave_idx_type nz = a.nnz ();
168 
169  MSparse<T> r (nr, nc, nz);
170 
171  for (octave_idx_type i = 0; i < nz; i++)
172  {
173  r.data (i) = op (a.data (i), s);
174  r.ridx (i) = a.ridx (i);
175  }
176  for (octave_idx_type i = 0; i < nc + 1; i++)
177  r.cidx (i) = a.cidx (i);
178  r.maybe_compress (true);
179  return r;
180 }
181 
182 template <typename T>
184 operator * (const MSparse<T>& a, const T& s)
185 {
186  return times_or_divide (a, s, std::multiplies<T> ());
187 }
188 
189 template <typename T>
191 operator / (const MSparse<T>& a, const T& s)
192 {
193  return times_or_divide (a, s, std::divides<T> ());
194 }
195 
196 
197 // Element by element scalar by MSparse ops.
198 
199 template <class T, class OP>
200 MArray<T>
201 plus_or_minus (const T& s, const MSparse<T>& a, OP op)
202 {
203  octave_idx_type nr = a.rows ();
204  octave_idx_type nc = a.cols ();
205 
206  MArray<T> r (dim_vector (nr, nc), op (s, 0.0));
207 
208  for (octave_idx_type j = 0; j < nc; j++)
209  for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++)
210  r.elem (a.ridx (i), j) = op (s, a.data (i));
211  return r;
212 }
213 
214 template <typename T>
215 MArray<T>
216 operator + (const T& s, const MSparse<T>& a)
217 {
218  return plus_or_minus (s, a, std::plus<T> ());
219 }
220 
221 template <typename T>
222 MArray<T>
223 operator - (const T& s, const MSparse<T>& a)
224 {
225  return plus_or_minus (s, a, std::minus<T> ());
226 }
227 
228 template <class T, class OP>
230 times_or_divides (const T& s, const MSparse<T>& a, OP op)
231 {
232  octave_idx_type nr = a.rows ();
233  octave_idx_type nc = a.cols ();
234  octave_idx_type nz = a.nnz ();
235 
236  MSparse<T> r (nr, nc, nz);
237 
238  for (octave_idx_type i = 0; i < nz; i++)
239  {
240  r.data (i) = op (s, a.data (i));
241  r.ridx (i) = a.ridx (i);
242  }
243  for (octave_idx_type i = 0; i < nc + 1; i++)
244  r.cidx (i) = a.cidx (i);
245  r.maybe_compress (true);
246  return r;
247 }
248 
249 template <class T>
251 operator * (const T& s, const MSparse<T>& a)
252 {
253  return times_or_divides (s, a, std::multiplies<T> ());
254 }
255 
256 template <class T>
258 operator / (const T& s, const MSparse<T>& a)
259 {
260  return times_or_divides (s, a, std::divides<T> ());
261 }
262 
263 
264 // Element by element MSparse by MSparse ops.
265 
266 template <class T, class OP>
268 plus_or_minus (const MSparse<T>& a, const MSparse<T>& b, OP op,
269  const char* op_name, bool negate)
270 {
271  MSparse<T> r;
272 
273  octave_idx_type a_nr = a.rows ();
274  octave_idx_type a_nc = a.cols ();
275 
276  octave_idx_type b_nr = b.rows ();
277  octave_idx_type b_nc = b.cols ();
278 
279  if (a_nr == 1 && a_nc == 1)
280  {
281  if (a.elem (0,0) == 0.)
282  if (negate)
283  r = -MSparse<T> (b);
284  else
285  r = MSparse<T> (b);
286  else
287  {
288  r = MSparse<T> (b_nr, b_nc, op (a.data (0), 0.));
289 
290  for (octave_idx_type j = 0 ; j < b_nc ; j++)
291  {
292  octave_quit ();
293  octave_idx_type idxj = j * b_nr;
294  for (octave_idx_type i = b.cidx (j) ; i < b.cidx (j+1) ; i++)
295  {
296  octave_quit ();
297  r.data (idxj + b.ridx (i)) = op (a.data (0), b.data (i));
298  }
299  }
300  r.maybe_compress ();
301  }
302  }
303  else if (b_nr == 1 && b_nc == 1)
304  {
305  if (b.elem (0,0) == 0.)
306  r = MSparse<T> (a);
307  else
308  {
309  r = MSparse<T> (a_nr, a_nc, op (0.0, b.data (0)));
310 
311  for (octave_idx_type j = 0 ; j < a_nc ; j++)
312  {
313  octave_quit ();
314  octave_idx_type idxj = j * a_nr;
315  for (octave_idx_type i = a.cidx (j) ; i < a.cidx (j+1) ; i++)
316  {
317  octave_quit ();
318  r.data (idxj + a.ridx (i)) = op (a.data (i), b.data (0));
319  }
320  }
321  r.maybe_compress ();
322  }
323  }
324  else if (a_nr != b_nr || a_nc != b_nc)
325  gripe_nonconformant (op_name, a_nr, a_nc, b_nr, b_nc);
326  else
327  {
328  r = MSparse<T> (a_nr, a_nc, (a.nnz () + b.nnz ()));
329 
330  octave_idx_type jx = 0;
331  r.cidx (0) = 0;
332  for (octave_idx_type i = 0 ; i < a_nc ; i++)
333  {
334  octave_idx_type ja = a.cidx (i);
335  octave_idx_type ja_max = a.cidx (i+1);
336  bool ja_lt_max= ja < ja_max;
337 
338  octave_idx_type jb = b.cidx (i);
339  octave_idx_type jb_max = b.cidx (i+1);
340  bool jb_lt_max = jb < jb_max;
341 
342  while (ja_lt_max || jb_lt_max)
343  {
344  octave_quit ();
345  if ((! jb_lt_max) || (ja_lt_max && (a.ridx (ja) < b.ridx (jb))))
346  {
347  r.ridx (jx) = a.ridx (ja);
348  r.data (jx) = op (a.data (ja), 0.);
349  jx++;
350  ja++;
351  ja_lt_max= ja < ja_max;
352  }
353  else if ((! ja_lt_max)
354  || (jb_lt_max && (b.ridx (jb) < a.ridx (ja))))
355  {
356  r.ridx (jx) = b.ridx (jb);
357  r.data (jx) = op (0., b.data (jb));
358  jx++;
359  jb++;
360  jb_lt_max= jb < jb_max;
361  }
362  else
363  {
364  if (op (a.data (ja), b.data (jb)) != 0.)
365  {
366  r.data (jx) = op (a.data (ja), b.data (jb));
367  r.ridx (jx) = a.ridx (ja);
368  jx++;
369  }
370  ja++;
371  ja_lt_max= ja < ja_max;
372  jb++;
373  jb_lt_max= jb < jb_max;
374  }
375  }
376  r.cidx (i+1) = jx;
377  }
378 
379  r.maybe_compress ();
380  }
381 
382  return r;
383 }
384 
385 template <class T>
387 operator+ (const MSparse<T>& a, const MSparse<T>& b)
388 {
389  return plus_or_minus (a, b, std::plus<T> (), "operator +", false);
390 }
391 
392 template <class T>
394 operator- (const MSparse<T>& a, const MSparse<T>& b)
395 {
396  return plus_or_minus (a, b, std::minus<T> (), "operator -", true);
397 }
398 
399 template <class T>
401 product (const MSparse<T>& a, const MSparse<T>& b)
402 {
403  MSparse<T> r;
404 
405  octave_idx_type a_nr = a.rows ();
406  octave_idx_type a_nc = a.cols ();
407 
408  octave_idx_type b_nr = b.rows ();
409  octave_idx_type b_nc = b.cols ();
410 
411  if (a_nr == 1 && a_nc == 1)
412  {
413  if (a.elem (0,0) == 0.)
414  r = MSparse<T> (b_nr, b_nc);
415  else
416  {
417  r = MSparse<T> (b);
418  octave_idx_type b_nnz = b.nnz ();
419 
420  for (octave_idx_type i = 0 ; i < b_nnz ; i++)
421  {
422  octave_quit ();
423  r.data (i) = a.data (0) * r.data (i);
424  }
425  r.maybe_compress ();
426  }
427  }
428  else if (b_nr == 1 && b_nc == 1)
429  {
430  if (b.elem (0,0) == 0.)
431  r = MSparse<T> (a_nr, a_nc);
432  else
433  {
434  r = MSparse<T> (a);
435  octave_idx_type a_nnz = a.nnz ();
436 
437  for (octave_idx_type i = 0 ; i < a_nnz ; i++)
438  {
439  octave_quit ();
440  r.data (i) = r.data (i) * b.data (0);
441  }
442  r.maybe_compress ();
443  }
444  }
445  else if (a_nr != b_nr || a_nc != b_nc)
446  gripe_nonconformant ("product", a_nr, a_nc, b_nr, b_nc);
447  else
448  {
449  r = MSparse<T> (a_nr, a_nc, (a.nnz () > b.nnz () ? a.nnz () : b.nnz ()));
450 
451  octave_idx_type jx = 0;
452  r.cidx (0) = 0;
453  for (octave_idx_type i = 0 ; i < a_nc ; i++)
454  {
455  octave_idx_type ja = a.cidx (i);
456  octave_idx_type ja_max = a.cidx (i+1);
457  bool ja_lt_max= ja < ja_max;
458 
459  octave_idx_type jb = b.cidx (i);
460  octave_idx_type jb_max = b.cidx (i+1);
461  bool jb_lt_max = jb < jb_max;
462 
463  while (ja_lt_max || jb_lt_max)
464  {
465  octave_quit ();
466  if ((! jb_lt_max) || (ja_lt_max && (a.ridx (ja) < b.ridx (jb))))
467  {
468  ja++; ja_lt_max= ja < ja_max;
469  }
470  else if ((! ja_lt_max)
471  || (jb_lt_max && (b.ridx (jb) < a.ridx (ja))))
472  {
473  jb++; jb_lt_max= jb < jb_max;
474  }
475  else
476  {
477  if ((a.data (ja) * b.data (jb)) != 0.)
478  {
479  r.data (jx) = a.data (ja) * b.data (jb);
480  r.ridx (jx) = a.ridx (ja);
481  jx++;
482  }
483  ja++; ja_lt_max= ja < ja_max;
484  jb++; jb_lt_max= jb < jb_max;
485  }
486  }
487  r.cidx (i+1) = jx;
488  }
489 
490  r.maybe_compress ();
491  }
492 
493  return r;
494 }
495 
496 template <class T>
498 quotient (const MSparse<T>& a, const MSparse<T>& b)
499 {
500  MSparse<T> r;
501  T Zero = T ();
502 
503  octave_idx_type a_nr = a.rows ();
504  octave_idx_type a_nc = a.cols ();
505 
506  octave_idx_type b_nr = b.rows ();
507  octave_idx_type b_nc = b.cols ();
508 
509  if (a_nr == 1 && a_nc == 1)
510  {
511  T val = a.elem (0,0);
512  T fill = val / T ();
513  if (fill == T ())
514  {
515  octave_idx_type b_nnz = b.nnz ();
516  r = MSparse<T> (b);
517  for (octave_idx_type i = 0 ; i < b_nnz ; i++)
518  r.data (i) = val / r.data (i);
519  r.maybe_compress ();
520  }
521  else
522  {
523  r = MSparse<T> (b_nr, b_nc, fill);
524  for (octave_idx_type j = 0 ; j < b_nc ; j++)
525  {
526  octave_quit ();
527  octave_idx_type idxj = j * b_nr;
528  for (octave_idx_type i = b.cidx (j) ; i < b.cidx (j+1) ; i++)
529  {
530  octave_quit ();
531  r.data (idxj + b.ridx (i)) = val / b.data (i);
532  }
533  }
534  r.maybe_compress ();
535  }
536  }
537  else if (b_nr == 1 && b_nc == 1)
538  {
539  T val = b.elem (0,0);
540  T fill = T () / val;
541  if (fill == T ())
542  {
543  octave_idx_type a_nnz = a.nnz ();
544  r = MSparse<T> (a);
545  for (octave_idx_type i = 0 ; i < a_nnz ; i++)
546  r.data (i) = r.data (i) / val;
547  r.maybe_compress ();
548  }
549  else
550  {
551  r = MSparse<T> (a_nr, a_nc, fill);
552  for (octave_idx_type j = 0 ; j < a_nc ; j++)
553  {
554  octave_quit ();
555  octave_idx_type idxj = j * a_nr;
556  for (octave_idx_type i = a.cidx (j) ; i < a.cidx (j+1) ; i++)
557  {
558  octave_quit ();
559  r.data (idxj + a.ridx (i)) = a.data (i) / val;
560  }
561  }
562  r.maybe_compress ();
563  }
564  }
565  else if (a_nr != b_nr || a_nc != b_nc)
566  gripe_nonconformant ("quotient", a_nr, a_nc, b_nr, b_nc);
567  else
568  {
569  r = MSparse<T> (a_nr, a_nc, (Zero / Zero));
570 
571  for (octave_idx_type i = 0 ; i < a_nc ; i++)
572  {
573  octave_idx_type ja = a.cidx (i);
574  octave_idx_type ja_max = a.cidx (i+1);
575  bool ja_lt_max= ja < ja_max;
576 
577  octave_idx_type jb = b.cidx (i);
578  octave_idx_type jb_max = b.cidx (i+1);
579  bool jb_lt_max = jb < jb_max;
580 
581  while (ja_lt_max || jb_lt_max)
582  {
583  octave_quit ();
584  if ((! jb_lt_max) || (ja_lt_max && (a.ridx (ja) < b.ridx (jb))))
585  {
586  r.elem (a.ridx (ja),i) = a.data (ja) / Zero;
587  ja++; ja_lt_max= ja < ja_max;
588  }
589  else if ((! ja_lt_max)
590  || (jb_lt_max && (b.ridx (jb) < a.ridx (ja))))
591  {
592  r.elem (b.ridx (jb),i) = Zero / b.data (jb);
593  jb++; jb_lt_max= jb < jb_max;
594  }
595  else
596  {
597  r.elem (a.ridx (ja),i) = a.data (ja) / b.data (jb);
598  ja++; ja_lt_max= ja < ja_max;
599  jb++; jb_lt_max= jb < jb_max;
600  }
601  }
602  }
603 
604  r.maybe_compress (true);
605  }
606 
607  return r;
608 }
609 
610 
611 
612 // Unary MSparse ops.
613 
614 template <class T>
617 {
618  return a;
619 }
620 
621 template <class T>
624 {
625  MSparse<T> retval (a);
626  octave_idx_type nz = a.nnz ();
627  for (octave_idx_type i = 0; i < nz; i++)
628  retval.data (i) = - retval.data (i);
629  return retval;
630 }
octave_idx_type cols(void) const
Definition: Sparse.h:264
void gripe_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
T * data(void)
Definition: Sparse.h:509
octave_idx_type rows(void) const
Definition: Sparse.h:263
T & elem(octave_idx_type n)
Definition: Sparse.h:362
MSparse< T > operator*(const MSparse< T > &a, const T &s)
Definition: MSparse.cc:184
MSparse< T > & operator-=(MSparse< T > &a, const MSparse< T > &b)
Definition: MSparse.cc:123
MSparse< T > product(const MSparse< T > &a, const MSparse< T > &b)
Definition: MSparse.cc:401
octave_idx_type * cidx(void)
Definition: Sparse.h:531
T & elem(octave_idx_type n)
Definition: Array.h:380
Definition: MArray.h:36
octave_idx_type nnz(void) const
Definition: Sparse.h:248
MSparse< T > times_or_divide(const MSparse< T > &a, const T &s, OP op)
Definition: MSparse.cc:163
MSparse< T > quotient(const MSparse< T > &a, const MSparse< T > &b)
Definition: MSparse.cc:498
Sparse< T > maybe_compress(bool remove_zeros=false)
Definition: Sparse.h:469
MSparse< T > times_or_divides(const T &s, const MSparse< T > &a, OP op)
Definition: MSparse.cc:230
MSparse< T > operator/(const MSparse< T > &a, const T &s)
Definition: MSparse.cc:191
octave_idx_type * ridx(void)
Definition: Sparse.h:518
MSparse< T > & operator+=(MSparse< T > &a, const MSparse< T > &b)
Definition: MSparse.cc:116
MArray< T > operator-(const MSparse< T > &a, const T &s)
Definition: MSparse.cc:155
MSparse< T > & plus_or_minus(MSparse< T > &a, const MSparse< T > &b, OP op, const char *op_name)
Definition: MSparse.cc:44
MArray< T > operator+(const MSparse< T > &a, const T &s)
Definition: MSparse.cc:148