GNU Octave  3.8.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-2013 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>
43 MSparse<T>&
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) ||
75  (ja_lt_max && (a.ridx (ja) < b.ridx (jb))))
76  {
77  r.ridx (jx) = a.ridx (ja);
78  r.data (jx) = op (a.data (ja), 0.);
79  jx++;
80  ja++;
81  ja_lt_max= ja < ja_max;
82  }
83  else if (( !ja_lt_max ) ||
84  (jb_lt_max && (b.ridx (jb) < a.ridx (ja)) ) )
85  {
86  r.ridx (jx) = b.ridx (jb);
87  r.data (jx) = op (0., b.data (jb));
88  jx++;
89  jb++;
90  jb_lt_max= jb < jb_max;
91  }
92  else
93  {
94  if (op (a.data (ja), b.data (jb)) != 0.)
95  {
96  r.data (jx) = op (a.data (ja), b.data (jb));
97  r.ridx (jx) = a.ridx (ja);
98  jx++;
99  }
100  ja++;
101  ja_lt_max= ja < ja_max;
102  jb++;
103  jb_lt_max= jb < jb_max;
104  }
105  }
106  r.cidx (i+1) = jx;
107  }
108 
109  a = r.maybe_compress ();
110  }
111 
112  return a;
113 }
114 
115 template <typename T>
116 MSparse<T>&
117 operator += (MSparse<T>& a, const MSparse<T>& b)
118 {
119  return plus_or_minus (a, b, std::plus<T> (), "operator +=");
120 }
121 
122 template <typename T>
123 MSparse<T>&
124 operator -= (MSparse<T>& a, const MSparse<T>& b)
125 {
126  return plus_or_minus (a, b, std::minus<T> (), "operator -=");
127 }
128 
129 
130 // Element by element MSparse by scalar ops.
131 
132 template <class T, class OP>
133 MArray<T>
134 plus_or_minus (const MSparse<T>& a, const T& s, OP op)
135 {
136  octave_idx_type nr = a.rows ();
137  octave_idx_type nc = a.cols ();
138 
139  MArray<T> r (dim_vector (nr, nc), op (0.0, s));
140 
141  for (octave_idx_type j = 0; j < nc; j++)
142  for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++)
143  r.elem (a.ridx (i), j) = op (a.data (i), s);
144  return r;
145 }
146 
147 template <typename T>
148 MArray<T>
149 operator + (const MSparse<T>& a, const T& s)
150 {
151  return plus_or_minus (a, s, std::plus<T> ());
152 }
153 
154 template <typename T>
155 MArray<T>
156 operator - (const MSparse<T>& a, const T& s)
157 {
158  return plus_or_minus (a, s, std::minus<T> ());
159 }
160 
161 
162 template <class T, class OP>
163 MSparse<T>
164 times_or_divide (const MSparse<T>& a, const T& s, OP op)
165 {
166  octave_idx_type nr = a.rows ();
167  octave_idx_type nc = a.cols ();
168  octave_idx_type nz = a.nnz ();
169 
170  MSparse<T> r (nr, nc, nz);
171 
172  for (octave_idx_type i = 0; i < nz; i++)
173  {
174  r.data (i) = op (a.data (i), s);
175  r.ridx (i) = a.ridx (i);
176  }
177  for (octave_idx_type i = 0; i < nc + 1; i++)
178  r.cidx (i) = a.cidx (i);
179  r.maybe_compress (true);
180  return r;
181 }
182 
183 template <typename T>
184 MSparse<T>
185 operator * (const MSparse<T>& a, const T& s)
186 {
187  return times_or_divide (a, s, std::multiplies<T> ());
188 }
189 
190 template <typename T>
191 MSparse<T>
192 operator / (const MSparse<T>& a, const T& s)
193 {
194  return times_or_divide (a, s, std::divides<T> ());
195 }
196 
197 
198 // Element by element scalar by MSparse ops.
199 
200 template <class T, class OP>
201 MArray<T>
202 plus_or_minus (const T& s, const MSparse<T>& a, OP op)
203 {
204  octave_idx_type nr = a.rows ();
205  octave_idx_type nc = a.cols ();
206 
207  MArray<T> r (dim_vector (nr, nc), op (s, 0.0));
208 
209  for (octave_idx_type j = 0; j < nc; j++)
210  for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++)
211  r.elem (a.ridx (i), j) = op (s, a.data (i));
212  return r;
213 }
214 
215 template <typename T>
216 MArray<T>
217 operator + (const T& s, const MSparse<T>& a)
218 {
219  return plus_or_minus (s, a, std::plus<T> ());
220 }
221 
222 template <typename T>
223 MArray<T>
224 operator - (const T& s, const MSparse<T>& a)
225 {
226  return plus_or_minus (s, a, std::minus<T> ());
227 }
228 
229 template <class T, class OP>
230 MSparse<T>
231 times_or_divides (const T& s, const MSparse<T>& a, OP op)
232 {
233  octave_idx_type nr = a.rows ();
234  octave_idx_type nc = a.cols ();
235  octave_idx_type nz = a.nnz ();
236 
237  MSparse<T> r (nr, nc, nz);
238 
239  for (octave_idx_type i = 0; i < nz; i++)
240  {
241  r.data (i) = op (s, a.data (i));
242  r.ridx (i) = a.ridx (i);
243  }
244  for (octave_idx_type i = 0; i < nc + 1; i++)
245  r.cidx (i) = a.cidx (i);
246  r.maybe_compress (true);
247  return r;
248 }
249 
250 template <class T>
251 MSparse<T>
252 operator * (const T& s, const MSparse<T>& a)
253 {
254  return times_or_divides (s, a, std::multiplies<T> ());
255 }
256 
257 template <class T>
258 MSparse<T>
259 operator / (const T& s, const MSparse<T>& a)
260 {
261  return times_or_divides (s, a, std::divides<T> ());
262 }
263 
264 
265 // Element by element MSparse by MSparse ops.
266 
267 template <class T, class OP>
268 MSparse<T>
269 plus_or_minus (const MSparse<T>& a, const MSparse<T>& b, OP op,
270  const char* op_name, bool negate)
271 {
272  MSparse<T> r;
273 
274  octave_idx_type a_nr = a.rows ();
275  octave_idx_type a_nc = a.cols ();
276 
277  octave_idx_type b_nr = b.rows ();
278  octave_idx_type b_nc = b.cols ();
279 
280  if (a_nr == 1 && a_nc == 1)
281  {
282  if (a.elem (0,0) == 0.)
283  if (negate)
284  r = -MSparse<T> (b);
285  else
286  r = MSparse<T> (b);
287  else
288  {
289  r = MSparse<T> (b_nr, b_nc, op (a.data (0), 0.));
290 
291  for (octave_idx_type j = 0 ; j < b_nc ; j++)
292  {
293  octave_quit ();
294  octave_idx_type idxj = j * b_nr;
295  for (octave_idx_type i = b.cidx (j) ; i < b.cidx (j+1) ; i++)
296  {
297  octave_quit ();
298  r.data (idxj + b.ridx (i)) = op (a.data (0), b.data (i));
299  }
300  }
301  r.maybe_compress ();
302  }
303  }
304  else if (b_nr == 1 && b_nc == 1)
305  {
306  if (b.elem (0,0) == 0.)
307  r = MSparse<T> (a);
308  else
309  {
310  r = MSparse<T> (a_nr, a_nc, op (0.0, b.data (0)));
311 
312  for (octave_idx_type j = 0 ; j < a_nc ; j++)
313  {
314  octave_quit ();
315  octave_idx_type idxj = j * a_nr;
316  for (octave_idx_type i = a.cidx (j) ; i < a.cidx (j+1) ; i++)
317  {
318  octave_quit ();
319  r.data (idxj + a.ridx (i)) = op (a.data (i), b.data (0));
320  }
321  }
322  r.maybe_compress ();
323  }
324  }
325  else if (a_nr != b_nr || a_nc != b_nc)
326  gripe_nonconformant (op_name, a_nr, a_nc, b_nr, b_nc);
327  else
328  {
329  r = MSparse<T> (a_nr, a_nc, (a.nnz () + b.nnz ()));
330 
331  octave_idx_type jx = 0;
332  r.cidx (0) = 0;
333  for (octave_idx_type i = 0 ; i < a_nc ; i++)
334  {
335  octave_idx_type ja = a.cidx (i);
336  octave_idx_type ja_max = a.cidx (i+1);
337  bool ja_lt_max= ja < ja_max;
338 
339  octave_idx_type jb = b.cidx (i);
340  octave_idx_type jb_max = b.cidx (i+1);
341  bool jb_lt_max = jb < jb_max;
342 
343  while (ja_lt_max || jb_lt_max )
344  {
345  octave_quit ();
346  if ((! jb_lt_max) ||
347  (ja_lt_max && (a.ridx (ja) < b.ridx (jb))))
348  {
349  r.ridx (jx) = a.ridx (ja);
350  r.data (jx) = op (a.data (ja), 0.);
351  jx++;
352  ja++;
353  ja_lt_max= ja < ja_max;
354  }
355  else if (( !ja_lt_max ) ||
356  (jb_lt_max && (b.ridx (jb) < a.ridx (ja)) ) )
357  {
358  r.ridx (jx) = b.ridx (jb);
359  r.data (jx) = op (0., b.data (jb));
360  jx++;
361  jb++;
362  jb_lt_max= jb < jb_max;
363  }
364  else
365  {
366  if (op (a.data (ja), b.data (jb)) != 0.)
367  {
368  r.data (jx) = op (a.data (ja), b.data (jb));
369  r.ridx (jx) = a.ridx (ja);
370  jx++;
371  }
372  ja++;
373  ja_lt_max= ja < ja_max;
374  jb++;
375  jb_lt_max= jb < jb_max;
376  }
377  }
378  r.cidx (i+1) = jx;
379  }
380 
381  r.maybe_compress ();
382  }
383 
384  return r;
385 }
386 
387 template <class T>
388 MSparse<T>
389 operator+ (const MSparse<T>& a, const MSparse<T>& b)
390 {
391  return plus_or_minus (a, b, std::plus<T> (), "operator +", false);
392 }
393 
394 template <class T>
395 MSparse<T>
396 operator- (const MSparse<T>& a, const MSparse<T>& b)
397 {
398  return plus_or_minus (a, b, std::minus<T> (), "operator -", true);
399 }
400 
401 template <class T>
402 MSparse<T>
403 product (const MSparse<T>& a, const MSparse<T>& b)
404 {
405  MSparse<T> r;
406 
407  octave_idx_type a_nr = a.rows ();
408  octave_idx_type a_nc = a.cols ();
409 
410  octave_idx_type b_nr = b.rows ();
411  octave_idx_type b_nc = b.cols ();
412 
413  if (a_nr == 1 && a_nc == 1)
414  {
415  if (a.elem (0,0) == 0.)
416  r = MSparse<T> (b_nr, b_nc);
417  else
418  {
419  r = MSparse<T> (b);
420  octave_idx_type b_nnz = b.nnz ();
421 
422  for (octave_idx_type i = 0 ; i < b_nnz ; i++)
423  {
424  octave_quit ();
425  r.data (i) = a.data (0) * r.data (i);
426  }
427  r.maybe_compress ();
428  }
429  }
430  else if (b_nr == 1 && b_nc == 1)
431  {
432  if (b.elem (0,0) == 0.)
433  r = MSparse<T> (a_nr, a_nc);
434  else
435  {
436  r = MSparse<T> (a);
437  octave_idx_type a_nnz = a.nnz ();
438 
439  for (octave_idx_type i = 0 ; i < a_nnz ; i++)
440  {
441  octave_quit ();
442  r.data (i) = r.data (i) * b.data (0);
443  }
444  r.maybe_compress ();
445  }
446  }
447  else if (a_nr != b_nr || a_nc != b_nc)
448  gripe_nonconformant ("product", a_nr, a_nc, b_nr, b_nc);
449  else
450  {
451  r = MSparse<T> (a_nr, a_nc, (a.nnz () > b.nnz () ? a.nnz () : b.nnz ()));
452 
453  octave_idx_type jx = 0;
454  r.cidx (0) = 0;
455  for (octave_idx_type i = 0 ; i < a_nc ; i++)
456  {
457  octave_idx_type ja = a.cidx (i);
458  octave_idx_type ja_max = a.cidx (i+1);
459  bool ja_lt_max= ja < ja_max;
460 
461  octave_idx_type jb = b.cidx (i);
462  octave_idx_type jb_max = b.cidx (i+1);
463  bool jb_lt_max = jb < jb_max;
464 
465  while (ja_lt_max || jb_lt_max )
466  {
467  octave_quit ();
468  if ((! jb_lt_max) ||
469  (ja_lt_max && (a.ridx (ja) < b.ridx (jb))))
470  {
471  ja++; ja_lt_max= ja < ja_max;
472  }
473  else if (( !ja_lt_max ) ||
474  (jb_lt_max && (b.ridx (jb) < a.ridx (ja)) ) )
475  {
476  jb++; jb_lt_max= jb < jb_max;
477  }
478  else
479  {
480  if ((a.data (ja) * b.data (jb)) != 0.)
481  {
482  r.data (jx) = a.data (ja) * b.data (jb);
483  r.ridx (jx) = a.ridx (ja);
484  jx++;
485  }
486  ja++; ja_lt_max= ja < ja_max;
487  jb++; jb_lt_max= jb < jb_max;
488  }
489  }
490  r.cidx (i+1) = jx;
491  }
492 
493  r.maybe_compress ();
494  }
495 
496  return r;
497 }
498 
499 template <class T>
500 MSparse<T>
501 quotient (const MSparse<T>& a, const MSparse<T>& b)
502 {
503  MSparse<T> r;
504  T Zero = T ();
505 
506  octave_idx_type a_nr = a.rows ();
507  octave_idx_type a_nc = a.cols ();
508 
509  octave_idx_type b_nr = b.rows ();
510  octave_idx_type b_nc = b.cols ();
511 
512  if (a_nr == 1 && a_nc == 1)
513  {
514  T val = a.elem (0,0);
515  T fill = val / T ();
516  if (fill == T ())
517  {
518  octave_idx_type b_nnz = b.nnz ();
519  r = MSparse<T> (b);
520  for (octave_idx_type i = 0 ; i < b_nnz ; i++)
521  r.data (i) = val / r.data (i);
522  r.maybe_compress ();
523  }
524  else
525  {
526  r = MSparse<T> (b_nr, b_nc, fill);
527  for (octave_idx_type j = 0 ; j < b_nc ; j++)
528  {
529  octave_quit ();
530  octave_idx_type idxj = j * b_nr;
531  for (octave_idx_type i = b.cidx (j) ; i < b.cidx (j+1) ; i++)
532  {
533  octave_quit ();
534  r.data (idxj + b.ridx (i)) = val / b.data (i);
535  }
536  }
537  r.maybe_compress ();
538  }
539  }
540  else if (b_nr == 1 && b_nc == 1)
541  {
542  T val = b.elem (0,0);
543  T fill = T () / val;
544  if (fill == T ())
545  {
546  octave_idx_type a_nnz = a.nnz ();
547  r = MSparse<T> (a);
548  for (octave_idx_type i = 0 ; i < a_nnz ; i++)
549  r.data (i) = r.data (i) / val;
550  r.maybe_compress ();
551  }
552  else
553  {
554  r = MSparse<T> (a_nr, a_nc, fill);
555  for (octave_idx_type j = 0 ; j < a_nc ; j++)
556  {
557  octave_quit ();
558  octave_idx_type idxj = j * a_nr;
559  for (octave_idx_type i = a.cidx (j) ; i < a.cidx (j+1) ; i++)
560  {
561  octave_quit ();
562  r.data (idxj + a.ridx (i)) = a.data (i) / val;
563  }
564  }
565  r.maybe_compress ();
566  }
567  }
568  else if (a_nr != b_nr || a_nc != b_nc)
569  gripe_nonconformant ("quotient", a_nr, a_nc, b_nr, b_nc);
570  else
571  {
572  r = MSparse<T>( a_nr, a_nc, (Zero / Zero));
573 
574  for (octave_idx_type i = 0 ; i < a_nc ; i++)
575  {
576  octave_idx_type ja = a.cidx (i);
577  octave_idx_type ja_max = a.cidx (i+1);
578  bool ja_lt_max= ja < ja_max;
579 
580  octave_idx_type jb = b.cidx (i);
581  octave_idx_type jb_max = b.cidx (i+1);
582  bool jb_lt_max = jb < jb_max;
583 
584  while (ja_lt_max || jb_lt_max )
585  {
586  octave_quit ();
587  if ((! jb_lt_max) ||
588  (ja_lt_max && (a.ridx (ja) < b.ridx (jb))))
589  {
590  r.elem (a.ridx (ja),i) = a.data (ja) / Zero;
591  ja++; ja_lt_max= ja < ja_max;
592  }
593  else if (( !ja_lt_max ) ||
594  (jb_lt_max && (b.ridx (jb) < a.ridx (ja)) ) )
595  {
596  r.elem (b.ridx (jb),i) = Zero / b.data (jb);
597  jb++; jb_lt_max= jb < jb_max;
598  }
599  else
600  {
601  r.elem (a.ridx (ja),i) = a.data (ja) / b.data (jb);
602  ja++; ja_lt_max= ja < ja_max;
603  jb++; jb_lt_max= jb < jb_max;
604  }
605  }
606  }
607 
608  r.maybe_compress (true);
609  }
610 
611  return r;
612 }
613 
614 
615 
616 // Unary MSparse ops.
617 
618 template <class T>
619 MSparse<T>
620 operator + (const MSparse<T>& a)
621 {
622  return a;
623 }
624 
625 template <class T>
626 MSparse<T>
627 operator - (const MSparse<T>& a)
628 {
629  MSparse<T> retval (a);
630  octave_idx_type nz = a.nnz ();
631  for (octave_idx_type i = 0; i < nz; i++)
632  retval.data (i) = - retval.data (i);
633  return retval;
634 }
Generated on Mon Dec 30 2013 03:04:44 for GNU Octave by   doxygen 1.8.1.2