Mercurial > octave
annotate libinterp/corefcn/sparse-xpow.cc @ 33608:5fba13104493 bytecode-interpreter tip
maint: merge default to bytecode-interpreter.
author | Nicholas R. Jankowski <jankowski.nicholas@gmail.com> |
---|---|
date | Sat, 18 May 2024 22:40:00 -0400 |
parents | 06a308cae32c |
children |
rev | line source |
---|---|
27923
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
1 //////////////////////////////////////////////////////////////////////// |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
2 // |
32632
2e484f9f1f18
maint: update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents:
31706
diff
changeset
|
3 // Copyright (C) 1998-2024 The Octave Project Developers |
27923
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
4 // |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
5 // See the file COPYRIGHT.md in the top-level directory of this |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
6 // distribution or <https://octave.org/copyright/>. |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
7 // |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
8 // This file is part of Octave. |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
9 // |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
10 // Octave is free software: you can redistribute it and/or modify it |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
11 // under the terms of the GNU General Public License as published by |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
12 // the Free Software Foundation, either version 3 of the License, or |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
13 // (at your option) any later version. |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
14 // |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
15 // Octave is distributed in the hope that it will be useful, but |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
16 // WITHOUT ANY WARRANTY; without even the implied warranty of |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
18 // GNU General Public License for more details. |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
19 // |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
20 // You should have received a copy of the GNU General Public License |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
21 // along with Octave; see the file COPYING. If not, see |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
22 // <https://www.gnu.org/licenses/>. |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
23 // |
bd51beb6205e
update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents:
27919
diff
changeset
|
24 //////////////////////////////////////////////////////////////////////// |
5164 | 25 |
21724
aba2e6293dd8
use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents:
21317
diff
changeset
|
26 #if defined (HAVE_CONFIG_H) |
21301
40de9f8f23a6
Use '#include "config.h"' rather than <config.h>.
Rik <rik@octave.org>
parents:
21200
diff
changeset
|
27 # include "config.h" |
5164 | 28 #endif |
29 | |
15215
9020dddc925a
use std::numeric_limits for integer max and min values
John W. Eaton <jwe@octave.org>
parents:
15195
diff
changeset
|
30 #include <limits> |
5164 | 31 |
32 #include "Array-util.h" | |
33 #include "oct-cmplx.h" | |
34 #include "quit.h" | |
35 | |
36 #include "error.h" | |
20940
48b2ad5ee801
maint: Rename oct-obj.[cc|h] to ovl.[cc|h] for clarity.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
37 #include "ovl.h" |
5164 | 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 | |
31605
e88a07dec498
maint: Use macros to begin/end C++ namespaces.
Rik <rik@octave.org>
parents:
31235
diff
changeset
|
46 OCTAVE_BEGIN_NAMESPACE(octave) |
29990
b839c36fd106
move sparse xdiv and xpow operator functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
29358
diff
changeset
|
47 |
31226
44cf6bbeeca9
sparse-xpow.cc: Change return type from int to bool for internal function
Arun Giridhar <arungiridhar@gmail.com>
parents:
31225
diff
changeset
|
48 static inline bool |
5164 | 49 xisint (double x) |
50 { | |
21782
2aef506f3fec
use namespace for lo-mappers.h functions
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
51 return (octave::math::x_nint (x) == x |
15215
9020dddc925a
use std::numeric_limits for integer max and min values
John W. Eaton <jwe@octave.org>
parents:
15195
diff
changeset
|
52 && ((x >= 0 && x < std::numeric_limits<int>::max ()) |
9020dddc925a
use std::numeric_limits for integer max and min values
John W. Eaton <jwe@octave.org>
parents:
15195
diff
changeset
|
53 || (x <= 0 && x > std::numeric_limits<int>::min ()))); |
5164 | 54 } |
55 | |
21751
b571fc85953f
maint: Use two spaces after period to indicate sentence break.
Rik <rik@octave.org>
parents:
21724
diff
changeset
|
56 // Safer pow functions. Only two make sense for sparse matrices, the |
5164 | 57 // others should all promote to full matrices. |
58 | |
59 octave_value | |
60 xpow (const SparseMatrix& a, double b) | |
61 { | |
62 octave_value retval; | |
63 | |
5275 | 64 octave_idx_type nr = a.rows (); |
65 octave_idx_type nc = a.cols (); | |
5164 | 66 |
31234
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
67 if (nr == 0 || nc == 0) |
31235
b542b88ad3b6
maint: Use space between function name and '(' in sparse-xpow.cc.
Rik <rik@octave.org>
parents:
31234
diff
changeset
|
68 return SparseMatrix (); |
31234
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
69 |
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
70 // If we are here, A is not empty ==> A needs to be square. |
31231
a026fb2be108
sparse-xpow.cc: Return empty matrix for empty input (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31227
diff
changeset
|
71 if (nr != nc) |
21751
b571fc85953f
maint: Use two spaces after period to indicate sentence break.
Rik <rik@octave.org>
parents:
21724
diff
changeset
|
72 error ("for A^b, A must be a square matrix. Use .^ for elementwise power."); |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
73 |
28144
c31c9eaa1f28
sparse-xpow.cc: use xisint instead of static_cast<int> to check int values
John W. Eaton <jwe@octave.org>
parents:
27923
diff
changeset
|
74 if (! xisint (b)) |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
75 error ("use full(a) ^ full(b)"); |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
76 |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
77 int btmp = static_cast<int> (b); |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
78 if (btmp == 0) |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
79 { |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
80 SparseMatrix tmp = SparseMatrix (nr, nr, nr); |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
81 for (octave_idx_type i = 0; i < nr; i++) |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
82 { |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
83 tmp.data (i) = 1.0; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
84 tmp.ridx (i) = i; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
85 } |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
86 for (octave_idx_type i = 0; i < nr + 1; i++) |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
87 tmp.cidx (i) = i; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
88 |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
89 retval = tmp; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
90 } |
5164 | 91 else |
92 { | |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
93 SparseMatrix atmp; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
94 if (btmp < 0) |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
95 { |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
96 btmp = -btmp; |
5164 | 97 |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
98 octave_idx_type info; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
99 double rcond = 0.0; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
100 MatrixType mattyp (a); |
5164 | 101 |
31234
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
102 // FIXME: This causes an error if the input sparse matrix is all-zeros. |
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
103 // That behavior is inconsistent with A ^ b when A is a full all-zeros |
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
104 // matrix, which just returns Inf of the same size with a warning. |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
105 atmp = a.inverse (mattyp, info, rcond, 1); |
5164 | 106 |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
107 if (info == -1) |
22197
e43d83253e28
refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents:
22022
diff
changeset
|
108 warning ("inverse: matrix singular to machine precision, rcond = %g", rcond); |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
109 } |
5164 | 110 else |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
111 atmp = a; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
112 |
31235
b542b88ad3b6
maint: Use space between function name and '(' in sparse-xpow.cc.
Rik <rik@octave.org>
parents:
31234
diff
changeset
|
113 if (atmp.nnz () == 0) // Fast return for all-zeros matrix |
31234
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
114 return atmp; |
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
115 |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
116 SparseMatrix result (atmp); |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
117 |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
118 btmp--; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
119 |
31224
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
120 // There are two approaches to the actual exponentiation. |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
121 // Exponentiation by squaring uses only a logarithmic number |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
122 // of multiplications but the matrices it multiplies tend to be dense |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
123 // towards the end. |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
124 // Linear multiplication uses a linear number of multiplications |
31235
b542b88ad3b6
maint: Use space between function name and '(' in sparse-xpow.cc.
Rik <rik@octave.org>
parents:
31234
diff
changeset
|
125 // but one of the matrices it uses will be as sparse as the original |
b542b88ad3b6
maint: Use space between function name and '(' in sparse-xpow.cc.
Rik <rik@octave.org>
parents:
31234
diff
changeset
|
126 // matrix. |
31224
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
127 // |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
128 // The time to multiply fixed-size matrices is strongly affected by their |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
129 // sparsity. Denser matrices take much longer to multiply together. |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
130 // See this URL for a worked-through example: |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
131 // https://octave.discourse.group/t/3216/4 |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
132 // |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
133 // The tradeoff is between many fast multiplications or a few slow ones. |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
134 // |
31235
b542b88ad3b6
maint: Use space between function name and '(' in sparse-xpow.cc.
Rik <rik@octave.org>
parents:
31234
diff
changeset
|
135 // Large exponents favor the squaring technique, and sparse matrices |
b542b88ad3b6
maint: Use space between function name and '(' in sparse-xpow.cc.
Rik <rik@octave.org>
parents:
31234
diff
changeset
|
136 // favor linear multiplication. |
31224
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
137 // |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
138 // We calculate a threshold based on the sparsity of the input |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
139 // and use squaring for exponents larger than that. |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
140 // |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
141 // FIXME: Improve this threshold calculation. |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
142 |
31235
b542b88ad3b6
maint: Use space between function name and '(' in sparse-xpow.cc.
Rik <rik@octave.org>
parents:
31234
diff
changeset
|
143 uint64_t sparsity = atmp.numel () / atmp.nnz (); // reciprocal of density |
31227
0dec459a4064
sparse-xpow.cc: Performance tweak for threshold selection
Arun Giridhar <arungiridhar@gmail.com>
parents:
31226
diff
changeset
|
144 int threshold = (sparsity >= 1000) ? 40 |
31607
aac27ad79be6
maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents:
31605
diff
changeset
|
145 : (sparsity >= 100) ? 20 |
aac27ad79be6
maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents:
31605
diff
changeset
|
146 : 3; |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
147 |
31224
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
148 if (btmp > threshold) // use squaring technique |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
149 { |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
150 while (btmp > 0) |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
151 { |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
152 if (btmp & 1) |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
153 result = result * atmp; |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
154 |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
155 btmp >>= 1; |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
156 |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
157 if (btmp > 0) |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
158 atmp = atmp * atmp; |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
159 } |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
160 } |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
161 else // use linear multiplication |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
162 { |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
163 for (int i = 0; i < btmp; i++) |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
164 result = result * atmp; |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
165 } |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
166 |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
167 retval = result; |
5164 | 168 } |
169 | |
170 return retval; | |
171 } | |
172 | |
173 octave_value | |
174 xpow (const SparseComplexMatrix& a, double b) | |
175 { | |
176 octave_value retval; | |
177 | |
5275 | 178 octave_idx_type nr = a.rows (); |
179 octave_idx_type nc = a.cols (); | |
5164 | 180 |
31234
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
181 if (nr == 0 || nc == 0) |
31235
b542b88ad3b6
maint: Use space between function name and '(' in sparse-xpow.cc.
Rik <rik@octave.org>
parents:
31234
diff
changeset
|
182 return SparseMatrix (); |
31234
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
183 |
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
184 // If we are here, A is not empty ==> A needs to be square. |
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
185 if (nr != nc) |
21751
b571fc85953f
maint: Use two spaces after period to indicate sentence break.
Rik <rik@octave.org>
parents:
21724
diff
changeset
|
186 error ("for A^b, A must be a square matrix. Use .^ for elementwise power."); |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
187 |
28144
c31c9eaa1f28
sparse-xpow.cc: use xisint instead of static_cast<int> to check int values
John W. Eaton <jwe@octave.org>
parents:
27923
diff
changeset
|
188 if (! xisint (b)) |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
189 error ("use full(a) ^ full(b)"); |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
190 |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
191 int btmp = static_cast<int> (b); |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
192 if (btmp == 0) |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
193 { |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
194 SparseMatrix tmp = SparseMatrix (nr, nr, nr); |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
195 for (octave_idx_type i = 0; i < nr; i++) |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
196 { |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
197 tmp.data (i) = 1.0; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
198 tmp.ridx (i) = i; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
199 } |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
200 for (octave_idx_type i = 0; i < nr + 1; i++) |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
201 tmp.cidx (i) = i; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
202 |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
203 retval = tmp; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
204 } |
5164 | 205 else |
206 { | |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
207 SparseComplexMatrix atmp; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
208 if (btmp < 0) |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
209 { |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
210 btmp = -btmp; |
5164 | 211 |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
212 octave_idx_type info; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
213 double rcond = 0.0; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
214 MatrixType mattyp (a); |
5164 | 215 |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
216 atmp = a.inverse (mattyp, info, rcond, 1); |
5164 | 217 |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
218 if (info == -1) |
22197
e43d83253e28
refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents:
22022
diff
changeset
|
219 warning ("inverse: matrix singular to machine precision, rcond = %g", rcond); |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
220 } |
5164 | 221 else |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
222 atmp = a; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
223 |
31235
b542b88ad3b6
maint: Use space between function name and '(' in sparse-xpow.cc.
Rik <rik@octave.org>
parents:
31234
diff
changeset
|
224 if (atmp.nnz () == 0) // Fast return for all-zeros matrix |
31234
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
225 return atmp; |
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
226 |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
227 SparseComplexMatrix result (atmp); |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
228 |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
229 btmp--; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
230 |
31225
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
231 // Select multiplication sequence based on sparsity of atmp. |
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
232 // See the long comment in xpow (const SparseMatrix& a, double b) |
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
233 // for more details. |
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
234 // |
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
235 // FIXME: Improve this threshold calculation. |
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
236 |
31235
b542b88ad3b6
maint: Use space between function name and '(' in sparse-xpow.cc.
Rik <rik@octave.org>
parents:
31234
diff
changeset
|
237 uint64_t sparsity = atmp.numel () / atmp.nnz (); // reciprocal of density |
31227
0dec459a4064
sparse-xpow.cc: Performance tweak for threshold selection
Arun Giridhar <arungiridhar@gmail.com>
parents:
31226
diff
changeset
|
238 int threshold = (sparsity >= 1000) ? 40 |
31607
aac27ad79be6
maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents:
31605
diff
changeset
|
239 : (sparsity >= 100) ? 20 |
aac27ad79be6
maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents:
31605
diff
changeset
|
240 : 3; |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
241 |
31225
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
242 if (btmp > threshold) // use squaring technique |
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
243 { |
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
244 while (btmp > 0) |
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
245 { |
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
246 if (btmp & 1) |
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
247 result = result * atmp; |
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
248 |
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
249 btmp >>= 1; |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
250 |
31225
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
251 if (btmp > 0) |
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
252 atmp = atmp * atmp; |
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
253 } |
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
254 } |
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
255 else // use linear multiplication |
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
256 { |
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
257 for (int i = 0; i < btmp; i++) |
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
258 result = result * atmp; |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
259 } |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
260 |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
261 retval = result; |
5164 | 262 } |
263 | |
264 return retval; | |
265 } | |
266 | |
267 // Safer pow functions that work elementwise for matrices. | |
268 // | |
269 // op2 \ op1: s m cs cm | |
270 // +-- +---+---+----+----+ | |
271 // scalar | | * | 3 | * | 9 | | |
272 // +---+---+----+----+ | |
273 // matrix | 1 | 4 | 7 | 10 | | |
274 // +---+---+----+----+ | |
275 // complex_scalar | * | 5 | * | 11 | | |
276 // +---+---+----+----+ | |
277 // complex_matrix | 2 | 6 | 8 | 12 | | |
278 // +---+---+----+----+ | |
279 // | |
280 // * -> not needed. | |
281 | |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17749
diff
changeset
|
282 // FIXME: these functions need to be fixed so that things |
5164 | 283 // like |
284 // | |
285 // a = -1; b = [ 0, 0.5, 1 ]; r = a .^ b | |
286 // | |
287 // and | |
288 // | |
289 // a = -1; b = [ 0, 0.5, 1 ]; for i = 1:3, r(i) = a .^ b(i), end | |
290 // | |
291 // produce identical results. Also, it would be nice if -1^0.5 | |
292 // produced a pure imaginary result instead of a complex number with a | |
293 // small real part. But perhaps that's really a problem with the math | |
294 // library... | |
295 | |
15282
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
296 // Handle special case of scalar-sparse-matrix .^ sparse-matrix. |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
297 // Forwarding to the scalar elem_xpow function and then converting the |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
298 // result back to a sparse matrix is a bit wasteful but it does not |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
299 // seem worth the effort to optimize -- how often does this case come up |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
300 // in practice? |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
301 |
21139
538b57866b90
consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents:
21121
diff
changeset
|
302 template <typename S, typename SM> |
15282
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
303 inline octave_value |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
304 scalar_xpow (const S& a, const SM& b) |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
305 { |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
306 octave_value val = elem_xpow (a, b); |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
307 |
23581
c3075ae020e1
maint: Deprecate is_complex_type and replace with iscomplex.
Rik <rik@octave.org>
parents:
23220
diff
changeset
|
308 if (val.iscomplex ()) |
15282
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
309 return SparseComplexMatrix (val.complex_matrix_value ()); |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
310 else |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
311 return SparseMatrix (val.matrix_value ()); |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
312 } |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
313 |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
314 /* |
21317
a4faec57f4c8
maint: remove semicolon after %!assert tests to follow Octave conventions.
Rik <rik@octave.org>
parents:
21301
diff
changeset
|
315 %!assert (sparse (2) .^ [3, 4], sparse ([8, 16])) |
22470
9d4cb0cf9cd2
maint: tag xtests and tests that fail on some systems with bug numbers
Mike Miller <mtmiller@octave.org>
parents:
22407
diff
changeset
|
316 %!assert <47775> (sparse (2i) .^ [3, 4], sparse ([-0-8i, 16])) |
31234
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
317 |
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
318 %!test <*63080> |
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
319 %! Z = sparse ([]); |
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
320 %! A = sparse (zeros (0, 2)); |
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
321 %! B = sparse (zeros (2, 0)); |
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
322 %! assert (Z ^ 1, Z); |
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
323 %! assert (Z ^ 0, Z); |
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
324 %! assert (Z ^ -1, Z); |
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
325 %! assert (A ^ 1, Z); |
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
326 %! assert (A ^ 0, Z); |
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
327 %! assert (A ^ -1, Z); |
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
328 %! assert (B ^ 1, Z); |
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
329 %! assert (B ^ 0, Z); |
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
330 %! assert (B ^ -1, Z); |
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
331 |
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
332 %!test <*63080> |
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
333 %! A = sparse (zeros (2, 2)); |
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
334 %! assert (A ^ 1, A); |
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
335 %! assert (A ^ 0, sparse (eye (2, 2))); |
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
336 |
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
337 %!test <63080> |
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
338 %! A = sparse (zeros (2, 2)); |
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
339 %! assert (A ^ -1, sparse (inf (2, 2))); |
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
340 |
15282
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
341 */ |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
342 |
5164 | 343 // -*- 1 -*- |
344 octave_value | |
345 elem_xpow (double a, const SparseMatrix& b) | |
346 { | |
347 octave_value retval; | |
348 | |
5275 | 349 octave_idx_type nr = b.rows (); |
350 octave_idx_type nc = b.cols (); | |
5164 | 351 |
352 double d1, d2; | |
353 | |
354 if (a < 0.0 && ! b.all_integers (d1, d2)) | |
355 { | |
356 Complex atmp (a); | |
357 ComplexMatrix result (nr, nc); | |
358 | |
5275 | 359 for (octave_idx_type j = 0; j < nc; j++) |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
360 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
361 for (octave_idx_type i = 0; i < nr; i++) |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
362 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
363 octave_quit (); |
30390
a61e1a0f6024
maint: style check C++ files in libinterp/ ahead of 7.1 release.
Rik <rik@octave.org>
parents:
29990
diff
changeset
|
364 result(i, j) = std::pow (atmp, b(i, j)); |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
365 } |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
366 } |
5164 | 367 |
368 retval = result; | |
369 } | |
370 else | |
371 { | |
372 Matrix result (nr, nc); | |
373 | |
5275 | 374 for (octave_idx_type j = 0; j < nc; j++) |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
375 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
376 for (octave_idx_type i = 0; i < nr; i++) |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
377 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
378 octave_quit (); |
30390
a61e1a0f6024
maint: style check C++ files in libinterp/ ahead of 7.1 release.
Rik <rik@octave.org>
parents:
29990
diff
changeset
|
379 result(i, j) = std::pow (a, b(i, j)); |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
380 } |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
381 } |
5164 | 382 |
383 retval = result; | |
384 } | |
385 | |
386 return retval; | |
387 } | |
388 | |
389 // -*- 2 -*- | |
390 octave_value | |
391 elem_xpow (double a, const SparseComplexMatrix& b) | |
392 { | |
5275 | 393 octave_idx_type nr = b.rows (); |
394 octave_idx_type nc = b.cols (); | |
5164 | 395 |
396 Complex atmp (a); | |
397 ComplexMatrix result (nr, nc); | |
398 | |
5275 | 399 for (octave_idx_type j = 0; j < nc; j++) |
5164 | 400 { |
5275 | 401 for (octave_idx_type i = 0; i < nr; i++) |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
402 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
403 octave_quit (); |
30390
a61e1a0f6024
maint: style check C++ files in libinterp/ ahead of 7.1 release.
Rik <rik@octave.org>
parents:
29990
diff
changeset
|
404 result(i, j) = std::pow (atmp, b(i, j)); |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
405 } |
5164 | 406 } |
407 | |
408 return result; | |
409 } | |
410 | |
411 // -*- 3 -*- | |
412 octave_value | |
413 elem_xpow (const SparseMatrix& a, double b) | |
414 { | |
17861
870f3e12e163
maint: Use phrase "FIXME:" for problem areas in code.
Rik <rik@octave.org>
parents:
17787
diff
changeset
|
415 // FIXME: What should a .^ 0 give? Matlab gives a |
5164 | 416 // sparse matrix with same structure as a, which is strictly |
21751
b571fc85953f
maint: Use two spaces after period to indicate sentence break.
Rik <rik@octave.org>
parents:
21724
diff
changeset
|
417 // incorrect. Keep compatibility. |
5164 | 418 |
419 octave_value retval; | |
420 | |
10527
b4d2080b6df7
Replace nzmax by nnz as needed
David Bateman <dbateman@free.fr>
parents:
10315
diff
changeset
|
421 octave_idx_type nz = a.nnz (); |
5164 | 422 |
423 if (b <= 0.0) | |
424 { | |
5275 | 425 octave_idx_type nr = a.rows (); |
426 octave_idx_type nc = a.cols (); | |
5164 | 427 |
28144
c31c9eaa1f28
sparse-xpow.cc: use xisint instead of static_cast<int> to check int values
John W. Eaton <jwe@octave.org>
parents:
27923
diff
changeset
|
428 if (! xisint (b) && a.any_element_is_negative ()) |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
429 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
430 ComplexMatrix result (nr, nc, Complex (std::pow (0.0, b))); |
5164 | 431 |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17749
diff
changeset
|
432 // FIXME: avoid apparent GNU libm bug by |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
433 // converting A and B to complex instead of just A. |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
434 Complex btmp (b); |
5164 | 435 |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
436 for (octave_idx_type j = 0; j < nc; j++) |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
437 for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++) |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
438 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
439 octave_quit (); |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
440 |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
441 Complex atmp (a.data (i)); |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
442 |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
443 result(a.ridx (i), j) = std::pow (atmp, btmp); |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
444 } |
5164 | 445 |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
446 retval = octave_value (result); |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
447 } |
5164 | 448 else |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
449 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
450 Matrix result (nr, nc, (std::pow (0.0, b))); |
5164 | 451 |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
452 for (octave_idx_type j = 0; j < nc; j++) |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
453 for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++) |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
454 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
455 octave_quit (); |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
456 result(a.ridx (i), j) = std::pow (a.data (i), b); |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
457 } |
5164 | 458 |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
459 retval = octave_value (result); |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
460 } |
5164 | 461 } |
28144
c31c9eaa1f28
sparse-xpow.cc: use xisint instead of static_cast<int> to check int values
John W. Eaton <jwe@octave.org>
parents:
27923
diff
changeset
|
462 else if (! xisint (b) && a.any_element_is_negative ()) |
5164 | 463 { |
464 SparseComplexMatrix result (a); | |
465 | |
5275 | 466 for (octave_idx_type i = 0; i < nz; i++) |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
467 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
468 octave_quit (); |
5164 | 469 |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17749
diff
changeset
|
470 // FIXME: avoid apparent GNU libm bug by |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
471 // converting A and B to complex instead of just A. |
5164 | 472 |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
473 Complex atmp (a.data (i)); |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
474 Complex btmp (b); |
5164 | 475 |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
476 result.data (i) = std::pow (atmp, btmp); |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
477 } |
5164 | 478 |
479 result.maybe_compress (true); | |
480 | |
481 retval = result; | |
482 } | |
483 else | |
484 { | |
485 SparseMatrix result (a); | |
486 | |
5275 | 487 for (octave_idx_type i = 0; i < nz; i++) |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
488 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
489 octave_quit (); |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
490 result.data (i) = std::pow (a.data (i), b); |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
491 } |
5164 | 492 |
493 result.maybe_compress (true); | |
494 | |
495 retval = result; | |
496 } | |
497 | |
498 return retval; | |
499 } | |
500 | |
501 // -*- 4 -*- | |
502 octave_value | |
503 elem_xpow (const SparseMatrix& a, const SparseMatrix& b) | |
504 { | |
505 octave_value retval; | |
506 | |
5275 | 507 octave_idx_type nr = a.rows (); |
508 octave_idx_type nc = a.cols (); | |
5164 | 509 |
5275 | 510 octave_idx_type b_nr = b.rows (); |
511 octave_idx_type b_nc = b.cols (); | |
5164 | 512 |
15282
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
513 if (a.numel () == 1 && b.numel () > 1) |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
514 return scalar_xpow (a(0), b); |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
515 |
5164 | 516 if (nr != b_nr || nc != b_nc) |
22327
d0562b3159c7
move more classes inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
22323
diff
changeset
|
517 octave::err_nonconformant ("operator .^", nr, nc, b_nr, b_nc); |
5164 | 518 |
519 int convert_to_complex = 0; | |
5275 | 520 for (octave_idx_type j = 0; j < nc; j++) |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
521 for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++) |
5164 | 522 { |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
523 if (a.data(i) < 0.0) |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
524 { |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
525 double btmp = b (a.ridx (i), j); |
28144
c31c9eaa1f28
sparse-xpow.cc: use xisint instead of static_cast<int> to check int values
John W. Eaton <jwe@octave.org>
parents:
27923
diff
changeset
|
526 if (! xisint (btmp)) |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
527 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
528 convert_to_complex = 1; |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
529 goto done; |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
530 } |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
531 } |
5164 | 532 } |
533 | |
534 done: | |
535 | |
5953 | 536 // This is a dumb operator for sparse matrices anyway, and there is |
21751
b571fc85953f
maint: Use two spaces after period to indicate sentence break.
Rik <rik@octave.org>
parents:
21724
diff
changeset
|
537 // no sensible way to handle the 0.^0 versus the 0.^x cases. Therefore |
5953 | 538 // allocate a full matrix filled for the 0.^0 case and shrink it later |
21751
b571fc85953f
maint: Use two spaces after period to indicate sentence break.
Rik <rik@octave.org>
parents:
21724
diff
changeset
|
539 // as needed. |
5164 | 540 |
541 if (convert_to_complex) | |
542 { | |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
543 SparseComplexMatrix complex_result (nr, nc, Complex (1.0, 0.0)); |
5164 | 544 |
5275 | 545 for (octave_idx_type j = 0; j < nc; j++) |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
546 { |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
547 for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++) |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
548 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
549 octave_quit (); |
27277
db687716fed6
style fixes: generally aim to break long lines before operators, not after
John W. Eaton <jwe@octave.org>
parents:
26376
diff
changeset
|
550 complex_result.xelem (a.ridx (i), j) |
db687716fed6
style fixes: generally aim to break long lines before operators, not after
John W. Eaton <jwe@octave.org>
parents:
26376
diff
changeset
|
551 = std::pow (Complex (a.data (i)), Complex (b(a.ridx (i), j))); |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
552 } |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
553 } |
5953 | 554 complex_result.maybe_compress (true); |
5164 | 555 retval = complex_result; |
556 } | |
557 else | |
558 { | |
5953 | 559 SparseMatrix result (nr, nc, 1.0); |
5164 | 560 |
5275 | 561 for (octave_idx_type j = 0; j < nc; j++) |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
562 { |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
563 for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++) |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
564 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
565 octave_quit (); |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
566 result.xelem (a.ridx (i), j) = std::pow (a.data (i), |
31607
aac27ad79be6
maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents:
31605
diff
changeset
|
567 b(a.ridx (i), j)); |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
568 } |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
569 } |
5953 | 570 result.maybe_compress (true); |
5164 | 571 retval = result; |
572 } | |
573 | |
574 return retval; | |
575 } | |
576 | |
577 // -*- 5 -*- | |
578 octave_value | |
579 elem_xpow (const SparseMatrix& a, const Complex& b) | |
580 { | |
581 octave_value retval; | |
582 | |
583 if (b == 0.0) | |
584 // Can this case ever happen, due to automatic retyping with maybe_mutate? | |
585 retval = octave_value (NDArray (a.dims (), 1)); | |
586 else | |
587 { | |
10527
b4d2080b6df7
Replace nzmax by nnz as needed
David Bateman <dbateman@free.fr>
parents:
10315
diff
changeset
|
588 octave_idx_type nz = a.nnz (); |
5164 | 589 SparseComplexMatrix result (a); |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
590 |
5275 | 591 for (octave_idx_type i = 0; i < nz; i++) |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
592 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
593 octave_quit (); |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
594 result.data (i) = std::pow (Complex (a.data (i)), b); |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
595 } |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
596 |
5164 | 597 result.maybe_compress (true); |
598 | |
599 retval = result; | |
600 } | |
601 | |
602 return retval; | |
603 } | |
604 | |
605 // -*- 6 -*- | |
606 octave_value | |
607 elem_xpow (const SparseMatrix& a, const SparseComplexMatrix& b) | |
608 { | |
5275 | 609 octave_idx_type nr = a.rows (); |
610 octave_idx_type nc = a.cols (); | |
5164 | 611 |
5275 | 612 octave_idx_type b_nr = b.rows (); |
613 octave_idx_type b_nc = b.cols (); | |
5164 | 614 |
15282
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
615 if (a.numel () == 1 && b.numel () > 1) |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
616 return scalar_xpow (a(0), b); |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
617 |
5164 | 618 if (nr != b_nr || nc != b_nc) |
22327
d0562b3159c7
move more classes inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
22323
diff
changeset
|
619 octave::err_nonconformant ("operator .^", nr, nc, b_nr, b_nc); |
5164 | 620 |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
621 SparseComplexMatrix result (nr, nc, Complex (1.0, 0.0)); |
5275 | 622 for (octave_idx_type j = 0; j < nc; j++) |
5164 | 623 { |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
624 for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++) |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
625 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
626 octave_quit (); |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
627 result.xelem (a.ridx(i), j) = std::pow (a.data (i), b(a.ridx (i), j)); |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
628 } |
5164 | 629 } |
630 | |
5953 | 631 result.maybe_compress (true); |
5164 | 632 |
633 return result; | |
634 } | |
635 | |
636 // -*- 7 -*- | |
637 octave_value | |
638 elem_xpow (const Complex& a, const SparseMatrix& b) | |
639 { | |
5275 | 640 octave_idx_type nr = b.rows (); |
641 octave_idx_type nc = b.cols (); | |
5164 | 642 |
643 ComplexMatrix result (nr, nc); | |
644 | |
5275 | 645 for (octave_idx_type j = 0; j < nc; j++) |
5164 | 646 { |
5275 | 647 for (octave_idx_type i = 0; i < nr; i++) |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
648 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
649 octave_quit (); |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
650 double btmp = b (i, j); |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
651 if (xisint (btmp)) |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
652 result (i, j) = std::pow (a, static_cast<int> (btmp)); |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
653 else |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
654 result (i, j) = std::pow (a, btmp); |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
655 } |
5164 | 656 } |
657 | |
658 return result; | |
659 } | |
660 | |
661 // -*- 8 -*- | |
662 octave_value | |
663 elem_xpow (const Complex& a, const SparseComplexMatrix& b) | |
664 { | |
5275 | 665 octave_idx_type nr = b.rows (); |
666 octave_idx_type nc = b.cols (); | |
5164 | 667 |
668 ComplexMatrix result (nr, nc); | |
5275 | 669 for (octave_idx_type j = 0; j < nc; j++) |
670 for (octave_idx_type i = 0; i < nr; i++) | |
5164 | 671 { |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
672 octave_quit (); |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
673 result (i, j) = std::pow (a, b (i, j)); |
5164 | 674 } |
675 | |
676 return result; | |
677 } | |
678 | |
679 // -*- 9 -*- | |
680 octave_value | |
681 elem_xpow (const SparseComplexMatrix& a, double b) | |
682 { | |
683 octave_value retval; | |
684 | |
685 if (b <= 0) | |
686 { | |
5275 | 687 octave_idx_type nr = a.rows (); |
688 octave_idx_type nc = a.cols (); | |
5164 | 689 |
5953 | 690 ComplexMatrix result (nr, nc, Complex (std::pow (0.0, b))); |
5164 | 691 |
692 if (xisint (b)) | |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
693 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
694 for (octave_idx_type j = 0; j < nc; j++) |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
695 for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++) |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
696 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
697 octave_quit (); |
27277
db687716fed6
style fixes: generally aim to break long lines before operators, not after
John W. Eaton <jwe@octave.org>
parents:
26376
diff
changeset
|
698 result (a.ridx (i), j) |
db687716fed6
style fixes: generally aim to break long lines before operators, not after
John W. Eaton <jwe@octave.org>
parents:
26376
diff
changeset
|
699 = std::pow (a.data (i), static_cast<int> (b)); |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
700 } |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
701 } |
5164 | 702 else |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
703 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
704 for (octave_idx_type j = 0; j < nc; j++) |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
705 for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++) |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
706 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
707 octave_quit (); |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
708 result (a.ridx (i), j) = std::pow (a.data (i), b); |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
709 } |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
710 } |
5164 | 711 |
712 retval = result; | |
713 } | |
714 else | |
715 { | |
10527
b4d2080b6df7
Replace nzmax by nnz as needed
David Bateman <dbateman@free.fr>
parents:
10315
diff
changeset
|
716 octave_idx_type nz = a.nnz (); |
5164 | 717 |
718 SparseComplexMatrix result (a); | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
719 |
5164 | 720 if (xisint (b)) |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
721 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
722 for (octave_idx_type i = 0; i < nz; i++) |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
723 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
724 octave_quit (); |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
725 result.data (i) = std::pow (a.data (i), static_cast<int> (b)); |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
726 } |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
727 } |
5164 | 728 else |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
729 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
730 for (octave_idx_type i = 0; i < nz; i++) |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
731 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
732 octave_quit (); |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
733 result.data (i) = std::pow (a.data (i), b); |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
734 } |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
735 } |
5164 | 736 |
737 result.maybe_compress (true); | |
738 | |
739 retval = result; | |
740 } | |
741 | |
742 return retval; | |
743 } | |
744 | |
745 // -*- 10 -*- | |
746 octave_value | |
747 elem_xpow (const SparseComplexMatrix& a, const SparseMatrix& b) | |
748 { | |
5275 | 749 octave_idx_type nr = a.rows (); |
750 octave_idx_type nc = a.cols (); | |
5164 | 751 |
5275 | 752 octave_idx_type b_nr = b.rows (); |
753 octave_idx_type b_nc = b.cols (); | |
5164 | 754 |
15282
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
755 if (a.numel () == 1 && b.numel () > 1) |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
756 return scalar_xpow (a(0), b); |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
757 |
5164 | 758 if (nr != b_nr || nc != b_nc) |
22327
d0562b3159c7
move more classes inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
22323
diff
changeset
|
759 octave::err_nonconformant ("operator .^", nr, nc, b_nr, b_nc); |
5164 | 760 |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
761 SparseComplexMatrix result (nr, nc, Complex (1.0, 0.0)); |
5275 | 762 for (octave_idx_type j = 0; j < nc; j++) |
5164 | 763 { |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
764 for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++) |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
765 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
766 octave_quit (); |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
767 double btmp = b(a.ridx (i), j); |
5164 | 768 |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
769 if (xisint (btmp)) |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
770 result.xelem (a.ridx (i), j) = std::pow (a.data (i), |
31607
aac27ad79be6
maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents:
31605
diff
changeset
|
771 static_cast<int> (btmp)); |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
772 else |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
773 result.xelem (a.ridx (i), j) = std::pow (a.data (i), btmp); |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
774 } |
5164 | 775 } |
776 | |
5953 | 777 result.maybe_compress (true); |
5164 | 778 |
779 return result; | |
780 } | |
781 | |
782 // -*- 11 -*- | |
783 octave_value | |
784 elem_xpow (const SparseComplexMatrix& a, const Complex& b) | |
785 { | |
786 octave_value retval; | |
787 | |
788 if (b == 0.0) | |
789 // Can this case ever happen, due to automatic retyping with maybe_mutate? | |
790 retval = octave_value (NDArray (a.dims (), 1)); | |
791 else | |
792 { | |
793 | |
10527
b4d2080b6df7
Replace nzmax by nnz as needed
David Bateman <dbateman@free.fr>
parents:
10315
diff
changeset
|
794 octave_idx_type nz = a.nnz (); |
5164 | 795 |
796 SparseComplexMatrix result (a); | |
797 | |
5275 | 798 for (octave_idx_type i = 0; i < nz; i++) |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
799 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
800 octave_quit (); |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
801 result.data (i) = std::pow (a.data (i), b); |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
802 } |
5164 | 803 |
804 result.maybe_compress (true); | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
805 |
5164 | 806 retval = result; |
807 } | |
808 | |
809 return retval; | |
810 } | |
811 | |
812 // -*- 12 -*- | |
813 octave_value | |
814 elem_xpow (const SparseComplexMatrix& a, const SparseComplexMatrix& b) | |
815 { | |
5275 | 816 octave_idx_type nr = a.rows (); |
817 octave_idx_type nc = a.cols (); | |
5164 | 818 |
5275 | 819 octave_idx_type b_nr = b.rows (); |
820 octave_idx_type b_nc = b.cols (); | |
5164 | 821 |
15282
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
822 if (a.numel () == 1 && b.numel () > 1) |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
823 return scalar_xpow (a(0), b); |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
824 |
5164 | 825 if (nr != b_nr || nc != b_nc) |
22327
d0562b3159c7
move more classes inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
22323
diff
changeset
|
826 octave::err_nonconformant ("operator .^", nr, nc, b_nr, b_nc); |
5164 | 827 |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
828 SparseComplexMatrix result (nr, nc, Complex (1.0, 0.0)); |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
829 for (octave_idx_type j = 0; j < nc; j++) |
5164 | 830 { |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
831 for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++) |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
832 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
833 octave_quit (); |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17749
diff
changeset
|
834 result.xelem (a.ridx (i), j) = std::pow (a.data (i), |
31607
aac27ad79be6
maint: Re-indent code after switch to using namespace macros.
Rik <rik@octave.org>
parents:
31605
diff
changeset
|
835 b(a.ridx (i), j)); |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
836 } |
5164 | 837 } |
838 result.maybe_compress (true); | |
839 | |
840 return result; | |
841 } | |
29990
b839c36fd106
move sparse xdiv and xpow operator functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
29358
diff
changeset
|
842 |
31234
a1318deb4584
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31231
diff
changeset
|
843 |
31605
e88a07dec498
maint: Use macros to begin/end C++ namespaces.
Rik <rik@octave.org>
parents:
31235
diff
changeset
|
844 OCTAVE_END_NAMESPACE(octave) |