home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / src / f-sort.cc < prev    next >
C/C++ Source or Header  |  1996-09-28  |  6KB  |  267 lines
  1. // f-sort.cc                                           -*- C++ -*-
  2. /*
  3.  
  4. Copyright (C) 1994, 1995 John W. Eaton
  5.  
  6. This file is part of Octave.
  7.  
  8. Octave is free software; you can redistribute it and/or modify it
  9. under the terms of the GNU General Public License as published by the
  10. Free Software Foundation; either version 2, or (at your option) any
  11. later version.
  12.  
  13. Octave is distributed in the hope that it will be useful, but WITHOUT
  14. ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  15. FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  16. for more details.
  17.  
  18. You should have received a copy of the GNU General Public License
  19. along with Octave; see the file COPYING.  If not, write to the Free
  20. Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
  21.  
  22. */
  23.  
  24. #ifdef HAVE_CONFIG_H
  25. #include "config.h"
  26. #endif
  27.  
  28. #include "tree-const.h"
  29. #include "error.h"
  30. #include "gripes.h"
  31. #include "help.h"
  32. #include "defun-dld.h"
  33.  
  34. static void
  35. mx_sort (Matrix& m, Matrix& idx, int return_idx)
  36. {
  37.   int nr = m.rows ();
  38.   int nc = m.columns ();
  39.   idx.resize (nr, nc);
  40.   int i, j;
  41.  
  42.   if (return_idx)
  43.     {
  44.       for (j = 0; j < nc; j++)
  45.     for (i = 0; i < nr; i++)
  46.       idx.elem (i, j) = i+1;
  47.     }
  48.  
  49.   for (j = 0; j < nc; j++)
  50.     {
  51.       for (int gap = nr/2; gap > 0; gap /= 2)
  52.     for (i = gap; i < nr; i++)
  53.       for (int k = i - gap;
  54.            k >= 0 && m.elem (k, j) > m.elem (k+gap, j);
  55.            k -= gap)
  56.         {
  57.           double tmp = m.elem (k, j);
  58.           m.elem (k, j) = m.elem (k+gap, j);
  59.           m.elem (k+gap, j) = tmp;
  60.  
  61.           if (return_idx)
  62.         {
  63.           double tmp = idx.elem (k, j);
  64.           idx.elem (k, j) = idx.elem (k+gap, j);
  65.           idx.elem (k+gap, j) = tmp;
  66.         }
  67.         }
  68.     }
  69. }
  70.  
  71. static void
  72. mx_sort (RowVector& v, RowVector& idx, int return_idx)
  73. {
  74.   int n = v.capacity ();
  75.   idx.resize (n);
  76.   int i;
  77.  
  78.   if (return_idx)
  79.     for (i = 0; i < n; i++)
  80.       idx.elem (i) = i+1;
  81.  
  82.   for (int gap = n/2; gap > 0; gap /= 2)
  83.     for (i = gap; i < n; i++)
  84.       for (int k = i - gap;
  85.        k >= 0 && v.elem (k) > v.elem (k+gap);
  86.        k -= gap)
  87.     {
  88.       double tmp = v.elem (k);
  89.       v.elem (k) = v.elem (k+gap);
  90.       v.elem (k+gap) = tmp;
  91.  
  92.       if (return_idx)
  93.         {
  94.           double tmp = idx.elem (k);
  95.           idx.elem (k) = idx.elem (k+gap);
  96.           idx.elem (k+gap) = tmp;
  97.         }
  98.     }
  99. }
  100.  
  101. static void
  102. mx_sort (ComplexMatrix& cm, Matrix& idx, int return_idx)
  103. {
  104.   int nr = cm.rows ();
  105.   int nc = cm.columns ();
  106.   idx.resize (nr, nc);
  107.   int i, j;
  108.  
  109.   if (return_idx)
  110.     {
  111.       for (j = 0; j < nc; j++)
  112.     for (i = 0; i < nr; i++)
  113.       idx.elem (i, j) = i+1;
  114.     }
  115.  
  116.   for (j = 0; j < nc; j++)
  117.     {
  118.       for (int gap = nr/2; gap > 0; gap /= 2)
  119.     for (i = gap; i < nr; i++)
  120.       for (int k = i - gap;
  121.            k >= 0 && abs (cm.elem (k, j)) > abs (cm.elem (k+gap, j));
  122.            k -= gap)
  123.         {
  124.           Complex ctmp = cm.elem (k, j);
  125.           cm.elem (k, j) = cm.elem (k+gap, j);
  126.           cm.elem (k+gap, j) = ctmp;
  127.  
  128.           if (return_idx)
  129.         {
  130.           double tmp = idx.elem (k, j);
  131.           idx.elem (k, j) = idx.elem (k+gap, j);
  132.           idx.elem (k+gap, j) = tmp;
  133.         }
  134.         }
  135.     }
  136. }
  137.  
  138. static void
  139. mx_sort (ComplexRowVector& cv, RowVector& idx, int return_idx)
  140. {
  141.   int n = cv.capacity ();
  142.   idx.resize (n);
  143.   int i;
  144.  
  145.   if (return_idx)
  146.     for (i = 0; i < n; i++)
  147.       idx.elem (i) = i+1;
  148.  
  149.   for (int gap = n/2; gap > 0; gap /= 2)
  150.     for (i = gap; i < n; i++)
  151.       for (int k = i - gap;
  152.        k >= 0 && abs (cv.elem (k)) > abs (cv.elem (k+gap));
  153.        k -= gap)
  154.     {
  155.       Complex tmp = cv.elem (k);
  156.       cv.elem (k) = cv.elem (k+gap);
  157.       cv.elem (k+gap) = tmp;
  158.  
  159.       if (return_idx)
  160.         {
  161.           double tmp = idx.elem (k);
  162.           idx.elem (k) = idx.elem (k+gap);
  163.           idx.elem (k+gap) = tmp;
  164.         }
  165.     }
  166. }
  167.  
  168. DEFUN_DLD_BUILTIN ("sort", Fsort, Ssort, 2, 2,
  169.   "[S, I] = sort (X)\n\
  170. \n\
  171. sort the columns of X, optionally return sort index")
  172. {
  173.   Octave_object retval;
  174.  
  175.   int nargin = args.length ();
  176.  
  177.   if (nargin != 1)
  178.     {
  179.       print_usage ("sort");
  180.       return retval;
  181.     }
  182.  
  183.   int return_idx = nargout > 1;
  184.   if (return_idx)
  185.     retval.resize (2);
  186.   else
  187.     retval.resize (1);
  188.  
  189.   tree_constant arg = args(0);
  190.  
  191.   if (arg.is_real_type ())
  192.     {
  193.       Matrix m = arg.matrix_value ();
  194.  
  195.       if (! error_state)
  196.     {
  197.       if (m.rows () == 1)
  198.         {
  199.           int nc = m.columns ();
  200.           RowVector v (nc);
  201.           for (int i = 0; i < nc; i++)
  202.         v.elem (i) = m.elem (0, i);
  203.           RowVector idx;
  204.           mx_sort (v, idx, return_idx);
  205.  
  206.           retval(0) = tree_constant (v, 0);
  207.           if (return_idx)
  208.         retval(1) = tree_constant (idx, 0);
  209.         }
  210.       else
  211.         {
  212. // Sorts m in place, optionally computes index Matrix.
  213.           Matrix idx;
  214.           mx_sort (m, idx, return_idx);
  215.  
  216.           retval(0) = m;
  217.           if (return_idx)
  218.         retval(1) = idx;
  219.         }
  220.     }
  221.     }
  222.   else if (arg.is_complex_type ())
  223.     {
  224.       ComplexMatrix cm = arg.complex_matrix_value ();
  225.  
  226.       if (! error_state)
  227.     {
  228.       if (cm.rows () == 1)
  229.         {
  230.           int nc = cm.columns ();
  231.           ComplexRowVector cv (nc);
  232.           for (int i = 0; i < nc; i++)
  233.         cv.elem (i) = cm.elem (0, i);
  234.           RowVector idx;
  235.           mx_sort (cv, idx, return_idx);
  236.  
  237.           retval(0) = tree_constant (cv, 0);
  238.           if (return_idx)
  239.         retval(1) = tree_constant (idx, 0);
  240.         }
  241.       else
  242.         {
  243. // Sorts cm in place, optionally computes index Matrix.
  244.           Matrix idx;
  245.           mx_sort (cm, idx, return_idx);
  246.  
  247.           retval(0) = cm;
  248.           if (return_idx)
  249.         retval(1) = idx;
  250.         }
  251.     }
  252.     }
  253.   else
  254.     {
  255.       gripe_wrong_type_arg ("sort", arg);
  256.     }
  257.  
  258.   return retval;
  259. }
  260.  
  261. /*
  262. ;;; Local Variables: ***
  263. ;;; mode: C++ ***
  264. ;;; page-delimiter: "^/\\*" ***
  265. ;;; End: ***
  266. */
  267.