Mercurial > octave-nkf

comparison liboctave/util/oct-binmap.h @ 18824:59975c3cea6b gui-release

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
maint: Periodic merge of stable to gui-release.
author Rik <rik@octave.org>
date Mon, 02 Jun 2014 09:00:22 -0700
parents 8e056300994b dccbc8bff5cb
children af41e41ad28e
comparison
equal deleted inserted replaced
18814:bd1fd4ed3d67 18824:59975c3cea6b
217 // scalar-Sparse 217 // scalar-Sparse
218 template <class U, class T, class R, class F> 218 template <class U, class T, class R, class F>
219 Sparse<U> 219 Sparse<U>
220 binmap (const T& x, const Sparse<R>& ys, F fcn) 220 binmap (const T& x, const Sparse<R>& ys, F fcn)
221 { 221 {
222 octave_idx_type nz = ys.nnz (); 222 R yzero = R ();
223 Sparse<U> retval (ys.rows (), ys.cols (), nz); 223 U fz = fcn (x, yzero);
224 for (octave_idx_type i = 0; i < nz; i++) 224
225 { 225 if (fz == U ()) // Sparsity preserving fcn
226 {
227 octave_idx_type nz = ys.nnz ();
228 Sparse<U> retval (ys.rows (), ys.cols (), nz);
229 copy_or_memcpy (nz, ys.ridx (), retval.ridx ());
230 copy_or_memcpy (ys.cols () + 1, ys.cidx (), retval.cidx ());
231
232 for (octave_idx_type i = 0; i < nz; i++)
233 {
234 octave_quit ();
235 // FIXME: Could keep track of whether fcn call results in a 0.
236 // If no zeroes are created could skip maybe_compress()
237 retval.xdata (i) = fcn (x, ys.data (i));
238 }
239
226 octave_quit (); 240 octave_quit ();
227 retval.xdata (i) = fcn (x, ys.data (i)); 241 retval.maybe_compress (true);
228 } 242 return retval;
229 243 }
230 octave_quit (); 244 else
231 retval.maybe_compress (); 245 return Sparse<U> (binmap<U, T, R, F> (x, ys.array_value (), fcn));
232 return retval;
233 } 246 }
234 247
235 // Sparse-scalar 248 // Sparse-scalar
236 template <class U, class T, class R, class F> 249 template <class U, class T, class R, class F>
237 Sparse<U> 250 Sparse<U>
238 binmap (const Sparse<T>& xs, const R& y, F fcn) 251 binmap (const Sparse<T>& xs, const R& y, F fcn)
239 { 252 {
240 octave_idx_type nz = xs.nnz (); 253 T xzero = T ();
241 Sparse<U> retval (xs.rows (), xs.cols (), nz); 254 U fz = fcn (xzero, y);
242 for (octave_idx_type i = 0; i < nz; i++) 255
243 { 256 if (fz == U ()) // Sparsity preserving fcn
257 {
258 octave_idx_type nz = xs.nnz ();
259 Sparse<U> retval (xs.rows (), xs.cols (), nz);
260 copy_or_memcpy (nz, xs.ridx (), retval.ridx ());
261 copy_or_memcpy (xs.cols () + 1, xs.cidx (), retval.cidx ());
262
263 for (octave_idx_type i = 0; i < nz; i++)
264 {
265 octave_quit ();
266 // FIXME: Could keep track of whether fcn call results in a 0.
267 // If no zeroes are created could skip maybe_compress()
268 retval.xdata (i) = fcn (xs.data (i), y);
269 }
270
244 octave_quit (); 271 octave_quit ();
245 retval.xdata (i) = fcn (xs.data (i), y); 272 retval.maybe_compress (true);
246 } 273 return retval;
247 274 }
248 octave_quit (); 275 else
249 retval.maybe_compress (); 276 return Sparse<U> (binmap<U, T, R, F> (xs.array_value (), y, fcn));
250 return retval;
251 } 277 }
252 278
253 // Sparse-Sparse (treats singletons as scalars) 279 // Sparse-Sparse (treats singletons as scalars)
254 template <class U, class T, class R, class F> 280 template <class U, class T, class R, class F>
255 Sparse<U> 281 Sparse<U>
262 else if (xs.dims () != ys.dims ()) 288 else if (xs.dims () != ys.dims ())
263 gripe_nonconformant (name, xs.dims (), ys.dims ()); 289 gripe_nonconformant (name, xs.dims (), ys.dims ());
264 290
265 T xzero = T (); 291 T xzero = T ();
266 R yzero = R (); 292 R yzero = R ();
267
268 U fz = fcn (xzero, yzero); 293 U fz = fcn (xzero, yzero);
294
269 if (fz == U ()) 295 if (fz == U ())
270 { 296 {
271 // Sparsity-preserving function. Do it efficiently. 297 // Sparsity-preserving function. Do it efficiently.
272 octave_idx_type nr = xs.rows (); 298 octave_idx_type nr = xs.rows ();
273 octave_idx_type nc = xs.cols (); 299 octave_idx_type nc = xs.cols ();
274 Sparse<T> retval (nr, nc); 300 Sparse<T> retval (nr, nc, xs.nnz () + ys.nnz ());
275 301
276 octave_idx_type nz = 0; 302 octave_idx_type nz = 0;
277 // Count nonzeros.
278 for (octave_idx_type j = 0; j < nc; j++) 303 for (octave_idx_type j = 0; j < nc; j++)
279 { 304 {
280 octave_quit (); 305 octave_quit ();
281 octave_idx_type ix = xs.cidx (j); 306
282 octave_idx_type iy = ys.cidx (j); 307 octave_idx_type jx = xs.cidx (j);
283 octave_idx_type ux = xs.cidx (j+1); 308 octave_idx_type jx_max = xs.cidx (j+1);
284 octave_idx_type uy = ys.cidx (j+1); 309 bool jx_lt_max = jx < jx_max;
285 while (ix != ux || iy != uy) 310
311 octave_idx_type jy = ys.cidx (j);
312 octave_idx_type jy_max = ys.cidx (j+1);
313 bool jy_lt_max = jy < jy_max;
314
315 while (jx_lt_max || jy_lt_max)
286 { 316 {
287 octave_idx_type rx = xs.ridx (ix); 317 if (! jy_lt_max
288 octave_idx_type ry = ys.ridx (ix); 318 || (jx_lt_max && (xs.ridx (jx) < ys.ridx (jy))))
289 ix += rx <= ry; 319 {
290 iy += ry <= rx; 320 retval.xridx (nz) = xs.ridx (jx);
321 retval.xdata (nz) = fcn (xs.data (jx), yzero);
322 jx++;
323 jx_lt_max = jx < jx_max;
324 }
325 else if (! jx_lt_max
326 || (jy_lt_max && (ys.ridx (jy) < xs.ridx (jx))))
327 {
328 retval.xridx (nz) = ys.ridx (jy);
329 retval.xdata (nz) = fcn (xzero, ys.data (jy));
330 jy++;
331 jy_lt_max = jy < jy_max;
332 }
333 else
334 {
335 retval.xridx (nz) = xs.ridx (jx);
336 retval.xdata (nz) = fcn (xs.data (jx), ys.data (jy));
337 jx++;
338 jx_lt_max = jx < jx_max;
339 jy++;
340 jy_lt_max = jy < jy_max;
341 }
291 nz++; 342 nz++;
292 } 343 }
293
294 retval.xcidx (j+1) = nz; 344 retval.xcidx (j+1) = nz;
295 } 345 }
296 346
297 // Allocate space. 347 retval.maybe_compress (true);
298 retval.change_capacity (retval.xcidx (nc));
299
300 // Fill.
301 nz = 0;
302 for (octave_idx_type j = 0; j < nc; j++)
303 {
304 octave_quit ();
305 octave_idx_type ix = xs.cidx (j);
306 octave_idx_type iy = ys.cidx (j);
307 octave_idx_type ux = xs.cidx (j+1);
308 octave_idx_type uy = ys.cidx (j+1);
309 while (ix != ux || iy != uy)
310 {
311 octave_idx_type rx = xs.ridx (ix);
312 octave_idx_type ry = ys.ridx (ix);
313 if (rx == ry)
314 {
315 retval.xridx (nz) = rx;
316 retval.xdata (nz) = fcn (xs.data (ix), ys.data (iy));
317 ix++;
318 iy++;
319 }
320 else if (rx < ry)
321 {
322 retval.xridx (nz) = rx;
323 retval.xdata (nz) = fcn (xs.data (ix), yzero);
324 ix++;
325 }
326 else if (ry < rx)
327 {
328 retval.xridx (nz) = ry;
329 retval.xdata (nz) = fcn (xzero, ys.data (iy));
330 iy++;
331 }
332
333 nz++;
334 }
335 }
336
337 retval.maybe_compress ();
338 return retval; 348 return retval;
339 } 349 }
340 else 350 else
341 return Sparse<U> (binmap<U, T, R, F> (xs.array_value (), ys.array_value (), 351 return Sparse<U> (binmap<U, T, R, F> (xs.array_value (), ys.array_value (),
342 fcn, name)); 352 fcn, name));