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.

Danijel
  • 8,198
  • 18
  • 69
  • 133
  • 7
    [An array of arrays is not the same as a pointer to a pointer](https://stackoverflow.com/questions/18440205/casting-void-to-2d-array-of-int-c/18440456#18440456). – Some programmer dude Feb 02 '22 at 13:58
  • 1
    Does this answer your question? [casting void\*\* to 2D array of int - C](https://stackoverflow.com/questions/18440205/casting-void-to-2d-array-of-int-c) – the busybee Feb 02 '22 at 14:00
  • 1
    That library is some seriously sloppy, bad code. It looks extremely data cache-unfriendly. – Lundin Feb 02 '22 at 14:15
  • 1
    Consider using a more recent SVD library, e.g. [this one](https://github.com/lucasmaystre/svdlibc). – einpoklum Feb 02 '22 at 17:02
  • @einpoklum I tried that one. Having some issues, check here: https://stackoverflow.com/questions/70970947/comparing-svd-in-c-with-matlab. – Danijel Feb 03 '22 at 12:10
  • @Danijel: Unfortunately, I'm not experienced with any SVD library, I actually just search GitHub for one... – einpoklum Feb 03 '22 at 12:12
  • @einpoklum Sure. I only meant to check the example main() C code, if input arrays are handled properly or similar. – Danijel Feb 03 '22 at 12:21
  • Your good to lose the sign on m and n ? You went from int to unsigned – Hogstrom Feb 03 '22 at 14:38

4 Answers4

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.)
einpoklum
  • 118,144
  • 57
  • 340
  • 684
  • I changed the function signature, and I got no warnings. The question is, does the function allow for this change. – Danijel Feb 02 '22 at 14:07
  • `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? – Aconcagua Feb 02 '22 at 14:15
  • 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. – Aconcagua Feb 02 '22 at 14:17
  • Both header and course have been changed, and recompiled. – Danijel Feb 02 '22 at 14:27
  • 1
    @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 Feb 02 '22 at 14:33
  • @Lundin That's presumably why the signature has been changed... – Aconcagua Feb 02 '22 at 14:50
  • Not sure if changing the sig is enough? I see the function tries to access the array out of bounds: there is a line `for (i=n;i>=1;i--) {` and n=4 in my case, so it breaks at runtime at `v[i][i]=1.0`. It should never be 4, but 0,1,2 or 3. Not sure how this can work at all... – Danijel Feb 02 '22 at 14:55
  • @Danijel: If you've changed both the header and the source of `svdcmp()`, you needed to clarify that in your question. I answered based on the assumption you only changed the header. Changing the signature of the source means rewriting the function. Without the function's source I can't tell whether that should work or not. – einpoklum Feb 02 '22 at 14:55
  • Sorry, of course I would change both. The func source is in the link. – Danijel Feb 02 '22 at 14:56
  • @Danijel It seems as if the function indeed is using one-based indices. You now need to peek into the documentation of the function what is actually expected to be provided as input. Should the arrays in both dimensions be larger by one? Do the indices 0 serve as sentinels then? Or has someone indeed been unaware of C's zero based indices??? If so, you might need to fix all those loops. – Aconcagua Feb 02 '22 at 15:16
  • @Danijel You might have mentioned more clearly that you are referring to *runtime* trouble, by the way, because so far I've been paying attention to compilation only... – Aconcagua Feb 02 '22 at 15:19
  • I've had compilation problems first, however, when I manage to compile without warnings, I get runtime err. Look at the line I mentioned above where the index goes out of bounds. – Danijel Feb 02 '22 at 15:23
  • This is getting out of hand, I will have to rephrase a question and provide a full minimal example. – Danijel Feb 02 '22 at 15:28
  • 2
    @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. – John Bollinger Feb 02 '22 at 16:30
  • @JohnBollinger This code ought to be trash in Fortran as well. To begin with it needs to be split up in several smaller functions. – Lundin Feb 03 '22 at 07:19
  • @Lundin, if "it needs to be split up" is the best you've got for a corresponding Fortran progenitor (which is not hypothetical, *Numerical Recipes in Fortran* having preceded *Numerical Recipes in C*), then I'm not impressed. Code style opinions and trends (both in Fortran and in C) were different when the code was written, to the extent that there is no justification for concluding that the style used reflects poorly on the professionalism or skill of the author of the original Fortran code. – John Bollinger Feb 03 '22 at 14:11
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.

Aconcagua
  • 24,880
  • 4
  • 34
  • 59
  • 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 Feb 02 '22 at 14:46
  • @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... – Aconcagua Feb 02 '22 at 15:02
  • ordering of arguments in the call to `svdcmp` is wrong – tstanisl Feb 02 '22 at 16:23
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.

John Bollinger
  • 160,171
  • 8
  • 81
  • 157
  • Thanks. My bet is that the original code somehow relies on the original data structure, therefore doesn't work with my changes. – Danijel Feb 03 '22 at 15:20
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 );
tstanisl
  • 13,520
  • 2
  • 25
  • 40