Mercurial > octave-nkf
annotate liboctave/NLEqn.cc @ 7961:a5d1e27ee1f4 ss-3-1-51
3.1.51 snapshot
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Tue, 22 Jul 2008 11:40:48 -0400 |
| parents | 29980c6b8604 |
| children |
| rev | line source |
|---|---|
| 3 | 1 /* |
| 2 | |
| 7017 | 3 Copyright (C) 1993, 1994, 1995, 1996, 1997, 2000, 2002, 2003, 2004, |
| 4 2005, 2007 John W. Eaton | |
| 3 | 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 | |
| 7016 | 10 Free Software Foundation; either version 3 of the License, or (at your |
| 11 option) any later version. | |
| 3 | 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 | |
| 7016 | 19 along with Octave; see the file COPYING. If not, see |
| 20 <http://www.gnu.org/licenses/>. | |
| 3 | 21 |
| 22 */ | |
| 23 | |
| 238 | 24 #ifdef HAVE_CONFIG_H |
| 1192 | 25 #include <config.h> |
| 3 | 26 #endif |
| 27 | |
| 28 #include "NLEqn.h" | |
| 465 | 29 #include "dMatrix.h" |
| 1847 | 30 #include "f77-fcn.h" |
| 227 | 31 #include "lo-error.h" |
| 4180 | 32 #include "quit.h" |
| 3 | 33 |
| 5275 | 34 typedef octave_idx_type (*hybrd1_fcn_ptr) (octave_idx_type*, double*, double*, octave_idx_type*); |
| 3507 | 35 |
| 5275 | 36 typedef octave_idx_type (*hybrj1_fcn_ptr) (octave_idx_type*, double*, double*, double*, octave_idx_type*, octave_idx_type*); |
| 3507 | 37 |
| 38 extern "C" | |
| 4552 | 39 { |
| 40 F77_RET_T | |
| 5275 | 41 F77_FUNC (hybrd1, HYBRD1) (hybrd1_fcn_ptr, const octave_idx_type&, double*, |
| 42 double*, const double&, octave_idx_type&, double*, | |
| 43 const octave_idx_type&); | |
| 4552 | 44 |
| 45 | |
| 46 F77_RET_T | |
| 5275 | 47 F77_FUNC (hybrj1, HYBRJ1) (hybrj1_fcn_ptr, const octave_idx_type&, double*, |
| 48 double*, double*, const octave_idx_type&, const | |
| 49 double&, octave_idx_type&, double*, const octave_idx_type&); | |
| 4552 | 50 } |
| 3 | 51 |
| 1536 | 52 static NLFunc::nonlinear_fcn user_fun; |
| 53 static NLFunc::jacobian_fcn user_jac; | |
| 3 | 54 |
| 55 // error handling | |
| 56 | |
| 57 void | |
| 58 NLEqn::error (const char* msg) | |
| 59 { | |
| 227 | 60 (*current_liboctave_error_handler) ("fatal NLEqn error: %s", msg); |
| 3 | 61 } |
| 62 | |
| 63 // Other operations | |
| 64 | |
| 5275 | 65 octave_idx_type |
| 66 hybrd1_fcn (octave_idx_type *n, double *x, double *fvec, octave_idx_type *iflag) | |
| 3 | 67 { |
| 4180 | 68 BEGIN_INTERRUPT_WITH_EXCEPTIONS; |
| 69 | |
| 5275 | 70 octave_idx_type nn = *n; |
| 1536 | 71 ColumnVector tmp_f (nn); |
| 72 ColumnVector tmp_x (nn); | |
| 3 | 73 |
| 5275 | 74 for (octave_idx_type i = 0; i < nn; i++) |
| 3 | 75 tmp_x.elem (i) = x[i]; |
| 76 | |
| 77 tmp_f = (*user_fun) (tmp_x); | |
| 78 | |
| 253 | 79 if (tmp_f.length () == 0) |
| 80 *iflag = -1; | |
| 81 else | |
| 82 { | |
| 5275 | 83 for (octave_idx_type i = 0; i < nn; i++) |
| 253 | 84 fvec[i] = tmp_f.elem (i); |
| 85 } | |
| 3 | 86 |
| 4180 | 87 END_INTERRUPT_WITH_EXCEPTIONS; |
| 88 | |
| 3 | 89 return 0; |
| 90 } | |
| 91 | |
| 5275 | 92 octave_idx_type |
| 93 hybrj1_fcn (octave_idx_type *n, double *x, double *fvec, double *fjac, | |
| 94 octave_idx_type *ldfjac, octave_idx_type *iflag) | |
| 3 | 95 { |
| 4180 | 96 BEGIN_INTERRUPT_WITH_EXCEPTIONS; |
| 97 | |
| 5275 | 98 octave_idx_type nn = *n; |
| 1536 | 99 ColumnVector tmp_x (nn); |
| 3 | 100 |
| 5275 | 101 for (octave_idx_type i = 0; i < nn; i++) |
| 3 | 102 tmp_x.elem (i) = x[i]; |
| 103 | |
| 5275 | 104 octave_idx_type flag = *iflag; |
| 3 | 105 if (flag == 1) |
| 106 { | |
| 1536 | 107 ColumnVector tmp_f (nn); |
| 3 | 108 |
| 109 tmp_f = (*user_fun) (tmp_x); | |
| 110 | |
| 253 | 111 if (tmp_f.length () == 0) |
| 112 *iflag = -1; | |
| 113 else | |
| 114 { | |
| 5275 | 115 for (octave_idx_type i = 0; i < nn; i++) |
| 253 | 116 fvec[i] = tmp_f.elem (i); |
| 117 } | |
| 3 | 118 } |
| 119 else | |
| 120 { | |
| 121 Matrix tmp_fj (nn, nn); | |
| 122 | |
| 123 tmp_fj = (*user_jac) (tmp_x); | |
| 124 | |
| 253 | 125 if (tmp_fj.rows () == 0 || tmp_fj.columns () == 0) |
| 126 *iflag = -1; | |
| 127 else | |
| 128 { | |
| 5275 | 129 octave_idx_type ld = *ldfjac; |
| 130 for (octave_idx_type j = 0; j < nn; j++) | |
| 131 for (octave_idx_type i = 0; i < nn; i++) | |
| 253 | 132 fjac[j*ld+i] = tmp_fj.elem (i, j); |
| 133 } | |
| 3 | 134 } |
| 135 | |
| 4180 | 136 END_INTERRUPT_WITH_EXCEPTIONS; |
| 137 | |
| 3 | 138 return 0; |
| 139 } | |
| 140 | |
| 1536 | 141 ColumnVector |
| 5275 | 142 NLEqn::solve (octave_idx_type& info) |
| 3 | 143 { |
| 1536 | 144 ColumnVector retval; |
| 145 | |
| 5275 | 146 octave_idx_type n = x.capacity (); |
| 1878 | 147 |
| 3 | 148 if (n == 0) |
| 227 | 149 { |
| 150 error ("equation set not initialized"); | |
| 1536 | 151 return retval; |
| 227 | 152 } |
| 3 | 153 |
| 289 | 154 double tol = tolerance (); |
| 3 | 155 |
| 1937 | 156 retval = x; |
| 157 double *px = retval.fortran_vec (); | |
| 3 | 158 |
| 159 user_fun = fun; | |
| 160 user_jac = jac; | |
| 161 | |
| 465 | 162 if (jac) |
| 3 | 163 { |
| 1937 | 164 Array<double> fvec (n); |
| 165 double *pfvec = fvec.fortran_vec (); | |
| 166 | |
| 5275 | 167 octave_idx_type lwa = (n*(n+13))/2; |
| 1937 | 168 Array<double> wa (lwa); |
| 169 double *pwa = wa.fortran_vec (); | |
| 3 | 170 |
| 1937 | 171 Array<double> fjac (n*n); |
| 172 double *pfjac = fjac.fortran_vec (); | |
| 3 | 173 |
| 1937 | 174 F77_XFCN (hybrj1, HYBRJ1, (hybrj1_fcn, n, px, pfvec, pfjac, n, |
| 175 tol, info, pwa, lwa)); | |
| 176 | |
| 3971 | 177 solution_status = info; |
| 178 | |
|
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7279
diff
changeset
|
179 fval = ColumnVector (fvec); |
| 3 | 180 } |
| 465 | 181 else |
| 182 { | |
| 1937 | 183 Array<double> fvec (n); |
| 184 double *pfvec = fvec.fortran_vec (); | |
| 465 | 185 |
| 5275 | 186 octave_idx_type lwa = (n*(3*n+13))/2; |
| 1937 | 187 Array<double> wa (lwa); |
| 188 double *pwa = wa.fortran_vec (); | |
| 465 | 189 |
| 1937 | 190 F77_XFCN (hybrd1, HYBRD1, (hybrd1_fcn, n, px, pfvec, tol, info, |
| 191 pwa, lwa)); | |
| 192 | |
| 3971 | 193 solution_status = info; |
| 194 | |
|
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7279
diff
changeset
|
195 fval = ColumnVector (fvec); |
| 465 | 196 } |
| 3 | 197 |
| 198 return retval; | |
| 199 } | |
| 200 | |
| 3971 | 201 std::string |
| 202 NLEqn::error_message (void) const | |
| 203 { | |
| 204 std::string retval; | |
| 205 | |
| 206 std::string prefix; | |
| 207 | |
| 5275 | 208 octave_idx_type info = solution_status; |
| 3971 | 209 if (info < 0) |
| 210 info = -info; | |
| 211 | |
| 212 switch (info) | |
| 213 { | |
| 214 case 0: | |
| 215 retval = "improper input parameters"; | |
| 216 break; | |
| 217 | |
| 218 case 1: | |
| 219 retval = "solution converged within specified tolerance"; | |
| 220 break; | |
| 221 | |
| 222 case 2: | |
| 223 retval = "number of function calls exceeded limit"; | |
| 224 break; | |
| 225 | |
| 226 case 3: | |
| 227 retval = "no further improvement possible (tolerance may be too small)"; | |
| 228 break; | |
| 229 | |
| 230 case 4: | |
| 231 retval = "iteration is not making good progress"; | |
| 232 break; | |
| 233 | |
| 234 default: | |
| 235 retval = "unknown error state"; | |
| 236 break; | |
| 237 } | |
| 238 | |
| 239 if (solution_status < 0) | |
| 240 retval = std::string ("user requested termination: ") + retval; | |
| 241 | |
| 242 return retval; | |
| 243 } | |
| 244 | |
| 3 | 245 /* |
| 246 ;;; Local Variables: *** | |
| 247 ;;; mode: C++ *** | |
| 248 ;;; End: *** | |
| 249 */ |