1

I am using SVD implementation from Numerical Recipes. The signature of the function was:
void svdcmp(float **a, int m, int n, float w[], float **v);
but I changed it to:
void svdcmp( float a[4][4], unsigned int m, unsigned int n, float *w, float v[4][4] );

My main() has static arrays, so I called svdcmp this way (the GCC compiler does not complain):

float A[4][4] = ...;
float W[4];
float V[4][4];

svdcmp( A, 4, 4, W, V );

Results I get in A seem to be wrong. Am I calling svdcmp properly?

I need the algo only for 4x4 case, so I need to figure out how to modify it and still get good results.

8

4 Answers 4

1

Am I calling svdcmp properly?

No, you are not. The C-style "literal" 2D array is actually a 1-D array in terms of its storage, but with the compiler setting things up so that double-square-bracketing it - works. There is no array-of-pointers, A[0], A1, A[2], A[3], which you index.

But before asking us, you should have listened to your compiler, which tells you:

source>: In function 'foo':
<source>:9:13: warning: passing argument 1 of 'svdcmp' from incompatible pointer type [-Wincompatible-pointer-types]
    9 |     svdcmp( A, 4, 4, W, V );
      |             ^
      |             |
      |             float (*)[4]
<source>:1:21: note: expected 'float **' but argument is of type 'float (*)[4]'
    1 | void svdcmp(float **a, int m, int n, float w[], float **v);
      |             ~~~~~~~~^
<source>:9:25: warning: passing argument 5 of 'svdcmp' from incompatible pointer type [-Wincompatible-pointer-types]
    9 |     svdcmp( A, 4, 4, W, V );
      |                         ^
      |                         |
      |                         float (*)[4]
<source>:1:57: note: expected 'float **' but argument is of type 'float (*)[4]'
    1 | void svdcmp(float **a, int m, int n, float w[], float **v);

It knows what it's talking about!

More generally:

  • You can't change functions' signatures and expect them to work. If you can't use the original signature, you have to investigate why that is.
  • Please read: Why should I always enable compiler warnings? (in this case the compiler warns you even without turning warnings on, but still.)
Sign up to request clarification or add additional context in comments.

16 Comments

I changed the function signature, and I got no warnings. The question is, does the function allow for this change.
int **a and int (*a)[n] use the same syntax to access the data, so no changes required inside the function, even though technically (under the hoods) the access changes pretty much, but that remains transparent. If the signature is changed, one needs to change the data being passed to as well, but in given case the arguments look fine, matching the function signature. So what's wrong with?
Ah, I assume the function signature is changed in both header and source... Of course if that's the header to a library not being re-compiled as well matter changes pretty much, I agree.
@Danijel The difference is that float** assumes a fragmented pointer-based lookup table, which is a plain awful thing to use for something as CPU-intensive as this horrible function, since the only thing it achieves is to lag everything down to oblivion. I don't think you can salvage this code overall, it is intrinsically badly written - they reference the inner index over the outer index in nested for loops all over the place. The code was clearly written by some math theory person and not a professional C programmer with knowledge of computers.
@Lundin, as an alternative to "The code was clearly written by some math theory person and not a professional C programmer", I submit that given the source, there is a high likelihood that the code was rather converted from Fortran to C by someone whose C programming expertise did not extend to computational methods.
|
1

Originally your function parameters have been int**, now they are int(*)[4]. While accessing single elements in both cases use the same syntax under the hoods the access differs pretty much (double pointer indirection vs. single indirection and implicit offset calculation)!

The change of the signature is only valid if you change it in both header and source and recompile the library as well, otherwise you get pretty much into trouble.

I'd recommend to define a constant for the dimension, though, replacing those magic numbers (4).

Assuming m and n are sizes of the matrices you might change the signature yet a bit more to remain more flexible:

void svdcmp(size_t m, size_t n, float a[m][n], float w[n /* or m??? */], float v[m][n]);

The use of variable length arrays (VLA) allows to still use matrices of different sizes even though you might not need.

In any case: The arguments you provide to the changed signature within main indeed match the signature, so the call is fine – recommending to avoid the magic numbers, though, would rather recommend

svdcmp(sizeof(A)/sizeof(*A), sizeof(*A)/sizeof(**A), A, W, V);

// or without re-ordered parameters:
svdcmp(A, sizeof(A)/sizeof(*A), sizeof(*A)/sizeof(**A), W, V);

or you use the constant recommended above.

3 Comments

I've eliminated m and n using a #define, recompiled it, but still no luck. Something seems to be terribly wrong with the code.
@Danijel Cannot reproduce any trouble – just changed the signature as you did and everything has been fine (as expected) – now if you drop m and n from the signature you need to replace them within the function as well, of course – or you make your life easier by simply adding size_t const n = YourDefine, m = YourDefine; to get around that work...
ordering of arguments in the call to svdcmp is wrong
1

I changed it to:

void svdcmp( float a[4][4], unsigned int m, unsigned int n, float *w, float v[4][4] );

My main() has static arrays, so I called svdcmp this way (the GCC compiler does not complain):

float A[4][4] = ...;
float W[4];
float V[4][4];

svdcmp( A, 4, 4, W, V );

Results I get in A seem to be wrong. Am I calling svdcmp properly?

The function call you present is consistent with the prototype you present. If you also made the corresponding changes in the function definition and used the so-modified version of the function then the call itself should be ok.

However, the data structures representing the matrices are only analogous, not equivalent, and it is possible that the svdcmp() code relies on the original data structure. It is also possible that the original code is buggy, or that there is something wrong with the data you are passing, or perhaps even that your expectations for the result are incorrect.

1 Comment

Thanks. My bet is that the original code somehow relies on the original data structure, therefore doesn't work with my changes.
0

You should not change the signature of the function if you are not able to recompile it.

Otherwise, expect crashes because the compiled function receives data in unexpected format.

In the perfect world, the VLA interfaces similar to answer should be used.

If the interface is fixed you have to introduce a proxy object, an array of pointers to rows of A.

float A[4][4] = { ... };
float* Aproxy[] = { A[0], A[1], A[2], A[3] };

... similar thing for V

svdcmp( Aproxy, 4, 4, W, Vproxy );

Comments

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.