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 / liboctave / CollocWt.cc < prev    next >
C/C++ Source or Header  |  1996-09-28  |  6KB  |  353 lines
  1. // CollocWt.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 <iostream.h>
  29.  
  30. #include "CollocWt.h"
  31. #include "f77-uscore.h"
  32. #include "lo-error.h"
  33.  
  34. extern "C"
  35. {
  36.   int F77_FCN (jcobi) (int*, int*, int*, int*, double*, double*,
  37.                double*, double*, double*, double*); 
  38.  
  39.   int F77_FCN (dfopr) (int*, int*, int*, int*, int*, int*,
  40.                double*, double*, double*, double*, double*);
  41. }
  42.  
  43. // Error handling.
  44.  
  45. void
  46. CollocWt::error (const char* msg)
  47. {
  48.   (*current_liboctave_error_handler) ("fatal CollocWt error: %s", msg);
  49. }
  50.  
  51. CollocWt::CollocWt (void)
  52. {
  53.   n = 0;
  54.   inc_left = 0;
  55.   inc_right = 0;
  56.   lb = 0.0;
  57.   rb = 1.0;
  58.  
  59.   Alpha = 0.0;
  60.   Beta = 0.0;
  61.  
  62.   initialized = 0;
  63. }
  64.  
  65. CollocWt::CollocWt (int nc, int il, int ir)
  66. {
  67.   n = nc;
  68.   inc_left = il;
  69.   inc_right = ir;
  70.   lb = 0.0;
  71.   rb = 1.0;
  72.  
  73.   Alpha = 0.0;
  74.   Beta = 0.0;
  75.  
  76.   initialized = 0;
  77. }
  78.  
  79. CollocWt::CollocWt (int nc, int ir, int il, double l, double r)
  80. {
  81.   n = nc;
  82.   inc_left = il;
  83.   inc_right = ir;
  84.   lb = l;
  85.   rb = r;
  86.  
  87.   Alpha = 0.0;
  88.   Beta = 0.0;
  89.  
  90.   initialized = 0;
  91. }
  92.  
  93. CollocWt::CollocWt (int nc, double a, double b, int il, int ir)
  94. {
  95.   n = nc;
  96.   inc_left = il;
  97.   inc_right = ir;
  98.   lb = 0.0;
  99.   rb = 1.0;
  100.  
  101.   Alpha = a;
  102.   Beta = b;
  103.  
  104.   initialized = 0;
  105. }
  106.  
  107. CollocWt::CollocWt (int nc, double a, double b, int ir, int il,
  108.             double l, double r)  
  109. {
  110.   n = nc;
  111.   inc_left = il;
  112.   inc_right = ir;
  113.   lb = l;
  114.   rb = r;
  115.  
  116.   Alpha = a;
  117.   Beta = b;
  118.  
  119.   initialized = 0;
  120. }
  121.  
  122. CollocWt::CollocWt (const CollocWt& a)
  123. {
  124.   n = a.n;
  125.   inc_left = a.inc_left;
  126.   inc_right = a.inc_right;
  127.   lb = a.lb;
  128.   rb = a.rb;
  129.   r = a.r;
  130.   q = a.q;
  131.   A = a.A;
  132.   B = a.B;
  133.  
  134.   nt = n + inc_left + inc_right;
  135.  
  136.   initialized = a.initialized;
  137. }
  138.  
  139. CollocWt&
  140. CollocWt::operator = (const CollocWt& a)
  141. {
  142.   n = a.n;
  143.   inc_left = a.inc_left;
  144.   inc_right = a.inc_right;
  145.   lb = a.lb;
  146.   rb = a.rb;
  147.   r = a.r;
  148.   q = a.q;
  149.   A = a.A;
  150.   B = a.B;
  151.  
  152.   nt = a.nt;
  153.  
  154.   initialized = a.initialized;
  155.  
  156.   return *this;
  157. }
  158.  
  159. CollocWt&
  160. CollocWt::resize (int ncol)
  161. {
  162.   n = ncol;
  163.   initialized = 0;
  164.   return *this;
  165. }
  166.  
  167. CollocWt&
  168. CollocWt::add_left (void)
  169. {
  170.   inc_left = 1;
  171.   initialized = 0;
  172.   return *this;
  173. }
  174.  
  175. CollocWt&
  176. CollocWt::delete_left (void)
  177. {
  178.   inc_left = 0;
  179.   initialized = 0;
  180.   return *this;
  181. }
  182.  
  183. CollocWt&
  184. CollocWt::set_left (double val)
  185. {
  186.   if (val >= rb)
  187.     {
  188.       error ("left bound greater than right bound");
  189.       return *this;
  190.     }
  191.  
  192.   lb = val;
  193.   initialized = 0;
  194.   return *this;
  195. }
  196.  
  197. CollocWt&
  198. CollocWt::add_right (void)
  199. {
  200.   inc_right = 1;
  201.   initialized = 0;
  202.   return *this;
  203. }
  204.  
  205. CollocWt&
  206. CollocWt::delete_right (void)
  207. {
  208.   inc_right = 0;
  209.   initialized = 0;
  210.   return *this;
  211. }
  212.  
  213. CollocWt&
  214. CollocWt::set_right (double val)
  215. {
  216.   if (val <= lb)
  217.     {
  218.       error ("right bound less than left bound");
  219.       return *this;
  220.     }
  221.  
  222.   rb = val;
  223.   initialized = 0;
  224.   return *this;
  225. }
  226.  
  227. CollocWt&
  228. CollocWt::set_alpha (double val)
  229. {
  230.   Alpha = val;
  231.   initialized = 0;
  232.   return *this;
  233. }
  234.  
  235. CollocWt&
  236. CollocWt::set_beta (double val)
  237. {
  238.   Beta = val;
  239.   initialized = 0;
  240.   return *this;
  241. }
  242.  
  243. void
  244. CollocWt::init (void)
  245. {
  246. // Check for possible errors.
  247.  
  248.   double wid = rb - lb;
  249.   if (wid <= 0.0)
  250.     {
  251.       error ("width less than or equal to zero");
  252.       return;
  253.     }
  254.  
  255.   nt = n + inc_left + inc_right;
  256.   if (nt < 0)
  257.     {
  258.       error ("total number of collocation points less than zero");
  259.       return;
  260.     }
  261.   else if (nt == 0)
  262.     return;
  263.  
  264.   double *dif1 = new double [nt];
  265.   double *dif2 = new double [nt];
  266.   double *dif3 = new double [nt];
  267.   double *vect = new double [nt];
  268.  
  269.   r.resize (nt);
  270.   q.resize (nt);
  271.   A.resize (nt, nt);
  272.   B.resize (nt, nt);
  273.  
  274.   double *pr = r.fortran_vec ();
  275.  
  276. // Compute roots.
  277.  
  278.   F77_FCN (jcobi) (&nt, &n, &inc_left, &inc_right, &Alpha, &Beta,
  279.            dif1, dif2, dif3, pr);
  280.  
  281.   int id;
  282.   int i, j;
  283.  
  284. // First derivative weights.
  285.  
  286.   id = 1;
  287.   for (i = 1; i <= nt; i++)
  288.     {
  289.       F77_FCN (dfopr) (&nt, &n, &inc_left, &inc_right, &i, &id, dif1,
  290.                dif2, dif3, pr, vect); 
  291.  
  292.       for (j = 0; j < nt; j++)
  293.     A (i-1, j) = vect[j];
  294.     }
  295.  
  296. // Second derivative weights.
  297.  
  298.   id = 2;
  299.   for (i = 1; i <= nt; i++)
  300.     {
  301.       F77_FCN (dfopr) (&nt, &n, &inc_left, &inc_right, &i, &id, dif1,
  302.                dif2, dif3, pr, vect); 
  303.  
  304.       for (j = 0; j < nt; j++)
  305.     B (i-1, j) = vect[j];
  306.     }
  307.  
  308. // Gaussian quadrature weights.
  309.  
  310.   id = 3;
  311.   double *pq = q.fortran_vec ();
  312.   F77_FCN (dfopr) (&nt, &n, &inc_left, &inc_right, &i, &id, dif1,
  313.            dif2, dif3, pr, pq);
  314.  
  315.   delete [] dif1;
  316.   delete [] dif2;
  317.   delete [] dif3;
  318.   delete [] vect;
  319.  
  320.   initialized = 1;
  321. }
  322.  
  323. ostream&
  324. operator << (ostream& os, const CollocWt& a)
  325. {
  326.   if (a.left_included ())
  327.     os << "left  boundary is included\n";
  328.   else
  329.     os << "left  boundary is not included\n";
  330.  
  331.   if (a.right_included ())
  332.     os << "right boundary is included\n";
  333.   else
  334.     os << "right boundary is not included\n";
  335.  
  336.   os << "\n";
  337.  
  338.   os << a.Alpha << " " << a.Beta << "\n\n"
  339.      << a.r << "\n\n"
  340.      << a.q << "\n\n"
  341.      << a.A << "\n"
  342.      << a.B << "\n";
  343.  
  344.   return os;
  345. }
  346.  
  347. /*
  348. ;;; Local Variables: ***
  349. ;;; mode: C++ ***
  350. ;;; page-delimiter: "^/\\*" ***
  351. ;;; End: ***
  352. */
  353.