Mercurial > hg > octave-lyh

comparison liboctave/MArray-defs.h @ 8303:b11c31849b44

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
improve norm computation capabilities
author Jaroslav Hajek <highegg@gmail.com>
date Fri, 31 Oct 2008 08:05:32 +0100
parents 82be108cc558
children b756ce0002db
comparison
equal deleted inserted replaced
8302:f2e050b62199 8303:b11c31849b44
341 (R, T, , T, dynamic_cast<const B<T>&>, R) \ 341 (R, T, , T, dynamic_cast<const B<T>&>, R) \
342 \ 342 \
343 MDIAGARRAY2_DADA_BINOP_FWD_DEFS \ 343 MDIAGARRAY2_DADA_BINOP_FWD_DEFS \
344 (R, T, dynamic_cast<const B<T>&>, R, dynamic_cast<const B<T>&>, R) 344 (R, T, dynamic_cast<const B<T>&>, R, dynamic_cast<const B<T>&>, R)
345 345
346 #define MARRAY_NORM_BODY(TYPE, RTYPE, blas_norm, BLAS_NORM, NAN_VALUE) \
347 \
348 RTYPE retval = NAN_VALUE; \
349 \
350 octave_idx_type len = length (); \
351 \
352 if (len > 0) \
353 { \
354 const TYPE *d = data (); \
355 \
356 if (p == -1) \
357 { \
358 /* Frobenius norm. */ \
359 retval = 0; \
360 \
361 /* precondition */ \
362 RTYPE inf_norm = 0.; \
363 for (octave_idx_type i = 0; i < len; i++) \
364 { \
365 RTYPE d_abs = std::abs (d[i]); \
366 if (d_abs > inf_norm) \
367 inf_norm = d_abs; \
368 } \
369 inf_norm = (inf_norm == octave_Inf || inf_norm == 0. ? 1.0 : \
370 inf_norm); \
371 RTYPE scale = 1. / inf_norm; \
372 \
373 for (octave_idx_type i = 0; i < len; i++) \
374 { \
375 RTYPE d_abs = std::abs (d[i]) * scale; \
376 retval += d_abs * d_abs; \
377 } \
378 \
379 retval = ::sqrt (retval) * inf_norm; \
380 } \
381 else if (p == 2) \
382 F77_FCN (blas_norm, BLAS_NORM) (len, d, 1, retval); \
383 else if (xisinf (p)) \
384 { \
385 octave_idx_type i = 0; \
386 \
387 while (i < len && xisnan (d[i])) \
388 i++; \
389 \
390 if (i < len) \
391 retval = std::abs (d[i]); \
392 \
393 if (p > 0) \
394 { \
395 while (i < len) \
396 { \
397 RTYPE d_abs = std::abs (d[i++]); \
398 \
399 if (d_abs > retval) \
400 retval = d_abs; \
401 } \
402 } \
403 else \
404 { \
405 while (i < len) \
406 { \
407 RTYPE d_abs = std::abs (d[i++]); \
408 \
409 if (d_abs < retval) \
410 retval = d_abs; \
411 } \
412 } \
413 } \
414 else \
415 { \
416 retval = 0; \
417 \
418 for (octave_idx_type i = 0; i < len; i++) \
419 { \
420 RTYPE d_abs = std::abs (d[i]); \
421 retval += pow (d_abs, p); \
422 } \
423 \
424 retval = pow (retval, 1/p); \
425 } \
426 } \
427 \
428 return retval
429
430 // Now we have all the definitions we need. 346 // Now we have all the definitions we need.
431 347
432 #endif 348 #endif