Problem passing array

[mrhegemon]

I am trying to pass the array "da" from the non_linear_min function to the lev_mrq_min_step function, have it filled with values, and then return it to add it into another array. The problem is, the "da" array appears to be given proper values from it's assignment in the lev_mrq_min_step function, as I can verify from the print_matrix function I have made, but when I try to access it after I am returned to the non_linear_min function, it appears empty. I am confused because it is a simple pointer.

Here is my code (excuse all the functions I defined but didn't include, just assume they work correctly):

double lev_mrq_min_step(double *datax,long p,double *datay,long n,double sig2,double *a,long m,long q,double *da,double (*func)(double*,double*,long,double*),double &lambda)
{
  long i=0,k=0,l=0,ratio=p/n,j=0;
  //beta(grad) is -0.5*dchisqdak, alpha(hess) is 0.5*d2chisq/dakdal
  double
  chisq=0,
  tmpy=0,
  arg=0,
  *xbuf=(double*)malloc(ratio*sizeof(double)),
  *dyda=(double*)malloc(3*sizeof(double)),
  *beta=(double*)calloc(q,sizeof(double)),
  *alpha=(double*)calloc(q*q,sizeof(double)),
  *tmpa=(double*)malloc(m*sizeof(double)),
  *tmpalpha=(double*)malloc(q*q*sizeof(double));

  i=0;while(i<m){tmpa[i]=a[i];++i;}
  if(lambda!=0){
    //matrix computation
    i=0;while(i<n){
      for(j=0;j<ratio;++j) *(xbuf+j)=*(datax+ratio*i+j);
      tmpy=func(xbuf,tmpa,m,dyda);
      arg=datay[i]-tmpy;
      k=0;while(k<q){
        //fill beta(gradient)
        beta[k]+=arg*dyda[k]/sig2;
        //fill alpha(hessian,linearized)
        l=0;while(l<q){*(alpha+k*q+l)+=dyda[k]*dyda[l]/sig2;++l;}
        ++k;}
      ++i;}

    printf(&quot;alpha:\n&quot;);
    print_matrix(alpha,q,q);
    printf(&quot;beta:\n&quot;);
    print_matrix(beta,1,q);

    //copy alpha and modify diagonal
    arg=1+lambda;i=0;while(i<q){k=0;while(k<q){*(tmpalpha+i*q+k)=(k==i)?(*(alpha+i*q+k))*arg:*(alpha+i*q+k);++k;}++i;}

    //solve system alpha*da=beta for da
    double *inv_alpha=matrix_invert(tmpalpha,q);
    da=matrix_mult(inv_alpha,q,q,beta,1);
    free(inv_alpha);

    printf(&quot;param step:\n&quot;);
    print_matrix(da,1,q);

    //create tmpa=a+da array for chisq
    i=0;while(i<q){tmpa[i]+=da[i];++i;}
  }
  //compute new chisq(a+da)
  chisq=0;i=0;while(i<n){for(j=0;j<ratio;++j) *(xbuf+j)=*(datax+ratio*i+j);tmpy=func(xbuf,tmpa,m,dyda);arg=datay[i]-tmpy;chisq+=arg*arg/sig2;++i;}

  free(xbuf);free(dyda);free(alpha);free(beta);free(tmpa);free(tmpalpha);
  return chisq;
}

double non_linear_min(double *datax,long p,double *datay,long n,double sig2,double *a,long m,long q,double (*func)(double*,double*,long,double*),double tolerance,long &iter)
{
  bool fail=false;
  long i=0,j=0;
  double lambda=0,oldchisq=0,chisq=0,*da;

  oldchisq=lev_mrq_min_step(datax,p,datay,n,sig2,a,m,q,da,func,lambda);

  free(da);

  lambda=0.001;

  while(i<iter){
    printf(&quot;chisq: % .8lf\nlamba: % .8lf\n----------\n&quot;,oldchisq,lambda);

    double *da;
    chisq=lev_mrq_min_step(datax,p,datay,n,sig2,a,m,q,da,func,lambda);
    ++i;

    print_matrix(da,1,q);

    //compare old and new chisq and update parameters if necessary
    if(chisq>=oldchisq){lambda*=10;fail=true;}
    else{j=0;while(j<q){a[j]+=da[j];++j;}lambda*=0.1;fail=false;}

    free(da);

    if(fail)printf(&quot;***fail***\n&quot;);

    if(((fabs(oldchisq-chisq)<tolerance)&&(!fail))||(chisq==0)) break;
    /*if(chisq<oldchisq)*/ oldchisq=chisq;}

  iter=i;
  return chisq;
}

Thanks for any help.

 

[Salem]

Oh dear, more C++ code posted on the C forum.

You're basically wanting to do this

void foo ( double **a ) {
  *a = malloc( sizeof **a ); 
}

int main ( void ) {
  double *ptr;
  foo( &ptr );
  free( ptr );
}

You need to pass a pointer to your pointer if you want changes made in your lev_mrq_min_step() function to be made available to your non_linear_min() function.

Or since this is C++ you're using, you can simply make a reference out of it (like you do with the iter parameter)

void foo ( double * &a ) {
  a = malloc( sizeof *a ); 
}

int main ( void ) {
  double *ptr;
  foo( ptr );
  free( ptr );
}

--