Mercurial > hg > octave-nkf

comparison src/DLD-FUNCTIONS/sparse.cc @ 7515:f3c00dc0912b

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Eliminate the rest of the dispatched sparse functions
author David Bateman <dbateman@free.fr>
date Fri, 22 Feb 2008 15:50:51 +0100
parents f5005d9510f4
children 8a42498edb30
comparison
equal deleted inserted replaced
7514:4f6a73fd8df9 7515:f3c00dc0912b
36 #include "quit.h" 36 #include "quit.h"
37 37
38 #include "ov-re-sparse.h" 38 #include "ov-re-sparse.h"
39 #include "ov-cx-sparse.h" 39 #include "ov-cx-sparse.h"
40 #include "ov-bool-sparse.h" 40 #include "ov-bool-sparse.h"
41
42 static bool
43 is_sparse (const octave_value& arg)
44 {
45 return (arg.is_sparse_type ());
46 }
47 41
48 DEFUN_DLD (issparse, args, , 42 DEFUN_DLD (issparse, args, ,
49 "-*- texinfo -*-\n\ 43 "-*- texinfo -*-\n\
50 @deftypefn {Loadable Function} {} issparse (@var{expr})\n\ 44 @deftypefn {Loadable Function} {} issparse (@var{expr})\n\
51 Return 1 if the value of the expression @var{expr} is a sparse matrix.\n\ 45 Return 1 if the value of the expression @var{expr} is a sparse matrix.\n\
55 { 49 {
56 print_usage (); 50 print_usage ();
57 return octave_value (); 51 return octave_value ();
58 } 52 }
59 else 53 else
60 return octave_value (is_sparse (args(0))); 54 return octave_value (args(0).is_sparse_type ());
61 } 55 }
62 56
63 DEFUN_DLD (sparse, args, , 57 DEFUN_DLD (sparse, args, ,
64 "-*- texinfo -*-\n\ 58 "-*- texinfo -*-\n\
65 @deftypefn {Loadable Function} {@var{s} =} sparse (@var{a})\n\ 59 @deftypefn {Loadable Function} {@var{s} =} sparse (@var{a})\n\
130 124
131 if (nargin == 1) 125 if (nargin == 1)
132 { 126 {
133 octave_value arg = args (0); 127 octave_value arg = args (0);
134 128
135 if (is_sparse (arg)) 129 if (arg.is_sparse_type ())
136 { 130 {
137 if (use_complex) 131 if (use_complex)
138 { 132 {
139 SparseComplexMatrix sm = arg.sparse_complex_matrix_value (); 133 SparseComplexMatrix sm = arg.sparse_complex_matrix_value ();
140 retval = new octave_sparse_complex_matrix (sm); 134 retval = new octave_sparse_complex_matrix (sm);
385 gripe_wrong_type_arg ("full", args(0)); 379 gripe_wrong_type_arg ("full", args(0));
386 380
387 return retval; 381 return retval;
388 } 382 }
389 383
390 #define SPARSE_DIM_ARG_BODY(NAME, FUNC) \
391 int nargin = args.length(); \
392 octave_value retval; \
393 if ((nargin != 1 ) && (nargin != 2)) \
394 print_usage (); \
395 else { \
396 int dim = (nargin == 1 ? -1 : args(1).int_value(true) - 1); \
397 if (error_state) return retval; \
398 if (dim < -1 || dim > 1) { \
399 error (#NAME ": invalid dimension argument = %d", dim + 1); \
400 return retval; \
401 } \
402 if (args(0).type_id () == \
403 octave_sparse_matrix::static_type_id () || args(0).type_id () == \
404 octave_sparse_bool_matrix::static_type_id ()) { \
405 retval = args(0).sparse_matrix_value () .FUNC (dim); \
406 } else if (args(0).type_id () == \
407 octave_sparse_complex_matrix::static_type_id ()) { \
408 retval = args(0).sparse_complex_matrix_value () .FUNC (dim); \
409 } else \
410 print_usage (); \
411 } \
412 return retval
413
414 // PKG_ADD: dispatch ("prod", "spprod", "sparse matrix");
415 // PKG_ADD: dispatch ("prod", "spprod", "sparse complex matrix");
416 // PKG_ADD: dispatch ("prod", "spprod", "sparse bool matrix");
417 DEFUN_DLD (spprod, args, ,
418 "-*- texinfo -*-\n\
419 @deftypefn {Loadable Function} {@var{y} =} spprod (@var{x},@var{dim})\n\
420 Product of elements along dimension @var{dim}. If @var{dim} is omitted,\n\
421 it defaults to 1 (column-wise products).\n\
422 @seealso{spsum, spsumsq}\n\
423 @end deftypefn")
424 {
425 SPARSE_DIM_ARG_BODY (spprod, prod);
426 }
427
428 // PKG_ADD: dispatch ("cumprod", "spcumprod", "sparse matrix");
429 // PKG_ADD: dispatch ("cumprod", "spcumprod", "sparse complex matrix");
430 // PKG_ADD: dispatch ("cumprod", "spcumprod", "sparse bool matrix");
431 DEFUN_DLD (spcumprod, args, ,
432 "-*- texinfo -*-\n\
433 @deftypefn {Loadable Function} {@var{y} =} spcumprod (@var{x},@var{dim})\n\
434 Cumulative product of elements along dimension @var{dim}. If @var{dim}\n\
435 is omitted, it defaults to 1 (column-wise cumulative products).\n\
436 @seealso{spcumsum}\n\
437 @end deftypefn")
438 {
439 SPARSE_DIM_ARG_BODY (spcumprod, cumprod);
440 }
441
442 // PKG_ADD: dispatch ("sum", "spsum", "sparse matrix");
443 // PKG_ADD: dispatch ("sum", "spsum", "sparse complex matrix");
444 // PKG_ADD: dispatch ("sum", "spsum", "sparse bool matrix");
445 DEFUN_DLD (spsum, args, ,
446 "-*- texinfo -*-\n\
447 @deftypefn {Loadable Function} {@var{y} =} spsum (@var{x},@var{dim})\n\
448 Sum of elements along dimension @var{dim}. If @var{dim} is omitted, it\n\
449 defaults to 1 (column-wise sum).\n\
450 @seealso{spprod, spsumsq}\n\
451 @end deftypefn")
452 {
453 SPARSE_DIM_ARG_BODY (spsum, sum);
454 }
455
456 // PKG_ADD: dispatch ("cumsum", "spcumsum", "sparse matrix");
457 // PKG_ADD: dispatch ("cumsum", "spcumsum", "sparse complex matrix");
458 // PKG_ADD: dispatch ("cumsum", "spcumsum", "sparse bool matrix");
459 DEFUN_DLD (spcumsum, args, ,
460 "-*- texinfo -*-\n\
461 @deftypefn {Loadable Function} {@var{y} =} spcumsum (@var{x},@var{dim})\n\
462 Cumulative sum of elements along dimension @var{dim}. If @var{dim}\n\
463 is omitted, it defaults to 1 (column-wise cumulative sums).\n\
464 @seealso{spcumprod}\n\
465 @end deftypefn")
466 {
467 SPARSE_DIM_ARG_BODY (spcumsum, cumsum);
468 }
469
470 // PKG_ADD: dispatch ("sumsq", "spsumsq", "sparse matrix");
471 // PKG_ADD: dispatch ("sumsq", "spsumsq", "sparse complex matrix");
472 // PKG_ADD: dispatch ("sumsq", "spsumsq", "sparse bool matrix");
473 DEFUN_DLD (spsumsq, args, ,
474 "-*- texinfo -*-\n\
475 @deftypefn {Loadable Function} {@var{y} =} spsumsq (@var{x},@var{dim})\n\
476 Sum of squares of elements along dimension @var{dim}. If @var{dim}\n\
477 is omitted, it defaults to 1 (column-wise sum of squares).\n\
478 This function is equivalent to computing\n\
479 @example\n\
480 spsum (x .* spconj (x), dim)\n\
481 @end example\n\
482 but it uses less memory and avoids calling @code{spconj} if @var{x} is\n\
483 real.\n\
484 @seealso{spprod, spsum}\n\
485 @end deftypefn")
486 {
487 SPARSE_DIM_ARG_BODY (spsumsq, sumsq);
488 }
489
490
491 static octave_value
492 make_spdiag (const octave_value& a, const octave_value& b)
493 {
494 octave_value retval;
495
496 if (a.is_complex_type ())
497 {
498 SparseComplexMatrix m = a.sparse_complex_matrix_value ();
499 octave_idx_type k = b.nint_value(true);
500
501 if (error_state)
502 return retval;
503
504 octave_idx_type nr = m.rows ();
505 octave_idx_type nc = m.columns ();
506
507 if (nr == 0 || nc == 0)
508 retval = m;
509 else if (nr == 1 || nc == 1)
510 {
511 octave_idx_type roff = 0;
512 octave_idx_type coff = 0;
513 if (k > 0)
514 {
515 roff = 0;
516 coff = k;
517 }
518 else if (k < 0)
519 {
520 k = -k;
521 roff = k;
522 coff = 0;
523 }
524
525 if (nr == 1)
526 {
527 octave_idx_type n = nc + k;
528 octave_idx_type nz = m.nzmax ();
529 SparseComplexMatrix r (n, n, nz);
530 for (octave_idx_type i = 0; i < coff+1; i++)
531 r.xcidx (i) = 0;
532 for (octave_idx_type j = 0; j < nc; j++)
533 {
534 for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++)
535 {
536 r.xdata (i) = m.data (i);
537 r.xridx (i) = j + roff;
538 }
539 r.xcidx (j+coff+1) = m.cidx(j+1);
540 }
541 for (octave_idx_type i = nc+coff+1; i < n+1; i++)
542 r.xcidx (i) = nz;
543 retval = r;
544 }
545 else
546 {
547 octave_idx_type n = nr + k;
548 octave_idx_type nz = m.nzmax ();
549 octave_idx_type ii = 0;
550 octave_idx_type ir = m.ridx(0);
551 SparseComplexMatrix r (n, n, nz);
552 for (octave_idx_type i = 0; i < coff+1; i++)
553 r.xcidx (i) = 0;
554 for (octave_idx_type i = 0; i < nr; i++)
555 {
556 if (ir == i)
557 {
558 r.xdata (ii) = m.data (ii);
559 r.xridx (ii++) = ir + roff;
560 if (ii != nz)
561 ir = m.ridx (ii);
562 }
563 r.xcidx (i+coff+1) = ii;
564 }
565 for (octave_idx_type i = nr+coff+1; i < n+1; i++)
566 r.xcidx (i) = nz;
567 retval = r;
568 }
569 }
570 else
571 {
572 SparseComplexMatrix r = m.diag (k);
573 // Don't use numel, since it can overflow for very large matrices
574 if (r.rows () > 0 && r.cols () > 0)
575 retval = r;
576 }
577 }
578 else if (a.is_real_type ())
579 {
580 SparseMatrix m = a.sparse_matrix_value ();
581
582 octave_idx_type k = b.nint_value(true);
583
584 if (error_state)
585 return retval;
586
587 octave_idx_type nr = m.rows ();
588 octave_idx_type nc = m.columns ();
589
590 if (nr == 0 || nc == 0)
591 retval = m;
592 else if (nr == 1 || nc == 1)
593 {
594 octave_idx_type roff = 0;
595 octave_idx_type coff = 0;
596 if (k > 0)
597 {
598 roff = 0;
599 coff = k;
600 }
601 else if (k < 0)
602 {
603 k = -k;
604 roff = k;
605 coff = 0;
606 }
607
608 if (nr == 1)
609 {
610 octave_idx_type n = nc + k;
611 octave_idx_type nz = m.nzmax ();
612 SparseMatrix r (n, n, nz);
613
614 for (octave_idx_type i = 0; i < coff+1; i++)
615 r.xcidx (i) = 0;
616 for (octave_idx_type j = 0; j < nc; j++)
617 {
618 for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++)
619 {
620 r.xdata (i) = m.data (i);
621 r.xridx (i) = j + roff;
622 }
623 r.xcidx (j+coff+1) = m.cidx(j+1);
624 }
625 for (octave_idx_type i = nc+coff+1; i < n+1; i++)
626 r.xcidx (i) = nz;
627 retval = r;
628 }
629 else
630 {
631 octave_idx_type n = nr + k;
632 octave_idx_type nz = m.nzmax ();
633 octave_idx_type ii = 0;
634 octave_idx_type ir = m.ridx(0);
635 SparseMatrix r (n, n, nz);
636 for (octave_idx_type i = 0; i < coff+1; i++)
637 r.xcidx (i) = 0;
638 for (octave_idx_type i = 0; i < nr; i++)
639 {
640 if (ir == i)
641 {
642 r.xdata (ii) = m.data (ii);
643 r.xridx (ii++) = ir + roff;
644 if (ii != nz)
645 ir = m.ridx (ii);
646 }
647 r.xcidx (i+coff+1) = ii;
648 }
649 for (octave_idx_type i = nr+coff+1; i < n+1; i++)
650 r.xcidx (i) = nz;
651 retval = r;
652 }
653 }
654 else
655 {
656 SparseMatrix r = m.diag (k);
657 if (r.rows () > 0 && r.cols () > 0)
658 retval = r;
659 }
660 }
661 else
662 gripe_wrong_type_arg ("spdiag", a);
663
664 return retval;
665 }
666
667 static octave_value
668 make_spdiag (const octave_value& a)
669 {
670 octave_value retval;
671 octave_idx_type nr = a.rows ();
672 octave_idx_type nc = a.columns ();
673
674 if (nr == 0 || nc == 0)
675 retval = SparseMatrix ();
676 else
677 retval = make_spdiag (a, octave_value (0.));
678
679 return retval;
680 }
681
682 // PKG_ADD: dispatch ("diag", "spdiag", "sparse matrix");
683 // PKG_ADD: dispatch ("diag", "spdiag", "sparse complex matrix");
684 // PKG_ADD: dispatch ("diag", "spdiag", "sparse bool matrix");
685 DEFUN_DLD (spdiag, args, ,
686 "-*- texinfo -*-\n\
687 @deftypefn {Loadable Function} {} spdiag (@var{v}, @var{k})\n\
688 Return a diagonal matrix with the sparse vector @var{v} on diagonal\n\
689 @var{k}. The second argument is optional. If it is positive, the vector is\n\
690 placed on the @var{k}-th super-diagonal. If it is negative, it is placed\n\
691 on the @var{-k}-th sub-diagonal. The default value of @var{k} is 0, and\n\
692 the vector is placed on the main diagonal. For example,\n\
693 \n\
694 @example\n\
695 @group\n\
696 spdiag ([1, 2, 3], 1)\n\
697 ans =\n\
698 \n\
699 Compressed Column Sparse (rows=4, cols=4, nnz=3)\n\
700 (1 , 2) -> 1\n\
701 (2 , 3) -> 2\n\
702 (3 , 4) -> 3\n\
703 @end group\n\
704 @end example\n\
705 \n\
706 @noindent\n\
707 Given a matrix argument, instead of a vector, @code{spdiag} extracts the\n\
708 @var{k}-th diagonal of the sparse matrix.\n\
709 @seealso{diag}\n\
710 @end deftypefn")
711 {
712 octave_value retval;
713
714 int nargin = args.length ();
715
716 if (nargin == 1 && args(0).is_defined ())
717 retval = make_spdiag (args(0));
718 else if (nargin == 2 && args(0).is_defined () && args(1).is_defined ())
719 retval = make_spdiag (args(0), args(1));
720 else
721 print_usage ();
722
723 return retval;
724 }
725
726 /* 384 /*
727 ;;; Local Variables: *** 385 ;;; Local Variables: ***
728 ;;; mode: C++ *** 386 ;;; mode: C++ ***
729 ;;; End: *** 387 ;;; End: ***
730 */ 388 */