GopherProxy

	tqrfunc.cpp - numeric - C++ library with numerical algorithms
	git clone git://src.adamsgaard.dk/numeric
	Log
	Files
	Refs
	LICENSE
	---
	tqrfunc.cpp (2033B)
	---
	1 #include <iostream>
	2 #include <armadillo>
	3 #include "header.h"
	4 #include "qrfunc.h"
	5
	6 // QR decomposition constructor
	7 QR::QR(arma::mat &A)
	8 : n(A.n_cols),
	9 A(A),
	10 Q(A)
	11 {
	12 // Initialize output structures
	13 R = arma::zeros<arma::mat> (n,n);
	14
	15 // Perform QR decomposition straight away
	16 decomp();
	17 }
	18
	19 // Class deconstructor (equivalent to compiler destructor)
	20 QR::~QR() { };
	21
	22 // Return system size
	23 Lengthtype QR::size()
	24 {
	25 return n;
	26 }
	27
	28 // QR decomposition function of Armadillo matrix.
	29 // Returns right matrix R, and modifies A into Q.
	30 // Uses Gram-Schmidt orthogonalization
	31 void QR::decomp()
	32 {
	33 Floattype r, s;
	34 Lengthtype j;
	35 for (Lengthtype i=0; i<n; ++i) {
	36 r = dot(Q.col(i), Q.col(i));
	37 R(i,i) = sqrt(r);
	38 Q.col(i) /= sqrt(r); // Normalization
	39 for (j=i+1; j<n; ++j) {
	40 s = dot(Q.col(i), Q.col(j));
	41 Q.col(j) -= s*Q.col(i); // Orthogonalization
	42 R(i,j) = s;
	43 }
	44 }
	45 }
	46
	47 // Solve the square system of linear equations with back
	48 // substitution. T is an upper triangular matrix.
	49 arma::vec QR::backsub(arma::vec &b)
	50 {
	51 Floattype tmpsum;
	52 arma::vec x = Q.t() * b;
	53 for (Lengthtype i=n-1; i>=0; --i) {
	54 tmpsum = 0.0f;
	55 for (Lengthtype k=i+1; k<n; ++k)
	56 tmpsum += R(i,k) * x(k);
	57
	58 x(i) = 1.0f/R(i,i) * (x(i) - tmpsum);
	59 }
	60 return x;
	61 }
	62
	63 // Calculate the (absolute value of the) determinant of
	64 // matrix A given the Q and R components.
	65 // det(A) = det(Q) * det(R), \|det(Q) = 1\|
	66 // => \|det(A)\| = \|det(R)\|
	67 Floattype QR::det()
	68 {
	69 Floattype determinant = 1.0f;
	70 for (Lengthtype i=0; i<n; ++i)
	71 determinant *= R(i,i);
	72 return fabs(determinant);
	73 }
	74
	75 // Calculate the inverse of matrix A given the Q and R
	76 // components.
	77 arma::mat QR::inverse()
	78 {
	79 arma::mat inv = arma::zeros<arma::mat> (n, n);
	80 // In vector z, all elements are equal to 0, except z(i) = 1
	81 arma::vec z = arma::zeros<arma::mat> (n);
	82
	83 for (Lengthtype i=0; i<n; ++i) {
	84 z(i) = 1.0f; // Element i changed to 1
	85 inv.col(i) = backsub(z);
	86 z(i) = 0.0f; // Element i changed back to 0
	87 }
	88
	89 return inv;
	90 }