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