GopherProxy

	tjacobi.cpp - numeric - C++ library with numerical algorithms
	git clone git://src.adamsgaard.dk/numeric
	Log
	Files
	Refs
	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 = cc app - 2.0f * s * c * apq + ss aqq;
	38 aqq1 = ss app + 2.0f * s * c * apq + cc 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 }