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 "NLEqn.h"
00032 #include "dMatrix.h"
00033 #include "f77-fcn.h"
00034 #include "lo-error.h"
00035 #include "quit.h"
00036
00037 typedef int (*hybrd1_fcn_ptr) (int*, double*, double*, int*);
00038
00039 typedef int (*hybrj1_fcn_ptr) (int*, double*, double*, double*, int*, int*);
00040
00041 extern "C"
00042 {
00043 F77_RET_T
00044 F77_FUNC (hybrd1, HYBRD1) (hybrd1_fcn_ptr, const int&, double*,
00045 double*, const double&, int&, double*,
00046 const int&);
00047
00048
00049 F77_RET_T
00050 F77_FUNC (hybrj1, HYBRJ1) (hybrj1_fcn_ptr, const int&, double*,
00051 double*, double*, const int&, const
00052 double&, int&, double*, const int&);
00053 }
00054
00055 static NLFunc::nonlinear_fcn user_fun;
00056 static NLFunc::jacobian_fcn user_jac;
00057
00058
00059
00060 void
00061 NLEqn::error (const char* msg)
00062 {
00063 (*current_liboctave_error_handler) ("fatal NLEqn error: %s", msg);
00064 }
00065
00066
00067
00068 int
00069 hybrd1_fcn (int *n, double *x, double *fvec, int *iflag)
00070 {
00071 BEGIN_INTERRUPT_WITH_EXCEPTIONS;
00072
00073 int nn = *n;
00074 ColumnVector tmp_f (nn);
00075 ColumnVector tmp_x (nn);
00076
00077 for (int i = 0; i < nn; i++)
00078 tmp_x.elem (i) = x[i];
00079
00080 tmp_f = (*user_fun) (tmp_x);
00081
00082 if (tmp_f.length () == 0)
00083 *iflag = -1;
00084 else
00085 {
00086 for (int i = 0; i < nn; i++)
00087 fvec[i] = tmp_f.elem (i);
00088 }
00089
00090 END_INTERRUPT_WITH_EXCEPTIONS;
00091
00092 return 0;
00093 }
00094
00095 int
00096 hybrj1_fcn (int *n, double *x, double *fvec, double *fjac,
00097 int *ldfjac, int *iflag)
00098 {
00099 BEGIN_INTERRUPT_WITH_EXCEPTIONS;
00100
00101 int nn = *n;
00102 ColumnVector tmp_x (nn);
00103
00104 for (int i = 0; i < nn; i++)
00105 tmp_x.elem (i) = x[i];
00106
00107 int flag = *iflag;
00108 if (flag == 1)
00109 {
00110 ColumnVector tmp_f (nn);
00111
00112 tmp_f = (*user_fun) (tmp_x);
00113
00114 if (tmp_f.length () == 0)
00115 *iflag = -1;
00116 else
00117 {
00118 for (int i = 0; i < nn; i++)
00119 fvec[i] = tmp_f.elem (i);
00120 }
00121 }
00122 else
00123 {
00124 Matrix tmp_fj (nn, nn);
00125
00126 tmp_fj = (*user_jac) (tmp_x);
00127
00128 if (tmp_fj.rows () == 0 || tmp_fj.columns () == 0)
00129 *iflag = -1;
00130 else
00131 {
00132 int ld = *ldfjac;
00133 for (int j = 0; j < nn; j++)
00134 for (int i = 0; i < nn; i++)
00135 fjac[j*ld+i] = tmp_fj.elem (i, j);
00136 }
00137 }
00138
00139 END_INTERRUPT_WITH_EXCEPTIONS;
00140
00141 return 0;
00142 }
00143
00144 ColumnVector
00145 NLEqn::solve (int& info)
00146 {
00147 ColumnVector retval;
00148
00149 int n = x.capacity ();
00150
00151 if (n == 0)
00152 {
00153 error ("equation set not initialized");
00154 return retval;
00155 }
00156
00157 double tol = tolerance ();
00158
00159 retval = x;
00160 double *px = retval.fortran_vec ();
00161
00162 user_fun = fun;
00163 user_jac = jac;
00164
00165 if (jac)
00166 {
00167 Array<double> fvec (n);
00168 double *pfvec = fvec.fortran_vec ();
00169
00170 int lwa = (n*(n+13))/2;
00171 Array<double> wa (lwa);
00172 double *pwa = wa.fortran_vec ();
00173
00174 Array<double> fjac (n*n);
00175 double *pfjac = fjac.fortran_vec ();
00176
00177 F77_XFCN (hybrj1, HYBRJ1, (hybrj1_fcn, n, px, pfvec, pfjac, n,
00178 tol, info, pwa, lwa));
00179
00180 solution_status = info;
00181
00182 if (f77_exception_encountered)
00183 (*current_liboctave_error_handler) ("unrecoverable error in hybrj1");
00184 }
00185 else
00186 {
00187 Array<double> fvec (n);
00188 double *pfvec = fvec.fortran_vec ();
00189
00190 int lwa = (n*(3*n+13))/2;
00191 Array<double> wa (lwa);
00192 double *pwa = wa.fortran_vec ();
00193
00194 F77_XFCN (hybrd1, HYBRD1, (hybrd1_fcn, n, px, pfvec, tol, info,
00195 pwa, lwa));
00196
00197 solution_status = info;
00198
00199 if (f77_exception_encountered)
00200 (*current_liboctave_error_handler) ("unrecoverable error in hybrd1");
00201 }
00202
00203 return retval;
00204 }
00205
00206 std::string
00207 NLEqn::error_message (void) const
00208 {
00209 std::string retval;
00210
00211 std::string prefix;
00212
00213 int info = solution_status;
00214 if (info < 0)
00215 info = -info;
00216
00217 switch (info)
00218 {
00219 case 0:
00220 retval = "improper input parameters";
00221 break;
00222
00223 case 1:
00224 retval = "solution converged within specified tolerance";
00225 break;
00226
00227 case 2:
00228 retval = "number of function calls exceeded limit";
00229 break;
00230
00231 case 3:
00232 retval = "no further improvement possible (tolerance may be too small)";
00233 break;
00234
00235 case 4:
00236 retval = "iteration is not making good progress";
00237 break;
00238
00239 default:
00240 retval = "unknown error state";
00241 break;
00242 }
00243
00244 if (solution_status < 0)
00245 retval = std::string ("user requested termination: ") + retval;
00246
00247 return retval;
00248 }
00249
00250
00251
00252
00253
00254