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