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