svd2x2.c

Go to the documentation of this file.

#include <math.h>


void Svd2x2(double *Theta, double *Sigma1, double *Sigma2,
    double *Phi, int *Sign1, int *Sign2, double A[2][2])
{
    double Temp, U[2][2], W[2][2];
    
    *Theta = 0.5*atan2(2*(A[0][0]*A[1][0] + A[0][1]*A[1][1]),
        (A[0][0]*A[0][0] + A[0][1]*A[0][1])
        - (A[1][0]*A[1][0] + A[1][1]*A[1][1]));
    *Phi = 0.5*atan2(2*(A[0][0]*A[0][1] + A[1][0]*A[1][1]),
        (A[0][0]*A[0][0] + A[1][0]*A[1][0])
        - (A[0][1]*A[0][1] + A[1][1]*A[1][1]));
    
    U[0][0] = cos(*Theta);
    U[1][0] = sin(*Theta);
    U[0][1] = -U[1][0];
    U[1][1] = U[0][0];
    
    W[0][0] = cos(*Phi);
    W[1][0] = sin(*Phi);
    W[0][1] = -W[1][0];
    W[1][1] = W[0][0];
    
    Temp = U[0][0]*(A[0][0]*W[0][0] + A[0][1]*W[1][0])
            + U[1][0]*(A[1][0]*W[0][0] + A[1][1]*W[1][0]);
    *Sigma1 = fabs(Temp);
    *Sign1 = (Temp < 0) ? -1 : 1;
    
    Temp = U[0][1]*(A[0][0]*W[0][1] + A[0][1]*W[1][1])
            + U[1][1]*(A[1][0]*W[0][1] + A[1][1]*W[1][1]);
    *Sigma2 = fabs(Temp);
    *Sign2 = (Temp < 0) ? -1 : 1;
}