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

Martin Jørgensen schrieb:
stathis gotsis wrote:

"Martin Jørgensen" <unoder.spam@xxxxxxxxxxxx> wrote in message
news:u4bmd3-v01.ln1@xxxxxxxxxxxxxx

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>

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.


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

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

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[i] = malloc(number_of_rows * sizeof *a[i]);
if(a[i] != NULL)
{
for(j = 0; j < number_of_rows; j++)
{
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);
}

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[i]==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[i] = 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
--
E-Mail: Mine is an /at/ gmx /dot/ de address.
.

Relevant Pages

  • Re: cpu type idea
    ... compiler related recently. ... int main ... float a; ... just to parallelize the vectror and tensor operations. ...
    (alt.lang.asm)
  • Re: getting avg of long and std::vector::size_type
    ... the compiler starts by evaluating a/b. ... int and b is an int, so it performs integer division on a and b. ... Now the compiler starts with the cast. ... first step is to convert a to a float. ...
    (comp.lang.cpp)
  • Re: problem with cast and unions
    ... Think of it as a compiler for Small programs. ... to a Small "cell" (actually an int). ... can anyone please tell me how is the structure assignment implemented ...
    (comp.lang.c)
  • Re: Abnormal program termination
    ... foo.c:4: warning: function declaration isn't a prototype ... Your version is legal under C90 rules, but int mainis more explicit ... float ave(), average; ...
    (comp.lang.c)
  • RTF Render with Scale
    ... float fHorzSizeInches = nHorzSize / 25.4f; ... IntPtr hdc = graphics.GetHdc; ... int nHorzSize = SafeNativeMethods.GetDeviceCaps(hdc, ... PHYSICALWIDTH = 110, ...
    (microsoft.public.dotnet.languages.csharp)