Mercurial > octave

comparison libinterp/corefcn/sparse-xpow.cc @ 31234:a1318deb4584

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
sparse-xpow.cc: Improve consistency between sparse and full matrices (bug #63080) sparse-xpow.cc: Return 0x0 matrix for empty inputs of all sizes. Return early for all-zeros matrix with positive power. Add FIXME for all-zeros matrix with negative power. Add BISTs.
author Arun Giridhar <arungiridhar@gmail.com>
date Tue, 20 Sep 2022 09:00:02 -0400
parents a026fb2be108
children b542b88ad3b6
comparison
equal deleted inserted replaced
31233:95fa3182f00c 31234:a1318deb4584
64 octave_value retval; 64 octave_value retval;
65 65
66 octave_idx_type nr = a.rows (); 66 octave_idx_type nr = a.rows ();
67 octave_idx_type nc = a.cols (); 67 octave_idx_type nc = a.cols ();
68 68
69 if (nr == 0 || nc == 0)
70 return SparseMatrix();
71
72 // If we are here, A is not empty ==> A needs to be square.
69 if (nr != nc) 73 if (nr != nc)
70 error ("for A^b, A must be a square matrix. Use .^ for elementwise power."); 74 error ("for A^b, A must be a square matrix. Use .^ for elementwise power.");
71
72 if (nr == 0 && nc == 0)
73 return a;
74 75
75 if (! xisint (b)) 76 if (! xisint (b))
76 error ("use full(a) ^ full(b)"); 77 error ("use full(a) ^ full(b)");
77 78
78 int btmp = static_cast<int> (b); 79 int btmp = static_cast<int> (b);
98 99
99 octave_idx_type info; 100 octave_idx_type info;
100 double rcond = 0.0; 101 double rcond = 0.0;
101 MatrixType mattyp (a); 102 MatrixType mattyp (a);
102 103
104 // FIXME: This causes an error if the input sparse matrix is all-zeros.
105 // That behavior is inconsistent with A ^ b when A is a full all-zeros
106 // matrix, which just returns Inf of the same size with a warning.
103 atmp = a.inverse (mattyp, info, rcond, 1); 107 atmp = a.inverse (mattyp, info, rcond, 1);
104 108
105 if (info == -1) 109 if (info == -1)
106 warning ("inverse: matrix singular to machine precision, rcond = %g", rcond); 110 warning ("inverse: matrix singular to machine precision, rcond = %g", rcond);
107 } 111 }
108 else 112 else
109 atmp = a; 113 atmp = a;
114
115 if (atmp.nnz() == 0) // Fast return for all-zeros matrix
116 return atmp;
110 117
111 SparseMatrix result (atmp); 118 SparseMatrix result (atmp);
112 119
113 btmp--; 120 btmp--;
114 121
170 octave_value retval; 177 octave_value retval;
171 178
172 octave_idx_type nr = a.rows (); 179 octave_idx_type nr = a.rows ();
173 octave_idx_type nc = a.cols (); 180 octave_idx_type nc = a.cols ();
174 181
175 if (nr == 0 || nc == 0 || nr != nc) 182 if (nr == 0 || nc == 0)
183 return SparseMatrix();
184
185 // If we are here, A is not empty ==> A needs to be square.
186 if (nr != nc)
176 error ("for A^b, A must be a square matrix. Use .^ for elementwise power."); 187 error ("for A^b, A must be a square matrix. Use .^ for elementwise power.");
177 188
178 if (! xisint (b)) 189 if (! xisint (b))
179 error ("use full(a) ^ full(b)"); 190 error ("use full(a) ^ full(b)");
180 191
208 if (info == -1) 219 if (info == -1)
209 warning ("inverse: matrix singular to machine precision, rcond = %g", rcond); 220 warning ("inverse: matrix singular to machine precision, rcond = %g", rcond);
210 } 221 }
211 else 222 else
212 atmp = a; 223 atmp = a;
224
225 if (atmp.nnz() == 0) // Fast return for all-zeros matrix
226 return atmp;
213 227
214 SparseComplexMatrix result (atmp); 228 SparseComplexMatrix result (atmp);
215 229
216 btmp--; 230 btmp--;
217 231
299 } 313 }
300 314
301 /* 315 /*
302 %!assert (sparse (2) .^ [3, 4], sparse ([8, 16])) 316 %!assert (sparse (2) .^ [3, 4], sparse ([8, 16]))
303 %!assert <47775> (sparse (2i) .^ [3, 4], sparse ([-0-8i, 16])) 317 %!assert <47775> (sparse (2i) .^ [3, 4], sparse ([-0-8i, 16]))
318
319 %!test <*63080>
320 %! Z = sparse ([]);
321 %! A = sparse (zeros (0, 2));
322 %! B = sparse (zeros (2, 0));
323 %! assert (Z ^ 1, Z);
324 %! assert (Z ^ 0, Z);
325 %! assert (Z ^ -1, Z);
326 %! assert (A ^ 1, Z);
327 %! assert (A ^ 0, Z);
328 %! assert (A ^ -1, Z);
329 %! assert (B ^ 1, Z);
330 %! assert (B ^ 0, Z);
331 %! assert (B ^ -1, Z);
332
333 %!test <*63080>
334 %! A = sparse (zeros (2, 2));
335 %! assert (A ^ 1, A);
336 %! assert (A ^ 0, sparse (eye (2, 2)));
337
338 %!test <63080>
339 %! A = sparse (zeros (2, 2));
340 %! assert (A ^ -1, sparse (inf (2, 2)));
341
304 */ 342 */
305 343
306 // -*- 1 -*- 344 // -*- 1 -*-
307 octave_value 345 octave_value
308 elem_xpow (double a, const SparseMatrix& b) 346 elem_xpow (double a, const SparseMatrix& b)
801 result.maybe_compress (true); 839 result.maybe_compress (true);
802 840
803 return result; 841 return result;
804 } 842 }
805 843
844
806 OCTAVE_NAMESPACE_END 845 OCTAVE_NAMESPACE_END