00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #if defined (__GNUG__) && defined (USE_PRAGMA_INTERFACE_IMPLEMENTATION)
00024 #pragma implementation
00025 #endif
00026
00027 #ifdef HAVE_CONFIG_H
00028 #include <config.h>
00029 #endif
00030
00031 #include <iostream>
00032
00033 #include "CollocWt.h"
00034 #include "f77-fcn.h"
00035 #include "lo-error.h"
00036
00037 extern "C"
00038 {
00039 F77_RET_T
00040 F77_FUNC (jcobi, JCOBI) (int&, int&, int&, int&, double&,
00041 double&, double*, double*, double*,
00042 double*);
00043
00044 F77_RET_T
00045 F77_FUNC (dfopr, DFOPR) (int&, int&, int&, int&, int&, int&,
00046 double*, double*, double*, double*,
00047 double*);
00048 }
00049
00050
00051
00052 void
00053 CollocWt::error (const char* msg)
00054 {
00055 (*current_liboctave_error_handler) ("fatal CollocWt error: %s", msg);
00056 }
00057
00058 CollocWt&
00059 CollocWt::set_left (double val)
00060 {
00061 if (val >= rb)
00062 {
00063 error ("left bound greater than right bound");
00064 return *this;
00065 }
00066
00067 lb = val;
00068 initialized = 0;
00069 return *this;
00070 }
00071
00072 CollocWt&
00073 CollocWt::set_right (double val)
00074 {
00075 if (val <= lb)
00076 {
00077 error ("right bound less than left bound");
00078 return *this;
00079 }
00080
00081 rb = val;
00082 initialized = 0;
00083 return *this;
00084 }
00085
00086 void
00087 CollocWt::init (void)
00088 {
00089
00090
00091 double wid = rb - lb;
00092 if (wid <= 0.0)
00093 {
00094 error ("width less than or equal to zero");
00095 return;
00096 }
00097
00098 int nt = n + inc_left + inc_right;
00099
00100 if (nt < 0)
00101 {
00102 error ("total number of collocation points less than zero");
00103 return;
00104 }
00105 else if (nt == 0)
00106 return;
00107
00108 Array<double> dif1 (nt);
00109 double *pdif1 = dif1.fortran_vec ();
00110
00111 Array<double> dif2 (nt);
00112 double *pdif2 = dif2.fortran_vec ();
00113
00114 Array<double> dif3 (nt);
00115 double *pdif3 = dif3.fortran_vec ();
00116
00117 Array<double> vect (nt);
00118 double *pvect = vect.fortran_vec ();
00119
00120 r.resize (nt);
00121 q.resize (nt);
00122 A.resize (nt, nt);
00123 B.resize (nt, nt);
00124
00125 double *pr = r.fortran_vec ();
00126
00127
00128
00129 F77_FUNC (jcobi, JCOBI) (nt, n, inc_left, inc_right, Alpha, Beta,
00130 pdif1, pdif2, pdif3, pr);
00131
00132 int id;
00133
00134
00135
00136 id = 1;
00137 for (int i = 1; i <= nt; i++)
00138 {
00139 F77_FUNC (dfopr, DFOPR) (nt, n, inc_left, inc_right, i, id, pdif1,
00140 pdif2, pdif3, pr, pvect);
00141
00142 for (int j = 0; j < nt; j++)
00143 A (i-1, j) = vect.elem (j);
00144 }
00145
00146
00147
00148 id = 2;
00149 for (int i = 1; i <= nt; i++)
00150 {
00151 F77_FUNC (dfopr, DFOPR) (nt, n, inc_left, inc_right, i, id, pdif1,
00152 pdif2, pdif3, pr, pvect);
00153
00154 for (int j = 0; j < nt; j++)
00155 B (i-1, j) = vect.elem (j);
00156 }
00157
00158
00159
00160 id = 3;
00161 double *pq = q.fortran_vec ();
00162 F77_FUNC (dfopr, DFOPR) (nt, n, inc_left, inc_right, id, id, pdif1,
00163 pdif2, pdif3, pr, pq);
00164
00165 initialized = 1;
00166 }
00167
00168 std::ostream&
00169 operator << (std::ostream& os, const CollocWt& a)
00170 {
00171 if (a.left_included ())
00172 os << "left boundary is included\n";
00173 else
00174 os << "left boundary is not included\n";
00175
00176 if (a.right_included ())
00177 os << "right boundary is included\n";
00178 else
00179 os << "right boundary is not included\n";
00180
00181 os << "\n";
00182
00183 os << a.Alpha << " " << a.Beta << "\n\n"
00184 << a.r << "\n\n"
00185 << a.q << "\n\n"
00186 << a.A << "\n"
00187 << a.B << "\n";
00188
00189 return os;
00190 }
00191
00192
00193
00194
00195
00196