b2singular_value_decomposition.H

Go to the documentation of this file.

//---------------------------------------------------------------------------
// b2singular_value_decomposition.H --
//
//     Singular value decomposition.
//
// written by Thomas Ludwig.
//
// Copyright (c) 2011-2012 SMR Engineering & Development SA
//               2502 Bienne, Switzerland
//
// All Rights Reserved.  Proprietary source code.  The contents of
// this file may not be disclosed to third parties, copied or
// duplicated in any form, in whole or in part, without the prior
// written permission of SMR.
//---------------------------------------------------------------------------
 
#ifndef __B2_SINGULAR_VALUE_DECOMPOSITION_H__
#define __B2_SINGULAR_VALUE_DECOMPOSITION_H__
 
#include <algorithm>
#include <cmath>
 
#include "b2eigen_decomposition.H"
#include "b2tensor_calculus.H"
 
namespace b2000 {
 
inline void svd_3_3(const double A[3][3], double U[3][3], double S[3], double V[3][3]) {
    // fprintf(stderr, "A=numpy.array([\n");
    // fprintf(stderr, "    [%.15g, %.15g, %.15g],\n",
    //         A[0][0], A[1][0], A[2][0]);
    // fprintf(stderr, "    [%.15g, %.15g, %.15g],\n",
    //         A[0][1], A[1][1], A[2][1]);
    // fprintf(stderr, "    [%.15g, %.15g, %.15g],\n",
    //         A[0][2], A[1][2], A[2][2]);
    // fprintf(stderr, "])\n");
 
    // Compute A^T * A.
    double AA[3][3];
    inner_product_3_3_TN(A, A, AA);
    // fprintf(stderr, "AA=numpy.array([\n");
    // fprintf(stderr, "    [%.15g, %.15g, %.15g],\n",
    //         AA[0][0], AA[1][0], AA[2][0]);
    // fprintf(stderr, "    [%.15g, %.15g, %.15g],\n",
    //         AA[0][1], AA[1][1], AA[2][1]);
    // fprintf(stderr, "    [%.15g, %.15g, %.15g],\n",
    //         AA[0][2], AA[1][2], AA[2][2]);
    // fprintf(stderr, "])\n");
 
    // Compute sorted eigenvalues and eigenvectors of AA.
    eigen_decomposition(AA, S);
 
    // Make sure that eigenvalues are non-negative.
    double s1 = std::max(0.0, S[0]);
    double s2 = std::max(0.0, S[1]);
    double s3 = std::max(0.0, S[2]);
 
    // Store the square root of the sorted eigenvalues in S.
    S[0] = std::sqrt(s1);
    S[1] = std::sqrt(s2);
    S[2] = std::sqrt(s3);
 
    // Store the sorted eigenvectors in V^T.
    std::copy(AA[0], AA[0] + 3 * 3, V[0]);
    // fprintf(stderr, "V_svd=numpy.array([\n");
    // fprintf(stderr, "    [%.15g, %.15g, %.15g],\n",
    //         V[0][0], V[1][0], V[2][0]);
    // fprintf(stderr, "    [%.15g, %.15g, %.15g],\n",
    //         V[0][1], V[1][1], V[2][1]);
    // fprintf(stderr, "    [%.15g, %.15g, %.15g],\n",
    //         V[0][2], V[1][2], V[2][2]);
    // fprintf(stderr, "])\n");
 
    // Compute U * S = A * V^T.
    inner_product_3_3_NT(A, V, U);
    // fprintf(stderr, "AV=numpy.array([\n");
    // fprintf(stderr, "    [%.15g, %.15g, %.15g],\n",
    //         U[0][0], U[1][0], U[2][0]);
    // fprintf(stderr, "    [%.15g, %.15g, %.15g],\n",
    //         U[0][1], U[1][1], U[2][1]);
    // fprintf(stderr, "    [%.15g, %.15g, %.15g],\n",
    //         U[0][2], U[1][2], U[2][2]);
    // fprintf(stderr, "])\n");
 
    // Normalize the columns of U.
    for (int i = 0; i != 3; ++i) {
        const double s = 1.0 / std::sqrt(U[i][0] * U[i][0] + U[i][1] * U[i][1] + U[i][2] * U[i][2]);
        // {
        //     const double ss =
        //         U[i][0] * U[i][0] + U[i][1] * U[i][1] + U[i][2] * U[i][2];
        //     fprintf(stderr, "i=%d ss=%.15g s=%.15g\n", i, ss, s);
        //     fprintf(stderr, "U[%d]=[%.15g,%.15g,%.15g]\n", i,
        //             U[i][0], U[i][1], U[i][2]);
        // }
        U[i][0] *= s;
        U[i][1] *= s;
        U[i][2] *= s;
    }
    // fprintf(stderr, "U_svd=numpy.array([\n");
    // fprintf(stderr, "    [%.15g, %.15g, %.15g],\n",
    //         U[0][0], U[1][0], U[2][0]);
    // fprintf(stderr, "    [%.15g, %.15g, %.15g],\n",
    //         U[0][1], U[1][1], U[2][1]);
    // fprintf(stderr, "    [%.15g, %.15g, %.15g],\n",
    //         U[0][2], U[1][2], U[2][2]);
    // fprintf(stderr, "])\n");
}
 
};  // namespace b2000
 
#endif /* __B2_SINGULAR_VALUE_DECOMPOSITION_H__ */