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

[Martin Jørgensen]

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("\nFinis hing 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

--
---------------------------------------------------------------------------
Home of Martin Jørgensen - http://www.martinjoergensen.dk

[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 *.

--
Richard Heathfield
"Usenet is a strange place" - dmr 29/7/1999
http://www.cpax.org.uk
email: rjh at above domain (but drop the www, obviously)

[Keith Thompson]

Martin Jørgensen <un*********@sp am.jay.net> writes:
[...]

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".

--
Keith Thompson (The_Other_Keit h) ks***@mib.org <http://www.ghoti.net/~kst>
San Diego Supercomputer Center <*> <http://users.sdsc.edu/~kst>
We must do something. This is something. Therefore, we must do this.

[Martin Jørgensen]

Richard Heathfield wrote:

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[i] = 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("\nFinis hing 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

--
---------------------------------------------------------------------------
Home of Martin Jørgensen - http://www.martinjoergensen.dk

[Martin Jørgensen]

Keith Thompson wrote:

Martin Jørgensen <un*********@sp am.jay.net> writes:
[...]

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

--
---------------------------------------------------------------------------
Home of Martin Jørgensen - http://www.martinjoergensen.dk

[stathis gotsis]

"Martin Jørgensen" <un*********@sp am.jay.net> wrote in message
news:q3******** ***@news.tdc.dk ...

I tried:
Hopefully you are including <stdlib.h>
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[i] = 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("\nFinis hing program now.\n\n");
Returning 0 would not be a bad idea.
}

[Richard Heathfield]

Martin Jørgensen said:

Richard Heathfield wrote:

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_o f_rows * sizeof *p);
if(p != NULL)
{
for(i = 0; i < number_of_rows; i++)
{
p[i] = malloc(number_o f_rows * sizeof *p[i]);
if(p[i] != NULL)
{
for(j = 0; j < number_of_rows; j++)
{
/* might as well use your other array to initialise with */
p[i][j] = a[i][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.

--
Richard Heathfield
"Usenet is a strange place" - dmr 29/7/1999
http://www.cpax.org.uk
email: rjh at above domain (but drop the www, obviously)

[Martin Jørgensen]

Richard Heathfield wrote:

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 :-)

This gives 3 problems:

1) warning C4047: 'initializing' : 'float **' differs in levels of
indirection from 'int' - this is the line: float **a =
malloc((numbe r_of_rows)*size of(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.

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_o f_rows * sizeof *a);
if(a != NULL)
{
for(i = 0; i < number_of_rows; i++)
{
a[i] = malloc(number_o f_rows * sizeof *a[i]);
if(a[i] != NULL)
{
for(j = 0; j < number_of_rows; j++)
{
/* might as well use your other array to initialise with */
a[i][j] = testvector[i][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("\nFinis hing 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_o f_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

--
---------------------------------------------------------------------------
Home of Martin Jørgensen - http://www.martinjoergensen.dk

[Martin Jørgensen]

stathis gotsis wrote:

"Martin Jørgensen" <un*********@sp am.jay.net> wrote in message
news:q3******** ***@news.tdc.dk ...

I tried:

Hopefully you are including <stdlib.h>

Damn... Thanks...

int main(void)
{
int i;

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[i] = 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...

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.

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

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

printf("\nFin ishing program now.\n\n");

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

--
---------------------------------------------------------------------------
Home of Martin Jørgensen - http://www.martinjoergensen.dk

[stathis gotsis]

"Martin Jørgensen" <un*********@sp am.jay.net> wrote in message
news:u4******** ****@news.tdc.d k...

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 */
You are missing <stdio.h>
#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_o f_rows * sizeof *a);
You are re-declaring a. Try: a = malloc(number_o f_rows * sizeof *a);
if(a != NULL)
{
for(i = 0; i < number_of_rows; i++)
{
a[i] = malloc(number_o f_rows * sizeof *a[i]);
if(a[i] != NULL)
{
for(j = 0; j < number_of_rows; j++)
{
/* might as well use your other array to initialise with */
a[i][j] = testvector[i][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("\nFinis hing 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_o f_rows * sizeof *a);

Therefore error 2 is probably connected to error 1.