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