GNU Octave  4.4.1
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
sparse-xpow.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 2004-2018 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
10 the Free Software Foundation, either version 3 of the License, or
11 (at your option) any later version.
12 
13 Octave is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License 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 <https://www.gnu.org/licenses/>.
21 
22 */
23 
24 #if defined (HAVE_CONFIG_H)
25 # include "config.h"
26 #endif
27 
28 #include <cassert>
29 
30 #include <limits>
31 
32 #include "Array-util.h"
33 #include "oct-cmplx.h"
34 #include "quit.h"
35 
36 #include "error.h"
37 #include "ovl.h"
38 #include "utils.h"
39 
40 #include "dSparse.h"
41 #include "CSparse.h"
42 #include "ov-re-sparse.h"
43 #include "ov-cx-sparse.h"
44 #include "sparse-xpow.h"
45 
46 static inline int
47 xisint (double x)
48 {
49  return (octave::math::x_nint (x) == x
50  && ((x >= 0 && x < std::numeric_limits<int>::max ())
51  || (x <= 0 && x > std::numeric_limits<int>::min ())));
52 }
53 
54 // Safer pow functions. Only two make sense for sparse matrices, the
55 // others should all promote to full matrices.
56 
58 xpow (const SparseMatrix& a, double b)
59 {
61 
62  octave_idx_type nr = a.rows ();
63  octave_idx_type nc = a.cols ();
64 
65  if (nr == 0 || nc == 0 || nr != nc)
66  error ("for A^b, A must be a square matrix. Use .^ for elementwise power.");
67 
68  if (static_cast<int> (b) != b)
69  error ("use full(a) ^ full(b)");
70 
71  int btmp = static_cast<int> (b);
72  if (btmp == 0)
73  {
74  SparseMatrix tmp = SparseMatrix (nr, nr, nr);
75  for (octave_idx_type i = 0; i < nr; i++)
76  {
77  tmp.data (i) = 1.0;
78  tmp.ridx (i) = i;
79  }
80  for (octave_idx_type i = 0; i < nr + 1; i++)
81  tmp.cidx (i) = i;
82 
83  retval = tmp;
84  }
85  else
86  {
87  SparseMatrix atmp;
88  if (btmp < 0)
89  {
90  btmp = -btmp;
91 
92  octave_idx_type info;
93  double rcond = 0.0;
94  MatrixType mattyp (a);
95 
96  atmp = a.inverse (mattyp, info, rcond, 1);
97 
98  if (info == -1)
99  warning ("inverse: matrix singular to machine precision, rcond = %g", rcond);
100  }
101  else
102  atmp = a;
103 
104  SparseMatrix result (atmp);
105 
106  btmp--;
107 
108  while (btmp > 0)
109  {
110  if (btmp & 1)
111  result = result * atmp;
112 
113  btmp >>= 1;
114 
115  if (btmp > 0)
116  atmp = atmp * atmp;
117  }
118 
119  retval = result;
120  }
121 
122  return retval;
123 }
124 
126 xpow (const SparseComplexMatrix& a, double b)
127 {
129 
130  octave_idx_type nr = a.rows ();
131  octave_idx_type nc = a.cols ();
132 
133  if (nr == 0 || nc == 0 || nr != nc)
134  error ("for A^b, A must be a square matrix. Use .^ for elementwise power.");
135 
136  if (static_cast<int> (b) != b)
137  error ("use full(a) ^ full(b)");
138 
139  int btmp = static_cast<int> (b);
140  if (btmp == 0)
141  {
142  SparseMatrix tmp = SparseMatrix (nr, nr, nr);
143  for (octave_idx_type i = 0; i < nr; i++)
144  {
145  tmp.data (i) = 1.0;
146  tmp.ridx (i) = i;
147  }
148  for (octave_idx_type i = 0; i < nr + 1; i++)
149  tmp.cidx (i) = i;
150 
151  retval = tmp;
152  }
153  else
154  {
155  SparseComplexMatrix atmp;
156  if (btmp < 0)
157  {
158  btmp = -btmp;
159 
160  octave_idx_type info;
161  double rcond = 0.0;
162  MatrixType mattyp (a);
163 
164  atmp = a.inverse (mattyp, info, rcond, 1);
165 
166  if (info == -1)
167  warning ("inverse: matrix singular to machine precision, rcond = %g", rcond);
168  }
169  else
170  atmp = a;
171 
173 
174  btmp--;
175 
176  while (btmp > 0)
177  {
178  if (btmp & 1)
179  result = result * atmp;
180 
181  btmp >>= 1;
182 
183  if (btmp > 0)
184  atmp = atmp * atmp;
185  }
186 
187  retval = result;
188  }
189 
190  return retval;
191 }
192 
193 // Safer pow functions that work elementwise for matrices.
194 //
195 // op2 \ op1: s m cs cm
196 // +-- +---+---+----+----+
197 // scalar | | * | 3 | * | 9 |
198 // +---+---+----+----+
199 // matrix | 1 | 4 | 7 | 10 |
200 // +---+---+----+----+
201 // complex_scalar | * | 5 | * | 11 |
202 // +---+---+----+----+
203 // complex_matrix | 2 | 6 | 8 | 12 |
204 // +---+---+----+----+
205 //
206 // * -> not needed.
207 
208 // FIXME: these functions need to be fixed so that things
209 // like
210 //
211 // a = -1; b = [ 0, 0.5, 1 ]; r = a .^ b
212 //
213 // and
214 //
215 // a = -1; b = [ 0, 0.5, 1 ]; for i = 1:3, r(i) = a .^ b(i), end
216 //
217 // produce identical results. Also, it would be nice if -1^0.5
218 // produced a pure imaginary result instead of a complex number with a
219 // small real part. But perhaps that's really a problem with the math
220 // library...
221 
222 // Handle special case of scalar-sparse-matrix .^ sparse-matrix.
223 // Forwarding to the scalar elem_xpow function and then converting the
224 // result back to a sparse matrix is a bit wasteful but it does not
225 // seem worth the effort to optimize -- how often does this case come up
226 // in practice?
227 
228 template <typename S, typename SM>
229 inline octave_value
230 scalar_xpow (const S& a, const SM& b)
231 {
233 
234  if (val.iscomplex ())
235  return SparseComplexMatrix (val.complex_matrix_value ());
236  else
237  return SparseMatrix (val.matrix_value ());
238 }
239 
240 /*
241 %!assert (sparse (2) .^ [3, 4], sparse ([8, 16]))
242 %!assert <47775> (sparse (2i) .^ [3, 4], sparse ([-0-8i, 16]))
243 */
244 
245 // -*- 1 -*-
247 elem_xpow (double a, const SparseMatrix& b)
248 {
250 
251  octave_idx_type nr = b.rows ();
252  octave_idx_type nc = b.cols ();
253 
254  double d1, d2;
255 
256  if (a < 0.0 && ! b.all_integers (d1, d2))
257  {
258  Complex atmp (a);
259  ComplexMatrix result (nr, nc);
260 
261  for (octave_idx_type j = 0; j < nc; j++)
262  {
263  for (octave_idx_type i = 0; i < nr; i++)
264  {
265  octave_quit ();
266  result(i, j) = std::pow (atmp, b(i,j));
267  }
268  }
269 
270  retval = result;
271  }
272  else
273  {
274  Matrix result (nr, nc);
275 
276  for (octave_idx_type j = 0; j < nc; j++)
277  {
278  for (octave_idx_type i = 0; i < nr; i++)
279  {
280  octave_quit ();
281  result(i, j) = std::pow (a, b(i,j));
282  }
283  }
284 
285  retval = result;
286  }
287 
288  return retval;
289 }
290 
291 // -*- 2 -*-
294 {
295  octave_idx_type nr = b.rows ();
296  octave_idx_type nc = b.cols ();
297 
298  Complex atmp (a);
299  ComplexMatrix result (nr, nc);
300 
301  for (octave_idx_type j = 0; j < nc; j++)
302  {
303  for (octave_idx_type i = 0; i < nr; i++)
304  {
305  octave_quit ();
306  result(i, j) = std::pow (atmp, b(i,j));
307  }
308  }
309 
310  return result;
311 }
312 
313 // -*- 3 -*-
315 elem_xpow (const SparseMatrix& a, double b)
316 {
317  // FIXME: What should a .^ 0 give? Matlab gives a
318  // sparse matrix with same structure as a, which is strictly
319  // incorrect. Keep compatibility.
320 
322 
323  octave_idx_type nz = a.nnz ();
324 
325  if (b <= 0.0)
326  {
327  octave_idx_type nr = a.rows ();
328  octave_idx_type nc = a.cols ();
329 
330  if (static_cast<int> (b) != b && a.any_element_is_negative ())
331  {
332  ComplexMatrix result (nr, nc, Complex (std::pow (0.0, b)));
333 
334  // FIXME: avoid apparent GNU libm bug by
335  // converting A and B to complex instead of just A.
336  Complex btmp (b);
337 
338  for (octave_idx_type j = 0; j < nc; j++)
339  for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++)
340  {
341  octave_quit ();
342 
343  Complex atmp (a.data (i));
344 
345  result(a.ridx (i), j) = std::pow (atmp, btmp);
346  }
347 
349  }
350  else
351  {
352  Matrix result (nr, nc, (std::pow (0.0, b)));
353 
354  for (octave_idx_type j = 0; j < nc; j++)
355  for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++)
356  {
357  octave_quit ();
358  result(a.ridx (i), j) = std::pow (a.data (i), b);
359  }
360 
362  }
363  }
364  else if (static_cast<int> (b) != b && a.any_element_is_negative ())
365  {
367 
368  for (octave_idx_type i = 0; i < nz; i++)
369  {
370  octave_quit ();
371 
372  // FIXME: avoid apparent GNU libm bug by
373  // converting A and B to complex instead of just A.
374 
375  Complex atmp (a.data (i));
376  Complex btmp (b);
377 
378  result.data (i) = std::pow (atmp, btmp);
379  }
380 
381  result.maybe_compress (true);
382 
383  retval = result;
384  }
385  else
386  {
388 
389  for (octave_idx_type i = 0; i < nz; i++)
390  {
391  octave_quit ();
392  result.data (i) = std::pow (a.data (i), b);
393  }
394 
395  result.maybe_compress (true);
396 
397  retval = result;
398  }
399 
400  return retval;
401 }
402 
403 // -*- 4 -*-
406 {
408 
409  octave_idx_type nr = a.rows ();
410  octave_idx_type nc = a.cols ();
411 
412  octave_idx_type b_nr = b.rows ();
413  octave_idx_type b_nc = b.cols ();
414 
415  if (a.numel () == 1 && b.numel () > 1)
416  return scalar_xpow (a(0), b);
417 
418  if (nr != b_nr || nc != b_nc)
419  octave::err_nonconformant ("operator .^", nr, nc, b_nr, b_nc);
420 
421  int convert_to_complex = 0;
422  for (octave_idx_type j = 0; j < nc; j++)
423  for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++)
424  {
425  if (a.data(i) < 0.0)
426  {
427  double btmp = b (a.ridx (i), j);
428  if (static_cast<int> (btmp) != btmp)
429  {
430  convert_to_complex = 1;
431  goto done;
432  }
433  }
434  }
435 
436 done:
437 
438  // This is a dumb operator for sparse matrices anyway, and there is
439  // no sensible way to handle the 0.^0 versus the 0.^x cases. Therefore
440  // allocate a full matrix filled for the 0.^0 case and shrink it later
441  // as needed.
442 
443  if (convert_to_complex)
444  {
445  SparseComplexMatrix complex_result (nr, nc, Complex (1.0, 0.0));
446 
447  for (octave_idx_type j = 0; j < nc; j++)
448  {
449  for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++)
450  {
451  octave_quit ();
452  complex_result.xelem (a.ridx (i), j) =
453  std::pow (Complex (a.data (i)), Complex (b(a.ridx (i), j)));
454  }
455  }
456  complex_result.maybe_compress (true);
457  retval = complex_result;
458  }
459  else
460  {
461  SparseMatrix result (nr, nc, 1.0);
462 
463  for (octave_idx_type j = 0; j < nc; j++)
464  {
465  for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++)
466  {
467  octave_quit ();
468  result.xelem (a.ridx (i), j) = std::pow (a.data (i),
469  b(a.ridx (i), j));
470  }
471  }
472  result.maybe_compress (true);
473  retval = result;
474  }
475 
476  return retval;
477 }
478 
479 // -*- 5 -*-
481 elem_xpow (const SparseMatrix& a, const Complex& b)
482 {
484 
485  if (b == 0.0)
486  // Can this case ever happen, due to automatic retyping with maybe_mutate?
487  retval = octave_value (NDArray (a.dims (), 1));
488  else
489  {
490  octave_idx_type nz = a.nnz ();
492 
493  for (octave_idx_type i = 0; i < nz; i++)
494  {
495  octave_quit ();
496  result.data (i) = std::pow (Complex (a.data (i)), b);
497  }
498 
499  result.maybe_compress (true);
500 
501  retval = result;
502  }
503 
504  return retval;
505 }
506 
507 // -*- 6 -*-
510 {
511  octave_idx_type nr = a.rows ();
512  octave_idx_type nc = a.cols ();
513 
514  octave_idx_type b_nr = b.rows ();
515  octave_idx_type b_nc = b.cols ();
516 
517  if (a.numel () == 1 && b.numel () > 1)
518  return scalar_xpow (a(0), b);
519 
520  if (nr != b_nr || nc != b_nc)
521  octave::err_nonconformant ("operator .^", nr, nc, b_nr, b_nc);
522 
523  SparseComplexMatrix result (nr, nc, Complex (1.0, 0.0));
524  for (octave_idx_type j = 0; j < nc; j++)
525  {
526  for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++)
527  {
528  octave_quit ();
529  result.xelem (a.ridx(i), j) = std::pow (a.data (i), b(a.ridx (i), j));
530  }
531  }
532 
533  result.maybe_compress (true);
534 
535  return result;
536 }
537 
538 // -*- 7 -*-
540 elem_xpow (const Complex& a, const SparseMatrix& b)
541 {
542  octave_idx_type nr = b.rows ();
543  octave_idx_type nc = b.cols ();
544 
545  ComplexMatrix result (nr, nc);
546 
547  for (octave_idx_type j = 0; j < nc; j++)
548  {
549  for (octave_idx_type i = 0; i < nr; i++)
550  {
551  octave_quit ();
552  double btmp = b (i, j);
553  if (xisint (btmp))
554  result (i, j) = std::pow (a, static_cast<int> (btmp));
555  else
556  result (i, j) = std::pow (a, btmp);
557  }
558  }
559 
560  return result;
561 }
562 
563 // -*- 8 -*-
566 {
567  octave_idx_type nr = b.rows ();
568  octave_idx_type nc = b.cols ();
569 
570  ComplexMatrix result (nr, nc);
571  for (octave_idx_type j = 0; j < nc; j++)
572  for (octave_idx_type i = 0; i < nr; i++)
573  {
574  octave_quit ();
575  result (i, j) = std::pow (a, b (i, j));
576  }
577 
578  return result;
579 }
580 
581 // -*- 9 -*-
584 {
586 
587  if (b <= 0)
588  {
589  octave_idx_type nr = a.rows ();
590  octave_idx_type nc = a.cols ();
591 
592  ComplexMatrix result (nr, nc, Complex (std::pow (0.0, b)));
593 
594  if (xisint (b))
595  {
596  for (octave_idx_type j = 0; j < nc; j++)
597  for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++)
598  {
599  octave_quit ();
600  result (a.ridx (i), j) =
601  std::pow (a.data (i), static_cast<int> (b));
602  }
603  }
604  else
605  {
606  for (octave_idx_type j = 0; j < nc; j++)
607  for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++)
608  {
609  octave_quit ();
610  result (a.ridx (i), j) = std::pow (a.data (i), b);
611  }
612  }
613 
614  retval = result;
615  }
616  else
617  {
618  octave_idx_type nz = a.nnz ();
619 
621 
622  if (xisint (b))
623  {
624  for (octave_idx_type i = 0; i < nz; i++)
625  {
626  octave_quit ();
627  result.data (i) = std::pow (a.data (i), static_cast<int> (b));
628  }
629  }
630  else
631  {
632  for (octave_idx_type i = 0; i < nz; i++)
633  {
634  octave_quit ();
635  result.data (i) = std::pow (a.data (i), b);
636  }
637  }
638 
639  result.maybe_compress (true);
640 
641  retval = result;
642  }
643 
644  return retval;
645 }
646 
647 // -*- 10 -*-
650 {
651  octave_idx_type nr = a.rows ();
652  octave_idx_type nc = a.cols ();
653 
654  octave_idx_type b_nr = b.rows ();
655  octave_idx_type b_nc = b.cols ();
656 
657  if (a.numel () == 1 && b.numel () > 1)
658  return scalar_xpow (a(0), b);
659 
660  if (nr != b_nr || nc != b_nc)
661  octave::err_nonconformant ("operator .^", nr, nc, b_nr, b_nc);
662 
663  SparseComplexMatrix result (nr, nc, Complex (1.0, 0.0));
664  for (octave_idx_type j = 0; j < nc; j++)
665  {
666  for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++)
667  {
668  octave_quit ();
669  double btmp = b(a.ridx (i), j);
670 
671  if (xisint (btmp))
672  result.xelem (a.ridx (i), j) = std::pow (a.data (i),
673  static_cast<int> (btmp));
674  else
675  result.xelem (a.ridx (i), j) = std::pow (a.data (i), btmp);
676  }
677  }
678 
679  result.maybe_compress (true);
680 
681  return result;
682 }
683 
684 // -*- 11 -*-
687 {
689 
690  if (b == 0.0)
691  // Can this case ever happen, due to automatic retyping with maybe_mutate?
692  retval = octave_value (NDArray (a.dims (), 1));
693  else
694  {
695 
696  octave_idx_type nz = a.nnz ();
697 
699 
700  for (octave_idx_type i = 0; i < nz; i++)
701  {
702  octave_quit ();
703  result.data (i) = std::pow (a.data (i), b);
704  }
705 
706  result.maybe_compress (true);
707 
708  retval = result;
709  }
710 
711  return retval;
712 }
713 
714 // -*- 12 -*-
717 {
718  octave_idx_type nr = a.rows ();
719  octave_idx_type nc = a.cols ();
720 
721  octave_idx_type b_nr = b.rows ();
722  octave_idx_type b_nc = b.cols ();
723 
724  if (a.numel () == 1 && b.numel () > 1)
725  return scalar_xpow (a(0), b);
726 
727  if (nr != b_nr || nc != b_nc)
728  octave::err_nonconformant ("operator .^", nr, nc, b_nr, b_nc);
729 
730  SparseComplexMatrix result (nr, nc, Complex (1.0, 0.0));
731  for (octave_idx_type j = 0; j < nc; j++)
732  {
733  for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++)
734  {
735  octave_quit ();
736  result.xelem (a.ridx (i), j) = std::pow (a.data (i),
737  b(a.ridx (i), j));
738  }
739  }
740  result.maybe_compress (true);
741 
742  return result;
743 }
octave_value xpow(const SparseMatrix &a, double b)
Definition: sparse-xpow.cc:58
octave_value scalar_xpow(const S &a, const SM &b)
Definition: sparse-xpow.cc:230
static int xisint(double x)
Definition: sparse-xpow.cc:47
identity matrix If supplied two scalar respectively For allows like xample val
Definition: data.cc:4986
T & xelem(octave_idx_type n)
Definition: Sparse.h:301
void error(const char *fmt,...)
Definition: error.cc:578
octave_idx_type b_nr
Definition: sylvester.cc:76
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:400
done
Definition: syscalls.cc:251
octave_value elem_xpow(double a, const SparseMatrix &b)
Definition: sparse-xpow.cc:247
octave_int< T > pow(const octave_int< T > &a, const octave_int< T > &b)
double tmp
Definition: data.cc:6252
void err_nonconformant(const char *op, octave_idx_type op1_len, octave_idx_type op2_len)
octave_value retval
Definition: data.cc:6246
Definition: dMatrix.h:36
With real return the complex result
Definition: data.cc:3260
void warning(const char *fmt,...)
Definition: error.cc:801
return octave_value(v1.char_array_value() . concat(v2.char_array_value(), ra_idx),((a1.is_sq_string()||a2.is_sq_string()) ? '\'' :'"'))
charNDArray max(char d, const charNDArray &m)
Definition: chNDArray.cc:227
octave_idx_type b_nc
Definition: sylvester.cc:77
b
Definition: cellfun.cc:400
Sparse< T > maybe_compress(bool remove_zeros=false)
Definition: Sparse.h:438
for i
Definition: data.cc:5264
std::complex< double > Complex
Definition: oct-cmplx.h:31
T x_nint(T x)
Definition: lo-mappers.h:284
F77_RET_T const F77_REAL const F77_REAL F77_REAL &F77_RET_T const F77_DBLE const F77_DBLE F77_DBLE &F77_RET_T const F77_DBLE F77_DBLE &F77_RET_T const F77_REAL F77_REAL &F77_RET_T const F77_DBLE * x
charNDArray min(char d, const charNDArray &m)
Definition: chNDArray.cc:204