Line data    Source code

       1             : #include <assert.h>
       2             : #include <float.h>
       3             : #include <stdio.h>
       4             : #include <stdlib.h>
       5             : #include <tgmath.h>
       6             : 
       7             : #include "qrfactorization.h"
       8             : 
       9             : #define map(i, j) ((j)-1) * lda + ((i)-1)
      10             : #define a(i, j) a[ map(i, j) ]
      11             : #define map1(i) (i) - 1
      12             : #define x(i) x[ map1(i) ]
      13             : #define y(i) y[ map1(i) ]
      14             : #define ipiv(i) ipiv[ map1(i) ]
      15             : #define tc(i) tc[ map1(i) ]
      16             : 
      17             : //  error code 1 if less than 10 coefficients to estimate Rc
      18             : 
      19           0 : int roc(int num_tc, double *tc, double scale, int kstart, double *rc,
      20             :         double *slope, double *intercept)
      21             : {
      22           0 :   if (num_tc - kstart < 11)
      23             :   {
      24           0 :     printf(
      25             :         "Message from ROC: Need at least 10 Taylor coefficients to estimate "
      26             :         "Rc\n");
      27           0 :     return 1;
      28             :   }
      29             : 
      30           0 :   int m = num_tc - kstart + 1;
      31           0 :   int lda = m, n = 2, min_m_n = 0;
      32           0 :   int res = 0, ier = 0;
      33           0 :   double *a = NULL;
      34           0 :   int *ipiv = NULL;
      35           0 :   double *tau = NULL;
      36           0 :   double *wrk = NULL;
      37           0 :   double *x = NULL;
      38           0 :   double *y = NULL;
      39           0 :   double safmin = DBL_MIN;
      40             : 
      41           0 :   min_m_n = (m < n) ? m : n;
      42             : 
      43           0 :   a = (double *)calloc((size_t)(lda * n), (size_t)sizeof(double));
      44           0 :   if (a == NULL)
      45             :   {
      46           0 :     printf("Failed to allocate matrix A memory!\n");
      47           0 :     assert(a);
      48             :   }
      49             : 
      50           0 :   tau = (double *)calloc((size_t)min_m_n, (size_t)sizeof(double));
      51           0 :   if (tau == NULL)
      52             :   {
      53           0 :     printf("Failed to allocate tau memory!\n");
      54           0 :     assert(tau);
      55             :   }
      56             : 
      57           0 :   wrk = (double *)calloc((size_t)n, (size_t)sizeof(double));
      58           0 :   if (wrk == NULL)
      59             :   {
      60           0 :     printf("Failed to allocate wrk memory!\n");
      61           0 :     assert(wrk);
      62             :   }
      63             : 
      64           0 :   ipiv = (int *)calloc((size_t)n, (size_t)sizeof(int));
      65           0 :   if (ipiv == NULL)
      66             :   {
      67           0 :     printf("Failed to allocate ipiv memory!\n");
      68           0 :     assert(ipiv);
      69             :   }
      70             : 
      71           0 :   x = (double *)calloc((size_t)n, (size_t)sizeof(double));
      72           0 :   if (x == NULL)
      73             :   {
      74           0 :     printf("Failed to allocate space for result x!\n");
      75           0 :     assert(x);
      76             :   }
      77             : 
      78           0 :   y = (double *)calloc((size_t)m, (size_t)sizeof(double));
      79           0 :   if (y == NULL)
      80             :   {
      81           0 :     printf("Failed to allocate space for result y!\n");
      82           0 :     assert(y);
      83             :   }
      84             : 
      85             :   // Construct the least squares linear system
      86             :   //
      87             :   // We want to fit the tail of the Taylor series, ignore the first kstart terms
      88             :   // of the TCs.
      89           0 :   for (int k = kstart; k <= num_tc; k++)
      90             :   {
      91           0 :     int row = k - kstart + 1;
      92           0 :     a(row, 1) = 1.0;
      93           0 :     a(row, 2) = (double)k;
      94           0 :     y(row) = log10(fabs(tc(k)));
      95             :   }
      96             : 
      97             :   // Find the QRFactorization of the least squares linear system
      98           0 :   res = qrf(m, n, a, lda, ipiv, tau, wrk, safmin, &ier);
      99           0 :   if (res != 0 || ier != 0)
     100             :   {
     101           0 :     printf("Error QR factoring A with res[%d]!=0 or ier[%d]!=0\n", res, ier);
     102             :   }
     103             : 
     104             :   // Find the least squares solution of the linear system via the
     105             :   // QRFactorization
     106             :   //
     107             :   // Solve on the first (two) n rows of R as in
     108             :   // transpose(Q)*b=transpose(Q)*A*E*E*x=transpose(Q)*Q*R*E*x=R*E*x
     109           0 :   res = qrs(m, n, a, lda, tau, y, x, &ier);
     110           0 :   if (res != 0 || ier != 0)
     111             :   {
     112           0 :     printf(
     113             :         "Error finding the Least Squares Solution with res[%d]!=0 or "
     114             :         "ier[%d]!=0\n",
     115             :         res, ier);
     116             :   }
     117             : 
     118             :   // If Taylor coefficients TC(i) are unscaled, then the scaled TCs are
     119             :   // T(i) = (1/scale)^i TC(i).
     120             :   //
     121             :   // Multiply Rc by scale to get Rc when computing with scaled TCs.
     122             :   //
     123             :   // See this from ratio test as Rc is 1/L where
     124             :   //
     125             :   //   L = (1/scale)*limit of abs(TC(i+1)/TC(i)) = limit of abs(T(i+1)/T(i)).
     126             :   //
     127             :   // To justify the formula for Rc, consider the problem of fitting m and b
     128             :   //
     129             :   // log10(abs(T(i))) = y(i) = b + m*i = [1 i](b,m)
     130             :   //
     131             :   // by linear least squares. Then
     132             :   //
     133             :   // Log10(L) = Log10(abs(T(i+1))/abs(T(i))) =
     134             :   // Log10(abs(T(i+1)))-Log10(abs(T(i)))
     135             :   //          = m*(i+1)+b - (m*i+b) = m
     136             :   //
     137             :   // It follows that
     138             :   //
     139             :   // rc = 1/L = scale/pow(10, m).
     140           0 :   *rc = scale / pow(10, x(ipiv(2)));
     141           0 :   *intercept = x(ipiv(1));
     142           0 :   *slope = x(ipiv(2));
     143             : 
     144           0 :   free(y);
     145           0 :   free(x);
     146           0 :   free(ipiv);
     147           0 :   free(wrk);
     148           0 :   free(tau);
     149           0 :   free(a);
     150             : 
     151           0 :   return 0;
     152             : }