Implement jacobi rotation method

include/jacobi.hpp

@@ -0,0 +1,48 @@
+/** @file jacobi.hpp
+ *  @brief Function prototypes for the jacobi rotation algorithm.
+ *
+ *  @author Cory Alexander Balaton (coryab)
+ *  @author Janita Ovidie Sandtrøen Willumsen (janitaws)
+ *  @bug No known bugs
+ */
+#ifndef __JACOBI__
+#define __JACOBI__
+
+#include <armadillo>
+
+/** @brief Computes a single rotation.
+ *
+ *  Description
+ *
+ *  @param A Matrix A<sup> (m) </sup>
+ *  @param R The rotation matrix R<sup> (m) </sup>
+ *  @param k Index of the row with the element of largest absolute value
+ *  @param l Index of the column with the element of largest absolute value
+ *
+ *  @return Void
+ */
+void jacobi_rotate(arma::mat& A, arma::mat& R, int k, int l);
+
+/** @ brief Solves the eigenvalue problem using the jacobi rotation method.
+ *
+ *  Description
+ *
+ *  @param A The initial matrix to be solved
+ *  @param eps Tolerance
+ *  @param eigenvalues A vector that will contain the computed eigenvalues
+ *  @param eigenvectors A Matrix that will contain the computed eigenvectors
+ *  @param maxiter The max number of iterations allowed
+ *  @param iterations To keep track of how many iterations the algorithm used
+ *  @param converged Tells the user if the algorithm has converged
+ *
+ *  @return Void
+ * */
+void jacobi_eigensolver(const arma::mat& A, 
+                        double eps, 
+                        arma::vec& eigenvalues, 
+                        arma::mat& eigenvectors, 
+                        const int maxiter, 
+                        int& iterations, 
+                        bool& converged);
+
+#endif

src/jacobi.cpp

@@ -0,0 +1,84 @@
+/** @file jacobi.cpp
+ *  @brief Implementation of the jacobi rotation method
+ *
+ *  @author Cory Alexander Balaton (coryab)
+ *  @author Janita Ovidie Sandtrøen Willumsen (janitaws)
+ *  @bug No known bugs
+ */
+#include <cmath>
+
+#include "utils.hpp"
+#include "jacobi.hpp"
+#include "matrix.hpp"
+
+void jacobi_rotate(arma::mat& A, arma::mat& R, int k, int l)
+{
+    double tau, t, c, s, sr, a_kk, a_ik, r_ik;
+
+    tau = (A(l,l) - A(k,k))/(2.*A(k,l));
+    
+    sr = sqrt(1 + tau*tau);
+    t = - tau + (tau > 0. ? sr : -sr);
+
+    c = 1. / std::sqrt(1. + t*t);
+
+    s = c*t;
+
+    // Calculating the diagonal elements at index k an l
+    a_kk = A(k,k);
+    A(k,k) = a_kk*c*c - 2.*A(k,l)*c*s + A(l,l)*s*s;
+    A(l,l) = A(l,l)*c*c + 2.*A(k,l)*c*s + a_kk*s*s;
+    A(k,l) = A(l,k) = 0.; // These are 0 by definition
+
+    // Calculating the other values on rows/cols k and l of A
+    for (int i=0; i < A.n_rows; i++) {
+        if (i == k || i == l) {
+            continue;
+        }
+        a_ik = A(i,k);
+        A(i,k) = a_ik*c - A(i,l)*s;
+        A(k,i) = A(i,k);
+        A(i,l) = A(i,l)*c + a_ik*s;
+        A(l,i) = A(i,l);
+    }
+
+    // Calculate the values for R
+    for (int i=0; i < R.n_rows; i++) {
+        r_ik = R(i,k);
+        R(i,k) = r_ik*c - R(i,l)*s;
+        R(i,l) = R(i,l)*c + r_ik*s;
+    }
+}
+
+void jacobi_eigensolver(const arma::mat& A, 
+                        double eps, 
+                        arma::vec& eigenvalues, 
+                        arma::mat& eigenvectors, 
+                        const int maxiter, 
+                        int& iterations, 
+                        bool& converged)
+{
+    int k, l;
+    double max_offdiag;
+
+    // Initialize matrices
+    arma::mat A_m = A;
+    arma::mat R = arma::mat(A.n_rows, A.n_cols, arma::fill::eye);
+    
+    iterations = 0;
+
+    // Use a do while so that max_offdiag gets a value before doing a comparison
+    do {
+    max_offdiag = max_offdiag_symmetric(A_m, k, l);
+        jacobi_rotate(A_m, R, k, l);
+
+    } while (max_offdiag >= eps && ++iterations < maxiter);
+
+    // The diagonal elements of A_m are the eigenvalues
+    eigenvalues = arma::diagvec(A_m);
+    eigenvectors = R;
+
+    converged = max_offdiag < eps;
+    A_m.print();
+    R.print();
+}