Line data    Source code

       1             : #include <cfloat>
       2             : #include <iostream>
       3             : #include <string>
       4             : #include <vector>
       5             : 
       6             : extern "C"
       7             : {
       8             : #include "qrfactorization.h"
       9             : }
      10             : #include "exceptions.hpp"
      11             : #include "matrix.hpp"
      12             : #include "qrfactorization.hpp"
      13             : #include "vectorf.hpp"
      14             : 
      15          38 : int qr(const int m, const int n, matrix<double> &W, vectorf<double> &b,
      16             :        vectorf<double> &x)
      17             : {
      18          38 :   int ier{0};
      19          76 :   vectorf<int> ipiv(n);
      20          40 :   vectorf<double> tau(n);
      21             : 
      22             :   // Factor
      23          38 :   factor(m, n, W, tau, ipiv);
      24             :   // Solve
      25          36 :   solve(m, n, W, tau, b, x);
      26             :   // Multiply solution x by permutation and store it in x
      27          36 :   permute(n, x, ipiv, tau);
      28             : 
      29          72 :   return ier;
      30             : }
      31             : 
      32          40 : int factor(const int m, const int n, matrix<double> &W, vectorf<double> &tau,
      33             :            vectorf<int> &ipiv)
      34             : {
      35          40 :   int ier{0};
      36          40 :   double safmin{DBL_MIN};
      37          44 :   vectorf<double> wrk(m);
      38          40 :   qrf(m, n, &W(1, 1), m, &ipiv(1), &tau(1), &wrk(1), safmin, &ier);
      39          40 :   if (ier != 0)
      40             :   {
      41             :     std::string message =
      42          12 :         "This QRFactorization error with ier= " + std::to_string(ier) + " \n";
      43           4 :     throw sayMessage(message);
      44             :   }
      45          72 :   return ier;
      46             : }
      47             : 
      48          38 : int solve(const int m, const int n, matrix<double> &W, vectorf<double> &tau,
      49             :           vectorf<double> &b, vectorf<double> &x)
      50             : {
      51          38 :   int ier{0};
      52          38 :   qrs(m, n, &W(1, 1), m, &tau(1), &b(1), &x(1), &ier);
      53          38 :   if (ier != 0)
      54             :   {
      55             :     std::string message =
      56           6 :         "QRSolver error with ier= " + std::to_string(ier) + " \n";
      57           2 :     throw sayMessage(message);
      58             :   }
      59          36 :   return ier;
      60             : }
      61             : 
      62          37 : int permute(const int n, vectorf<double> &x, vectorf<int> &ipiv,
      63             :             vectorf<double> &wrk)
      64             : {
      65          37 :   int ier{0};
      66         139 :   for (int i{1}; i <= n; i++)
      67             :   {
      68         102 :     wrk(i) = x(ipiv(i));
      69             :   }
      70         139 :   for (int i{1}; i <= n; i++)
      71             :   {
      72         102 :     x(i) = wrk(i);
      73             :   }
      74          37 :   return ier;
      75             : }