I have written some code in MatLab (R2017b), and I am now trying to get the same results by writing the script's equivalent in C. For the code's purposes, it is important that some quantities are equal up to numerical precision (using double precision).
The quantity to process is a [900x522] matrix matrix_in. Using MatLab, I have written the contents of this matrix --preceded by its number of dimensions and the extent in each dimension-- as doubles (8 bytes per value) to the binary file data.bin, i.e:
dim = size(matrix_in);
N = length(size(dim));
fileID = fopen(file_out,'w');
fwrite(fileID,N,'double'); % number of dimensions, usually 2
fwrite(fileID,dim,'double'); % [nrows, ncols, ..]
fwrite(fileID,matrix_in,'double');
fclose(fileID);
Next, I want to read in the contents of this binary file in c-code. To this end, I have written the following code in visual studio 15:
#include <cstdio>
#include <cstdlib>
int main(int argc, char *argv[])
{
const char *fname = argv[1];
FILE * myPtr;
myPtr = fopen(fname, "r");
// Test if input file exists, abort if it does not:
if (!myPtr) {
printf("\nUnable to open file: > %s <, are you sure it exists? Aborting program.\n", fname);
return -1;
}
/* Read in values from binary file here
Format:
- First byte: N = number of dimensions of matrix (bytes 1:8)
- N*8 next bytes: dim = extent of matrix in all N dimensions (bytes 9:(n+1)*8)
- prod(dim)*8 next bytes: A = values of matrix (column-wise) (8 bytes per value)
*/
// Read N:
double n;
fread(&n, sizeof(double), 1, myPtr);
int N = (int)n;
printf("The value of N = %d\n", N); // This prints 2, as it should
// Read dim:
double *dim = (double *)malloc(N * sizeof(double));
fread(dim, sizeof(double), N, myPtr);
printf("The values of dim are %f,%f.\n", dim[0],dim[1]); // This prints 900.000000,522.000000, as it should
// Read A:
int Nel = (int)(dim[0] * dim[1]); // Number of elements in A;
double *A = (double *)malloc(Nel * sizeof(double)); // Allocate memory to store A;
fread(A, sizeof(double), 150, myPtr); // Read in values to A;
printf("There are %d elements in array A.\n", Nel); // This prints 469800, as it should
// This is where the incorrectness shows:
for (int i = 0; i < 136; i++) {
printf("%d: %.8f: %p \n", i, A[i],*A+i);
}
// Close file, deallocate memory, and exit
fclose(myPtr);
free(dim);
free(A);
printf("End of program reached.\n");
return 0;
}
So what goes wrong? Well, when I build the project (using Visual Studio 15), and run the executable with the .bin file, the output is as follows:
The value of N = 2
The values of dim are 900.000000,522.000000.
There are 469800 elements in array A.
0: -0.19632467: BFC9212AADC25959
1: -0.05302308: 3FE9B7B5548F69AA
2: 0.12151127: 3FFCDBDAAA47B4D5
3: -0.11154920: 40066DED5523DA6A
// and so forth, until..
132: 0.11697594: 406079B7B5548F6A
133:-6277436239688796889934943738416506716724917685420501794480371793920.00: 406099B7B5548F6A
134:-6277438562204192487878988888393020692503707483087375482269988814848.00: 4060B9B7B5548F6A
135:-6277438562204192487878988888393020692503707483087375482269988814848.00: 4060D9B7B5548F6A
This continues until the last value, i.e. A[Nel-1]. The first 133 values read in are exactly the same as in MatLab up to double precision, as you would expect. For some reason, it goes wrong afterwards. Now, the .bin file is exactly the right size for its contents, and when I read it into matlab instead of C, the results are fine. Does anybody have any clue as to what is going wrong in my C code?
