Changeset 389 for cpp/frams/model/similarity

Timestamp:

05/24/15 21:03:46 (10 years ago)

Author:

Maciej Komosinski

Message:

Code formatting

Location:

cpp/frams/model/similarity/SVD

Files:

• lapack.cpp (modified) (13 diffs)
• lapack.h (modified) (1 diff)
• matrix_tools.cpp (modified) (2 diffs)
• matrix_tools.h (modified) (1 diff)

cpp/frams/model/similarity/SVD/lapack.cpp

@@ -65 +65 @@
 #endif
 
-template<typename _Tp, size_t fixed_size = 1024 / sizeof (_Tp) + 8 > class AutoBuffer
+template<typename _Tp, size_t fixed_size = 1024 / sizeof(_Tp) + 8 > class AutoBuffer
 {
 public:
-    typedef _Tp value_type;
-
-    //! the default constructor
-    AutoBuffer();
-    //! constructor taking the real buffer size
-    AutoBuffer(size_t _size);
-
-    //! the copy constructor
-    AutoBuffer(const AutoBuffer<_Tp, fixed_size>& buf);
-    //! the assignment operator
-    AutoBuffer<_Tp, fixed_size>& operator=(const AutoBuffer<_Tp, fixed_size>& buf);
-
-    //! destructor. calls deallocate()
-    ~AutoBuffer();
-
-    //! allocates the new buffer of size _size. if the _size is small enough, stack-allocated buffer is used
-    void allocate(size_t _size);
-    //! deallocates the buffer if it was dynamically allocated
-    void deallocate();
-    //! resizes the buffer and preserves the content
-    void resize(size_t _size);
-    //! returns the current buffer size
-    size_t size() const;
-    //! returns pointer to the real buffer, stack-allocated or head-allocated
-    operator _Tp* ();
-    //! returns read-only pointer to the real buffer, stack-allocated or head-allocated
-    operator const _Tp* () const;
+        typedef _Tp value_type;
+
+        //! the default constructor
+        AutoBuffer();
+        //! constructor taking the real buffer size
+        AutoBuffer(size_t _size);
+
+        //! the copy constructor
+        AutoBuffer(const AutoBuffer<_Tp, fixed_size>& buf);
+        //! the assignment operator
+        AutoBuffer<_Tp, fixed_size>& operator=(const AutoBuffer<_Tp, fixed_size>& buf);
+
+        //! destructor. calls deallocate()
+        ~AutoBuffer();
+
+        //! allocates the new buffer of size _size. if the _size is small enough, stack-allocated buffer is used
+        void allocate(size_t _size);
+        //! deallocates the buffer if it was dynamically allocated
+        void deallocate();
+        //! resizes the buffer and preserves the content
+        void resize(size_t _size);
+        //! returns the current buffer size
+        size_t size() const;
+        //! returns pointer to the real buffer, stack-allocated or head-allocated
+        operator _Tp* ();
+        //! returns read-only pointer to the real buffer, stack-allocated or head-allocated
+        operator const _Tp* () const;
 
 protected:
-    //! pointer to the real buffer, can point to buf if the buffer is small enough
-    _Tp* ptr;
-    //! size of the real buffer
-    size_t sz;
-    //! pre-allocated buffer. At least 1 element to confirm C++ standard reqirements
-    _Tp buf[(fixed_size > 0) ? fixed_size : 1];
+        //! pointer to the real buffer, can point to buf if the buffer is small enough
+        _Tp* ptr;
+        //! size of the real buffer
+        size_t sz;
+        //! pre-allocated buffer. At least 1 element to confirm C++ standard reqirements
+        _Tp buf[(fixed_size > 0) ? fixed_size : 1];
 };
 
@@ -110 +110 @@
 AutoBuffer<_Tp, fixed_size>::AutoBuffer()
 {
-    ptr = buf;
-    sz = fixed_size;
+        ptr = buf;
+        sz = fixed_size;
 }
 
@@ -117 +117 @@
 AutoBuffer<_Tp, fixed_size>::AutoBuffer(size_t _size)
 {
-    ptr = buf;
-    sz = fixed_size;
-    allocate(_size);
+        ptr = buf;
+        sz = fixed_size;
+        allocate(_size);
 }
 
@@ -125 +125 @@
 AutoBuffer<_Tp, fixed_size>::AutoBuffer(const AutoBuffer<_Tp, fixed_size>& abuf)
 {
-    ptr = buf;
-    sz = fixed_size;
-    allocate(abuf.size());
-    for (size_t i = 0; i < sz; i++)
-        ptr[i] = abuf.ptr[i];
+        ptr = buf;
+        sz = fixed_size;
+        allocate(abuf.size());
+        for (size_t i = 0; i < sz; i++)
+                ptr[i] = abuf.ptr[i];
 }
 
@@ -135 +135 @@
 AutoBuffer<_Tp, fixed_size>::operator=(const AutoBuffer<_Tp, fixed_size>& abuf)
 {
-    if (this != &abuf)
-    {
-        deallocate();
-        allocate(abuf.size());
-        for (size_t i = 0; i < sz; i++)
-            ptr[i] = abuf.ptr[i];
-    }
-    return *this;
+        if (this != &abuf)
+        {
+                deallocate();
+                allocate(abuf.size());
+                for (size_t i = 0; i < sz; i++)
+                        ptr[i] = abuf.ptr[i];
+        }
+        return *this;
 }
 
@@ -148 +148 @@
 AutoBuffer<_Tp, fixed_size>::~AutoBuffer()
 {
-    deallocate();
+        deallocate();
 }
 
@@ -154 +154 @@
 AutoBuffer<_Tp, fixed_size>::allocate(size_t _size)
 {
-    if (_size <= sz)
-    {
-        sz = _size;
-        return;
-    }
-    deallocate();
-    if (_size > fixed_size)
-    {
-        ptr = new _Tp[_size];
-        sz = _size;
-    }
+        if (_size <= sz)
+        {
+                sz = _size;
+                return;
+        }
+        deallocate();
+        if (_size > fixed_size)
+        {
+                ptr = new _Tp[_size];
+                sz = _size;
+        }
 }
 
@@ -170 +170 @@
 AutoBuffer<_Tp, fixed_size>::deallocate()
 {
-    if (ptr != buf)
-    {
-        delete[] ptr;
-        ptr = buf;
-        sz = fixed_size;
-    }
+        if (ptr != buf)
+        {
+                delete[] ptr;
+                ptr = buf;
+                sz = fixed_size;
+        }
 }
 
@@ -181 +181 @@
 AutoBuffer<_Tp, fixed_size>::resize(size_t _size)
 {
-    if (_size <= sz)
-    {
-        sz = _size;
-        return;
-    }
-    //size_t i, prevsize = sz, minsize = MIN(prevsize, _size);
-    size_t i, prevsize = sz, minsize = prevsize < _size ? prevsize : _size;
-    _Tp* prevptr = ptr;
-
-    ptr = _size > fixed_size ? new _Tp[_size] : buf;
-    sz = _size;
-
-    if (ptr != prevptr)
-        for (i = 0; i < minsize; i++)
-            ptr[i] = prevptr[i];
-    for (i = prevsize; i < _size; i++)
-        ptr[i] = _Tp();
-
-    if (prevptr != buf)
-        delete[] prevptr;
+        if (_size <= sz)
+        {
+                sz = _size;
+                return;
+        }
+        //size_t i, prevsize = sz, minsize = MIN(prevsize, _size);
+        size_t i, prevsize = sz, minsize = prevsize < _size ? prevsize : _size;
+        _Tp* prevptr = ptr;
+
+        ptr = _size > fixed_size ? new _Tp[_size] : buf;
+        sz = _size;
+
+        if (ptr != prevptr)
+                for (i = 0; i < minsize; i++)
+                        ptr[i] = prevptr[i];
+        for (i = prevsize; i < _size; i++)
+                ptr[i] = _Tp();
+
+        if (prevptr != buf)
+                delete[] prevptr;
 }
 
@@ -206 +206 @@
 AutoBuffer<_Tp, fixed_size>::size() const
 {
-    return sz;
+        return sz;
 }
 
@@ -212 +212 @@
 AutoBuffer<_Tp, fixed_size>::operator _Tp* ()
 {
-    return ptr;
+        return ptr;
 }
 
@@ -218 +218 @@
 AutoBuffer<_Tp, fixed_size>::operator const _Tp* () const
 {
-    return ptr;
+        return ptr;
 }
 
@@ -224 +224 @@
 {
 
-    int dot(const T*, const T*, int, T*) const
-    {
-        return 0;
-    }
-
-    int givens(T*, T*, int, T, T) const
-    {
-        return 0;
-    }
-
-    int givensx(T*, T*, int, T, T, T*, T*) const
-    {
-        return 0;
-    }
+        int dot(const T*, const T*, int, T*) const
+        {
+                return 0;
+        }
+
+        int givens(T*, T*, int, T, T) const
+        {
+                return 0;
+        }
+
+        int givensx(T*, T*, int, T, T, T*, T*) const
+        {
+                return 0;
+        }
 };
 
 template<typename _Tp> void
 JacobiSVDImpl_(_Tp* At, size_t astep, _Tp* _W, _Tp* Vt, size_t vstep,
-               int m, int n, int n1, double minval, _Tp eps)
-{
-    VBLAS<_Tp> vblas;
-    AutoBuffer<double> Wbuf(n);
-    double* W = Wbuf;
-    int i, j, k, iter, max_iter = std::max(m, 30);
-    _Tp c, s;
-    double sd;
-    astep /= sizeof (At[0]);
-    vstep /= sizeof (Vt[0]);
-
-    for (i = 0; i < n; i++)
-    {
-        for (k = 0, sd = 0; k < m; k++)
-        {
-            _Tp t = At[i * astep + k];
-            sd += (double) t*t;
-        }
-        W[i] = sd;
-
-        if (Vt)
-        {
-            for (k = 0; k < n; k++)
-                Vt[i * vstep + k] = 0;
-            Vt[i * vstep + i] = 1;
-        }
-    }
-
-    for (iter = 0; iter < max_iter; iter++)
-    {
-        bool changed = false;
-
-        for (i = 0; i < n - 1; i++)
-            for (j = i + 1; j < n; j++)
-            {
-                _Tp *Ai = At + i*astep, *Aj = At + j*astep;
-                double a = W[i], p = 0, b = W[j];
-
-                for (k = 0; k < m; k++)
-                    p += (double) Ai[k] * Aj[k];
-
-                if (std::abs(p) <= eps * std::sqrt((double) a * b))
-                    continue;
-
-                p *= 2;
-                double beta = a - b, gamma = hypot((double) p, beta);
-                if (beta < 0)
-                {
-                    double delta = (gamma - beta)*0.5;
-                    s = (_Tp) std::sqrt(delta / gamma);
-                    c = (_Tp) (p / (gamma * s * 2));
-                }
-                else
-                {
-                    c = (_Tp) std::sqrt((gamma + beta) / (gamma * 2));
-                    s = (_Tp) (p / (gamma * c * 2));
-                }
-
-                a = b = 0;
-                for (k = 0; k < m; k++)
-                {
-                    _Tp t0 = c * Ai[k] + s * Aj[k];
-                    _Tp t1 = -s * Ai[k] + c * Aj[k];
-                    Ai[k] = t0;
-                    Aj[k] = t1;
-
-                    a += (double) t0*t0;
-                    b += (double) t1*t1;
-                }
-                W[i] = a;
-                W[j] = b;
-
-                changed = true;
-
-                if (Vt)
-                {
-                    _Tp *Vi = Vt + i*vstep, *Vj = Vt + j*vstep;
-                    k = vblas.givens(Vi, Vj, n, c, s);
-
-                    for (; k < n; k++)
-                    {
-                        _Tp t0 = c * Vi[k] + s * Vj[k];
-                        _Tp t1 = -s * Vi[k] + c * Vj[k];
-                        Vi[k] = t0;
-                        Vj[k] = t1;
-                    }
-                }
-            }
-        if (!changed)
-            break;
-    }
-
-    for (i = 0; i < n; i++)
-    {
-        for (k = 0, sd = 0; k < m; k++)
-        {
-            _Tp t = At[i * astep + k];
-            sd += (double) t*t;
-        }
-        W[i] = std::sqrt(sd);
-    }
-
-    for (i = 0; i < n - 1; i++)
-    {
-        j = i;
-        for (k = i + 1; k < n; k++)
-        {
-            if (W[j] < W[k])
-                j = k;
-        }
-        if (i != j)
-        {
-            std::swap(W[i], W[j]);
-            if (Vt)
-            {
-                for (k = 0; k < m; k++)
-                    std::swap(At[i * astep + k], At[j * astep + k]);
-
-                for (k = 0; k < n; k++)
-                    std::swap(Vt[i * vstep + k], Vt[j * vstep + k]);
-            }
-        }
-    }
-
-    for (i = 0; i < n; i++)
-        _W[i] = (_Tp) W[i];
-
-    if (!Vt)
-        return;
+int m, int n, int n1, double minval, _Tp eps)
+{
+        VBLAS<_Tp> vblas;
+        AutoBuffer<double> Wbuf(n);
+        double* W = Wbuf;
+        int i, j, k, iter, max_iter = std::max(m, 30);
+        _Tp c, s;
+        double sd;
+        astep /= sizeof(At[0]);
+        vstep /= sizeof(Vt[0]);
+
+        for (i = 0; i < n; i++)
+        {
+                for (k = 0, sd = 0; k < m; k++)
+                {
+                        _Tp t = At[i * astep + k];
+                        sd += (double)t*t;
+                }
+                W[i] = sd;
+
+                if (Vt)
+                {
+                        for (k = 0; k < n; k++)
+                                Vt[i * vstep + k] = 0;
+                        Vt[i * vstep + i] = 1;
+                }
+        }
+
+        for (iter = 0; iter < max_iter; iter++)
+        {
+                bool changed = false;
+
+                for (i = 0; i < n - 1; i++)
+                        for (j = i + 1; j < n; j++)
+                        {
+                        _Tp *Ai = At + i*astep, *Aj = At + j*astep;
+                        double a = W[i], p = 0, b = W[j];
+
+                        for (k = 0; k < m; k++)
+                                p += (double)Ai[k] * Aj[k];
+
+                        if (std::abs(p) <= eps * std::sqrt((double)a * b))
+                                continue;
+
+                        p *= 2;
+                        double beta = a - b, gamma = hypot((double)p, beta);
+                        if (beta < 0)
+                        {
+                                double delta = (gamma - beta)*0.5;
+                                s = (_Tp)std::sqrt(delta / gamma);
+                                c = (_Tp)(p / (gamma * s * 2));
+                        }
+                        else
+                        {
+                                c = (_Tp)std::sqrt((gamma + beta) / (gamma * 2));
+                                s = (_Tp)(p / (gamma * c * 2));
+                        }
+
+                        a = b = 0;
+                        for (k = 0; k < m; k++)
+                        {
+                                _Tp t0 = c * Ai[k] + s * Aj[k];
+                                _Tp t1 = -s * Ai[k] + c * Aj[k];
+                                Ai[k] = t0;
+                                Aj[k] = t1;
+
+                                a += (double)t0*t0;
+                                b += (double)t1*t1;
+                        }
+                        W[i] = a;
+                        W[j] = b;
+
+                        changed = true;
+
+                        if (Vt)
+                        {
+                                _Tp *Vi = Vt + i*vstep, *Vj = Vt + j*vstep;
+                                k = vblas.givens(Vi, Vj, n, c, s);
+
+                                for (; k < n; k++)
+                                {
+                                        _Tp t0 = c * Vi[k] + s * Vj[k];
+                                        _Tp t1 = -s * Vi[k] + c * Vj[k];
+                                        Vi[k] = t0;
+                                        Vj[k] = t1;
+                                }
+                        }
+                        }
+                if (!changed)
+                        break;
+        }
+
+        for (i = 0; i < n; i++)
+        {
+                for (k = 0, sd = 0; k < m; k++)
+                {
+                        _Tp t = At[i * astep + k];
+                        sd += (double)t*t;
+                }
+                W[i] = std::sqrt(sd);
+        }
+
+        for (i = 0; i < n - 1; i++)
+        {
+                j = i;
+                for (k = i + 1; k < n; k++)
+                {
+                        if (W[j] < W[k])
+                                j = k;
+                }
+                if (i != j)
+                {
+                        std::swap(W[i], W[j]);
+                        if (Vt)
+                        {
+                                for (k = 0; k < m; k++)
+                                        std::swap(At[i * astep + k], At[j * astep + k]);
+
+                                for (k = 0; k < n; k++)
+                                        std::swap(Vt[i * vstep + k], Vt[j * vstep + k]);
+                        }
+                }
+        }
+
+        for (i = 0; i < n; i++)
+                _W[i] = (_Tp)W[i];
+
+        if (!Vt)
+                return;
 
         uint8_t rndstate = 0x02; //PRBS-7 from http://en.wikipedia.org/wiki/Pseudorandom_binary_sequence
         for (i = 0; i < n1; i++)
-    {
-        sd = i < n ? W[i] : 0;
-
-        while (sd <= minval)
-        {
-            // if we got a zero singular value, then in order to get the corresponding left singular vector
-            // we generate a random vector, project it to the previously computed left singular vectors,
-            // subtract the projection and normalize the difference.
-            const _Tp val0 = (_Tp) (1. / m);
-            for (k = 0; k < m; k++)
-            {
+        {
+                sd = i < n ? W[i] : 0;
+
+                while (sd <= minval)
+                {
+                        // if we got a zero singular value, then in order to get the corresponding left singular vector
+                        // we generate a random vector, project it to the previously computed left singular vectors,
+                        // subtract the projection and normalize the difference.
+                        const _Tp val0 = (_Tp)(1. / m);
+                        for (k = 0; k < m; k++)
+                        {
                                 int rndbit = (((rndstate >> 6) ^ (rndstate >> 5)) & 1);
                                 rndstate = ((rndstate << 1) | rndbit) & 0x7f;
                                 _Tp val = rndbit == 0 ? val0 : -val0;
-                At[i * astep + k] = val;
-            }
-            for (iter = 0; iter < 2; iter++)
-            {
-                for (j = 0; j < i; j++)
-                {
-                    sd = 0;
-                    for (k = 0; k < m; k++)
-                        sd += At[i * astep + k] * At[j * astep + k];
-                    _Tp asum = 0;
-                    for (k = 0; k < m; k++)
-                    {
-                        _Tp t = (_Tp) (At[i * astep + k] - sd * At[j * astep + k]);
-                        At[i * astep + k] = t;
-                        asum += std::abs(t);
-                    }
-                    asum = asum ? 1 / asum : 0;
-                    for (k = 0; k < m; k++)
-                        At[i * astep + k] *= asum;
-                }
-            }
-            sd = 0;
-            for (k = 0; k < m; k++)
-            {
-                _Tp t = At[i * astep + k];
-                sd += (double) t*t;
-            }
-            sd = std::sqrt(sd);
-        }
-
-        s = (_Tp) (1 / sd);
-        for (k = 0; k < m; k++)
-            At[i * astep + k] *= s;
-    }
+                                At[i * astep + k] = val;
+                        }
+                        for (iter = 0; iter < 2; iter++)
+                        {
+                                for (j = 0; j < i; j++)
+                                {
+                                        sd = 0;
+                                        for (k = 0; k < m; k++)
+                                                sd += At[i * astep + k] * At[j * astep + k];
+                                        _Tp asum = 0;
+                                        for (k = 0; k < m; k++)
+                                        {
+                                                _Tp t = (_Tp)(At[i * astep + k] - sd * At[j * astep + k]);
+                                                At[i * astep + k] = t;
+                                                asum += std::abs(t);
+                                        }
+                                        asum = asum ? 1 / asum : 0;
+                                        for (k = 0; k < m; k++)
+                                                At[i * astep + k] *= asum;
+                                }
+                        }
+                        sd = 0;
+                        for (k = 0; k < m; k++)
+                        {
+                                _Tp t = At[i * astep + k];
+                                sd += (double)t*t;
+                        }
+                        sd = std::sqrt(sd);
+                }
+
+                s = (_Tp)(1 / sd);
+                for (k = 0; k < m; k++)
+                        At[i * astep + k] *= s;
+        }
 }
 
 void Lapack::JacobiSVD(double* At, size_t astep, double* W, double* Vt, size_t vstep, int m, int n, int n1 = -1)
 {
-    JacobiSVDImpl_(At, astep, W, Vt, vstep, m, n, !Vt ? 0 : n1 < 0 ? n : n1, DBL_MIN, DBL_EPSILON * 10);
-}
+        JacobiSVDImpl_(At, astep, W, Vt, vstep, m, n, !Vt ? 0 : n1 < 0 ? n : n1, DBL_MIN, DBL_EPSILON * 10);
+}

cpp/frams/model/similarity/SVD/lapack.h

@@ -5 +5 @@
 {
 public:
-    static void JacobiSVD(double* At, size_t astep, double* W, double* Vt, size_t vstep, int m, int n, int n1);
+        static void JacobiSVD(double* At, size_t astep, double* W, double* Vt, size_t vstep, int m, int n, int n1);
 };
 

cpp/frams/model/similarity/SVD/matrix_tools.cpp

@@ -15 +15 @@
 double *Create(int nSize)
 {
-        double *matrix = (double *) malloc(nSize * sizeof (double));
-
-    for (int i = 0; i < nSize; i++)
-    {
-        matrix[i] = 0;
-    }
-
-    return matrix;
+        double *matrix = (double *)malloc(nSize * sizeof(double));
+
+        for (int i = 0; i < nSize; i++)
+        {
+                matrix[i] = 0;
+        }
+
+        return matrix;
 }
 
 double *Multiply(double *&a, double *&b, int nrow, int ncol, double ncol2, double *&toDel, int delSize)
 {
-    double *c = Create(nrow * ncol2);
-    int i = 0, j = 0, k = 0;
-
-    for (i = 0; i < nrow; i++)
-    {
-        for (j = 0; j < ncol2; j++)
-        {
-            for (k = 0; k < ncol; k++)
-                c[i * nrow + j] += a[i * nrow + k] * b[k * ncol + j];
-        }
-    }
-
-    if (delSize != 0)
-        free(toDel);
-    return c;
+        double *c = Create(nrow * ncol2);
+        int i = 0, j = 0, k = 0;
+
+        for (i = 0; i < nrow; i++)
+        {
+                for (j = 0; j < ncol2; j++)
+                {
+                        for (k = 0; k < ncol; k++)
+                                c[i * nrow + j] += a[i * nrow + k] * b[k * ncol + j];
+                }
+        }
+
+        if (delSize != 0)
+                free(toDel);
+        return c;
 }
 
 double *Power(double *&array, int nrow, int ncol, double pow, double *&toDel, int delSize)
 {
-    double *m_Power = Create(nrow * ncol);
-    if (pow == 2)
-    {
-        for (int i = 0; i < nrow; i++)
-        {
-            for (int j = 0; j < ncol; j++)
-            {
-                m_Power[i * nrow + j] = array[i * nrow + j] * array[i * nrow + j];
-            }
-
-        }
-    }
-    else
-    {
-        for (int i = 0; i < nrow; i++)
-        {
-            for (int j = 0; j < ncol; j++)
-            {
-                m_Power[i * nrow + j] = sqrt(array[i * nrow + j]);
-            }
-
-        }
-    }
-
-    if (delSize != 0)
-        free(toDel);
-
-    return m_Power;
+        double *m_Power = Create(nrow * ncol);
+        if (pow == 2)
+        {
+                for (int i = 0; i < nrow; i++)
+                {
+                        for (int j = 0; j < ncol; j++)
+                        {
+                                m_Power[i * nrow + j] = array[i * nrow + j] * array[i * nrow + j];
+                        }
+
+                }
+        }
+        else
+        {
+                for (int i = 0; i < nrow; i++)
+                {
+                        for (int j = 0; j < ncol; j++)
+                        {
+                                m_Power[i * nrow + j] = sqrt(array[i * nrow + j]);
+                        }
+
+                }
+        }
+
+        if (delSize != 0)
+                free(toDel);
+
+        return m_Power;
 }
 
 void Print(double *&mat, int nelems)
 {
-    for (int i = 0; i < nelems; i++)
-        printf("%6.2f ", mat[i]);
-    printf("\n");
+        for (int i = 0; i < nelems; i++)
+                printf("%6.2f ", mat[i]);
+        printf("\n");
 
 }
@@ -86 +86 @@
 double *Transpose(double *&A, int arow, int acol)
 {
-    double *result = Create(acol * arow);
-
-    for (int i = 0; i < acol; i++)
-        for (int j = 0; j < arow; j++)
-        {
-            result[i * arow + j] = A[j * acol + i];
-        }
-
-    return result;
+        double *result = Create(acol * arow);
+
+        for (int i = 0; i < acol; i++)
+                for (int j = 0; j < arow; j++)
+                {
+                result[i * arow + j] = A[j * acol + i];
+                }
+
+        return result;
 
 }
 
 /** Computes the SVD of the nSize x nSize distance matrix
-        @param vdEigenvalues [OUT] Vector of doubles. On return holds the eigenvalues of the 
-        decomposed distance matrix (or rather, to be strict, of the matrix of scalar products 
-        created from the matrix of distances). The vector is assumed to be empty before the function call and
-        all variance percentages are pushed at the end of it.
-        @param nSize size of the matrix of distances.
-        @param pDistances [IN] matrix of distances between parts.
-        @param Coordinates [OUT] array of three dimensional coordinates obtained from SVD of pDistances matrix.
- */
+                @param vdEigenvalues [OUT] Vector of doubles. On return holds the eigenvalues of the
+                decomposed distance matrix (or rather, to be strict, of the matrix of scalar products
+                created from the matrix of distances). The vector is assumed to be empty before the function call and
+                all variance percentages are pushed at the end of it.
+                @param nSize size of the matrix of distances.
+                @param pDistances [IN] matrix of distances between parts.
+                @param Coordinates [OUT] array of three dimensional coordinates obtained from SVD of pDistances matrix.
+                */
 void MatrixTools::SVD(std::vector<double> &vdEigenvalues, int nSize, double *pDistances, Pt3D *&Coordinates)
 {
-    // compute squares of elements of this array
-    // compute the matrix B that is the parameter of SVD
-    double *B = Create(nSize * nSize);
-    {
-        // use additional scope to delete temporary matrices
-        double *Ones, *Eye, *Z, *D;
-
-        D = Create(nSize * nSize);
-        D = Power(pDistances, nSize, nSize, 2.0, D, nSize);
-
-        Ones = Create(nSize * nSize);
-        for (int i = 0; i < nSize; i++)
-            for (int j = 0; j < nSize; j++)
-            {
-                Ones[i * nSize + j] = 1;
-            }
-
-        Eye = Create(nSize * nSize);
-        for (int i = 0; i < nSize; i++)
-        {
-            for (int j = 0; j < nSize; j++)
-            {
-                if (i == j)
-                {
-                    Eye[i * nSize + j] = 1;
-                }
-                else
-                {
-                    Eye[i * nSize + j] = 0;
-                }
-            }
-        }
-
-        Z = Create(nSize * nSize);
-        for (int i = 0; i < nSize; i++)
-        {
-            for (int j = 0; j < nSize; j++)
-            {
-                Z[i * nSize + j] = 1.0 / ((double) nSize) * Ones[i * nSize + j];
-            }
-        }
-
-        for (int i = 0; i < nSize; i++)
-        {
-            for (int j = 0; j < nSize; j++)
-            {
-                Z[i * nSize + j] = Eye[i * nSize + j] - Z[i * nSize + j];
-            }
-        }
-
-        for (int i = 0; i < nSize; i++)
-        {
-            for (int j = 0; j < nSize; j++)
-            {
-                B[i * nSize + j] = Z[i * nSize + j] * -0.5;
-            }
-        }
-
-        B = Multiply(B, D, nSize, nSize, nSize, B, nSize);
-        B = Multiply(B, Z, nSize, nSize, nSize, B, nSize);
-
-        free(Ones);
-        free(Eye);
-        free(Z);
-        free(D);
-    }
-
-    double *Eigenvalues = Create(nSize);
-    double *S = Create(nSize * nSize);
-
-    // call SVD function
-    double *Vt = Create(nSize * nSize);
-    size_t astep = nSize * sizeof (double);
-    Lapack::JacobiSVD(B, astep, Eigenvalues, Vt, astep, nSize, nSize, nSize);
-
-    double *W = Transpose(Vt, nSize, nSize);
-
-    free(B);
-    free(Vt);
-
-    for (int i = 0; i < nSize; i++)
-        for (int j = 0; j < nSize; j++)
-        {
-            if (i == j)
-                S[i * nSize + j] = Eigenvalues[i];
-            else
-                S[i * nSize + j] = 0;
-        }
-
-    // compute coordinates of points 
-    double *sqS, *dCoordinates;
-    sqS = Power(S, nSize, nSize, 0.5, S, nSize);
-    dCoordinates = Multiply(W, sqS, nSize, nSize, nSize, W, nSize);
-    free(sqS);
-
-    for (int i = 0; i < nSize; i++)
-    {
-        // set coordinate from the SVD solution
-        Coordinates[ i ].x = dCoordinates[i * nSize];
-        Coordinates[ i ].y = dCoordinates[i * nSize + 1 ];
-        if (nSize > 2)
-            Coordinates[ i ].z = dCoordinates[i * nSize + 2 ];
-        else
-            Coordinates[ i ].z = 0;
-    }
-
-    // store the eigenvalues in the output vector
-    for (int i = 0; i < nSize; i++)
-    {
-        double dElement = Eigenvalues[i];
-        vdEigenvalues.push_back(dElement);
-    }
-
-    free(Eigenvalues);
-    free(dCoordinates);
-}
+        // compute squares of elements of this array
+        // compute the matrix B that is the parameter of SVD
+        double *B = Create(nSize * nSize);
+        {
+                // use additional scope to delete temporary matrices
+                double *Ones, *Eye, *Z, *D;
+
+                D = Create(nSize * nSize);
+                D = Power(pDistances, nSize, nSize, 2.0, D, nSize);
+
+                Ones = Create(nSize * nSize);
+                for (int i = 0; i < nSize; i++)
+                        for (int j = 0; j < nSize; j++)
+                        {
+                        Ones[i * nSize + j] = 1;
+                        }
+
+                Eye = Create(nSize * nSize);
+                for (int i = 0; i < nSize; i++)
+                {
+                        for (int j = 0; j < nSize; j++)
+                        {
+                                if (i == j)
+                                {
+                                        Eye[i * nSize + j] = 1;
+                                }
+                                else
+                                {
+                                        Eye[i * nSize + j] = 0;
+                                }
+                        }
+                }
+
+                Z = Create(nSize * nSize);
+                for (int i = 0; i < nSize; i++)
+                {
+                        for (int j = 0; j < nSize; j++)
+                        {
+                                Z[i * nSize + j] = 1.0 / ((double)nSize) * Ones[i * nSize + j];
+                        }
+                }
+
+                for (int i = 0; i < nSize; i++)
+                {
+                        for (int j = 0; j < nSize; j++)
+                        {
+                                Z[i * nSize + j] = Eye[i * nSize + j] - Z[i * nSize + j];
+                        }
+                }
+
+                for (int i = 0; i < nSize; i++)
+                {
+                        for (int j = 0; j < nSize; j++)
+                        {
+                                B[i * nSize + j] = Z[i * nSize + j] * -0.5;
+                        }
+                }
+
+                B = Multiply(B, D, nSize, nSize, nSize, B, nSize);
+                B = Multiply(B, Z, nSize, nSize, nSize, B, nSize);
+
+                free(Ones);
+                free(Eye);
+                free(Z);
+                free(D);
+        }
+
+        double *Eigenvalues = Create(nSize);
+        double *S = Create(nSize * nSize);
+
+        // call SVD function
+        double *Vt = Create(nSize * nSize);
+        size_t astep = nSize * sizeof(double);
+        Lapack::JacobiSVD(B, astep, Eigenvalues, Vt, astep, nSize, nSize, nSize);
+
+        double *W = Transpose(Vt, nSize, nSize);
+
+        free(B);
+        free(Vt);
+
+        for (int i = 0; i < nSize; i++)
+                for (int j = 0; j < nSize; j++)
+                {
+                if (i == j)
+                        S[i * nSize + j] = Eigenvalues[i];
+                else
+                        S[i * nSize + j] = 0;
+                }
+
+        // compute coordinates of points 
+        double *sqS, *dCoordinates;
+        sqS = Power(S, nSize, nSize, 0.5, S, nSize);
+        dCoordinates = Multiply(W, sqS, nSize, nSize, nSize, W, nSize);
+        free(sqS);
+
+        for (int i = 0; i < nSize; i++)
+        {
+                // set coordinate from the SVD solution
+                Coordinates[i].x = dCoordinates[i * nSize];
+                Coordinates[i].y = dCoordinates[i * nSize + 1];
+                if (nSize > 2)
+                        Coordinates[i].z = dCoordinates[i * nSize + 2];
+                else
+                        Coordinates[i].z = 0;
+        }
+
+        // store the eigenvalues in the output vector
+        for (int i = 0; i < nSize; i++)
+        {
+                double dElement = Eigenvalues[i];
+                vdEigenvalues.push_back(dElement);
+        }
+
+        free(Eigenvalues);
+        free(dCoordinates);
+}

cpp/frams/model/similarity/SVD/matrix_tools.h

@@ -12 +12 @@
 {
 public:
-    static void SVD(std::vector<double> &vdEigenvalues, int nSize, double *pDistances, Pt3D *&Coordinates);
+        static void SVD(std::vector<double> &vdEigenvalues, int nSize, double *pDistances, Pt3D *&Coordinates);
 };