The middle matrix is not diagonal

[serpheroth]

I try to run svd on input 1.630981, -0.325926, 0.367789 2.747989, -2.976772, 2.575702 18.267998, -26.747061, 18.224825

The cuda version does not produce a diagonal matrix for the the matrix 'w' in the code. Any way to enforce the svd to produce the matrix 'w' as a diagonal matrix?

[ericjang]

are you getting this behavior on the CPU version or just the CUDA version? unfortunately I don't have access to CUDA-enabled computers until I get back to school to test this.

On Thu, May 29, 2014 at 5:20 PM, serpheroth notifications@github.com wrote:

I try to run svd on input 1.630981, -0.325926, 0.367789 2.747989, -2.976772, 2.575702 18.267998, -26.747061, 18.224825

The cuda version does not produce a diagonal matrix for the the matrix 'w' in the code. Any way to enforce the svd to produce the matrix 'w' as a diagonal matrix?

— Reply to this email directly or view it on GitHub https://github.com/ericjang/svd3/issues/2.

[serpheroth]

Both versions have this problem. Here is the cpu result I get by replacing rsqrt(x) with 1.0f / sqrt(x) for higher precision U: 0.032541 0.859952 0.509336 0.126920 -0.509036 0.851336 0.991379 0.036942 -0.125710 S: 37.488434 0.000000 0.000007 0.000000 0.432530 0.661659 -0.000001 0.000000 1.099874 V: 0.493815 0.353704 0.794380 -0.717685 -0.350040 0.601997 0.490994 -0.867390 0.080993 USV* : 1.630981 -0.325926 0.367789 2.747989 -2.976772 2.575702 18.267998 -26.747051 18.224823

[ericjang]

Very strange. I'll take a look at it tomorrow evening. Thanks for catching this !

On Thursday, May 29, 2014, serpheroth notifications@github.com wrote:

Both versions have this problem. Here is the result I get by replacing rsqrt(x) with 1.0f / sqrt(x) for higher precision U: 0.032541 0.859952 0.509336 0.126920 -0.509036 0.851336 0.991379 0.036942 -0.125710 S: 37.488434 0.000000 0.000007 0.000000 0.432530 0.661659 -0.000001 0.000000 1.099874 V: 0.493815 0.353704 0.794380 -0.717685 -0.350040 0.601997 0.490994 -0.867390 0.080993 USV* : 1.630981 -0.325926 0.367789 2.747989 -2.976772 2.575702 18.267998 -26.747051 18.224823

— Reply to this email directly or view it on GitHub https://github.com/ericjang/svd3/issues/2#issuecomment-44617431.

[serpheroth]

The problem is in the iteration number in jacobiEigenanlysis method. 4 iterations is not enough for make S close to diagonal for some matrix.

[ericjang]

thanks a lot for looking into this - ah, that does make sense. according to the paper,

10−15 iterations are typically sufficient to diagonalize S to within single-precision roundoff error.

4 iterations seemed to be sufficient for achieving within 0.0005 for 99.9% of the input matrices. so that was the number i went with. Out of curiosity, what are you using the 3x3 SVD for?

On Thu, Jun 19, 2014 at 11:45 PM, serpheroth notifications@github.com wrote:

The problem is in the iteration number in jacobiEigenanlysis method. 4 iterations is not enough for make S close to diagonal for some matrix.

— Reply to this email directly or view it on GitHub https://github.com/ericjang/svd3/issues/2#issuecomment-46649339.

[serpheroth]

I am working on a psysically based simulation project using strain limiting. I need to use 3x3 SVD to analyze the deformation of a tetrahedron. I also found your snow project interesting, I will play with it once I have access to a CUDA machine.

I also compare your 3x3 svd code with the code from book "numeric recipes in C", which can do any mxn svd, I found your code has similar performance as their code. If I increase the iteration count, your code is slower. Is it not optimized enough?

[sohale]

My S is not diagonal too. I use the same example code on CPU (see below). Increasing the number of iterations from 4 to 14 didn't help.

#include <iostream>
#include <svd3/svd3.h>

int main() {
    float a11, a12, a13, a21, a22, a23, a31, a32, a33;

    a11= -0.558253; a12 = -0.0461681; a13 = -0.505735;
    a21 = -0.411397; a22 = 0.0365854; a23 = 0.199707;
    a31 = 0.285389; a32 =-0.313789; a33 = 0.200189;

    float   u11, u12, u13,
            u21, u22, u23,
            u31, u32, u33;

    float   s11, s12, s13,
            s21, s22, s23,
            s31, s32, s33;

    float   v11, v12, v13,
            v21, v22, v23,
            v31, v32, v33;

    svd(
        a11, a12, a13, a21, a22, a23, a31, a32, a33,
        u11, u12, u13, u21, u22, u23, u31, u32, u33,
        s11, s12, s13, s21, s22, s23, s31, s32, s33,
        v11, v12, v13, v21, v22, v23, v31, v32, v33);

    std::cout << "result: A= ";
    std::cout << a11 <<" , "<< a12 <<" , "<< a13 <<" , "<< a21 <<" , "<< a22 <<" , "<< a23 <<" , "<< a31 <<" , "<< a32 <<" , "<< a33;
    std::cout << std::endl << "U= ";
    std::cout << u11 <<" , "<< u12 <<" , "<< u13 <<" , "<< u21 <<" , "<< u22 <<" , "<< u23 <<" , "<< u31 <<" , "<< u32 <<" , "<< u33;
    std::cout << std::endl << "S= ";
    std::cout << s11 <<" , "<< s12 <<" , "<< s13 <<" , "<< s21 <<" , "<< s22 <<" , "<< s23 <<" , "<< s31 <<" , "<< s32 <<" , "<< s33;
    std::cout << std::endl << "V= ";
    std::cout << v11 <<" , "<< v12 <<" , "<< v13 <<" , "<< v21 <<" , "<< v22 <<" , "<< v23 <<" , "<< v31 <<" , "<< v32 <<" , "<< v33;
    std::cout << std::endl;
}