Mercurial > octave
annotate libinterp/corefcn/sparse-xpow.cc @ 31231:a026fb2be108
sparse-xpow.cc: Return empty matrix for empty input (bug #63080)
author | Arun Giridhar <arungiridhar@gmail.com> |
---|---|
date | Mon, 19 Sep 2022 07:05:31 -0400 |
parents | 0dec459a4064 |
children | a1318deb4584 |
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 // |
30564
796f54d4ddbf
update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents:
30390
diff
changeset
|
3 // Copyright (C) 1998-2022 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 | |
30 #include <cassert> | |
15215
9020dddc925a
use std::numeric_limits for integer max and min values
John W. Eaton <jwe@octave.org>
parents:
15195
diff
changeset
|
31 |
9020dddc925a
use std::numeric_limits for integer max and min values
John W. Eaton <jwe@octave.org>
parents:
15195
diff
changeset
|
32 #include <limits> |
5164 | 33 |
34 #include "Array-util.h" | |
35 #include "oct-cmplx.h" | |
36 #include "quit.h" | |
37 | |
38 #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
|
39 #include "ovl.h" |
5164 | 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 | |
29990
b839c36fd106
move sparse xdiv and xpow operator functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
29358
diff
changeset
|
48 OCTAVE_NAMESPACE_BEGIN |
b839c36fd106
move sparse xdiv and xpow operator functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
29358
diff
changeset
|
49 |
31226
44cf6bbeeca9
sparse-xpow.cc: Change return type from int to bool for internal function
Arun Giridhar <arungiridhar@gmail.com>
parents:
31225
diff
changeset
|
50 static inline bool |
5164 | 51 xisint (double x) |
52 { | |
21782
2aef506f3fec
use namespace for lo-mappers.h functions
John W. Eaton <jwe@octave.org>
parents:
21751
diff
changeset
|
53 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
|
54 && ((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
|
55 || (x <= 0 && x > std::numeric_limits<int>::min ()))); |
5164 | 56 } |
57 | |
21751
b571fc85953f
maint: Use two spaces after period to indicate sentence break.
Rik <rik@octave.org>
parents:
21724
diff
changeset
|
58 // Safer pow functions. Only two make sense for sparse matrices, the |
5164 | 59 // others should all promote to full matrices. |
60 | |
61 octave_value | |
62 xpow (const SparseMatrix& a, double b) | |
63 { | |
64 octave_value retval; | |
65 | |
5275 | 66 octave_idx_type nr = a.rows (); |
67 octave_idx_type nc = a.cols (); | |
5164 | 68 |
31231
a026fb2be108
sparse-xpow.cc: Return empty matrix for empty input (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31227
diff
changeset
|
69 if (nr != nc) |
21751
b571fc85953f
maint: Use two spaces after period to indicate sentence break.
Rik <rik@octave.org>
parents:
21724
diff
changeset
|
70 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
|
71 |
31231
a026fb2be108
sparse-xpow.cc: Return empty matrix for empty input (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31227
diff
changeset
|
72 if (nr == 0 && nc == 0) |
a026fb2be108
sparse-xpow.cc: Return empty matrix for empty input (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31227
diff
changeset
|
73 return a; |
a026fb2be108
sparse-xpow.cc: Return empty matrix for empty input (bug #63080)
Arun Giridhar <arungiridhar@gmail.com>
parents:
31227
diff
changeset
|
74 |
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
|
75 if (! xisint (b)) |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
76 error ("use full(a) ^ full(b)"); |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
77 |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
78 int btmp = static_cast<int> (b); |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
79 if (btmp == 0) |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
80 { |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
81 SparseMatrix tmp = SparseMatrix (nr, nr, nr); |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
82 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
|
83 { |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
84 tmp.data (i) = 1.0; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
85 tmp.ridx (i) = i; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
86 } |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
87 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
|
88 tmp.cidx (i) = i; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
89 |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
90 retval = tmp; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
91 } |
5164 | 92 else |
93 { | |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
94 SparseMatrix atmp; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
95 if (btmp < 0) |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
96 { |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
97 btmp = -btmp; |
5164 | 98 |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
99 octave_idx_type info; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
100 double rcond = 0.0; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
101 MatrixType mattyp (a); |
5164 | 102 |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
103 atmp = a.inverse (mattyp, info, rcond, 1); |
5164 | 104 |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
105 if (info == -1) |
22197
e43d83253e28
refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents:
22022
diff
changeset
|
106 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
|
107 } |
5164 | 108 else |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
109 atmp = a; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
110 |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
111 SparseMatrix result (atmp); |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
112 |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
113 btmp--; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
114 |
31224
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
115 // 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
|
116 // 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
|
117 // 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
|
118 // 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
|
119 // Linear multiplication uses a linear number of multiplications |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
120 // but one of the matrices it uses will be as sparse as the original matrix. |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
121 // |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
122 // 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
|
123 // 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
|
124 // 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
|
125 // 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
|
126 // |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
127 // 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
|
128 // |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
129 // Large exponents favor the squaring technique, and sparse matrices favor |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
130 // linear multiplication. |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
131 // |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
132 // 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
|
133 // 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
|
134 // |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
135 // FIXME: Improve this threshold calculation. |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
136 |
31224
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
137 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
|
138 int threshold = (sparsity >= 1000) ? 40 |
0dec459a4064
sparse-xpow.cc: Performance tweak for threshold selection
Arun Giridhar <arungiridhar@gmail.com>
parents:
31226
diff
changeset
|
139 : (sparsity >= 100) ? 20 |
31224
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
140 : 3; |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
141 |
31224
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
142 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
|
143 { |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
144 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
|
145 { |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
146 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
|
147 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
|
148 |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
149 btmp >>= 1; |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
150 |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
151 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
|
152 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
|
153 } |
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 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
|
156 { |
45984c799215
sparse-xpow.cc: Use faster multiplication technique based on input matrix sparsity
Arun Giridhar <arungiridhar@gmail.com>
parents:
30564
diff
changeset
|
157 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
|
158 result = result * atmp; |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
159 } |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
160 |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
161 retval = result; |
5164 | 162 } |
163 | |
164 return retval; | |
165 } | |
166 | |
167 octave_value | |
168 xpow (const SparseComplexMatrix& a, double b) | |
169 { | |
170 octave_value retval; | |
171 | |
5275 | 172 octave_idx_type nr = a.rows (); |
173 octave_idx_type nc = a.cols (); | |
5164 | 174 |
175 if (nr == 0 || nc == 0 || nr != nc) | |
21751
b571fc85953f
maint: Use two spaces after period to indicate sentence break.
Rik <rik@octave.org>
parents:
21724
diff
changeset
|
176 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
|
177 |
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
|
178 if (! xisint (b)) |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
179 error ("use full(a) ^ full(b)"); |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
180 |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
181 int btmp = static_cast<int> (b); |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
182 if (btmp == 0) |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
183 { |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
184 SparseMatrix tmp = SparseMatrix (nr, nr, nr); |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
185 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
|
186 { |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
187 tmp.data (i) = 1.0; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
188 tmp.ridx (i) = i; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
189 } |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
190 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
|
191 tmp.cidx (i) = i; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
192 |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
193 retval = tmp; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
194 } |
5164 | 195 else |
196 { | |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
197 SparseComplexMatrix atmp; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
198 if (btmp < 0) |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
199 { |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
200 btmp = -btmp; |
5164 | 201 |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
202 octave_idx_type info; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
203 double rcond = 0.0; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
204 MatrixType mattyp (a); |
5164 | 205 |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
206 atmp = a.inverse (mattyp, info, rcond, 1); |
5164 | 207 |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
208 if (info == -1) |
22197
e43d83253e28
refill multi-line macro definitions
John W. Eaton <jwe@octave.org>
parents:
22022
diff
changeset
|
209 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
|
210 } |
5164 | 211 else |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
212 atmp = a; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
213 |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
214 SparseComplexMatrix result (atmp); |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
215 |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
216 btmp--; |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
217 |
31225
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
218 // 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
|
219 // 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
|
220 // for more details. |
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
221 // |
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
222 // 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
|
223 |
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
224 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
|
225 int threshold = (sparsity >= 1000) ? 40 |
0dec459a4064
sparse-xpow.cc: Performance tweak for threshold selection
Arun Giridhar <arungiridhar@gmail.com>
parents:
31226
diff
changeset
|
226 : (sparsity >= 100) ? 20 |
31225
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
227 : 3; |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
228 |
31225
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
229 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
|
230 { |
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
231 while (btmp > 0) |
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
232 { |
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
233 if (btmp & 1) |
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
234 result = result * atmp; |
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
235 |
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
236 btmp >>= 1; |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
237 |
31225
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
238 if (btmp > 0) |
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
239 atmp = atmp * atmp; |
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
240 } |
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
241 } |
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
242 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
|
243 { |
3eab70385569
sparse-xpow.cc: Use faster multiplication technique, this time for complex
Arun Giridhar <arungiridhar@gmail.com>
parents:
31224
diff
changeset
|
244 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
|
245 result = result * atmp; |
20982
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
246 } |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
247 |
d27f66b4b8e6
maint: invert if/else/error instances.
John W. Eaton <jwe@octave.org>
parents:
20940
diff
changeset
|
248 retval = result; |
5164 | 249 } |
250 | |
251 return retval; | |
252 } | |
253 | |
254 // Safer pow functions that work elementwise for matrices. | |
255 // | |
256 // op2 \ op1: s m cs cm | |
257 // +-- +---+---+----+----+ | |
258 // scalar | | * | 3 | * | 9 | | |
259 // +---+---+----+----+ | |
260 // matrix | 1 | 4 | 7 | 10 | | |
261 // +---+---+----+----+ | |
262 // complex_scalar | * | 5 | * | 11 | | |
263 // +---+---+----+----+ | |
264 // complex_matrix | 2 | 6 | 8 | 12 | | |
265 // +---+---+----+----+ | |
266 // | |
267 // * -> not needed. | |
268 | |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17749
diff
changeset
|
269 // FIXME: these functions need to be fixed so that things |
5164 | 270 // like |
271 // | |
272 // a = -1; b = [ 0, 0.5, 1 ]; r = a .^ b | |
273 // | |
274 // and | |
275 // | |
276 // a = -1; b = [ 0, 0.5, 1 ]; for i = 1:3, r(i) = a .^ b(i), end | |
277 // | |
278 // produce identical results. Also, it would be nice if -1^0.5 | |
279 // produced a pure imaginary result instead of a complex number with a | |
280 // small real part. But perhaps that's really a problem with the math | |
281 // library... | |
282 | |
15282
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
283 // 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
|
284 // 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
|
285 // 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
|
286 // 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
|
287 // in practice? |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
288 |
21139
538b57866b90
consistently use "typename" intead of "class" in template declarations
John W. Eaton <jwe@octave.org>
parents:
21121
diff
changeset
|
289 template <typename S, typename SM> |
15282
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
290 inline octave_value |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
291 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
|
292 { |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
293 octave_value val = elem_xpow (a, b); |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
294 |
23581
c3075ae020e1
maint: Deprecate is_complex_type and replace with iscomplex.
Rik <rik@octave.org>
parents:
23220
diff
changeset
|
295 if (val.iscomplex ()) |
15282
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
296 return SparseComplexMatrix (val.complex_matrix_value ()); |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
297 else |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
298 return SparseMatrix (val.matrix_value ()); |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
299 } |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
300 |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
301 /* |
21317
a4faec57f4c8
maint: remove semicolon after %!assert tests to follow Octave conventions.
Rik <rik@octave.org>
parents:
21301
diff
changeset
|
302 %!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
|
303 %!assert <47775> (sparse (2i) .^ [3, 4], sparse ([-0-8i, 16])) |
15282
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
304 */ |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
305 |
5164 | 306 // -*- 1 -*- |
307 octave_value | |
308 elem_xpow (double a, const SparseMatrix& b) | |
309 { | |
310 octave_value retval; | |
311 | |
5275 | 312 octave_idx_type nr = b.rows (); |
313 octave_idx_type nc = b.cols (); | |
5164 | 314 |
315 double d1, d2; | |
316 | |
317 if (a < 0.0 && ! b.all_integers (d1, d2)) | |
318 { | |
319 Complex atmp (a); | |
320 ComplexMatrix result (nr, nc); | |
321 | |
5275 | 322 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
|
323 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
324 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
|
325 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
326 octave_quit (); |
30390
a61e1a0f6024
maint: style check C++ files in libinterp/ ahead of 7.1 release.
Rik <rik@octave.org>
parents:
29990
diff
changeset
|
327 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
|
328 } |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
329 } |
5164 | 330 |
331 retval = result; | |
332 } | |
333 else | |
334 { | |
335 Matrix result (nr, nc); | |
336 | |
5275 | 337 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
|
338 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
339 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
|
340 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
341 octave_quit (); |
30390
a61e1a0f6024
maint: style check C++ files in libinterp/ ahead of 7.1 release.
Rik <rik@octave.org>
parents:
29990
diff
changeset
|
342 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
|
343 } |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
344 } |
5164 | 345 |
346 retval = result; | |
347 } | |
348 | |
349 return retval; | |
350 } | |
351 | |
352 // -*- 2 -*- | |
353 octave_value | |
354 elem_xpow (double a, const SparseComplexMatrix& b) | |
355 { | |
5275 | 356 octave_idx_type nr = b.rows (); |
357 octave_idx_type nc = b.cols (); | |
5164 | 358 |
359 Complex atmp (a); | |
360 ComplexMatrix result (nr, nc); | |
361 | |
5275 | 362 for (octave_idx_type j = 0; j < nc; j++) |
5164 | 363 { |
5275 | 364 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
|
365 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
366 octave_quit (); |
30390
a61e1a0f6024
maint: style check C++ files in libinterp/ ahead of 7.1 release.
Rik <rik@octave.org>
parents:
29990
diff
changeset
|
367 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
|
368 } |
5164 | 369 } |
370 | |
371 return result; | |
372 } | |
373 | |
374 // -*- 3 -*- | |
375 octave_value | |
376 elem_xpow (const SparseMatrix& a, double b) | |
377 { | |
17861
870f3e12e163
maint: Use phrase "FIXME:" for problem areas in code.
Rik <rik@octave.org>
parents:
17787
diff
changeset
|
378 // FIXME: What should a .^ 0 give? Matlab gives a |
5164 | 379 // 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
|
380 // incorrect. Keep compatibility. |
5164 | 381 |
382 octave_value retval; | |
383 | |
10527
b4d2080b6df7
Replace nzmax by nnz as needed
David Bateman <dbateman@free.fr>
parents:
10315
diff
changeset
|
384 octave_idx_type nz = a.nnz (); |
5164 | 385 |
386 if (b <= 0.0) | |
387 { | |
5275 | 388 octave_idx_type nr = a.rows (); |
389 octave_idx_type nc = a.cols (); | |
5164 | 390 |
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
|
391 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
|
392 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
393 ComplexMatrix result (nr, nc, Complex (std::pow (0.0, b))); |
5164 | 394 |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17749
diff
changeset
|
395 // FIXME: avoid apparent GNU libm bug by |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
396 // 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
|
397 Complex btmp (b); |
5164 | 398 |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
399 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
|
400 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
|
401 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
402 octave_quit (); |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
403 |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
404 Complex atmp (a.data (i)); |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
405 |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
406 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
|
407 } |
5164 | 408 |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
409 retval = octave_value (result); |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
410 } |
5164 | 411 else |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
412 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
413 Matrix result (nr, nc, (std::pow (0.0, b))); |
5164 | 414 |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
415 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
|
416 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
|
417 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
418 octave_quit (); |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
419 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
|
420 } |
5164 | 421 |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
422 retval = octave_value (result); |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
423 } |
5164 | 424 } |
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
|
425 else if (! xisint (b) && a.any_element_is_negative ()) |
5164 | 426 { |
427 SparseComplexMatrix result (a); | |
428 | |
5275 | 429 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
|
430 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
431 octave_quit (); |
5164 | 432 |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17749
diff
changeset
|
433 // FIXME: avoid apparent GNU libm bug by |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
434 // converting A and B to complex instead of just A. |
5164 | 435 |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
436 Complex atmp (a.data (i)); |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
437 Complex btmp (b); |
5164 | 438 |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
439 result.data (i) = std::pow (atmp, btmp); |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
440 } |
5164 | 441 |
442 result.maybe_compress (true); | |
443 | |
444 retval = result; | |
445 } | |
446 else | |
447 { | |
448 SparseMatrix result (a); | |
449 | |
5275 | 450 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
|
451 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
452 octave_quit (); |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
453 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
|
454 } |
5164 | 455 |
456 result.maybe_compress (true); | |
457 | |
458 retval = result; | |
459 } | |
460 | |
461 return retval; | |
462 } | |
463 | |
464 // -*- 4 -*- | |
465 octave_value | |
466 elem_xpow (const SparseMatrix& a, const SparseMatrix& b) | |
467 { | |
468 octave_value retval; | |
469 | |
5275 | 470 octave_idx_type nr = a.rows (); |
471 octave_idx_type nc = a.cols (); | |
5164 | 472 |
5275 | 473 octave_idx_type b_nr = b.rows (); |
474 octave_idx_type b_nc = b.cols (); | |
5164 | 475 |
15282
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
476 if (a.numel () == 1 && b.numel () > 1) |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
477 return scalar_xpow (a(0), b); |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
478 |
5164 | 479 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
|
480 octave::err_nonconformant ("operator .^", nr, nc, b_nr, b_nc); |
5164 | 481 |
482 int convert_to_complex = 0; | |
5275 | 483 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
|
484 for (octave_idx_type i = a.cidx (j); i < a.cidx (j+1); i++) |
5164 | 485 { |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
486 if (a.data(i) < 0.0) |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
487 { |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
488 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
|
489 if (! xisint (btmp)) |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
490 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
491 convert_to_complex = 1; |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
492 goto done; |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
493 } |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
494 } |
5164 | 495 } |
496 | |
497 done: | |
498 | |
5953 | 499 // 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
|
500 // no sensible way to handle the 0.^0 versus the 0.^x cases. Therefore |
5953 | 501 // 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
|
502 // as needed. |
5164 | 503 |
504 if (convert_to_complex) | |
505 { | |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
506 SparseComplexMatrix complex_result (nr, nc, Complex (1.0, 0.0)); |
5164 | 507 |
5275 | 508 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
|
509 { |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
510 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
|
511 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
512 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
|
513 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
|
514 = 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
|
515 } |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
516 } |
5953 | 517 complex_result.maybe_compress (true); |
5164 | 518 retval = complex_result; |
519 } | |
520 else | |
521 { | |
5953 | 522 SparseMatrix result (nr, nc, 1.0); |
5164 | 523 |
5275 | 524 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
|
525 { |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
526 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
|
527 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
528 octave_quit (); |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
529 result.xelem (a.ridx (i), j) = std::pow (a.data (i), |
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
530 b(a.ridx (i), j)); |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
531 } |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
532 } |
5953 | 533 result.maybe_compress (true); |
5164 | 534 retval = result; |
535 } | |
536 | |
537 return retval; | |
538 } | |
539 | |
540 // -*- 5 -*- | |
541 octave_value | |
542 elem_xpow (const SparseMatrix& a, const Complex& b) | |
543 { | |
544 octave_value retval; | |
545 | |
546 if (b == 0.0) | |
547 // Can this case ever happen, due to automatic retyping with maybe_mutate? | |
548 retval = octave_value (NDArray (a.dims (), 1)); | |
549 else | |
550 { | |
10527
b4d2080b6df7
Replace nzmax by nnz as needed
David Bateman <dbateman@free.fr>
parents:
10315
diff
changeset
|
551 octave_idx_type nz = a.nnz (); |
5164 | 552 SparseComplexMatrix result (a); |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
553 |
5275 | 554 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
|
555 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
556 octave_quit (); |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
557 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
|
558 } |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
559 |
5164 | 560 result.maybe_compress (true); |
561 | |
562 retval = result; | |
563 } | |
564 | |
565 return retval; | |
566 } | |
567 | |
568 // -*- 6 -*- | |
569 octave_value | |
570 elem_xpow (const SparseMatrix& a, const SparseComplexMatrix& b) | |
571 { | |
5275 | 572 octave_idx_type nr = a.rows (); |
573 octave_idx_type nc = a.cols (); | |
5164 | 574 |
5275 | 575 octave_idx_type b_nr = b.rows (); |
576 octave_idx_type b_nc = b.cols (); | |
5164 | 577 |
15282
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
578 if (a.numel () == 1 && b.numel () > 1) |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
579 return scalar_xpow (a(0), b); |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
580 |
5164 | 581 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
|
582 octave::err_nonconformant ("operator .^", nr, nc, b_nr, b_nc); |
5164 | 583 |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
584 SparseComplexMatrix result (nr, nc, Complex (1.0, 0.0)); |
5275 | 585 for (octave_idx_type j = 0; j < nc; j++) |
5164 | 586 { |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
587 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
|
588 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
589 octave_quit (); |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
590 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
|
591 } |
5164 | 592 } |
593 | |
5953 | 594 result.maybe_compress (true); |
5164 | 595 |
596 return result; | |
597 } | |
598 | |
599 // -*- 7 -*- | |
600 octave_value | |
601 elem_xpow (const Complex& a, const SparseMatrix& b) | |
602 { | |
5275 | 603 octave_idx_type nr = b.rows (); |
604 octave_idx_type nc = b.cols (); | |
5164 | 605 |
606 ComplexMatrix result (nr, nc); | |
607 | |
5275 | 608 for (octave_idx_type j = 0; j < nc; j++) |
5164 | 609 { |
5275 | 610 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
|
611 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
612 octave_quit (); |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
613 double btmp = b (i, j); |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
614 if (xisint (btmp)) |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
615 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
|
616 else |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
617 result (i, j) = std::pow (a, btmp); |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
618 } |
5164 | 619 } |
620 | |
621 return result; | |
622 } | |
623 | |
624 // -*- 8 -*- | |
625 octave_value | |
626 elem_xpow (const Complex& a, const SparseComplexMatrix& b) | |
627 { | |
5275 | 628 octave_idx_type nr = b.rows (); |
629 octave_idx_type nc = b.cols (); | |
5164 | 630 |
631 ComplexMatrix result (nr, nc); | |
5275 | 632 for (octave_idx_type j = 0; j < nc; j++) |
633 for (octave_idx_type i = 0; i < nr; i++) | |
5164 | 634 { |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
635 octave_quit (); |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
636 result (i, j) = std::pow (a, b (i, j)); |
5164 | 637 } |
638 | |
639 return result; | |
640 } | |
641 | |
642 // -*- 9 -*- | |
643 octave_value | |
644 elem_xpow (const SparseComplexMatrix& a, double b) | |
645 { | |
646 octave_value retval; | |
647 | |
648 if (b <= 0) | |
649 { | |
5275 | 650 octave_idx_type nr = a.rows (); |
651 octave_idx_type nc = a.cols (); | |
5164 | 652 |
5953 | 653 ComplexMatrix result (nr, nc, Complex (std::pow (0.0, b))); |
5164 | 654 |
655 if (xisint (b)) | |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
656 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
657 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
|
658 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
|
659 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
660 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
|
661 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
|
662 = 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
|
663 } |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
664 } |
5164 | 665 else |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
666 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
667 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
|
668 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
|
669 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
670 octave_quit (); |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
671 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
|
672 } |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
673 } |
5164 | 674 |
675 retval = result; | |
676 } | |
677 else | |
678 { | |
10527
b4d2080b6df7
Replace nzmax by nnz as needed
David Bateman <dbateman@free.fr>
parents:
10315
diff
changeset
|
679 octave_idx_type nz = a.nnz (); |
5164 | 680 |
681 SparseComplexMatrix result (a); | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
682 |
5164 | 683 if (xisint (b)) |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
684 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
685 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
|
686 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
687 octave_quit (); |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
688 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
|
689 } |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
690 } |
5164 | 691 else |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
692 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
693 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
|
694 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
695 octave_quit (); |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
696 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
|
697 } |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
698 } |
5164 | 699 |
700 result.maybe_compress (true); | |
701 | |
702 retval = result; | |
703 } | |
704 | |
705 return retval; | |
706 } | |
707 | |
708 // -*- 10 -*- | |
709 octave_value | |
710 elem_xpow (const SparseComplexMatrix& a, const SparseMatrix& b) | |
711 { | |
5275 | 712 octave_idx_type nr = a.rows (); |
713 octave_idx_type nc = a.cols (); | |
5164 | 714 |
5275 | 715 octave_idx_type b_nr = b.rows (); |
716 octave_idx_type b_nc = b.cols (); | |
5164 | 717 |
15282
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
718 if (a.numel () == 1 && b.numel () > 1) |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
719 return scalar_xpow (a(0), b); |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
720 |
5164 | 721 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
|
722 octave::err_nonconformant ("operator .^", nr, nc, b_nr, b_nc); |
5164 | 723 |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
724 SparseComplexMatrix result (nr, nc, Complex (1.0, 0.0)); |
5275 | 725 for (octave_idx_type j = 0; j < nc; j++) |
5164 | 726 { |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
727 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
|
728 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
729 octave_quit (); |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
730 double btmp = b(a.ridx (i), j); |
5164 | 731 |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
732 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
|
733 result.xelem (a.ridx (i), j) = std::pow (a.data (i), |
25103
078b795c5219
maint: style check C++ ahead of 4.4 release.
Rik <rik@octave.org>
parents:
25054
diff
changeset
|
734 static_cast<int> (btmp)); |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
735 else |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
736 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
|
737 } |
5164 | 738 } |
739 | |
5953 | 740 result.maybe_compress (true); |
5164 | 741 |
742 return result; | |
743 } | |
744 | |
745 // -*- 11 -*- | |
746 octave_value | |
747 elem_xpow (const SparseComplexMatrix& a, const Complex& b) | |
748 { | |
749 octave_value retval; | |
750 | |
751 if (b == 0.0) | |
752 // Can this case ever happen, due to automatic retyping with maybe_mutate? | |
753 retval = octave_value (NDArray (a.dims (), 1)); | |
754 else | |
755 { | |
756 | |
10527
b4d2080b6df7
Replace nzmax by nnz as needed
David Bateman <dbateman@free.fr>
parents:
10315
diff
changeset
|
757 octave_idx_type nz = a.nnz (); |
5164 | 758 |
759 SparseComplexMatrix result (a); | |
760 | |
5275 | 761 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
|
762 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
763 octave_quit (); |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
764 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
|
765 } |
5164 | 766 |
767 result.maybe_compress (true); | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
768 |
5164 | 769 retval = result; |
770 } | |
771 | |
772 return retval; | |
773 } | |
774 | |
775 // -*- 12 -*- | |
776 octave_value | |
777 elem_xpow (const SparseComplexMatrix& a, const SparseComplexMatrix& b) | |
778 { | |
5275 | 779 octave_idx_type nr = a.rows (); |
780 octave_idx_type nc = a.cols (); | |
5164 | 781 |
5275 | 782 octave_idx_type b_nr = b.rows (); |
783 octave_idx_type b_nc = b.cols (); | |
5164 | 784 |
15282
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
785 if (a.numel () == 1 && b.numel () > 1) |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
786 return scalar_xpow (a(0), b); |
06ce57277bfb
handle scalar-sparse-matrix .^ matrix ops
John W. Eaton <jwe@octave.org>
parents:
14138
diff
changeset
|
787 |
5164 | 788 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
|
789 octave::err_nonconformant ("operator .^", nr, nc, b_nr, b_nc); |
5164 | 790 |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
791 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
|
792 for (octave_idx_type j = 0; j < nc; j++) |
5164 | 793 { |
14861
f7afecdd87ef
maint: Use Octave coding conventions for cuddling parentheses in src/ directory
Rik <octave@nomad.inbox5.com>
parents:
14138
diff
changeset
|
794 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
|
795 { |
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
796 octave_quit (); |
17787
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17749
diff
changeset
|
797 result.xelem (a.ridx (i), j) = std::pow (a.data (i), |
175b392e91fe
Use GNU style coding conventions for code in libinterp/
Rik <rik@octave.org>
parents:
17749
diff
changeset
|
798 b(a.ridx (i), j)); |
10315
57a59eae83cc
untabify src C++ source files
John W. Eaton <jwe@octave.org>
parents:
10160
diff
changeset
|
799 } |
5164 | 800 } |
801 result.maybe_compress (true); | |
802 | |
803 return result; | |
804 } | |
29990
b839c36fd106
move sparse xdiv and xpow operator functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
29358
diff
changeset
|
805 |
b839c36fd106
move sparse xdiv and xpow operator functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents:
29358
diff
changeset
|
806 OCTAVE_NAMESPACE_END |