1 /* ////////////////////////////////////////////////////////////////////////////////
2 //
3 // Matlab MEX file for the Levenberg - Marquardt minimization algorithm
4 // Copyright (C) 2007 Manolis Lourakis (lourakis at ics forth gr)
5 // Institute of Computer Science, Foundation for Research & Technology - Hellas
6 // Heraklion, Crete, Greece.
7 //
8 // This program is free software; you can redistribute it and/or modify
9 // it under the terms of the GNU General Public License as published by
10 // the Free Software Foundation; either version 2 of the License, or
11 // (at your option) any later version.
12 //
13 // This program is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 // GNU General Public License for more details.
17 //
18 //////////////////////////////////////////////////////////////////////////////// */
19
20 #include <stdio.h>
21 #include <stdlib.h>
22 #include <stdarg.h>
23 #include <math.h>
24 #include <string.h>
25 #include <ctype.h>
26
27 #include <lm.h>
28
29 #include <mex.h>
30
31 /**
32 #define DEBUG
33 **/
34
35 #ifndef HAVE_LAPACK
36 #ifdef _MSC_VER
37 #pragma message("LAPACK not available, certain functionalities cannot be compiled!")
38 #else
39 #warning LAPACK not available, certain functionalities cannot be compiled
40 #endif /* _MSC_VER */
41 #endif /* HAVE_LAPACK */
42
43 #define __MAX__(A, B) ((A)>=(B)? (A) : (B))
44
45 #define MIN_UNCONSTRAINED 0
46 #define MIN_CONSTRAINED_BC 1
47 #define MIN_CONSTRAINED_LEC 2
48 #define MIN_CONSTRAINED_BLEC 3
49
50 struct mexdata {
51 /* matlab names of the fitting function & its Jacobian */
52 char *fname, *jacname;
53
54 /* binary flags specifying if input p0 is a row or column vector */
55 int isrow_p0;
56
57 /* rhs args to be passed to matlab. rhs[0] is reserved for
58 * passing the parameter vector. If present, problem-specific
59 * data are passed in rhs[1], rhs[2], etc
60 */
61 mxArray **rhs;
62 int nrhs; /* >= 1 */
63 };
64
65 /* display printf-style error messages in matlab */
matlabFmtdErrMsgTxt(char * fmt,...)66 static void matlabFmtdErrMsgTxt(char *fmt, ...)
67 {
68 char buf[256];
69 va_list args;
70
71 va_start(args, fmt);
72 vsprintf(buf, fmt, args);
73 va_end(args);
74
75 mexErrMsgTxt(buf);
76 }
77
78 /* display printf-style warning messages in matlab */
matlabFmtdWarnMsgTxt(char * fmt,...)79 static void matlabFmtdWarnMsgTxt(char *fmt, ...)
80 {
81 char buf[256];
82 va_list args;
83
84 va_start(args, fmt);
85 vsprintf(buf, fmt, args);
86 va_end(args);
87
88 mexWarnMsgTxt(buf);
89 }
90
func(double * p,double * hx,int m,int n,void * adata)91 static void func(double *p, double *hx, int m, int n, void *adata)
92 {
93 mxArray *lhs[1];
94 double *mp, *mx;
95 register int i;
96 struct mexdata *dat=(struct mexdata *)adata;
97
98 /* prepare to call matlab */
99 mp=mxGetPr(dat->rhs[0]);
100 for(i=0; i<m; ++i)
101 mp[i]=p[i];
102
103 /* invoke matlab */
104 mexCallMATLAB(1, lhs, dat->nrhs, dat->rhs, dat->fname);
105
106 /* copy back results & cleanup */
107 mx=mxGetPr(lhs[0]);
108 for(i=0; i<n; ++i)
109 hx[i]=mx[i];
110
111 /* delete the matrix created by matlab */
112 mxDestroyArray(lhs[0]);
113 }
114
jacfunc(double * p,double * j,int m,int n,void * adata)115 static void jacfunc(double *p, double *j, int m, int n, void *adata)
116 {
117 mxArray *lhs[1];
118 double *mp;
119 double *mj;
120 register int i, k;
121 struct mexdata *dat=(struct mexdata *)adata;
122
123 /* prepare to call matlab */
124 mp=mxGetPr(dat->rhs[0]);
125 for(i=0; i<m; ++i)
126 mp[i]=p[i];
127
128 /* invoke matlab */
129 mexCallMATLAB(1, lhs, dat->nrhs, dat->rhs, dat->jacname);
130
131 /* copy back results & cleanup. Note that the nxm Jacobian
132 * computed by matlab should be transposed so that
133 * levmar gets it in row major, as expected
134 */
135 mj=mxGetPr(lhs[0]);
136 for(i=0; i<n; ++i)
137 for(k=0; k<m; ++k)
138 j[i*m+k]=mj[i+k*n];
139
140 /* delete the matrix created by matlab */
141 mxDestroyArray(lhs[0]);
142 }
143
144 /* matlab matrices are in column-major, this routine converts them to row major for levmar */
getTranspose(mxArray * Am)145 static double *getTranspose(mxArray *Am)
146 {
147 int m, n;
148 double *At, *A;
149 register int i, j;
150
151 m=mxGetM(Am);
152 n=mxGetN(Am);
153 A=mxGetPr(Am);
154 At=mxMalloc(m*n*sizeof(double));
155
156 for(i=0; i<m; i++)
157 for(j=0; j<n; j++)
158 At[i*n+j]=A[i+j*m];
159
160 return At;
161 }
162
163 /* check the supplied matlab function and its Jacobian. Returns 1 on error, 0 otherwise */
checkFuncAndJacobian(double * p,int m,int n,int havejac,struct mexdata * dat)164 static int checkFuncAndJacobian(double *p, int m, int n, int havejac, struct mexdata *dat)
165 {
166 mxArray *lhs[1];
167 register int i;
168 int ret=0;
169 double *mp;
170
171 mexSetTrapFlag(1); /* handle errors in the MEX-file */
172
173 mp=mxGetPr(dat->rhs[0]);
174 for(i=0; i<m; ++i)
175 mp[i]=p[i];
176
177 /* attempt to call the supplied func */
178 i=mexCallMATLAB(1, lhs, dat->nrhs, dat->rhs, dat->fname);
179 if(i){
180 fprintf(stderr, "levmar: error calling '%s'.\n", dat->fname);
181 ret=1;
182 }
183 else if(!mxIsDouble(lhs[0]) || mxIsComplex(lhs[0]) || !(mxGetM(lhs[0])==1 || mxGetN(lhs[0])==1) ||
184 __MAX__(mxGetM(lhs[0]), mxGetN(lhs[0]))!=n){
185 fprintf(stderr, "levmar: '%s' should produce a real vector with %d elements (got %d).\n",
186 dat->fname, m, __MAX__(mxGetM(lhs[0]), mxGetN(lhs[0])));
187 ret=1;
188 }
189 /* delete the matrix created by matlab */
190 mxDestroyArray(lhs[0]);
191
192 if(havejac){
193 /* attempt to call the supplied jac */
194 i=mexCallMATLAB(1, lhs, dat->nrhs, dat->rhs, dat->jacname);
195 if(i){
196 fprintf(stderr, "levmar: error calling '%s'.\n", dat->jacname);
197 ret=1;
198 }
199 else if(!mxIsDouble(lhs[0]) || mxIsComplex(lhs[0]) || mxGetM(lhs[0])!=n || mxGetN(lhs[0])!=m){
200 fprintf(stderr, "levmar: '%s' should produce a real %dx%d matrix (got %dx%d).\n",
201 dat->jacname, n, m, mxGetM(lhs[0]), mxGetN(lhs[0]));
202 ret=1;
203 }
204 else if(mxIsSparse(lhs[0])){
205 fprintf(stderr, "levmar: '%s' should produce a real dense matrix (got a sparse one).\n", dat->jacname);
206 ret=1;
207 }
208 /* delete the matrix created by matlab */
209 mxDestroyArray(lhs[0]);
210 }
211
212 mexSetTrapFlag(0); /* on error terminate the MEX-file and return control to the MATLAB prompt */
213
214 return ret;
215 }
216
217
218 /*
219 [ret, p, info, covar]=levmar_der (f, j, p0, x, itmax, opts, 'unc' ...)
220 [ret, p, info, covar]=levmar_bc (f, j, p0, x, itmax, opts, 'bc', lb, ub, ...)
221 [ret, p, info, covar]=levmar_lec (f, j, p0, x, itmax, opts, 'lec', A, b, ...)
222 [ret, p, info, covar]=levmar_blec(f, j, p0, x, itmax, opts, 'blec', lb, ub, A, b, wghts, ...)
223 */
224
mexFunction(int nlhs,mxArray * plhs[],int nrhs,const mxArray * Prhs[])225 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *Prhs[])
226 {
227 register int i;
228 register double *pdbl;
229 mxArray **prhs=(mxArray **)&Prhs[0], *At;
230 struct mexdata mdata;
231 int len, status;
232 double *p, *p0, *ret, *x;
233 int m, n, havejac, Arows, itmax, nopts, mintype, nextra;
234 double opts[LM_OPTS_SZ]={LM_INIT_MU, LM_STOP_THRESH, LM_STOP_THRESH, LM_STOP_THRESH, LM_DIFF_DELTA};
235 double info[LM_INFO_SZ];
236 double *lb=NULL, *ub=NULL, *A=NULL, *b=NULL, *wghts=NULL, *covar=NULL;
237
238 /* parse input args; start by checking their number */
239 if((nrhs<5))
240 matlabFmtdErrMsgTxt("levmar: at least 5 input arguments required (got %d).", nrhs);
241 if(nlhs>4)
242 matlabFmtdErrMsgTxt("levmar: too many output arguments (max. 4, got %d).", nlhs);
243 else if(nlhs<2)
244 matlabFmtdErrMsgTxt("levmar: too few output arguments (min. 2, got %d).", nlhs);
245
246 /* note that in order to accommodate optional args, prhs & nrhs are adjusted accordingly below */
247
248 /** func **/
249 /* first argument must be a string , i.e. a char row vector */
250 if(mxIsChar(prhs[0])!=1)
251 mexErrMsgTxt("levmar: first argument must be a string.");
252 if(mxGetM(prhs[0])!=1)
253 mexErrMsgTxt("levmar: first argument must be a string (i.e. char row vector).");
254 /* store supplied name */
255 len=mxGetN(prhs[0])+1;
256 mdata.fname=mxCalloc(len, sizeof(char));
257 status=mxGetString(prhs[0], mdata.fname, len);
258 if(status!=0)
259 mexErrMsgTxt("levmar: not enough space. String is truncated.");
260
261 /** jac (optional) **/
262 /* check whether second argument is a string */
263 if(mxIsChar(prhs[1])==1){
264 if(mxGetM(prhs[1])!=1)
265 mexErrMsgTxt("levmar: second argument must be a string (i.e. row vector).");
266 /* store supplied name */
267 len=mxGetN(prhs[1])+1;
268 mdata.jacname=mxCalloc(len, sizeof(char));
269 status=mxGetString(prhs[1], mdata.jacname, len);
270 if(status!=0)
271 mexErrMsgTxt("levmar: not enough space. String is truncated.");
272 havejac=1;
273
274 ++prhs;
275 --nrhs;
276 }
277 else{
278 mdata.jacname=NULL;
279 havejac=0;
280 }
281
282 #ifdef DEBUG
283 fflush(stderr);
284 fprintf(stderr, "LEVMAR: %s analytic Jacobian\n", havejac? "with" : "no");
285 #endif /* DEBUG */
286
287 /* CHECK
288 if(!mxIsDouble(prhs[1]) || mxIsComplex(prhs[1]) || !(mxGetM(prhs[1])==1 && mxGetN(prhs[1])==1))
289 */
290
291 /** p0 **/
292 /* the second required argument must be a real row or column vector */
293 if(!mxIsDouble(prhs[1]) || mxIsComplex(prhs[1]) || !(mxGetM(prhs[1])==1 || mxGetN(prhs[1])==1))
294 mexErrMsgTxt("levmar: p0 must be a real vector.");
295 p0=mxGetPr(prhs[1]);
296 /* determine if we have a row or column vector and retrieve its
297 * size, i.e. the number of parameters
298 */
299 if(mxGetM(prhs[1])==1){
300 m=mxGetN(prhs[1]);
301 mdata.isrow_p0=1;
302 }
303 else{
304 m=mxGetM(prhs[1]);
305 mdata.isrow_p0=0;
306 }
307 /* copy input parameter vector to avoid destroying it */
308 p=mxMalloc(m*sizeof(double));
309 for(i=0; i<m; ++i)
310 p[i]=p0[i];
311
312 /** x **/
313 /* the third required argument must be a real row or column vector */
314 if(!mxIsDouble(prhs[2]) || mxIsComplex(prhs[2]) || !(mxGetM(prhs[2])==1 || mxGetN(prhs[2])==1))
315 mexErrMsgTxt("levmar: x must be a real vector.");
316 x=mxGetPr(prhs[2]);
317 n=__MAX__(mxGetM(prhs[2]), mxGetN(prhs[2]));
318
319 /** itmax **/
320 /* the fourth required argument must be a scalar */
321 if(!mxIsDouble(prhs[3]) || mxIsComplex(prhs[3]) || mxGetM(prhs[3])!=1 || mxGetN(prhs[3])!=1)
322 mexErrMsgTxt("levmar: itmax must be a scalar.");
323 itmax=(int)mxGetScalar(prhs[3]);
324
325 /** opts **/
326 /* the fifth required argument must be a real row or column vector */
327 if(!mxIsDouble(prhs[4]) || mxIsComplex(prhs[4]) || (!(mxGetM(prhs[4])==1 || mxGetN(prhs[4])==1) &&
328 !(mxGetM(prhs[4])==0 && mxGetN(prhs[4])==0)))
329 mexErrMsgTxt("levmar: opts must be a real vector.");
330 pdbl=mxGetPr(prhs[4]);
331 nopts=__MAX__(mxGetM(prhs[4]), mxGetN(prhs[4]));
332 if(nopts!=0){ /* if opts==[], nothing needs to be done and the defaults are used */
333 if(nopts>LM_OPTS_SZ)
334 matlabFmtdErrMsgTxt("levmar: opts must have at most %d elements, got %d.", LM_OPTS_SZ, nopts);
335 else if(nopts<((havejac)? LM_OPTS_SZ-1 : LM_OPTS_SZ))
336 matlabFmtdWarnMsgTxt("levmar: only the %d first elements of opts specified, remaining set to defaults.", nopts);
337 for(i=0; i<nopts; ++i)
338 opts[i]=pdbl[i];
339 }
340 #ifdef DEBUG
341 else{
342 fflush(stderr);
343 fprintf(stderr, "LEVMAR: empty options vector, using defaults\n");
344 }
345 #endif /* DEBUG */
346
347 /** mintype (optional) **/
348 /* check whether sixth argument is a string */
349 if(nrhs>=6 && mxIsChar(prhs[5])==1 && mxGetM(prhs[5])==1){
350 char *minhowto;
351
352 /* examine supplied name */
353 len=mxGetN(prhs[5])+1;
354 minhowto=mxCalloc(len, sizeof(char));
355 status=mxGetString(prhs[5], minhowto, len);
356 if(status!=0)
357 mexErrMsgTxt("levmar: not enough space. String is truncated.");
358
359 for(i=0; minhowto[i]; ++i)
360 minhowto[i]=tolower(minhowto[i]);
361 if(!strncmp(minhowto, "unc", 3)) mintype=MIN_UNCONSTRAINED;
362 else if(!strncmp(minhowto, "bc", 2)) mintype=MIN_CONSTRAINED_BC;
363 else if(!strncmp(minhowto, "lec", 3)) mintype=MIN_CONSTRAINED_LEC;
364 else if(!strncmp(minhowto, "blec", 4)) mintype=MIN_CONSTRAINED_BLEC;
365 else matlabFmtdErrMsgTxt("levmar: unknown minimization type '%s'.", minhowto);
366
367 mxFree(minhowto);
368
369 ++prhs;
370 --nrhs;
371 }
372 else
373 mintype=MIN_UNCONSTRAINED;
374
375 if(mintype==MIN_UNCONSTRAINED) goto extraargs;
376
377 /* arguments below this point are optional and their presence depends
378 * upon the minimization type determined above
379 */
380 /** lb, ub **/
381 if(nrhs>=7 && (mintype==MIN_CONSTRAINED_BC || mintype==MIN_CONSTRAINED_BLEC)){
382 /* check if the next two arguments are real row or column vectors */
383 if(mxIsDouble(prhs[5]) && !mxIsComplex(prhs[5]) && (mxGetM(prhs[5])==1 || mxGetN(prhs[5])==1)){
384 if(mxIsDouble(prhs[6]) && !mxIsComplex(prhs[6]) && (mxGetM(prhs[6])==1 || mxGetN(prhs[6])==1)){
385 if((i=__MAX__(mxGetM(prhs[5]), mxGetN(prhs[5])))!=m)
386 matlabFmtdErrMsgTxt("levmar: lb must have %d elements, got %d.", m, i);
387 if((i=__MAX__(mxGetM(prhs[6]), mxGetN(prhs[6])))!=m)
388 matlabFmtdErrMsgTxt("levmar: ub must have %d elements, got %d.", m, i);
389
390 lb=mxGetPr(prhs[5]);
391 ub=mxGetPr(prhs[6]);
392
393 prhs+=2;
394 nrhs-=2;
395 }
396 }
397 }
398
399 /** A, b **/
400 if(nrhs>=7 && (mintype==MIN_CONSTRAINED_LEC || mintype==MIN_CONSTRAINED_BLEC)){
401 /* check if the next two arguments are a real matrix and a real row or column vector */
402 if(mxIsDouble(prhs[5]) && !mxIsComplex(prhs[5]) && mxGetM(prhs[5])>=1 && mxGetN(prhs[5])>=1){
403 if(mxIsDouble(prhs[6]) && !mxIsComplex(prhs[6]) && (mxGetM(prhs[6])==1 || mxGetN(prhs[6])==1)){
404 if((i=mxGetN(prhs[5]))!=m)
405 matlabFmtdErrMsgTxt("levmar: A must have %d columns, got %d.", m, i);
406 if((i=__MAX__(mxGetM(prhs[6]), mxGetN(prhs[6])))!=(Arows=mxGetM(prhs[5])))
407 matlabFmtdErrMsgTxt("levmar: b must have %d elements, got %d.", Arows, i);
408
409 At=prhs[5];
410 b=mxGetPr(prhs[6]);
411 A=getTranspose(At);
412
413 prhs+=2;
414 nrhs-=2;
415 }
416 }
417 }
418
419 /* wghts */
420 /* check if we have a weights vector */
421 if(nrhs>=6 && mintype==MIN_CONSTRAINED_BLEC){ /* only check if we have seen both box & linear constraints */
422 if(mxIsDouble(prhs[5]) && !mxIsComplex(prhs[5]) && (mxGetM(prhs[5])==1 || mxGetN(prhs[5])==1)){
423 if(__MAX__(mxGetM(prhs[5]), mxGetN(prhs[5]))==m){
424 wghts=mxGetPr(prhs[5]);
425
426 ++prhs;
427 --nrhs;
428 }
429 }
430 }
431 /* arguments below this point are assumed to be extra arguments passed
432 * to every invocation of the fitting function and its Jacobian
433 */
434
435 extraargs:
436 /* handle any extra args and allocate memory for
437 * passing the current parameter estimate to matlab
438 */
439 nextra=nrhs-5;
440 mdata.nrhs=nextra+1;
441 mdata.rhs=(mxArray **)mxMalloc(mdata.nrhs*sizeof(mxArray *));
442 for(i=0; i<nextra; ++i)
443 mdata.rhs[i+1]=(mxArray *)prhs[nrhs-nextra+i]; /* discard 'const' modifier */
444 #ifdef DEBUG
445 fflush(stderr);
446 fprintf(stderr, "LEVMAR: %d extra args\n", nextra);
447 #endif /* DEBUG */
448
449 if(mdata.isrow_p0){ /* row vector */
450 mdata.rhs[0]=mxCreateDoubleMatrix(1, m, mxREAL);
451 /*
452 mxSetM(mdata.rhs[0], 1);
453 mxSetN(mdata.rhs[0], m);
454 */
455 }
456 else{ /* column vector */
457 mdata.rhs[0]=mxCreateDoubleMatrix(m, 1, mxREAL);
458 /*
459 mxSetM(mdata.rhs[0], m);
460 mxSetN(mdata.rhs[0], 1);
461 */
462 }
463
464 /* ensure that the supplied function & Jacobian are as expected */
465 if(checkFuncAndJacobian(p, m, n, havejac, &mdata)){
466 status=LM_ERROR;
467 goto cleanup;
468 }
469
470 if(nlhs>3) /* covariance output required */
471 covar=mxMalloc(m*m*sizeof(double));
472
473 /* invoke levmar */
474 if(!lb && !ub){
475 if(!A && !b){ /* no constraints */
476 if(havejac)
477 status=dlevmar_der(func, jacfunc, p, x, m, n, itmax, opts, info, NULL, covar, (void *)&mdata);
478 else
479 status=dlevmar_dif(func, p, x, m, n, itmax, opts, info, NULL, covar, (void *)&mdata);
480 #ifdef DEBUG
481 fflush(stderr);
482 fprintf(stderr, "LEVMAR: calling dlevmar_der()/dlevmar_dif()\n");
483 #endif /* DEBUG */
484 }
485 else{ /* linear constraints */
486 #ifdef HAVE_LAPACK
487 if(havejac)
488 status=dlevmar_lec_der(func, jacfunc, p, x, m, n, A, b, Arows, itmax, opts, info, NULL, covar, (void *)&mdata);
489 else
490 status=dlevmar_lec_dif(func, p, x, m, n, A, b, Arows, itmax, opts, info, NULL, covar, (void *)&mdata);
491 #else
492 mexErrMsgTxt("levmar: no linear constraints support, HAVE_LAPACK was not defined during MEX-file compilation.");
493 #endif /* HAVE_LAPACK */
494
495 #ifdef DEBUG
496 fflush(stderr);
497 fprintf(stderr, "LEVMAR: calling dlevmar_lec_der()/dlevmar_lec_dif()\n");
498 #endif /* DEBUG */
499 }
500 }
501 else{
502 if(!A && !b){ /* box constraints */
503 if(havejac)
504 status=dlevmar_bc_der(func, jacfunc, p, x, m, n, lb, ub, itmax, opts, info, NULL, covar, (void *)&mdata);
505 else
506 status=dlevmar_bc_dif(func, p, x, m, n, lb, ub, itmax, opts, info, NULL, covar, (void *)&mdata);
507 #ifdef DEBUG
508 fflush(stderr);
509 fprintf(stderr, "LEVMAR: calling dlevmar_bc_der()/dlevmar_bc_dif()\n");
510 #endif /* DEBUG */
511 }
512 else{ /* box & linear constraints */
513 #ifdef HAVE_LAPACK
514 if(havejac)
515 status=dlevmar_blec_der(func, jacfunc, p, x, m, n, lb, ub, A, b, Arows, wghts, itmax, opts, info, NULL, covar, (void *)&mdata);
516 else
517 status=dlevmar_blec_dif(func, p, x, m, n, lb, ub, A, b, Arows, wghts, itmax, opts, info, NULL, covar, (void *)&mdata);
518 #else
519 mexErrMsgTxt("levmar: no box & linear constraints support, HAVE_LAPACK was not defined during MEX-file compilation.");
520 #endif /* HAVE_LAPACK */
521
522 #ifdef DEBUG
523 fflush(stderr);
524 fprintf(stderr, "LEVMAR: calling dlevmar_blec_der()/dlevmar_blec_dif()\n");
525 #endif /* DEBUG */
526 }
527 }
528 #ifdef DEBUG
529 fflush(stderr);
530 printf("LEVMAR: minimization returned %d in %g iter, reason %g\n\tSolution: ", status, info[5], info[6]);
531 for(i=0; i<m; ++i)
532 printf("%.7g ", p[i]);
533 printf("\n\n\tMinimization info:\n\t");
534 for(i=0; i<LM_INFO_SZ; ++i)
535 printf("%g ", info[i]);
536 printf("\n");
537 #endif /* DEBUG */
538
539 /* copy back return results */
540 /** ret **/
541 plhs[0]=mxCreateDoubleMatrix(1, 1, mxREAL);
542 ret=mxGetPr(plhs[0]);
543 ret[0]=(double)status;
544
545 /** popt **/
546 plhs[1]=(mdata.isrow_p0==1)? mxCreateDoubleMatrix(1, m, mxREAL) : mxCreateDoubleMatrix(m, 1, mxREAL);
547 pdbl=mxGetPr(plhs[1]);
548 for(i=0; i<m; ++i)
549 pdbl[i]=p[i];
550
551 /** info **/
552 if(nlhs>2){
553 plhs[2]=mxCreateDoubleMatrix(1, LM_INFO_SZ, mxREAL);
554 pdbl=mxGetPr(plhs[2]);
555 for(i=0; i<LM_INFO_SZ; ++i)
556 pdbl[i]=info[i];
557 }
558
559 /** covar **/
560 if(nlhs>3){
561 plhs[3]=mxCreateDoubleMatrix(m, m, mxREAL);
562 pdbl=mxGetPr(plhs[3]);
563 for(i=0; i<m*m; ++i) /* covariance matrices are symmetric, thus no need to transpose! */
564 pdbl[i]=covar[i];
565 }
566
567 cleanup:
568 /* cleanup */
569 mxDestroyArray(mdata.rhs[0]);
570 if(A) mxFree(A);
571
572 mxFree(mdata.fname);
573 if(havejac) mxFree(mdata.jacname);
574 mxFree(p);
575 mxFree(mdata.rhs);
576 if(covar) mxFree(covar);
577
578 if(status==LM_ERROR)
579 mexWarnMsgTxt("levmar: optimization returned with an error!");
580 }
581