home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-base.tgz / octave-1.1.1p1-base.tar / fsf / octave / src / xdiv.cc < prev    next >
C/C++ Source or Header  |  1995-01-03  |  8KB  |  360 lines
  1. // xdiv.cc                                               -*- C++ -*-
  2. /*
  3.  
  4. Copyright (C) 1992, 1993, 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 <assert.h>
  29. #include <Complex.h>
  30.  
  31. #include "xdiv.h"
  32. #include "dMatrix.h"
  33. #include "CMatrix.h"
  34. #include "tree-const.h"
  35. #include "error.h"
  36.  
  37. static inline int
  38. result_ok (int info, double rcond, int warn = 1)
  39. {
  40.   assert (info != -1);
  41.  
  42.   if (info == -2)
  43.     {
  44.       if (warn)
  45.     warning ("matrix singular to machine precision, rcond = %g", rcond);
  46.       else
  47.     error ("matrix singular to machine precision, rcond = %g", rcond);
  48.  
  49.       return 0;
  50.     }
  51.   else
  52.     return 1;
  53. }
  54.  
  55. static inline int
  56. mx_leftdiv_conform (int a_nr, int a_nc, int b_nr, int warn = 1)
  57. {
  58.   if (a_nr != b_nr)
  59.     {
  60.       error ("number of rows must be the same for left division");
  61.       return 0;
  62.     }
  63.  
  64.   return 1;
  65. }
  66.  
  67. static inline int
  68. mx_div_conform (int b_nr, int b_nc, int a_nc, int warn = 1)
  69. {
  70.   if (a_nc != b_nc)
  71.     {
  72.       error ("number of columns must be the same for right division");
  73.       return 0;
  74.     }
  75.  
  76.   return 1;
  77. }
  78.  
  79. // Right division functions.
  80. //
  81. //       op2 / op1:   m   cm
  82. //            +--   +---+----+
  83. //   matrix         | 1 |  3 |
  84. //                  +---+----+
  85. //   complex_matrix | 2 |  4 |
  86. //                  +---+----+
  87.  
  88. // -*- 1 -*-
  89. tree_constant
  90. xdiv (const Matrix& a, const Matrix& b)
  91. {
  92.   if (! mx_div_conform (b.rows (), b.columns (), a.columns ()))
  93.     return tree_constant ();
  94.  
  95.   Matrix atmp = a.transpose ();
  96.   Matrix btmp = b.transpose ();
  97.  
  98.   int info;
  99.   if (btmp.rows () == btmp.columns ())
  100.     {
  101.       double rcond = 0.0;
  102.       Matrix result = btmp.solve (atmp, info, rcond);
  103.       if (result_ok (info, rcond))
  104.     return tree_constant (result.transpose ());
  105.     }
  106.  
  107.   int rank;
  108.   Matrix result = btmp.lssolve (atmp, info, rank);
  109.  
  110.   return tree_constant (result.transpose ());
  111. }
  112.  
  113. // -*- 2 -*-
  114. tree_constant
  115. xdiv (const Matrix& a, const ComplexMatrix& b)
  116. {
  117.   if (! mx_div_conform (b.rows (), b.columns (), a.columns ()))
  118.     return tree_constant ();
  119.  
  120.   Matrix atmp = a.transpose ();
  121.   ComplexMatrix btmp = b.hermitian ();
  122.  
  123.   int info;
  124.   if (btmp.rows () == btmp.columns ())
  125.     {
  126.       double rcond = 0.0;
  127.       ComplexMatrix result = btmp.solve (atmp, info, rcond);
  128.       if (result_ok (info, rcond))
  129.     return tree_constant (result.hermitian ());
  130.     }
  131.  
  132.   int rank;
  133.   ComplexMatrix result = btmp.lssolve (atmp, info, rank);
  134.  
  135.   return tree_constant (result.hermitian ());
  136. }
  137.  
  138. // -*- 3 -*-
  139. tree_constant
  140. xdiv (const ComplexMatrix& a, const Matrix& b)
  141. {
  142.   if (! mx_div_conform (b.rows (), b.columns (), a.columns ()))
  143.     return tree_constant ();
  144.  
  145.   ComplexMatrix atmp = a.hermitian ();
  146.   Matrix btmp = b.transpose ();
  147.  
  148.   int info;
  149.   if (btmp.rows () == btmp.columns ())
  150.     {
  151.       double rcond = 0.0;
  152.       ComplexMatrix result = btmp.solve (atmp, info, rcond);
  153.       if (result_ok (info, rcond))
  154.     return tree_constant (result.hermitian ());
  155.     }
  156.  
  157.   int rank;
  158.   ComplexMatrix result = btmp.lssolve (atmp, info, rank);
  159.  
  160.   return tree_constant (result.hermitian ());
  161. }
  162.  
  163. // -*- 4 -*-
  164. tree_constant
  165. xdiv (const ComplexMatrix& a, const ComplexMatrix& b)
  166. {
  167.   if (! mx_div_conform (b.rows (), b.columns (), a.columns ()))
  168.     return tree_constant ();
  169.  
  170.   ComplexMatrix atmp = a.hermitian ();
  171.   ComplexMatrix btmp = b.hermitian ();
  172.  
  173.   int info;
  174.   if (btmp.rows () == btmp.columns ())
  175.     {
  176.       double rcond = 0.0;
  177.       ComplexMatrix result = btmp.solve (atmp, info, rcond);
  178.       if (result_ok (info, rcond))
  179.     return tree_constant (result.hermitian ());
  180.     }
  181.  
  182.   int rank;
  183.   ComplexMatrix result = btmp.lssolve (atmp, info, rank);
  184.  
  185.   return tree_constant (result.hermitian ());
  186. }
  187.  
  188. // Funny element by element division operations.
  189. //
  190. //       op2 \ op1:   s   cs
  191. //            +--   +---+----+
  192. //   matrix         | 1 |  3 |
  193. //                  +---+----+
  194. //   complex_matrix | 2 |  4 |
  195. //                  +---+----+
  196.  
  197. tree_constant
  198. x_el_div (double a, const Matrix& b)
  199. {
  200.   int nr = b.rows ();
  201.   int nc = b.columns ();
  202.  
  203.   Matrix result (nr, nc);
  204.  
  205.   for (int j = 0; j < nc; j++)
  206.     for (int i = 0; i < nr; i++)
  207.       result.elem (i, j) = a / b.elem (i, j);
  208.  
  209.   return tree_constant (result);
  210. }
  211.  
  212. tree_constant
  213. x_el_div (double a, const ComplexMatrix& b)
  214. {
  215.   int nr = b.rows ();
  216.   int nc = b.columns ();
  217.  
  218.   ComplexMatrix result (nr, nc);
  219.  
  220.   for (int j = 0; j < nc; j++)
  221.     for (int i = 0; i < nr; i++)
  222.       result.elem (i, j) = a / b.elem (i, j);
  223.  
  224.   return tree_constant (result);
  225. }
  226.  
  227. tree_constant
  228. x_el_div (const Complex a, const Matrix& b)
  229. {
  230.   int nr = b.rows ();
  231.   int nc = b.columns ();
  232.  
  233.   ComplexMatrix result (nr, nc);
  234.  
  235.   for (int j = 0; j < nc; j++)
  236.     for (int i = 0; i < nr; i++)
  237.       result.elem (i, j) = a / b.elem (i, j);
  238.  
  239.   return tree_constant (result);
  240. }
  241.  
  242. tree_constant
  243. x_el_div (const Complex a, const ComplexMatrix& b)
  244. {
  245.   int nr = b.rows ();
  246.   int nc = b.columns ();
  247.  
  248.   ComplexMatrix result (nr, nc);
  249.  
  250.   for (int j = 0; j < nc; j++)
  251.     for (int i = 0; i < nr; i++)
  252.       result.elem (i, j) = a / b.elem (i, j);
  253.  
  254.   return tree_constant (result);
  255. }
  256.  
  257. // Left division functions.
  258. //
  259. //       op2 \ op1:   m   cm
  260. //            +--   +---+----+
  261. //   matrix         | 1 |  3 |
  262. //                  +---+----+
  263. //   complex_matrix | 2 |  4 |
  264. //                  +---+----+
  265.  
  266. // -*- 1 -*-
  267. tree_constant
  268. xleftdiv (const Matrix& a, const Matrix& b)
  269. {
  270.   if (! mx_leftdiv_conform (a.rows (), a.columns (), b.rows ()))
  271.     return tree_constant ();
  272.  
  273.   int info;
  274.   if (a.rows () == a.columns ())
  275.     {
  276.       double rcond = 0.0;
  277.       Matrix result = a.solve (b, info, rcond);
  278.       if (result_ok (info, rcond))
  279.     return tree_constant (result);
  280.     }
  281.  
  282.   int rank;
  283.   Matrix result = a.lssolve (b, info, rank);
  284.  
  285.   return tree_constant (result);
  286. }
  287.  
  288. // -*- 2 -*-
  289. tree_constant
  290. xleftdiv (const Matrix& a, const ComplexMatrix& b)
  291. {
  292.   if (! mx_leftdiv_conform (a.rows (), a.columns (), b.rows ()))
  293.     return tree_constant ();
  294.  
  295.   int info;
  296.   if (a.rows () == a.columns ())
  297.     {
  298.       double rcond = 0.0;
  299.       ComplexMatrix result = a.solve (b, info, rcond);
  300.       if (result_ok (info, rcond))
  301.     return tree_constant (result);
  302.     }
  303.  
  304.   int rank;
  305.   ComplexMatrix result = a.lssolve (b, info, rank);
  306.  
  307.   return tree_constant (result);
  308. }
  309.  
  310. // -*- 3 -*-
  311. tree_constant
  312. xleftdiv (const ComplexMatrix& a, const Matrix& b)
  313. {
  314.   if (! mx_leftdiv_conform (a.rows (), a.columns (), b.rows ()))
  315.     return tree_constant ();
  316.  
  317.   int info;
  318.   if (a.rows () == a.columns ())
  319.     {
  320.       double rcond = 0.0;
  321.       ComplexMatrix result = a.solve (b, info, rcond);
  322.       if (result_ok (info, rcond))
  323.     return tree_constant (result);
  324.     }
  325.  
  326.   int rank;
  327.   ComplexMatrix result = a.lssolve (b, info, rank);
  328.  
  329.   return tree_constant (result);
  330. }
  331.  
  332. // -*- 4 -*-
  333. tree_constant
  334. xleftdiv (const ComplexMatrix& a, const ComplexMatrix& b)
  335. {
  336.   if (! mx_leftdiv_conform (a.rows (), a.columns (), b.rows ()))
  337.     return tree_constant ();
  338.  
  339.   int info;
  340.   if (a.rows () == a.columns ())
  341.     {
  342.       double rcond = 0.0;
  343.       ComplexMatrix result = a.solve (b, info, rcond);
  344.       if (result_ok (info, rcond))
  345.     return tree_constant (result);
  346.     }
  347.  
  348.   int rank;
  349.   ComplexMatrix result = a.lssolve (b, info, rank);
  350.  
  351.   return tree_constant (result);
  352. }
  353.  
  354. /*
  355. ;;; Local Variables: ***
  356. ;;; mode: C++ ***
  357. ;;; page-delimiter: "^/\\*" ***
  358. ;;; End: ***
  359. */
  360.