tjacobi.cpp - numeric - C++ library with numerical algorithms
 (HTM) git clone git://src.adamsgaard.dk/numeric
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) LICENSE
       ---
       tjacobi.cpp (2253B)
       ---
            1 #include <iostream>
            2 #include <cmath>
            3 #include <armadillo>
            4 #include "header.h"
            5 #include "jacobi.h"
            6 
            7 // Constructor: Perform rotation straight away
            8 Jacobi::Jacobi(const arma::Mat<Floattype> &A)
            9                 : n(A.n_rows), At(A), e(n), V(n,n)
           10 {
           11   // Initialize diagonal vector to be filled with eigenvalues
           12   e = A.diag();
           13 
           14   // Set main diagonal values to 1, off-diagonal values to 0
           15   V.eye();
           16 
           17   Lengthtype p, q, i; // Iterator vars
           18   int changed;
           19   Lengthtype rotations = 0; // Number of rotations performed
           20 
           21   // Cell value variables, used as local storage for mat/vec vals.
           22   Floattype app, aqq, apq, phi, c, s, app1, aqq1;
           23   Floattype aip, api, aiq, aqi, vip, viq;
           24 
           25   do {
           26     changed = 0;
           27     for (p=0; p<n; ++p) { // Row iteration
           28       for (q=p+1; q<n; ++q) { // Cols right of diagonal
           29         
           30         // Initialize cell-relevant data
           31         app = e(p);
           32         aqq = e(q);
           33         apq = At(p,q);
           34         phi = 0.5f * atan2(2.0f * apq, aqq-app);
           35         c = cos(phi);
           36         s = sin(phi);
           37         app1 = c*c * app - 2.0f * s * c * apq + s*s * aqq;
           38         aqq1 = s*s * app + 2.0f * s * c * apq + c*c * aqq;
           39 
           40         if (app1 != app || aqq1 != aqq) {
           41           changed = 1;
           42           ++rotations;
           43           e(p)    = app1;
           44           e(q)    = aqq1;
           45           At(p,q) = 0.0f;
           46 
           47           for (i = 0; i<p; ++i) {
           48             aip = At(i,p);
           49             aiq = At(i,q);
           50             At(i,p) = c * aip - s * aiq;
           51             At(i,q) = c * aiq + s * aip;
           52           }
           53 
           54           for (i=p+1; i<q; ++i) {
           55             api = At(p,i);
           56             aiq = At(i,q);
           57             At(p,i) = c * api - s * aiq;
           58             At(i,q) = c * aiq + s * api;
           59           }
           60 
           61           for (i=q+1; i<n; ++i) {
           62             api = At(p,i);
           63             aqi = At(q,i);
           64             At(p,i) = c * api - s * aqi;
           65             At(q,i) = c * aqi + s * api;
           66           }
           67 
           68           for (i=0; i<n; ++i) {
           69             vip = V(i,p);
           70             viq = V(i,q);
           71             V(i,p) = c * vip - s * viq;
           72             V(i,q) = c * viq + s * vip;
           73           }
           74         
           75         } // if end
           76       } // q loop end
           77     } // p loop end
           78   } while (changed != 0); // do-while end
           79 
           80   // Make transformed matrix symmetric 
           81   At = arma::symmatu(At);
           82 
           83   if (verbose == true)
           84     std::cout << "\nPerformed " << rotations << " Jacobi rotations.\n";
           85 }
           86 
           87 // Return transformed matrix
           88 arma::Mat<Floattype> Jacobi::trans()
           89 {
           90   return At;
           91 }
           92 
           93 // Return matrix of eigenvectors
           94 arma::Mat<Floattype> Jacobi::eigenvectors()
           95 {
           96   return V;
           97 }
           98 
           99 // Return vector of eigenvalues
          100 arma::Col<Floattype> Jacobi::eigenvalues()
          101 {
          102   return e;
          103 }