Mercurial > octave-nkf

comparison liboctave/Sparse.cc @ 6988:c7484dcadd4d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-10-09 19:58:32 by dbateman]
author dbateman
date Tue, 09 Oct 2007 19:59:51 +0000
parents 2e7f62e52c13
children 93c65f2a5668
comparison
equal deleted inserted replaced
6987:deb175b6e4a1 6988:c7484dcadd4d
1809 } 1809 }
1810 1810
1811 return retval; 1811 return retval;
1812 } 1812 }
1813 1813
1814 struct
1815 idx_node
1816 {
1817 octave_idx_type i;
1818 struct idx_node *next;
1819 };
1820
1814 template <class T> 1821 template <class T>
1815 Sparse<T> 1822 Sparse<T>
1816 Sparse<T>::index (idx_vector& idx_i, idx_vector& idx_j, int resize_ok) const 1823 Sparse<T>::index (idx_vector& idx_i, idx_vector& idx_j, int resize_ok) const
1817 { 1824 {
1818 Sparse<T> retval; 1825 Sparse<T> retval;
1901 for (octave_idx_type j = 0; j < m; j++) 1908 for (octave_idx_type j = 0; j < m; j++)
1902 { 1909 {
1903 octave_idx_type jj = idx_j.elem (j); 1910 octave_idx_type jj = idx_j.elem (j);
1904 for (octave_idx_type i = cidx(jj); i < cidx(jj+1); i++) 1911 for (octave_idx_type i = cidx(jj); i < cidx(jj+1); i++)
1905 { 1912 {
1913 OCTAVE_QUIT;
1914
1906 octave_idx_type ii = itmp [ridx(i)]; 1915 octave_idx_type ii = itmp [ridx(i)];
1907 if (ii >= 0) 1916 if (ii >= 0)
1908 { 1917 {
1909 X [ii] = data (i); 1918 X [ii] = data (i);
1910 retval.xridx (kk++) = ii; 1919 retval.xridx (kk++) = ii;
1917 } 1926 }
1918 retval.maybe_compress (); 1927 retval.maybe_compress ();
1919 } 1928 }
1920 else 1929 else
1921 { 1930 {
1931 OCTAVE_LOCAL_BUFFER (struct idx_node, nodes, n);
1932 OCTAVE_LOCAL_BUFFER (octave_idx_type, start_nodes, nr);
1933
1934 for (octave_idx_type i = 0; i < nr; i++)
1935 start_nodes[i] = -1;
1936
1937 for (octave_idx_type i = 0; i < n; i++)
1938 {
1939 octave_idx_type ii = idx_i.elem (i);
1940 nodes[i].i = i;
1941 nodes[i].next = 0;
1942
1943 octave_idx_type node = start_nodes[ii];
1944 if (node == -1)
1945 start_nodes[ii] = i;
1946 else
1947 {
1948 struct idx_node inode = nodes[node];
1949 while (inode.next)
1950 inode = *inode.next;
1951 inode.next = nodes + i;
1952 }
1953 }
1954
1922 // First count the number of non-zero elements 1955 // First count the number of non-zero elements
1923 octave_idx_type new_nzmx = 0; 1956 octave_idx_type new_nzmx = 0;
1924 for (octave_idx_type j = 0; j < m; j++) 1957 for (octave_idx_type j = 0; j < m; j++)
1925 { 1958 {
1926 octave_idx_type jj = idx_j.elem (j); 1959 octave_idx_type jj = idx_j.elem (j);
1927 for (octave_idx_type i = 0; i < n; i++) 1960
1961 if (jj < nc)
1928 { 1962 {
1929 OCTAVE_QUIT; 1963 for (octave_idx_type i = cidx(jj);
1930 1964 i < cidx(jj+1); i++)
1931 octave_idx_type ii = idx_i.elem (i);
1932 if (ii < nr && jj < nc)
1933 { 1965 {
1934 for (octave_idx_type k = cidx(jj); k < cidx(jj+1); k++) 1966 OCTAVE_QUIT;
1967
1968 octave_idx_type ii = start_nodes [ridx(i)];
1969
1970 if (ii >= 0)
1935 { 1971 {
1936 if (ridx(k) == ii) 1972 struct idx_node inode = nodes[ii];
1937 new_nzmx++; 1973
1938 if (ridx(k) >= ii) 1974 while (true)
1939 break; 1975 {
1976 if (inode.i >= 0 &&
1977 idx_i.elem (inode.i) < nc)
1978 new_nzmx ++;
1979 if (inode.next == 0)
1980 break;
1981 else
1982 inode = *inode.next;
1983 }
1940 } 1984 }
1941 } 1985 }
1942 } 1986 }
1943 } 1987 }
1944 1988
1989 std::vector<T> X (n);
1945 retval = Sparse<T> (n, m, new_nzmx); 1990 retval = Sparse<T> (n, m, new_nzmx);
1991 octave_idx_type *ri = retval.xridx ();
1946 1992
1947 octave_idx_type kk = 0; 1993 octave_idx_type kk = 0;
1948 retval.xcidx(0) = 0; 1994 retval.xcidx(0) = 0;
1949 for (octave_idx_type j = 0; j < m; j++) 1995 for (octave_idx_type j = 0; j < m; j++)
1950 { 1996 {
1951 octave_idx_type jj = idx_j.elem (j); 1997 octave_idx_type jj = idx_j.elem (j);
1952 for (octave_idx_type i = 0; i < n; i++) 1998 if (jj < nc)
1953 { 1999 {
1954 OCTAVE_QUIT; 2000 for (octave_idx_type i = cidx(jj);
1955 2001 i < cidx(jj+1); i++)
1956 octave_idx_type ii = idx_i.elem (i);
1957 if (ii < nr && jj < nc)
1958 { 2002 {
1959 for (octave_idx_type k = cidx(jj); k < cidx(jj+1); k++) 2003 OCTAVE_QUIT;
2004
2005 octave_idx_type ii = start_nodes [ridx(i)];
2006
2007 if (ii >= 0)
1960 { 2008 {
1961 if (ridx(k) == ii) 2009 struct idx_node inode = nodes[ii];
2010
2011 while (true)
1962 { 2012 {
1963 retval.xdata(kk) = data(k); 2013 if (inode.i >= 0 &&
1964 retval.xridx(kk++) = i; 2014 idx_i.elem (inode.i) < nc)
2015 {
2016 X [inode.i] = data (i);
2017 retval.xridx (kk++) = inode.i;
2018 }
2019
2020 if (inode.next == 0)
2021 break;
2022 else
2023 inode = *inode.next;
1965 } 2024 }
1966 if (ridx(k) >= ii)
1967 break;
1968 } 2025 }
1969 } 2026 }
2027 sort.sort (ri + retval.xcidx (j),
2028 kk - retval.xcidx (j));
2029 for (octave_idx_type p = retval.xcidx (j);
2030 p < kk; p++)
2031 retval.xdata (p) = X [retval.xridx (p)];
2032 retval.xcidx(j+1) = kk;
1970 } 2033 }
1971 retval.xcidx(j+1) = kk;
1972 } 2034 }
1973 } 2035 }
1974 } 2036 }
1975 } 2037 }
1976 } 2038 }
2394 idx_j = tmp[1]; 2456 idx_j = tmp[1];
2395 2457
2396 if (n_idx > 0) 2458 if (n_idx > 0)
2397 idx_i = tmp[0]; 2459 idx_i = tmp[0];
2398 2460
2461 // Take a constant copy of lhs. This means that ridx and family won't
2462 // call make_unique.
2463 const Sparse<LT> c_lhs (lhs);
2464
2399 if (n_idx == 2) 2465 if (n_idx == 2)
2400 { 2466 {
2401 octave_idx_type n = idx_i.freeze (lhs_nr, "row", true); 2467 octave_idx_type n = idx_i.freeze (lhs_nr, "row", true);
2402 octave_idx_type m = idx_j.freeze (lhs_nc, "column", true); 2468 octave_idx_type m = idx_j.freeze (lhs_nc, "column", true);
2403 2469
2451 2517
2452 octave_idx_type ii = idx_i.elem (i); 2518 octave_idx_type ii = idx_i.elem (i);
2453 2519
2454 if (ii < lhs_nr) 2520 if (ii < lhs_nr)
2455 { 2521 {
2456 for (octave_idx_type k = lhs.cidx(jj); 2522 for (octave_idx_type k = c_lhs.cidx(jj);
2457 k < lhs.cidx(jj+1); k++) 2523 k < c_lhs.cidx(jj+1); k++)
2458 { 2524 {
2459 if (lhs.ridx(k) == ii) 2525 if (c_lhs.ridx(k) == ii)
2460 new_nzmx--; 2526 new_nzmx--;
2461 if (lhs.ridx(k) >= ii) 2527 if (c_lhs.ridx(k) >= ii)
2462 break; 2528 break;
2463 } 2529 }
2464 } 2530 }
2465 } 2531 }
2466 } 2532 }
2481 { 2547 {
2482 octave_idx_type iii = 0; 2548 octave_idx_type iii = 0;
2483 octave_idx_type ii = idx_i.elem (iii); 2549 octave_idx_type ii = idx_i.elem (iii);
2484 octave_idx_type ppp = 0; 2550 octave_idx_type ppp = 0;
2485 octave_idx_type ppi = (j >= lhs_nc ? 0 : 2551 octave_idx_type ppi = (j >= lhs_nc ? 0 :
2486 lhs.cidx(j+1) - 2552 c_lhs.cidx(j+1) -
2487 lhs.cidx(j)); 2553 c_lhs.cidx(j));
2488 octave_idx_type pp = (ppp < ppi ? 2554 octave_idx_type pp = (ppp < ppi ?
2489 lhs.ridx(lhs.cidx(j)+ppp) : 2555 c_lhs.ridx(c_lhs.cidx(j)+ppp) :
2490 new_nr); 2556 new_nr);
2491 while (ppp < ppi || iii < n) 2557 while (ppp < ppi || iii < n)
2492 { 2558 {
2493 if (iii < n && ii <= pp) 2559 if (iii < n && ii <= pp)
2494 { 2560 {
2495 if (scalar != RT ()) 2561 if (scalar != RT ())
2496 { 2562 {
2497 stmp.data(kk) = scalar; 2563 stmp.data(kk) = scalar;
2498 stmp.ridx(kk++) = ii; 2564 stmp.ridx(kk++) = ii;
2499 } 2565 }
2500 if (ii == pp) 2566 if (ii == pp)
2501 pp = (++ppp < ppi ? lhs.ridx(lhs.cidx(j)+ppp) : new_nr); 2567 pp = (++ppp < ppi ? c_lhs.ridx(c_lhs.cidx(j)+ppp) : new_nr);
2502 if (++iii < n) 2568 if (++iii < n)
2503 ii = idx_i.elem(iii); 2569 ii = idx_i.elem(iii);
2504 } 2570 }
2505 else 2571 else
2506 { 2572 {
2507 stmp.data(kk) = 2573 stmp.data(kk) =
2508 lhs.data(lhs.cidx(j)+ppp); 2574 c_lhs.data(c_lhs.cidx(j)+ppp);
2509 stmp.ridx(kk++) = pp; 2575 stmp.ridx(kk++) = pp;
2510 pp = (++ppp < ppi ? lhs.ridx(lhs.cidx(j)+ppp) : new_nr); 2576 pp = (++ppp < ppi ? c_lhs.ridx(c_lhs.cidx(j)+ppp) : new_nr);
2511 } 2577 }
2512 } 2578 }
2513 if (++jji < m) 2579 if (++jji < m)
2514 jj = idx_j.elem(jji); 2580 jj = idx_j.elem(jji);
2515 } 2581 }
2516 else if (j < lhs.cols()) 2582 else if (j < lhs_nc)
2517 { 2583 {
2518 for (octave_idx_type i = lhs.cidx(j); 2584 for (octave_idx_type i = c_lhs.cidx(j);
2519 i < lhs.cidx(j+1); i++) 2585 i < c_lhs.cidx(j+1); i++)
2520 { 2586 {
2521 stmp.data(kk) = lhs.data(i); 2587 stmp.data(kk) = c_lhs.data(i);
2522 stmp.ridx(kk++) = lhs.ridx(i); 2588 stmp.ridx(kk++) = c_lhs.ridx(i);
2523 } 2589 }
2524 } 2590 }
2525 stmp.cidx(j+1) = kk; 2591 stmp.cidx(j+1) = kk;
2526 } 2592 }
2527 2593
2634 { 2700 {
2635 octave_idx_type iii = 0; 2701 octave_idx_type iii = 0;
2636 octave_idx_type ii = idx_i.elem (iii); 2702 octave_idx_type ii = idx_i.elem (iii);
2637 octave_idx_type ppp = 0; 2703 octave_idx_type ppp = 0;
2638 octave_idx_type ppi = (j >= lhs_nc ? 0 : 2704 octave_idx_type ppi = (j >= lhs_nc ? 0 :
2639 lhs.cidx(j+1) - 2705 c_lhs.cidx(j+1) -
2640 lhs.cidx(j)); 2706 c_lhs.cidx(j));
2641 octave_idx_type pp = (ppp < ppi ? 2707 octave_idx_type pp = (ppp < ppi ?
2642 lhs.ridx(lhs.cidx(j)+ppp) : 2708 c_lhs.ridx(c_lhs.cidx(j)+ppp) :
2643 new_nr); 2709 new_nr);
2644 while (ppp < ppi || iii < n) 2710 while (ppp < ppi || iii < n)
2645 { 2711 {
2646 if (iii < n && ii <= pp) 2712 if (iii < n && ii <= pp)
2647 { 2713 {
2648 RT rtmp = rhs.elem (rhs_idx_i[iii], 2714 RT rtmp = rhs.elem (rhs_idx_i[iii],
2651 { 2717 {
2652 stmp.data(kk) = rtmp; 2718 stmp.data(kk) = rtmp;
2653 stmp.ridx(kk++) = ii; 2719 stmp.ridx(kk++) = ii;
2654 } 2720 }
2655 if (ii == pp) 2721 if (ii == pp)
2656 pp = (++ppp < ppi ? lhs.ridx(lhs.cidx(j)+ppp) : new_nr); 2722 pp = (++ppp < ppi ? c_lhs.ridx(c_lhs.cidx(j)+ppp) : new_nr);
2657 if (++iii < n) 2723 if (++iii < n)
2658 ii = idx_i.elem(iii); 2724 ii = idx_i.elem(iii);
2659 } 2725 }
2660 else 2726 else
2661 { 2727 {
2662 stmp.data(kk) = 2728 stmp.data(kk) =
2663 lhs.data(lhs.cidx(j)+ppp); 2729 c_lhs.data(c_lhs.cidx(j)+ppp);
2664 stmp.ridx(kk++) = pp; 2730 stmp.ridx(kk++) = pp;
2665 pp = (++ppp < ppi ? lhs.ridx(lhs.cidx(j)+ppp) : new_nr); 2731 pp = (++ppp < ppi ? c_lhs.ridx(c_lhs.cidx(j)+ppp) : new_nr);
2666 } 2732 }
2667 } 2733 }
2668 if (++jji < m) 2734 if (++jji < m)
2669 jj = idx_j.elem(jji); 2735 jj = idx_j.elem(jji);
2670 } 2736 }
2671 else if (j < lhs.cols()) 2737 else if (j < lhs_nc)
2672 { 2738 {
2673 for (octave_idx_type i = lhs.cidx(j); 2739 for (octave_idx_type i = c_lhs.cidx(j);
2674 i < lhs.cidx(j+1); i++) 2740 i < c_lhs.cidx(j+1); i++)
2675 { 2741 {
2676 stmp.data(kk) = lhs.data(i); 2742 stmp.data(kk) = c_lhs.data(i);
2677 stmp.ridx(kk++) = lhs.ridx(i); 2743 stmp.ridx(kk++) = c_lhs.ridx(i);
2678 } 2744 }
2679 } 2745 }
2680 stmp.cidx(j+1) = kk; 2746 stmp.cidx(j+1) = kk;
2681 } 2747 }
2682 2748
2793 2859
2794 octave_idx_type len = idx_i.freeze (lhs_nr * lhs_nc, "matrix"); 2860 octave_idx_type len = idx_i.freeze (lhs_nr * lhs_nc, "matrix");
2795 2861
2796 if (idx_i) 2862 if (idx_i)
2797 { 2863 {
2798 // Take a constant copy of lhs. This means that elem won't
2799 // create missing elements.
2800 const Sparse<LT> c_lhs (lhs);
2801
2802 if (rhs_nr == 0 && rhs_nc == 0) 2864 if (rhs_nr == 0 && rhs_nc == 0)
2803 lhs.maybe_delete_elements (idx_i); 2865 lhs.maybe_delete_elements (idx_i);
2804 else if (len == 0) 2866 else if (len == 0)
2805 { 2867 {
2806 if (! ((rhs_nr == 1 && rhs_nc == 1) 2868 if (! ((rhs_nr == 1 && rhs_nc == 1)