tqrfunc.cpp - numeric - C++ library with numerical algorithms
 (HTM) git clone git://src.adamsgaard.dk/numeric
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) 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 }