469,631 Members | 1,302 Online
Bytes | Developer Community
New Post

Home Posts Topics Members FAQ

Post your question to a community of 469,631 developers. It's quick & easy.

allocating memory for multi dimensional arrays

Hi,

After my question today, I am not sure if this is the right group to ask but
I have a simple question I guess.

I have multi dimensional array called cover what I do is this

iT GIVES AN ERROR that the memory cannot be read.

So I believe I do not how how to handle multi dimensional arrays.

Can anyone help me please?

#define SWAP(a,b) {swap=(a);(a)=(b);(b)=swap;}

void main()
{
Covar = (float**)malloc(sizeof(float)*ma);
covsrt(covar, ma, ia, mfit);

}

void covsrt(float **covar, int ma, int ia[], int mfit)
{
int i,j,k;
float swap;

for (i=mfit+1;i<=ma;i++)
for (j=1;j<=i;j++) covar[i][j]=covar[j][i]=0.0;
k=mfit;
for (j=ma;j>=1;j--) {
if (ia[j]) {
for (i=1;i<=ma;i++) SWAP(covar[i][k],covar[i][j])
for (i=1;i<=ma;i++) SWAP(covar[k][i],covar[j][i])
k--;
}
}
}
Jul 22 '05 #1
4 2139

"kak3012" <s0*****@student.dtu.dk> wrote in message
news:ct**********@gnd.k-net.dk...
Hi,

After my question today, I am not sure if this is the right group to ask but I have a simple question I guess.
You have taken a C approach to this, and not a C++ (you use malloc instead
of new []). But anyway, malloc wants to know how many bytes you want to
allocate. Normally you would do MatrixWidth * MatrixHeight * sizeof (
DataType ). I am not sure if your variable ma covers MatrixWidth *
MatrixHeight .

Instead of trying to figure that out I would do something like...

/// Simple Matrix

template <typename T,unsigned int W,unsigned int H>
struct Matrix {
struct MatrixColumn {
T data_[H];

unsigned int size() { return H; };

T& operator [] (unsigned int x) {
return data_[x];
}
};
MatrixColumn data_[W];

unsigned int size() { return W; };

MatrixColumn& operator [] (unsigned int x) {
return data_[x];
}
void initializeTo( T initialValue ) {
for (unsigned long x = 0; x<W; ++x)
for (unsigned long y = 0; y<H; ++y)
(*this)[x][y] = initialValue;
}

};

And declare my matrix as:

{
Matrix< float,MatrixWidth,MatrixHeight > myMatrix;
myMatrix.initializeTo( 0.0 );
}
The rest of your code I have a hard time deciphering, since the purpose of
the code is not stated anywhere.
But to me it looks like you are allocating only one row or column of data
with the malloc...
I have multi dimensional array called cover what I do is this

iT GIVES AN ERROR that the memory cannot be read.

So I believe I do not how how to handle multi dimensional arrays.

Can anyone help me please?

#define SWAP(a,b) {swap=(a);(a)=(b);(b)=swap;}

void main()
{
Covar = (float**)malloc(sizeof(float)*ma);
covsrt(covar, ma, ia, mfit);

}

void covsrt(float **covar, int ma, int ia[], int mfit)
{
int i,j,k;
float swap;

for (i=mfit+1;i<=ma;i++)
for (j=1;j<=i;j++) covar[i][j]=covar[j][i]=0.0;
k=mfit;
for (j=ma;j>=1;j--) {
if (ia[j]) {
for (i=1;i<=ma;i++) SWAP(covar[i][k],covar[i][j])
for (i=1;i<=ma;i++) SWAP(covar[k][i],covar[j][i])
k--;
}
}
}

Jul 22 '05 #2
Thanks very much for your reply Jesper.

I am trying to use Numerical Reciies to fit my data into a sinusoidal curve
and I have read the NR that I can do it with some function in it.

I do not want to change anything in NR codes to let the next user to be able
to use it as a standart.

Your solution includes a change in the 1 of the NR codes.

I can send the whole code if you are intersted...

Jul 22 '05 #3

"kak3012" <s0*****@student.dtu.dk> wrote in message
news:ct**********@gnd.k-net.dk...
Thanks very much for your reply Jesper.

I am trying to use Numerical Reciies to fit my data into a sinusoidal curve and I have read the NR that I can do it with some function in it.

I do not want to change anything in NR codes to let the next user to be able to use it as a standart.

Your solution includes a change in the 1 of the NR codes.

I can send the whole code if you are intersted...


Please do - I might get it up and running if I have something compilable..
Jul 22 '05 #4
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "nrutil.h"
#define NR_END 1
#define FREE_ARG char*
#define SWAP(a,b) {swap=(a);(a)=(b);(b)=swap;}

char Buffer[1024];

typedef void (*MYFUNCTYPE)(float, float [], float *, float [], int);
void funcs(float x, float a[], float *y, float dyda[], int na);
void covsrt(float **covar, int ma, int ia[], int mfit);
void free_vector(float *v, long nl, long nh);
float *vector(long nl, long nh);
int *ivector(long nl, long nh);
void free_ivector(int *v, long nl, long nh);
void free_matrix(float **m, long nrl, long nrh, long ncl, long nch);
float **matrix(long nrl, long nrh, long ncl, long nch);
void nrerror(char error_text[]);
void gaussj(float **a, int n, float **b, int m);
void mrqmin(float x[], float y[], float sig[], int ndata, float a[], int
ia[],int ma, float **covar, float **alpha, float *chisq,MYFUNCTYPE funcs,
float *alamda);


int main()
{

int i;
int ma=9,mfit=0; // number of parameters
float *a; // parameters array
int *ia;
float **alpha; // curvature matrix
float *beta; // beta vector
float *chisq,ochisq;
float *lam;
float *x,*y,*sig,*yfit;
float **covar;
int N=100;

// Start the process

a = (float*)malloc(sizeof(float)*ma);
ia = (int*)malloc(sizeof(int)*ma);
beta = (float*)malloc(sizeof(float)*ma);
x = (float*)malloc(sizeof(float)*N);
y = (float*)malloc(sizeof(float)*N);
yfit = (float*)malloc(sizeof(float)*N);
sig = (float*)malloc(sizeof(float)*N);
alpha = (float**)malloc(sizeof(float)*ma);
lam = (float*)malloc(sizeof(float)*ma);
//covar = (float**)malloc(sizeof(float)*ma*ma);
covar = new float*[ma];

//covar[0][0]=(float)1.2;
//sprintf(Buffer, " covar[1][1]=%f\n ", covar[0][0]);
// printf(Buffer);

for (i=0;i<N;i++)
{
x[i]=(float)rand();
y[i]=x[i]*sin(x[i]);
sig[i]=1;
sprintf(Buffer, " x_in=%f y_in=%f\n", x[i],y[i]);
printf(Buffer);
}

for (i=0;i<ma;i++)
{
a[i]=0;
ia[i]=1;
beta[i]=0;
}

*lam=(float)0.01;
chisq=0;

covsrt(covar,ma,ia,mfit);

mrqmin(x,y,sig,N,a,ia,ma,covar,alpha,chisq,funcs,l am);

for (i=0;i<N;i++)
{
sprintf(Buffer, " x=%f y=%f ", x[i],y[i]);
printf(Buffer);
}

ochisq=0;

}

void covsrt(float **covar, int ma, int ia[], int mfit)
{
int i,j,k;
float swap;

covar = new float*[ma];

for (i=mfit+1;i<=ma;i++)
covar[i] = new float[ma];
covar[i][1]=(float)0;
exit(0);
for (j=1;j<=i;j++) {covar[i][j]=covar[j][i]=0.0;}

k=mfit;
for (j=ma;j>=1;j--) {
if (ia[j]) {
for (i=1;i<=ma;i++) SWAP(covar[i][k],covar[i][j])
for (i=1;i<=ma;i++) SWAP(covar[k][i],covar[j][i])
k--;
}
}
}

void funcs(float x, float a[], float *y, float dyda[], int na)
{
float fac,Mcos,Msin;

*y=0.0;
for(int i=1;i<na-1;i+=3)
{
Mcos=cos(x+a[i+2]);
Msin=sin(x+a[i+2]);
fac=-a[i+1]*sin(x+a[i+2]);
*y+=a[i]+(a[i+1]*Mcos);
dyda[i]=1;
dyda[i+1]=Mcos-(a[i+1]*Msin);
dyda[i+2]=-a[i+1]*Msin;
}
}
typedef void (*MYFUNCTYPE)(float, float [], float *, float [], int);
void mrqmin(float x[], float y[], float sig[], int ndata, float a[], int
ia[],
int ma, float **covar, float **alpha, float *chisq,MYFUNCTYPE funcs, float
*alamda)
{
//void covsrt(float **covar, int ma, int ia[], int mfit);
//void gaussj(float **a, int n, float **b, int m);
void mrqcof(float x[], float y[], float sig[], int ndata, float a[],int
ia[], int ma, float **alpha, float beta[], float *chisq,MYFUNCTYPE funcs);

int j,k,l;
static int mfit;
static float ochisq,*atry,*beta,*da,**oneda;
if (*alamda < 0.0) {
atry=vector(1,ma);
beta=vector(1,ma);
da=vector(1,ma);
for (mfit=0,j=1;j<=ma;j++)
if (ia[j]) mfit++;
oneda=matrix(1,mfit,1,1);
*alamda=(float)0.001;
mrqcof(x,y,sig,ndata,a,ia,ma,alpha,beta,chisq,func s);
ochisq=(*chisq);
for (j=1;j<=ma;j++) atry[j]=a[j];
}
for (j=1;j<=mfit;j++) {
for (k=1;k<=mfit;k++) covar[j][k]=alpha[j][k];
covar[j][j]=(float)(alpha[j][j]*(1.0+(*alamda)));
oneda[j][1]=beta[j];
}
gaussj(covar,mfit,oneda,1);
for (j=1;j<=mfit;j++) da[j]=oneda[j][1];
if (*alamda == 0.0) {
covsrt(covar,ma,ia,mfit);
covsrt(alpha,ma,ia,mfit);
free_matrix(oneda,1,mfit,1,1);
free_vector(da,1,ma);
free_vector(beta,1,ma);
free_vector(atry,1,ma);
return;
}
for (j=0,l=1;l<=ma;l++)
if (ia[l]) atry[l]=a[l]+da[++j];
mrqcof(x,y,sig,ndata,atry,ia,ma,covar,da,chisq,fun cs);
if (*chisq < ochisq) {
*alamda *= (float)0.1;
ochisq=(*chisq);
for (j=1;j<=mfit;j++) {
for (k=1;k<=mfit;k++) alpha[j][k]=covar[j][k];
beta[j]=da[j];
}
for (l=1;l<=ma;l++) a[l]=atry[l];
} else {
*alamda *= 10.0;
*chisq=ochisq;
}
}

float *vector(long nl, long nh)
/* allocate a float vector with subscript range v[nl..nh] */
{
float *v;

v=(float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(float)));
if (!v) nrerror("allocation failure in vector()");
return v-nl+NR_END;
}

void nrerror(char error_text[])
/* Numerical Recipes standard error handler */
{
fprintf(stderr,"Numerical Recipes run-time error...\n");
fprintf(stderr,"%s\n",error_text);
fprintf(stderr,"...now exiting to system...\n");
exit(1);
}

float **matrix(long nrl, long nrh, long ncl, long nch)
/* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
float **m;

/* allocate pointers to rows */
m=(float **) malloc((size_t)((nrow+NR_END)*sizeof(float*)));
if (!m) nrerror("allocation failure 1 in matrix()");
m += NR_END;
m -= nrl;

/* allocate rows and set pointers to them */
m[nrl]=(float *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float))) ;
if (!m[nrl]) nrerror("allocation failure 2 in matrix()");
m[nrl] += NR_END;
m[nrl] -= ncl;

for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;

/* return pointer to array of pointers to rows */
return m;
}
void free_matrix(float **m, long nrl, long nrh, long ncl, long nch)
/* free a float matrix allocated by matrix() */
{
free((FREE_ARG) (m[nrl]+ncl-NR_END));
free((FREE_ARG) (m+nrl-NR_END));
}
void free_vector(float *v, long nl, long nh)
/* free a float vector allocated with vector() */
{
free((FREE_ARG) (v+nl-NR_END));
}

void free_ivector(int *v, long nl, long nh)
/* free an int vector allocated with ivector() */
{
free((FREE_ARG) (v+nl-NR_END));
}

int *ivector(long nl, long nh)
/* allocate an int vector with subscript range v[nl..nh] */
{
int *v;

v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
if (!v) nrerror("allocation failure in ivector()");
return v-nl+NR_END;
}

void mrqcof(float x[], float y[], float sig[], int ndata, float a[], int
ia[],int ma, float **alpha, float beta[], float *chisq,MYFUNCTYPE funcs)

{
int i,j,k,l,m,mfit=0;
float ymod,wt,sig2i,dy,*dyda;
dyda=vector(1,ma);
for (j=1;j<=ma;j++)
if (ia[j]) mfit++;
for (j=1;j<=mfit;j++) {
for (k=1;k<=j;k++) alpha[j][k]=0.0;
beta[j]=0.0;
}
*chisq=0.0;
for (i=1;i<=ndata;i++) {
(*funcs)(x[i],a,&ymod,dyda,ma);
sig2i=(float)(1.0/(sig[i]*sig[i]));
dy=y[i]-ymod;
for (j=0,l=1;l<=ma;l++) {
if (ia[l]) {
wt=dyda[l]*sig2i;
for (j++,k=0,m=1;m<=l;m++)
if (ia[m]) alpha[j][++k] += wt*dyda[m];
beta[j] += dy*wt;
}
}
*chisq += dy*dy*sig2i;
}
for (j=2;j<=mfit;j++)
for (k=1;k<j;k++) alpha[k][j]=alpha[j][k];
free_vector(dyda,1,ma);
}
void gaussj(float **a, int n, float **b, int m)
{

int *indxc,*indxr,*ipiv;
int i,icol,irow,j,k,l,ll;
float big,dum,pivinv;//,temp;
float swap;

indxc=ivector(1,n);
indxr=ivector(1,n);
ipiv=ivector(1,n);
for (j=1;j<=n;j++) ipiv[j]=0;
for (i=1;i<=n;i++) {
big=0.0;
for (j=1;j<=n;j++)
if (ipiv[j] != 1)
for (k=1;k<=n;k++) {
if (ipiv[k] == 0) {
if (fabs(a[j][k]) >= big) {
big=(float)(fabs(a[j][k]));
irow=j;
icol=k;
}
} else if (ipiv[k] > 1) nrerror("gaussj: Singular Matrix-1");
}
++(ipiv[icol]);
if (irow != icol) {
for (l=1;l<=n;l++) SWAP(a[irow][l],a[icol][l])
for (l=1;l<=m;l++) SWAP(b[irow][l],b[icol][l])
}
indxr[i]=irow;
indxc[i]=icol;
if (a[icol][icol] == 0.0) nrerror("gaussj: Singular Matrix-2");
pivinv=(float)(1.0/a[icol][icol]);
a[icol][icol]=1.0;
for (l=1;l<=n;l++) a[icol][l] *= pivinv;
for (l=1;l<=m;l++) b[icol][l] *= pivinv;
for (ll=1;ll<=n;ll++)
if (ll != icol) {
dum=a[ll][icol];
a[ll][icol]=0.0;
for (l=1;l<=n;l++) a[ll][l] -= a[icol][l]*dum;
for (l=1;l<=m;l++) b[ll][l] -= b[icol][l]*dum;
}
}
for (l=n;l>=1;l--) {
if (indxr[l] != indxc[l])
for (k=1;k<=n;k++)
SWAP(a[k][indxr[l]],a[k][indxc[l]]);
}
free_ivector(ipiv,1,n);
free_ivector(indxr,1,n);
free_ivector(indxc,1,n);

}
#undef SWAP

Jul 22 '05 #5

This discussion thread is closed

Replies have been disabled for this discussion.

Similar topics

11 posts views Thread by fivelitermustang | last post: by
5 posts views Thread by Cant Think Today | last post: by
4 posts views Thread by Richard Hayden | last post: by
4 posts views Thread by Balaskas Evaggelos | last post: by
4 posts views Thread by =?Utf-8?B?SGVucmlrIFNjaG1pZA==?= | last post: by
reply views Thread by gheharukoh7 | last post: by
By using this site, you agree to our Privacy Policy and Terms of Use.