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-log.cc < prev    next >
C/C++ Source or Header  |  1996-09-28  |  6KB  |  270 lines
  1. // f-log.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 "EIG.h"
  29.  
  30. #include "tree-const.h"
  31. #include "user-prefs.h"
  32. #include "error.h"
  33. #include "gripes.h"
  34. #include "utils.h"
  35. #include "help.h"
  36. #include "defun-dld.h"
  37.  
  38. // XXX FIXME XXX -- the next two functions should really be just one...
  39.  
  40. DEFUN_DLD_BUILTIN ("logm", Flogm, Slogm, 2, 1,
  41.   "logm (X): matrix logarithm")
  42. {
  43.   Octave_object retval;
  44.  
  45.   int nargin = args.length ();
  46.  
  47.   if (nargin != 1)
  48.     {
  49.       print_usage ("logm");
  50.       return retval;
  51.     }
  52.  
  53.   tree_constant arg = args(0);
  54.  
  55.   int arg_is_empty = empty_arg ("logm", arg.rows (), arg.columns ());
  56.  
  57.   if (arg_is_empty < 0)
  58.     return retval;
  59.   else if (arg_is_empty > 0)
  60.     return Matrix ();
  61.  
  62.   if (arg.is_real_scalar ())
  63.     {
  64.       double d = arg.double_value ();
  65.       if (d > 0.0)
  66.     retval(0) = log (d);
  67.       else
  68.     {
  69.       Complex dtmp (d);
  70.       retval(0) = log (dtmp);
  71.     }
  72.     }
  73.   else if (arg.is_complex_scalar ())
  74.     {
  75.       Complex c = arg.complex_value ();
  76.       retval(0) = log (c);
  77.     }
  78.   else if (arg.is_real_type ())
  79.     {
  80.       Matrix m = arg.matrix_value ();
  81.  
  82.       if (! error_state)
  83.     {
  84.       int nr = m.rows ();
  85.       int nc = m.columns ();
  86.  
  87.       if (nr == 0 || nc == 0 || nr != nc)
  88.         gripe_square_matrix_required ("logm");
  89.       else
  90.         {
  91.           EIG m_eig (m);
  92.           ComplexColumnVector lambda (m_eig.eigenvalues ());
  93.           ComplexMatrix Q (m_eig.eigenvectors ());
  94.  
  95.           for (int i = 0; i < nr; i++)
  96.         {
  97.           Complex elt = lambda.elem (i);
  98.           if (imag (elt) == 0.0 && real (elt) > 0.0)
  99.             lambda.elem (i) = log (real (elt));
  100.           else
  101.             lambda.elem (i) = log (elt);
  102.         }
  103.  
  104.           ComplexDiagMatrix D (lambda);
  105.           ComplexMatrix result = Q * D * Q.inverse ();
  106.  
  107.           retval(0) = result;
  108.         }
  109.     }
  110.     }
  111.   else if (arg.is_complex_type ())
  112.     {
  113.       ComplexMatrix m = arg.complex_matrix_value ();
  114.  
  115.       if (! error_state)
  116.     {
  117.       int nr = m.rows ();
  118.       int nc = m.columns ();
  119.  
  120.       if (nr == 0 || nc == 0 || nr != nc)
  121.         gripe_square_matrix_required ("logm");
  122.       else
  123.         {
  124.           EIG m_eig (m);
  125.           ComplexColumnVector lambda (m_eig.eigenvalues ());
  126.           ComplexMatrix Q (m_eig.eigenvectors ());
  127.  
  128.           for (int i = 0; i < nr; i++)
  129.         {
  130.           Complex elt = lambda.elem (i);
  131.           if (imag (elt) == 0.0 && real (elt) > 0.0)
  132.             lambda.elem (i) = log (real (elt));
  133.           else
  134.             lambda.elem (i) = log (elt);
  135.         }
  136.  
  137.           ComplexDiagMatrix D (lambda);
  138.           ComplexMatrix result = Q * D * Q.inverse ();
  139.  
  140.           retval(0) = result;
  141.         }
  142.     }
  143.     }
  144.   else
  145.     {
  146.       gripe_wrong_type_arg ("logm", arg);
  147.     }
  148.  
  149.   return retval;
  150. }
  151.  
  152. DEFUN_DLD_BUILTIN ("sqrtm", Fsqrtm, Ssqrtm, 2, 1,
  153.  "sqrtm (X): matrix sqrt")
  154. {
  155.   Octave_object retval;
  156.  
  157.   int nargin = args.length ();
  158.  
  159.   if (nargin != 1)
  160.     {
  161.       print_usage ("sqrtm");
  162.       return retval;
  163.     }
  164.  
  165.   tree_constant arg = args(0);
  166.  
  167.   int arg_is_empty = empty_arg ("sqrtm", arg.rows (), arg.columns ());
  168.  
  169.   if (arg_is_empty < 0)
  170.     return retval;
  171.   else if (arg_is_empty > 0)
  172.     return Matrix ();
  173.  
  174.   if (arg.is_real_scalar ())
  175.     {
  176.       double d = arg.double_value ();
  177.       if (d > 0.0)
  178.     retval(0) = sqrt (d);
  179.       else
  180.     {
  181.       Complex dtmp (d);
  182.       retval(0) = sqrt (dtmp);
  183.     }
  184.     }
  185.   else if (arg.is_complex_scalar ())
  186.     {
  187.       Complex c = arg.complex_value ();
  188.       retval(0) = log (c);
  189.     }
  190.   else if (arg.is_real_type ())
  191.     {
  192.       Matrix m = arg.matrix_value ();
  193.  
  194.       if (! error_state)
  195.     {
  196.       int nr = m.rows ();
  197.       int nc = m.columns ();
  198.  
  199.       if (nr == 0 || nc == 0 || nr != nc)
  200.         gripe_square_matrix_required ("sqrtm");
  201.       else
  202.         {
  203.           EIG m_eig (m);
  204.           ComplexColumnVector lambda (m_eig.eigenvalues ());
  205.           ComplexMatrix Q (m_eig.eigenvectors ());
  206.  
  207.           for (int i = 0; i < nr; i++)
  208.         {
  209.           Complex elt = lambda.elem (i);
  210.           if (imag (elt) == 0.0 && real (elt) > 0.0)
  211.             lambda.elem (i) = sqrt (real (elt));
  212.           else
  213.             lambda.elem (i) = sqrt (elt);
  214.         }
  215.  
  216.           ComplexDiagMatrix D (lambda);
  217.           ComplexMatrix result = Q * D * Q.inverse ();
  218.  
  219.           retval(0) = result;
  220.         }
  221.     }
  222.     }
  223.   else if (arg.is_complex_type ())
  224.     {
  225.       ComplexMatrix m = arg.complex_matrix_value ();
  226.  
  227.       if (! error_state)
  228.     {
  229.       int nr = m.rows ();
  230.       int nc = m.columns ();
  231.  
  232.       if (nr == 0 || nc == 0 || nr != nc)
  233.         gripe_square_matrix_required ("sqrtm");
  234.       else
  235.         {
  236.           EIG m_eig (m);
  237.           ComplexColumnVector lambda (m_eig.eigenvalues ());
  238.           ComplexMatrix Q (m_eig.eigenvectors ());
  239.  
  240.           for (int i = 0; i < nr; i++)
  241.         {
  242.           Complex elt = lambda.elem (i);
  243.           if (imag (elt) == 0.0 && real (elt) > 0.0)
  244.             lambda.elem (i) = sqrt (real (elt));
  245.           else
  246.             lambda.elem (i) = sqrt (elt);
  247.         }
  248.  
  249.           ComplexDiagMatrix D (lambda);
  250.           ComplexMatrix result = Q * D * Q.inverse ();
  251.  
  252.           retval(0) = result;
  253.         }
  254.     }
  255.     }
  256.   else
  257.     {
  258.       gripe_wrong_type_arg ("sqrtm", arg);
  259.     }
  260.  
  261.   return retval;
  262. }
  263.  
  264. /*
  265. ;;; Local Variables: ***
  266. ;;; mode: C++ ***
  267. ;;; page-delimiter: "^/\\*" ***
  268. ;;; End: ***
  269. */
  270.