matrix stuff (solving b = A*x) --> using numerical recipes

[=?ISO-8859-1?Q?Martin_J=F8rgensen?=]

Hi,

I've made a program from numerical recipes. Looks like I'm not allowed
to distribute the source code from numerical recipes but it shouldn't
even be necessary to do that.

My problem is that I'm not very experienced with pointers, pointers to
pointers and the like and I got 4 compiler warnings + I don't completely
understand how to build this "compact matrix" (see later).

I'm also not completely sure if I understand the text correct... So far
my program looks like this:



#include <math.h>
#include "nrutil.h"

/* definitions */
#define number_of_rows 7
#define TINY 1.0e-20
#define SWAP(a,b) {dum=(a);(a)=(b);(b)=dum;}

/* prototypes */
void banmul(float **a, unsigned long n, int left, int right, float x[],
float b[]);
void bandec(float **a, unsigned long n, int left, int right, float **al,
unsigned long indx[], float *d);
void banbks(float **a, unsigned long n, int left, int right, float **al,
unsigned long indx[], float b[]);

/* Start of main program */

int main(void)
{
/* declaration */
double x[number_of_rows], b[number_of_rows]; // remember that we're
solving "b = A * x"

/* initialisation */
double a[][number_of_rows] = // number_of_cols = number_of_rows
{
{3.,1.,0.,0.,0.,0.,0.},
{4.,1.,5.,0.,0.,0.,0.},
{9.,2.,6.,5.,0.,0.,0.},
{0.,3.,5.,8.,9.,0.,0.},
{0.,0.,7.,9.,3.,2.,0.},
{0.,0.,0.,3.,8.,4.,6.},
{0.,0.,0.,0.,2.,4.,4.}
};

banmul(a, number_of_rows, 2, 1, x, b);

printf("\nFinishing program now.\n\n");
}


If there's anything you need to know you can read it here:

http://www.library.cornell.edu/nr/bookcpdf/c2-4.pdf (start from the
chapter with "Band diagonal system, p.52".

You can see my 2D-matrix is the same... The text talks about transposing
the matrix to a more compact form... I must admit I need some help
because I didn't understand the text :-(

My compiler warnings (don't know what to do):

----
warning C4047: 'function' : 'float **' differs in levels of indirection
from 'double [7][7]'

warning C4024: 'banmul' : different types for formal and actual parameter 1

warning C4133: 'function' : incompatible types - from 'double [7]' to
'float *'

warning C4133: 'function' : incompatible types - from 'double [7]' to
'float *'
----

The problem is the line: "banmul(a, ....)";


So, generally, I don't really understand what I'm doing :-(

The x-array (vector) is not defined. If somebody could show me an
example of how to solve this matrix system, I would be very happy (you
just make x whatever you like: x=[2,5,2,1,6,2,1] or whatever...

AFAIK the text from numerical recipes and my code should be enough for
you to understand my code and I hope somebody can give me some hints
that can get me closer to a solution.

Thanks in advance for any (hopefully good) inputs...


Best regards / Med venlig hilsen
Martin Jørgensen

 

[Richard Heathfield]

Martin Jørgensen said:

My problem is that I'm not very experienced with pointers, pointers to
pointers and the like and I got 4 compiler warnings + I don't completely
understand how to build this "compact matrix" (see later).

doubles are not floats, and arrays are not pointers.

void banmul(float **a, unsigned long n, int left, int right, float x[],
float b[]);

/* declaration */
double x[number_of_rows], b[number_of_rows]; // remember that we're
solving "b = A * x"

/* initialisation */
double a[][number_of_rows] = // number_of_cols = number_of_rows
{ [...]
};

banmul(a, number_of_rows, 2, 1, x, b);

banmul takes float ** for parameter 1, but you're actually passing it a
double (*)[7], which is not the same thing by several rows of apple trees.

banmul takes float * for parameters 5 and 6, and you are trying to pass it
double *.

 

[Keith Thompson]

warning C4047: 'function' : 'float **' differs in levels of
indirection from 'double [7][7]'

Arrays are not pointers, and 2-dimensional arrays are not pointers to
pointers. The comp.lang.c FAQ is at <http://www.c-faq.com/>; start
with section 6, "Arrays and Pointers".

 

[=?ISO-8859-1?Q?Martin_J=F8rgensen?=]

Martin Jørgensen said: -snip-

banmul(a, number_of_rows, 2, 1, x, b);

banmul takes float ** for parameter 1, but you're actually passing it a
double (*)[7], which is not the same thing by several rows of apple trees.

Damn... I understand what you're writing. I just don't know how to fix
the problem... Could you/somebody please help?

I tried:

int main(void)
{
int i;

/* declaration */
float x[number_of_rows], b[number_of_rows]; // remember that we're
solving "b = A * x"
float **a = malloc((number_of_rows)*sizeof(float*));
a[0] = malloc((number_of_rows*number_of_rows)*sizeof(float));

for (i = 1; i < number_of_rows ;i++)
a = a[i-1] + number_of_rows+1;

a = //
{
{3.,1.,0.,0.,0.,0.,0.},
{4.,1.,5.,0.,0.,0.,0.},
{9.,2.,6.,5.,0.,0.,0.},
{0.,3.,5.,8.,9.,0.,0.},
{0.,0.,7.,9.,3.,2.,0.},
{0.,0.,0.,3.,8.,4.,6.},
{0.,0.,0.,0.,2.,4.,4.}
};

banmul(a, number_of_rows, 2, 1, x, b);

printf("\nFinishing program now.\n\n");
}

This gives 3 problems:

1) warning C4047: 'initializing' : 'float **' differs in levels of
indirection from 'int' - this is the line: float **a =
malloc((number_of_rows)*sizeof(float*));

2) warning C4047: '=' : 'float *' differs in levels of indirection from
'int' - this is the line: a[0] =
malloc((number_of_rows*number_of_rows)*sizeof(float));

3) syntax error : '{' - I now have a problem with initializing the 2D
"a"-matrix with this notation...

banmul takes float * for parameters 5 and 6, and you are trying to pass it
double *.

Damn... Stupid mistake. I fixed that just by changing the definition.


Best regards / Med venlig hilsen
Martin Jørgensen

 

[=?ISO-8859-1?Q?Martin_J=F8rgensen?=]

warning C4047: 'function' : 'float **' differs in levels of
indirection from 'double [7][7]'

Arrays are not pointers, and 2-dimensional arrays are not pointers to
pointers. The comp.lang.c FAQ is at <http://www.c-faq.com/>; start
with section 6, "Arrays and Pointers".

I read it and felt I understood it... But I need more help yet... See my
other post for a more thorough description of what I now tried...


Best regards / Med venlig hilsen
Martin Jørgensen

 

[stathis gotsis]

I tried:

int main(void)
{
int i;

Where do you declare number_of_rows?

/* declaration */
float x[number_of_rows], b[number_of_rows]; // remember that we're
solving "b = A * x"
float **a = malloc((number_of_rows)*sizeof(float*));
a[0] = malloc((number_of_rows*number_of_rows)*sizeof(float));

for (i = 1; i < number_of_rows ;i++)
a = a[i-1] + number_of_rows+1;

Why "+1"?

a = //
{
{3.,1.,0.,0.,0.,0.,0.},
{4.,1.,5.,0.,0.,0.,0.},
{9.,2.,6.,5.,0.,0.,0.},
{0.,3.,5.,8.,9.,0.,0.},
{0.,0.,7.,9.,3.,2.,0.},
{0.,0.,0.,3.,8.,4.,6.},
{0.,0.,0.,0.,2.,4.,4.}
};

You would better initialize some other 2D-array this way and copy its
contents into a.

banmul(a, number_of_rows, 2, 1, x, b);

printf("\nFinishing program now.\n\n");

Returning 0 would not be a bad idea.

 

[Richard Heathfield]

Martin Jørgensen said:

Martin Jørgensen said: -snip-

banmul(a, number_of_rows, 2, 1, x, b);

banmul takes float ** for parameter 1, but you're actually passing it a
double (*)[7], which is not the same thing by several rows of apple
trees.

Damn... I understand what you're writing. I just don't know how to fix
the problem...

Presuming you want a square:

float **p = malloc(number_of_rows * sizeof *p);
if(p != NULL)
{
for(i = 0; i < number_of_rows; i++)
{
p = malloc(number_of_rows * sizeof *p);
if(p != NULL)
{
for(j = 0; j < number_of_rows; j++)
{
/* might as well use your other array to initialise with */
p[j] = a[j];
}
}
else
{
you're running low on memory, and you might want to do something
about that
}
}
}
else
{
you're running low on memory, and you might want to do something
about that
}

Now, this method is a bit wasteful in the sense that it uses your original
array *as well* - for initialising. But I can't help wondering whether,
right now, you'd rather have something than nothing.

This gives 3 problems:

1) warning C4047: 'initializing' : 'float **' differs in levels of
indirection from 'int' - this is the line: float **a =
malloc((number_of_rows)*sizeof(float*));

That's because you forgot to prototype malloc, so the compiler is forced to
assume it returns int, whereas in fact it returns void *. Simply add:

#include <stdlib.h>

to fix that problem.

 

[=?ISO-8859-1?Q?Martin_J=F8rgensen?=]

Martin Jørgensen said: -snip-

Presuming you want a square: -snip-

Now, this method is a bit wasteful in the sense that it uses your original
array *as well* - for initialising. But I can't help wondering whether,
right now, you'd rather have something than nothing.

Completely correct! You read my mind, oh great mindreader

That's because you forgot to prototype malloc, so the compiler is forced to
assume it returns int, whereas in fact it returns void *. Simply add:

#include <stdlib.h>

to fix that problem.

Dough! Fixed that now.

Since I'm a newbie, I still need some comments on my code. For avoiding
confusion, I post all of the code I got except the functions from
Numerical recipes (can be seen on
http://www.library.cornell.edu/nr/bookcpdf/c2-4.pdf p.52-54 AFAIR).

-------------
/* Including header files */
#include <stdlib.h>
#include <math.h>
#include "nrutil.h"

/* definitions */
#define number_of_rows 7
#define TINY 1.0e-20
#define SWAP(a,b) {dum=(a);(a)=(b);(b)=dum;}

/* prototypes */
void banmul(float **a, unsigned long n, int left, int right, float x[],
float b[]);
void bandec(float **a, unsigned long n, int left, int right, float **al,
unsigned long indx[], float *d);
void banbks(float **a, unsigned long n, int left, int right, float **al,
unsigned long indx[], float b[]);

/* Start of main program */

int main(void)
{
/* declaration */
float **a;
float testvector[number_of_rows][number_of_rows];
int i, j;
float x[number_of_rows], b[number_of_rows]; // remember that we're
solving "b = A * x"

testvector[][number_of_rows] = // number_of_cols = number_of_rows
{
{3.,1.,0.,0.,0.,0.,0.},
{4.,1.,5.,0.,0.,0.,0.},
{9.,2.,6.,5.,0.,0.,0.},
{0.,3.,5.,8.,9.,0.,0.},
{0.,0.,7.,9.,3.,2.,0.},
{0.,0.,0.,3.,8.,4.,6.},
{0.,0.,0.,0.,2.,4.,4.}
};


float **a = malloc(number_of_rows * sizeof *a);
if(a != NULL)
{
for(i = 0; i < number_of_rows; i++)
{
a = malloc(number_of_rows * sizeof *a);
if(a != NULL)
{
for(j = 0; j < number_of_rows; j++)
{
/* might as well use your other array to initialise with */
a[j] = testvector[j];
}
}
else
{
printf("You're running low on memory - you should do something.\n");
exit(1);
}
}
}
else
{
printf("You're running low on memory - you should do something.\n");
exit(1);
}

banmul(a, number_of_rows, 2, 1, x, b);

printf("\nFinishing program now.\n\n");
return(0);
}

/* Here comes the NR functions */
.....
-------------

Compiler errors:
1) error C2059: syntax error : ']'

This is really strange! I don't understand that - it used to work before
I think.... Error 1) happen in the line: testvector[][number_of_rows] =
..... { {3.,1.,0.,0.,......},{4.,1.,..} etc};


2) error C2143: syntax error : missing ';' before 'type'
Error 2) Happens just after the testvector (matrix actually) is defined
- in the line: float **a = malloc(number_of_rows * sizeof *a);

Therefore error 2 is probably connected to error 1.

Thanks in advance for a solution that will fix my problems...


Best regards / Med venlig hilsen
Martin Jørgensen

 

[=?ISO-8859-1?Q?Martin_J=F8rgensen?=]

Hopefully you are including <stdlib.h>

Damn... Thanks...

Where do you declare number_of_rows?

Yeah, sorry. It's in the top and I forgot to include that part in my
posting. You can see my new code which is still causing me some trouble
in my other reply to Richard.

/* declaration */
float x[number_of_rows], b[number_of_rows]; // remember that we're
solving "b = A * x"
float **a = malloc((number_of_rows)*sizeof(float*));
a[0] = malloc((number_of_rows*number_of_rows)*sizeof(float));

for (i = 1; i < number_of_rows ;i++)
a = a[i-1] + number_of_rows+1;

Why "+1"?

Good question... I think it's because I took it from some other code
where I had a mesh consisting of nx=6 elements (or whatever) but they
didn't start at index 0 - but at index 1... Using +1 was probably a
mistake in this new code, so I changed it...

You would better initialize some other 2D-array this way and copy its
contents into a.

Understood... But I now get a stupid error, which I can't figure out
of... Probably something *really* simple - see me other reply...

Returning 0 would not be a bad idea.

You're right

Code changed a bit. Thanks for the comments.


Best regards / Med venlig hilsen
Martin Jørgensen

 

[stathis gotsis]

Since I'm a newbie, I still need some comments on my code. For avoiding
confusion, I post all of the code I got except the functions from
Numerical recipes (can be seen on
http://www.library.cornell.edu/nr/bookcpdf/c2-4.pdf p.52-54 AFAIR).

#include <stdlib.h>
#include <math.h>
#include "nrutil.h"

/* definitions */
#define number_of_rows 7
#define TINY 1.0e-20
#define SWAP(a,b) {dum=(a);(a)=(b);(b)=dum;}

/* prototypes */
void banmul(float **a, unsigned long n, int left, int right, float x[],
float b[]);
void bandec(float **a, unsigned long n, int left, int right, float **al,
unsigned long indx[], float *d);
void banbks(float **a, unsigned long n, int left, int right, float **al,
unsigned long indx[], float b[]);

/* Start of main program */

int main(void)
{
/* declaration */
float **a;
float testvector[number_of_rows][number_of_rows];

You intend to initialize this array, you should consider doing that upon
declaration.

int i, j;
float x[number_of_rows], b[number_of_rows]; // remember that we're
solving "b = A * x"
testvector[][number_of_rows] = // number_of_cols = number_of_rows
{
{3.,1.,0.,0.,0.,0.,0.},
{4.,1.,5.,0.,0.,0.,0.},
{9.,2.,6.,5.,0.,0.,0.},
{0.,3.,5.,8.,9.,0.,0.},
{0.,0.,7.,9.,3.,2.,0.},
{0.,0.,0.,3.,8.,4.,6.},
{0.,0.,0.,0.,2.,4.,4.}
};

This is not a proper initialization, that is why the compiler complained.

float **a = malloc(number_of_rows * sizeof *a);

You are re-declaring a. Try: a = malloc(number_of_rows * sizeof *a);

if(a != NULL)
{
for(i = 0; i < number_of_rows; i++)
{
a = malloc(number_of_rows * sizeof *a);
if(a != NULL)
{
for(j = 0; j < number_of_rows; j++)
{
/* might as well use your other array to initialise with */
a[j] = testvector[j];
}
}
else
{
printf("You're running low on memory - you should do something.\n");
exit(1);
}
}
}
else
{
printf("You're running low on memory - you should do something.\n");
exit(1);
}

banmul(a, number_of_rows, 2, 1, x, b);

printf("\nFinishing program now.\n\n");
return(0);
}

/* Here comes the NR functions */
....
-------------

Compiler errors:
1) error C2059: syntax error : ']'

This is really strange! I don't understand that - it used to work before
I think.... Error 1) happen in the line: testvector[][number_of_rows] =
.... { {3.,1.,0.,0.,......},{4.,1.,..} etc};


2) error C2143: syntax error : missing ';' before 'type'
Error 2) Happens just after the testvector (matrix actually) is defined
- in the line: float **a = malloc(number_of_rows * sizeof *a);

Therefore error 2 is probably connected to error 1.

 

[=?ISO-8859-1?Q?Martin_J=F8rgensen?=]

You are missing <stdio.h>

Ok, I added it now (actually I also had it at some earlier moment, but
then I removed it since I couldn't see any compiler warnings/errors
without stdio.h).
-snip-

/* Start of main program */

int main(void)
{
/* declaration */
float **a;
float testvector[number_of_rows][number_of_rows];

You intend to initialize this array, you should consider doing that upon
declaration.

You're right... The reason I didn't was because I had a problem with
doing malloc earlier if I had some code before... I don't know why...
Why isn't it possible to declare the array and later initialize it as
with other variables?

int i, j;
float x[number_of_rows], b[number_of_rows]; // remember that we're
solving "b = A * x"
testvector[][number_of_rows] = // number_of_cols = number_of_rows
{
{3.,1.,0.,0.,0.,0.,0.},
{4.,1.,5.,0.,0.,0.,0.},
{9.,2.,6.,5.,0.,0.,0.},
{0.,3.,5.,8.,9.,0.,0.},
{0.,0.,7.,9.,3.,2.,0.},
{0.,0.,0.,3.,8.,4.,6.},
{0.,0.,0.,0.,2.,4.,4.}
};

This is not a proper initialization, that is why the compiler complained.

Yeah, but the compiler error message was really lousy. It complained
about a "]"-sign... I don't understand why...

You are re-declaring a. Try: a = malloc(number_of_rows * sizeof *a);

Ok. Works. And a is a pointer to pointer, right?

-snip the rest-

This is really weird. My program now compiles but stops execution
because of an "unhandled win32 exception":


Unhandled exception at 0x004155a3 in
band_diag_system_(numerical_recipes).exe: 0xC0000005: Access violation
reading location 0xabababaf.

I think I found something wrong in the code from "numerical recipes in
C"... But it's strange... My code (just copy/paste):

-----
/* Including header files */
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include "nrutil.h"

/* definitions */
#define number_of_rows 7
#define TINY 1.0e-20
#define SWAP(a,b) {dum=(a);(a)=(b);(b)=dum;}

/* prototypes */
void banmul(float **a, unsigned long n, int left, int right, float x[],
float b[]);
void bandec(float **a, unsigned long n, int left, int right, float **al,
unsigned long indx[], float *d);
void banbks(float **a, unsigned long n, int left, int right, float **al,
unsigned long indx[], float b[]);

/* Start of main program */

int main(void)
{
/* declaration */
float **a;
int i, j;
float x[number_of_rows], b[number_of_rows]; // remember that we're
solving "b = A * x"
float testvector[number_of_rows][number_of_rows] = // number_of_cols =
number_of_rows
{
{3.,1.,0.,0.,0.,0.,0.},
{4.,1.,5.,0.,0.,0.,0.},
{9.,2.,6.,5.,0.,0.,0.},
{0.,3.,5.,8.,9.,0.,0.},
{0.,0.,7.,9.,3.,2.,0.},
{0.,0.,0.,3.,8.,4.,6.},
{0.,0.,0.,0.,2.,4.,4.}
};


a = malloc(number_of_rows * sizeof *a);
if(a != NULL)
{
for(i = 0; i < number_of_rows; i++)
{
a = malloc(number_of_rows * sizeof *a);
if(a != NULL)
{
for(j = 0; j < number_of_rows; j++)
{
a[j] = testvector[j];
}
}
else
{
printf("You're running low on memory - you should do something.\n");
exit(1);
}
}
}
else
{
printf("You're running low on memory - you should do something.\n");
exit(1);
}

banmul(a, number_of_rows, 2, 1, x, b);

printf("\nFinishing program now.\n\n");
return(0);
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////
// m1= left, m2= right.
////////////////////////////////////////////////////////////////////////////////////////////////////////////

void banmul(float **a, unsigned long n, int left, int right, float x[],
float b[])
{
unsigned long i,j,k,tmploop;
for (i=1;i<=n;i++) {
k=i-left-1;
// 2 lines extra code from numerical recipes in C
// http://www.library.cornell.edu/nr/bookcpdf/c2-4.pdf
}
}
----

When I debug and get to the line in banmul that says: k=i-left-1, then i
is 1, left is 2 and from that we subtract 1. The total sum is zero so k
should be = 0. But it doesn't become 0 - it becomes 4294967294 according
to my debugger. I don't understand that... Any explanation?

I assume I just need to cast, like writing: "k = (unsigned
long)i-left-1"; but I would like to know why this exception error happens...


Best regards / Med venlig hilsen
Martin Jørgensen

 

[Michael Mair]

Ok, I added it now (actually I also had it at some earlier moment, but
then I removed it since I couldn't see any compiler warnings/errors
without stdio.h).

Then
- invoke your compiler in strict C89 or C99 mode
- turn the warning level up to the maximum
If this does not help or if the former is not possible, I'd
feel a certain amount of distrust towards the compiler.

-snip-

/* Start of main program */

int main(void)
{
/* declaration */
float **a;
float testvector[number_of_rows][number_of_rows];

You intend to initialize this array, you should consider doing that upon
declaration.

You're right... The reason I didn't was because I had a problem with
doing malloc earlier if I had some code before... I don't know why...
Why isn't it possible to declare the array and later initialize it as
with other variables?

It just is not.
One of the problems: Using an array name yields a pointer to the
array's first element in most contexts. Assigning an
"array literal" to a pointer obviously may lead to problems...

int i, j;
float x[number_of_rows], b[number_of_rows]; // remember that we're
solving "b = A * x"
testvector[][number_of_rows] = // number_of_cols = number_of_rows
{
{3.,1.,0.,0.,0.,0.,0.},
{4.,1.,5.,0.,0.,0.,0.},
{9.,2.,6.,5.,0.,0.,0.},
{0.,3.,5.,8.,9.,0.,0.},
{0.,0.,7.,9.,3.,2.,0.},
{0.,0.,0.,3.,8.,4.,6.},
{0.,0.,0.,0.,2.,4.,4.}
};

This is not a proper initialization, that is why the compiler complained.

Yeah, but the compiler error message was really lousy. It complained
about a "]"-sign... I don't understand why...

Because [] is only allowed in _declarations_.
From a parser's point of view, the closing bracket is not allowed
to appear after the opening bracket -- you need something which can
be interpreted as array index in between.

Ok. Works. And a is a pointer to pointer, right?
Yes.

-snip the rest-

This is really weird. My program now compiles but stops execution
because of an "unhandled win32 exception":


Unhandled exception at 0x004155a3 in
band_diag_system_(numerical_recipes).exe: 0xC0000005: Access violation
reading location 0xabababaf.

And this location is part of ??

I think I found something wrong in the code from "numerical recipes in
C"... But it's strange... My code (just copy/paste):

Perfect. Copy/paste, I mean

-----
/* Including header files */
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include "nrutil.h"

Compiler error; I do not have nrutil.h -- it obviously is not
used here.

/* definitions */
#define number_of_rows 7
#define TINY 1.0e-20
#define SWAP(a,b) {dum=(a);(a)=(b);(b)=dum;}

Note that this swap is not exactly typesafe.

/* prototypes */
void banmul(float **a, unsigned long n, int left, int right, float x[],
float b[]);
void bandec(float **a, unsigned long n, int left, int right, float **al,
unsigned long indx[], float *d);
void banbks(float **a, unsigned long n, int left, int right, float **al,
unsigned long indx[], float b[]);

/* Start of main program */

int main(void)
{
/* declaration */
float **a;
int i, j;
float x[number_of_rows], b[number_of_rows]; // remember that we're
solving "b = A * x"

Please do not use // comments when posting code to usenet -- line
breaks may change semantics.

float testvector[number_of_rows][number_of_rows] = // number_of_cols
= number_of_rows
{
{3.,1.,0.,0.,0.,0.,0.},
{4.,1.,5.,0.,0.,0.,0.},
{9.,2.,6.,5.,0.,0.,0.},
{0.,3.,5.,8.,9.,0.,0.},
{0.,0.,7.,9.,3.,2.,0.},
{0.,0.,0.,3.,8.,4.,6.},
{0.,0.,0.,0.,2.,4.,4.}
};


a = malloc(number_of_rows * sizeof *a);
if(a != NULL)
{
for(i = 0; i < number_of_rows; i++)
{
a = malloc(number_of_rows * sizeof *a);
if(a != NULL)
{
for(j = 0; j < number_of_rows; j++)
{
a[j] = testvector[j];
}
}
else
{
printf("You're running low on memory - you should do
something.\n");
exit(1);
}
}
}
else
{
printf("You're running low on memory - you should do
something.\n");
exit(1);
}

1) exit(EXIT_FAILURE) is portable, exit(1) is not.
2) If you change the logic, the nesting depth decreases:
if (a==NULL) {/*bleargh*/}
for ....
....
if (a==NULL) {/*bleargh*/}
for ....
}

banmul(a, number_of_rows, 2, 1, x, b);

Note: x and b are still uninitialized here.

printf("\nFinishing program now.\n\n");
return(0);
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////

// m1= left, m2= right.
////////////////////////////////////////////////////////////////////////////////////////////////////////////


void banmul(float **a, unsigned long n, int left, int right, float x[],
float b[])
{
unsigned long i,j,k,tmploop;
for (i=1;i<=n;i++) {
k=i-left-1;
// 2 lines extra code from numerical recipes in C
// http://www.library.cornell.edu/nr/bookcpdf/c2-4.pdf
}
}
----

When I debug and get to the line in banmul that says: k=i-left-1, then i
is 1, left is 2 and from that we subtract 1. The total sum is zero so k
should be = 0. But it doesn't become 0 - it becomes 4294967294 according
to my debugger. I don't understand that... Any explanation?

Yes. 1-2-1 is -2 which wraps around to 4294967296-2.
Everything fine there.

I assume I just need to cast, like writing: "k = (unsigned
long)i-left-1"; but I would like to know why this exception error
happens...

Won't help.

Are you sure that the above is the faulty code?
I cannot see any problem with it.

Q: Does the error also occur for the following variant?

/* Including header files */
#include <stdlib.h>
#include <stdio.h>
#include <string.h>

/* definitions */
#define number_of_rows 7

/* prototypes */
void banmul(float **a, unsigned long n, int left, int right, float x[],
float b[]);

/* Start of main program */

int main(void)
{
/* declaration */
float **a;
int i;
float x[number_of_rows] = {0.f}, b[number_of_rows] = {0.f}; /*
remember that we're solving "b = A * x"*/
float testvector[number_of_rows][number_of_rows] = /*
number_of_cols = number_of_rows*/
{
{3.,1.,0.,0.,0.,0.,0.},
{4.,1.,5.,0.,0.,0.,0.},
{9.,2.,6.,5.,0.,0.,0.},
{0.,3.,5.,8.,9.,0.,0.},
{0.,0.,7.,9.,3.,2.,0.},
{0.,0.,0.,3.,8.,4.,6.},
{0.,0.,0.,0.,2.,4.,4.}
};


a = malloc(number_of_rows * sizeof *a);
if (a == NULL)
{
fprintf(stderr, "%s, %d: You're running low on memory -"
" you should do something.\n",
__FILE__, __LINE__);
exit(EXIT_FAILURE);
}

a[0] = malloc(number_of_rows * number_of_rows * sizeof *a[0]);
if (a[0] == NULL)
{
fprintf(stderr, "%s, %d: You're running low on memory -"
" you should do something.\n",
__FILE__, __LINE__);
exit(EXIT_FAILURE);
}
memcpy(&a[0][0], &testvector[0][0], sizeof testvector);

for(i = 1; i < number_of_rows; i++)
{
a = a[0] + i*number_of_rows;
}

banmul(a, number_of_rows, 2, 1, x, b);

printf("\nFinishing program now.\n\n");

free(a[0]);
free(a);
return(0);
}


void banmul(float **a, unsigned long n, int left, int right, float x[],
float b[])
{
unsigned long i,j,k,tmploop;
for (i=1;i<=n;i++) {
k=i-left-1;
/* 2 lines extra code from numerical recipes in C
** http://www.library.cornell.edu/nr/bookcpdf/c2-4.pdf
*/
}
}


Cheers
Michael

 

[=?ISO-8859-15?Q?Martin_J=F8rgensen?=]

Martin Jørgensen schrieb: -snip-



Then
- invoke your compiler in strict C89 or C99 mode
- turn the warning level up to the maximum
If this does not help or if the former is not possible, I'd
feel a certain amount of distrust towards the compiler.

Hmmm.. Strange. I tried it again and it doesn't complain. Which
functions are needed by stdio.h?

I don't know how to turn up the warning levels or put the compiler in
strict C89 or C99 mode... The compiler is Microsoft visual studio 2005
c/c++. I don't even know how to invoke it with the command-line or what
the file-name for the compiler is. I assume I could just type
(compiler-name.exe) /help in a dos cmd-prompt window to learn about the
options...

-snip-

It just is not.
One of the problems: Using an array name yields a pointer to the
array's first element in most contexts. Assigning an
"array literal" to a pointer obviously may lead to problems...

Ok, makes sense.

-snip-

Yeah, but the compiler error message was really lousy. It complained
about a "]"-sign... I don't understand why...

Because [] is only allowed in _declarations_.
From a parser's point of view, the closing bracket is not allowed
to appear after the opening bracket -- you need something which can
be interpreted as array index in between.

Ok, understood now. Ofcourse I now see that it needs to reserve space
for the array, logically...

-snip-

And this location is part of ??

I think (at least hope) you should be able to get the same message with
the new code I will be posting now...

-snip-

Compiler error; I do not have nrutil.h -- it obviously is not
used here.

It is public domain: http://www.nr.com/pubdom/nrutil_nr.h.txt

But you won't need it so far I think...

Note that this swap is not exactly typesafe.
Why?

-snip-



Yes. 1-2-1 is -2 which wraps around to 4294967296-2.
Everything fine there.

Ok, but actually there's something wrong with the "numerical recipes"
code, I think then - it seems to try to access a pointer with index k
which is far away from the malloc'ed memory... Explanation follows later...

Won't help.

Are you sure that the above is the faulty code?
I cannot see any problem with it.

Q: Does the error also occur for the following variant?

Yep, but you didn't have a chance to see that before now (I hoped it
wasn't necessary to make it more complicated but I now post about 3
extra lines to the banmul function) + 2 definitions from the (public
domain) header file which is included to make things easy for you...

I took your proposal and modified my program. I made a small error with
the matrix I think - but it's corrected now (I think it should be a
"compact matrix" form for the bandmul function to work properly).

Try this program (just copy/paste) and I hope you'll see my error (I
also included the original matrix outcommented for clarity reasons):
----------

/* Including header files */
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
//#include "nrutil.h"

/* definitions */
#define number_of_rows 7
#define number_of_cols 4
#define TINY 1.0e-20
#define SWAP(a,b) {dum=(a);(a)=(b);(b)=dum;}

static long lminarg1,lminarg2;
#define LMIN(a,b) (lminarg1=(a),lminarg2=(b),(lminarg1) < (lminarg2) ?\
(lminarg1) : (lminarg2))

static long lmaxarg1,lmaxarg2;
#define LMAX(a,b) (lmaxarg1=(a),lmaxarg2=(b),(lmaxarg1) > (lmaxarg2) ?\
(lmaxarg1) : (lmaxarg2))

/* prototypes */
void banmul(float **a, unsigned long n, int left, int right, float x[],
float b[]);


/* Start of main program */
int main(void)
{
/* declaration */
float **a;
int i;

/* remember that we're solving "b = A * x" */
float x[number_of_rows] = {1., 2., 3., 7., 6., 5., 4.};
float b[number_of_rows] = {0.f};
float testmatrix[][number_of_rows] =
{
{0.,0.,3.,1.},
{0.,4.,1.,5.},
{9.,2.,6.,5.},
{3.,5.,8.,9.},
{7.,9.,3.,2.},
{3.,8.,4.,6.},
{2.,4.,4.,0.}
};

//{
// {3.,1.,0.,0.,0.,0.,0.},
// {4.,1.,5.,0.,0.,0.,0.},
// {9.,2.,6.,5.,0.,0.,0.},
// {0.,3.,5.,8.,9.,0.,0.},
// {0.,0.,7.,9.,3.,2.,0.},
// {0.,0.,0.,3.,8.,4.,6.},
// {0.,0.,0.,0.,2.,4.,4.}
//};

a = malloc(number_of_rows * sizeof *a);
if (a == NULL)
{
fprintf(stderr, "%s, %d: You're running low on memory -"
" you should do something.\n",
__FILE__, __LINE__);
exit(EXIT_FAILURE);
}

a[0] = malloc(number_of_rows * number_of_rows * sizeof *a[0]);
if (a[0] == NULL)
{
fprintf(stderr, "%s, %d: You're running low on memory -"
" you should do something.\n",
__FILE__, __LINE__);
exit(EXIT_FAILURE);
}
memcpy(&a[0][0], &testmatrix[0][0], sizeof testmatrix);

for(i = 1; i < number_of_rows; i++)
{
a = a[0] + i*number_of_rows;
}

banmul(a, number_of_rows, 2, 1, x, b);
/* the number 2 = number of rows to the left of diagonal (see matrix)
and the number 1 = number rows to the right of the diagonal (see matrix) */

printf("\nFinishing program now.\n\n");

free(a[0]);
free(a);
return(0);
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////
// m1= left, m2= right.
////////////////////////////////////////////////////////////////////////////////////////////////////////////


void banmul(float **a, unsigned long n, int left, int right, float x[],
float b[])
{
unsigned long i,j,k,tmploop;
for (i=1;i<=n;i++) {
k=i-left-1;
tmploop=LMIN(left+right+1,n-k);
b=0.0;
for (j=LMAX(1,1-k); j<=tmploop;j++) b += a[j]*x[j+k];
}
}


----------

It put a breakpoint in the banmul function. I get that k = 4294967294
(first time the computer enters the function, so i=1). tmploop = 4. b[1]
= 0 (ofcourse). Now I wanted to try my visual studio compiler function
which prints something like the following out if I change something and
continue debugging:

-------- Edit and Continue build started --------

--------------------- Done ----------------------

My problem is that I wanted to insert lines that calculates 1-k =
1-4294967294 such as:

tmploop = 1-k; // = 1-4294967294

With the debugger I then want to see tmploop (I don't know what the
result is because I guess I would have to think "binary" then). Then I
would like to continue and put tmploop back to its original value:

tmploop=LMIN(left+right+1,n-k);

Is that possible with this "edit and continue build"-stuff with my compiler?

Ok, the real problem AFAICS is last line: b += a[j]*x[j+k],
because at this moment k is 4294967294 ! So it's trying to multiply
a[j] with x[a very big number outside the real data]...

What do you think?

If you or somebody else want there's a brief and okay description at
http://www.library.cornell.edu/nr/bookcpdf/c2-4.pdf (banmul described on
page 3 (p.52) ). I hope I didn't miss anything from there...

I hope there are some clever guys in this forum that can solve the
problem


Best regards / Med venlig hilsen
Martin Jørgensen

 

[Michael Mair]

Hmmm.. Strange. I tried it again and it doesn't complain. Which
functions are needed by stdio.h?

I don't know how to turn up the warning levels or put the compiler in
strict C89 or C99 mode... The compiler is Microsoft visual studio 2005
c/c++. I don't even know how to invoke it with the command-line or what
the file-name for the compiler is. I assume I could just type
(compiler-name.exe) /help in a dos cmd-prompt window to learn about the
options...

<OT>Right-click on your project, go to Properties.
Within the properties window, click on "C/C++"/General, turn up the
warning level to W4; consider switching on Treat Warnings As Errors
as long as you are still learning.
Go on to "C/C++"/Language. Disable language extensions. Consider
making char defaulting to unsigned char.
Go on to "C/C++"/Advanced. Set to Compile As C.

Note:<stdio.h> is often included by "stdafx.h".
Try to keep your includes where they belong.

I did not actually try it out because I do not have VC 2005 at hand.

I think (at least hope) you should be able to get the same message with
the new code I will be posting now...

-snip-


It is public domain: http://www.nr.com/pubdom/nrutil_nr.h.txt

But you won't need it so far I think...


Why?

Let us make it a game:
- If dum is of your largest integer type, I swap double variables
containing values too large to be represented by dum.
- If dum is of your largest floating point type and this type has
less mantissa bits than your largest integer type has bits, I
use variables of that type containing a sufficiently large number
with enough bits set to be not representable.
<OT>With MS VC++, you have only 64 bit long doubles, so I can
leave it at using 64 bit integers</OT>.

If the above conversions can take place implicitly, you do not
get the least warning that your swap throws away everything.
In addition: Where do we get "dum"? Is it a global variable?

-snip-


Ok, but actually there's something wrong with the "numerical recipes"
code, I think then - it seems to try to access a pointer with index k
which is far away from the malloc'ed memory... Explanation follows later...

Okay. Throw in a check that your indices are in the valid range.
Provide enough context and let the programme die gracefully.

Won't help.

Are you sure that the above is the faulty code?
I cannot see any problem with it.

Q: Does the error also occur for the following variant?

Yep, but you didn't have a chance to see that before now (I hoped it
wasn't necessary to make it more complicated but I now post about 3
extra lines to the banmul function) + 2 definitions from the (public
domain) header file which is included to make things easy for you...

I took your proposal and modified my program. I made a small error with
the matrix I think - but it's corrected now (I think it should be a
"compact matrix" form for the bandmul function to work properly).
Okay.

Try this program (just copy/paste) and I hope you'll see my error (I
also included the original matrix outcommented for clarity reasons):
----------

/* Including header files */
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
//#include "nrutil.h"

/* definitions */
#define number_of_rows 7
#define number_of_cols 4
#define TINY 1.0e-20
#define SWAP(a,b) {dum=(a);(a)=(b);(b)=dum;}

static long lminarg1,lminarg2;
#define LMIN(a,b) (lminarg1=(a),lminarg2=(b),(lminarg1) < (lminarg2) ?\
(lminarg1) : (lminarg2))

static long lmaxarg1,lmaxarg2;
#define LMAX(a,b) (lmaxarg1=(a),lmaxarg2=(b),(lmaxarg1) > (lmaxarg2) ?\
(lmaxarg1) : (lmaxarg2))

/* prototypes */
void banmul(float **a, unsigned long n, int left, int right, float x[],
float b[]);


/* Start of main program */
int main(void)
{
/* declaration */
float **a;
int i;

/* remember that we're solving "b = A * x" */
float x[number_of_rows] = {1., 2., 3., 7., 6., 5., 4.};
float b[number_of_rows] = {0.f};
float testmatrix[][number_of_rows] =
{
{0.,0.,3.,1.},
{0.,4.,1.,5.},
{9.,2.,6.,5.},
{3.,5.,8.,9.},
{7.,9.,3.,2.},
{3.,8.,4.,6.},
{2.,4.,4.,0.}
};

Are you sure that this is the memory layout you want to have?

//{
// {3.,1.,0.,0.,0.,0.,0.},
// {4.,1.,5.,0.,0.,0.,0.},
// {9.,2.,6.,5.,0.,0.,0.},
// {0.,3.,5.,8.,9.,0.,0.},
// {0.,0.,7.,9.,3.,2.,0.},
// {0.,0.,0.,3.,8.,4.,6.},
// {0.,0.,0.,0.,2.,4.,4.}
//};

a = malloc(number_of_rows * sizeof *a);
if (a == NULL)
{
fprintf(stderr, "%s, %d: You're running low on memory -"
" you should do something.\n",
__FILE__, __LINE__);
exit(EXIT_FAILURE);
}

a[0] = malloc(number_of_rows * number_of_rows * sizeof *a[0]);

Okay, this is more or less my mistake: If you change the size of
testmatrix, this will foul up the whole thing.
As you are only copying testmatrix, replace the line by
a[0] = malloc(sizeof testmatrix);

if (a[0] == NULL)
{
fprintf(stderr, "%s, %d: You're running low on memory -"
" you should do something.\n",
__FILE__, __LINE__);
exit(EXIT_FAILURE);
}
memcpy(&a[0][0], &testmatrix[0][0], sizeof testmatrix);

for(i = 1; i < number_of_rows; i++)
{
a = a[0] + i*number_of_rows;

Depending on what kind of matrix you are trying to arrive at,
this is not necessarily correct. Replace number_of_rows by the actual
number of columns (which is
sizeof testmatrix[0]/ sizeof testmatrix[0][0]).

If you do not want to keep testmatrix as it is at initialization,
then you also can skip the second allocation and the memcpy() and
just set up a to point to the start of the "i'th row".

}


banmul(a, number_of_rows, 2, 1, x, b);
/* the number 2 = number of rows to the left of diagonal (see matrix)
and the number 1 = number rows to the right of the diagonal (see
matrix) */

printf("\nFinishing program now.\n\n");

free(a[0]);
free(a);
return(0);
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////

// m1= left, m2= right.
////////////////////////////////////////////////////////////////////////////////////////////////////////////



void banmul(float **a, unsigned long n, int left, int right, float x[],
float b[])
{
unsigned long i,j,k,tmploop;
for (i=1;i<=n;i++) {
k=i-left-1;
tmploop=LMIN(left+right+1,n-k);
b=0.0;
for (j=LMAX(1,1-k); j<=tmploop;j++) b += a[j]*x[j+k];
}

Question: Why don't you just output the value of i, k, tmploop
within the outer loop and the value of j and b within the
inner loop?
Consider calculating and outputting the values as signed values and
unsigned values.
<OT>As you have a 64 bit integer type, you can introduce auxiliary
variables of that type, do the computations in these variables,
assign the value to the intended variables and output the auxiliary

}


----------

It put a breakpoint in the banmul function. I get that k = 4294967294
(first time the computer enters the function, so i=1). tmploop = 4. b[1]
= 0 (ofcourse). Now I wanted to try my visual studio compiler function
which prints something like the following out if I change something and
continue debugging:

-------- Edit and Continue build started --------

--------------------- Done ----------------------

My problem is that I wanted to insert lines that calculates 1-k =
1-4294967294 such as:

tmploop = 1-k; // = 1-4294967294

With the debugger I then want to see tmploop (I don't know what the
result is because I guess I would have to think "binary" then). Then I
would like to continue and put tmploop back to its original value:

tmploop=LMIN(left+right+1,n-k);

Is that possible with this "edit and continue build"-stuff with my
compiler?

Sorry, this is a compiler/debugger specific thing. I do not
work with the VC debugger that much, so I cannot help you there.
An appropriate newsgroup might be better for that.

Ok, the real problem AFAICS is last line: b += a[j]*x[j+k],
because at this moment k is 4294967294 ! So it's trying to multiply
a[j] with x[a very big number outside the real data]...

What do you think?

Not necessarily. The authors may rely on defined overflow
behavior. If j is large enough that the unsigned long representation
of j+k is >= 0, everything is alright. This seems to be the case
on first glance.
If your compiler implicitly promotes the result to __int64 or
something like that, you are in trouble; then you may want to
cast j+k to unsigned long to be on the safe side. However, first
try to find out what really happens before introducing casts.

If you or somebody else want there's a brief and okay description at
http://www.library.cornell.edu/nr/bookcpdf/c2-4.pdf (banmul described on
page 3 (p.52) ). I hope I didn't miss anything from there...

I hope there are some clever guys in this forum that can solve the
problem

Work with the above. If I have the time later on, I will have
a closer look at what actually happens (if you did not report
any results by then).


Cheers
Michael

 

[Martin Joergensen]

<OT>Right-click on your project, go to Properties.
Within the properties window, click on "C/C++"/General, turn up the
warning level to W4; consider switching on Treat Warnings As Errors
as long as you are still learning.

Nice... But gives plenty of warnings now. I guess this code is pretty old -
I hate those deprecation warnings... But for programs I make myself, this
looks very nice.

Go on to "C/C++"/Language. Disable language extensions. Consider
making char defaulting to unsigned char.

What does that setting (language extension) mean?

Go on to "C/C++"/Advanced. Set to Compile As C.

Done.

Note:<stdio.h> is often included by "stdafx.h".
Try to keep your includes where they belong.

Ok, but stdafx.h doesn't say anything to me. Where do you mean I should put
the include?

I did not actually try it out because I do not have VC 2005 at hand.
</OT>

Ok, thanks - detailed enough for me, because I could find all those
settings.

-snip-

Are you sure that this is the memory layout you want to have?

OMFG! It says in the text that the matrix should go from [1..n][1...4], but
my matrix goes from [0..n-1][0...3]! Might be that they just took the code
from somewhere else because this is a stupid way (I think). Problem solved
by changing a bit (implementing your suggesions from this posting) - I now
get the correct results out which is a major breakthrough!

But I have a problem with the free() statements you suggested - I'll post my
code in the end. I don't know what's wrong but you can probably very easy
tell (no compiler problems).


-snip (implementing your suggestions in the code)

void banmul(float **a, unsigned long n, int left, int right, float x[],
float b[])
{
unsigned long i,j,k,tmploop;
for (i=1;i<=n;i++) {
k=i-left-1;
tmploop=LMIN(left+right+1,n-k);
b=0.0;
for (j=LMAX(1,1-k); j<=tmploop;j++) b += a[j]*x[j+k];
}

Question: Why don't you just output the value of i, k, tmploop
within the outer loop and the value of j and b within the
inner loop?

You mean with printf(" text %format specifier", variable)-statements? I
hoped I could do it without having to recompile each time I modify such
small things... But ok.

Consider calculating and outputting the values as signed values and
unsigned values.

Don't need to. I get the correct results now.

<OT>As you have a 64 bit integer type, you can introduce auxiliary
variables of that type, do the computations in these variables,
assign the value to the intended variables and output the auxiliary

Sounds very interesting... Can you show an example or just briefly tell how
that can be done? Is it a compiler/debugger issue or a C-programming issue?
Sound to me like a compiler/debugger issue...

It put a breakpoint in the banmul function. I get that k = 4294967294
(first time the computer enters the function, so i=1). tmploop = 4. b[1]
= 0 (ofcourse). Now I wanted to try my visual studio compiler function
which prints something like the following out if I change something and
continue debugging:

-------- Edit and Continue build started --------

--------------------- Done ----------------------

My problem is that I wanted to insert lines that calculates 1-k =
1-4294967294 such as:

tmploop = 1-k; // = 1-4294967294

With the debugger I then want to see tmploop (I don't know what the
result is because I guess I would have to think "binary" then). Then I
would like to continue and put tmploop back to its original value:

tmploop=LMIN(left+right+1,n-k);

Is that possible with this "edit and continue build"-stuff with my
compiler?

Sorry, this is a compiler/debugger specific thing. I do not
work with the VC debugger that much, so I cannot help you there.
An appropriate newsgroup might be better for that.

Ok, I'll find out later...

Ok, the real problem AFAICS is last line: b += a[j]*x[j+k], because
at this moment k is 4294967294 ! So it's trying to multiply a[j] with
x[a very big number outside the real data]...

What do you think?

Not necessarily. The authors may rely on defined overflow
behavior. If j is large enough that the unsigned long representation
of j+k is >= 0, everything is alright. This seems to be the case
on first glance.

You're right... It seems to work that way...

If your compiler implicitly promotes the result to __int64 or
something like that, you are in trouble; then you may want to
cast j+k to unsigned long to be on the safe side. However, first
try to find out what really happens before introducing casts.


Work with the above. If I have the time later on, I will have
a closer look at what actually happens (if you did not report
any results by then).

As said: The program (almost) works now... I appreciate the comments a lot.
I just have one last thing that bugs me...... I run the program and without
debugging (ctrl+F5) I get "an unhandled win32 exception" due to the
free()-function in the end of the program. If I debug the program, setting a
breakpoint in the last command in main() - that is the return (success)
function, I get a window popping up saying that: "Windows has triggered a
breakpoint in (file).exe. This may be due to a corruption of the heap and
indicates a bug in (file).exe or any of the DLL's it has loaded."

The debugger then opens "osfinfo.c" which I have no idea about what does...
But it's something with the memory allocation because the problem happens at
the lines:

printf("\nFinishing program now.\n\n");
free(a[0]);
free(a); // here the error pops up
return(0); // here is the breakpoint

The problem is something with the malloc's I guess.... Just copy/paste my
code:
-------
/* Including header files */
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include "nrutil.h"
/* definitions */
#define number_of_rows 7+1
#define number_of_cols 4+1
#define number_of_eqs number_of_rows-1
#define TINY 1.0e-20
#define SWAP(a,b) {dum=(a);(a)=(b);(b)=dum;}

static long lminarg1,lminarg2;
#define LMIN(a,b) (lminarg1=(a),lminarg2=(b),(lminarg1) < (lminarg2) ?\
(lminarg1) : (lminarg2))
static long lmaxarg1,lmaxarg2;
#define LMAX(a,b) (lmaxarg1=(a),lmaxarg2=(b),(lmaxarg1) > (lmaxarg2) ?\
(lmaxarg1) : (lmaxarg2))

/* prototypes */
void banmul(float **a, unsigned long n, int left, int right, float x[],
float b[]);

/* Start of main program */
int main(void)
{
/* declaration */
float **a;
int i;
/* remember that we're solving "b = A * x" */
float x[number_of_rows] = {0., 1., 1., 1., 1., 1., 1., 1.}; // x[0] = dummy
float b[number_of_rows] = {0.f}; // b[0] = dummy
float testmatrix[][number_of_cols] =
{
{0., 0.,0.,0.,0.},
/* dummy row - testmatrix starts from index 1->n */
{0., 0.,0.,3.,1.},
{0., 0.,4.,1.,5.},
{0., 9.,2.,6.,5.},
{0., 3.,5.,8.,9.},
{0., 7.,9.,3.,2.},
{0., 3.,8.,4.,6.},
{0., 2.,4.,4.,0.}
};
a = malloc(number_of_rows * sizeof *a);
if (a == NULL)
{
fprintf(stderr, "%s, %d: You're running low on memory -"
" you should do something.\n",
__FILE__, __LINE__);
exit(EXIT_FAILURE);
}
a[0] = malloc(sizeof testmatrix); /* a[0] = malloc(number_of_rows *
number_of_rows * sizeof *a[0]); */
if (a[0] == NULL)
{
fprintf(stderr, "%s, %d: You're running low on memory -"
" you should do something.\n",
__FILE__, __LINE__);
exit(EXIT_FAILURE);
}
for(i = 1; i < number_of_rows; i++)
{
a = a[0] + i*(sizeof testmatrix[0]/ sizeof testmatrix[0][0]); /* a =
a[0] + i*number_of_columns; */
}
memcpy(&a[0][0], &testmatrix[0][0], sizeof testmatrix);
banmul(a, number_of_eqs, 2, 1, x, b); /* result of b = A*x is: 4, 10, 22,
25, 21, 21, 10.

printf("\nFinishing program now.\n\n");
free(a[0]);
free(a);
return(0);
}
////////////////////////////////////////////////
void banmul(float **a, unsigned long n, int left, int right, float x[],
float b[])
{
unsigned long i,j,k,tmploop;
for (i=1;i<=n;i++) {
k=i-left-1;
tmploop=LMIN(left+right+1,n-k);
b=0.0;
for (j=LMAX(1,1-k); j<=tmploop;j++) b += a[j]*x[j+k];
}
}
------

I was thinking about first I do this:

a = malloc(number_of_rows * sizeof *a);

Then I do this:

a[0] = malloc(sizeof testmatrix);

Then this:

free(a[0]);

And last:

free(a);


Perhaps the order is wrong?



Best regards / Med venlig hilsen
Martin Jørgensen

 

[Pedro Graca]

The problem is something with the malloc's I guess.... Just copy/paste my
code:
-------
/* Including header files */
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include "nrutil.h"
/* definitions */
#define number_of_rows 7+1

Oops, I expect something bad out of this.
Instead, try

#define NUMBER_OF_ROWS (7+1)

#define number_of_cols 4+1

/* Start of main program */
int main(void)
{
/* declaration */
float **a;
int i;
/* remember that we're solving "b = A * x" */
float x[number_of_rows] = {0., 1., 1., 1., 1., 1., 1., 1.}; // x[0] = dummy

a = malloc(number_of_rows * sizeof *a);

Aha!
This line, after preprocessing becomes

a = malloc(7+1 * sizeof *a);

<snip the rest (I'm at my brother's computer and it's really tough to
work with the portuguese keyboard, especially his version)>

 

[=?ISO-8859-1?Q?Martin_J=F8rgensen?=]

Pedro Graca wrote:
-snip-

Aha!
This line, after preprocessing becomes

a = malloc(7+1 * sizeof *a);

<snip the rest (I'm at my brother's computer and it's really tough to
work with the portuguese keyboard, especially his version)>

Thanks... I should have seen that...


Best regards / Med venlig hilsen
Martin Jørgensen