MATLAB Answers

Problem in converting SVD written in C to Mex

2 views (last 30 days)
Hariprasad
Hariprasad on 11 Dec 2014
Commented: Geoff Hayes on 16 Dec 2014
I'm trying to modify the code of SVD which is in C to Mex way of indexing. I have changed the indexing at all the necessary places but the results are not matching with the Matlab svd function results for the same set of input data. I have the attached codes here. SVD is the C code and SVD1 is the one I modified. Any help is greatly appreciated. Updated the SVD1 file.

  0 Comments

Sign in to comment.

Answers (1)

Geoff Hayes
Geoff Hayes on 12 Dec 2014
Edited: Geoff Hayes on 12 Dec 2014
Hariprasad - why have you changed the code to access elements in the 2D arrays (matrices) in a different manner? You should be able to use the original SVD code in MEX (if that is what you are trying to do?).
That being said, one of the first changes that you have made is to the following code (from SVD.txt)
for (i = 0; i < n; i++)
{
/* left-hand reduction */
l = i + 1;
rv1[i] = scale * g;
g = s = scale = 0.0;
if (i < m)
{
for (k = i; k < m; k++)
scale += fabs((double)a[k][i]);
% etc
where m is the number of rows, and n is the number of columns. The above code iterates over each column i and then calculates the sum scale of the absolute value of each element in column i starting from row i (so the main diagonal). You have replaced the sum with
for (k = i; k < m; k++)
scale += fabs(a[k + i*m]);
In the above code, a[k + i*m], let us suppose that i is zero, m is 10, and n is 10. Then the indices into a are
0 + 0*10 = 0
1 + 0*10 = 1
2 + 0*10 = 2
.
.
9 + 0*10 = 9
So we are not accessing the first column of a but rather the first row! What you want to do instead, is the following
for (k = i; k < m; k++)
scale += fabs(a[i + k*n]);
Now if we try this again, the indices into a are
0 + 0*10 = 0
0 + 1*10 = 10
0 + 2*10 = 20
.
.
0 + 9*10 = 90
So we are grabbing the first element which would a[0][0] and every tenth element after that which corresponds then to the first column.
Try modifying your code using this indexing scheme and see what happens. Also ensure that the original SVD code produces similar results to MATLAB's svd.
EDIT
You have to decide which indexing scheme to use given the original code. If you are accessing all elements in a single column, then use the method described above. If accessing all elements in a row, then you will need to use a different scheme.

  2 Comments

Hariprasad
Hariprasad on 15 Dec 2014
Thank You Geoff Hayes.
1. When I try to use the original SVD code, I get error in all the lines where the indexing is of the form [i][j]. The error says "subscript requires array or pointer type". Honestly I thought I'm sending in an array in a and hence thought this way of indexing is not compatible in MEX.
2. [k + i*m] doesn't traverse through each row of a column if m is no.of rows and n is no.of cols? I thought it that way and in this manner I got the results where most of the data match except for a few columns.
I have attached a .jpg where I have shown the indexing which I thought I was doing in this case(considering a is 11x11 matrix). As you have pointed out
0 + 0*10 = 0
1 + 0*10 = 1
2 + 0*10 = 2
.
.
9 + 0*10 = 9
As per the code,
for (k = i; k < m; k++)
scale += fabs((double)a[k][i]);
are we not traversing all the rows (k) of 'a' for the first column (i) in the loop??
Will also try the method as you suggested and get back.
Geoff Hayes
Geoff Hayes on 16 Dec 2014
For your first point - you would have to show how you are using this function from MEX (is that what you are doing?) or however you are calling dsvd. Please include the errors so we can figure out what is happening.
As for the second point, you have to be careful of what k and i refer to, which, in this case, are the rows and columns respectively as can be seen from a[k][i].

Sign in to comment.

Sign in to answer this question.

Tags