Re: Octave/C++ matrix Inv() comparison

[ionone]

i thought i had sent the Matlab/Octave's code to the mailing list, i screwed
up somewhere.

here is the Octave's code i used

x = [1/1024:1/1024:1];

p = 20
h = 128;
w = 2*h;
ov = 1;

if (size(x,2) == 1)
  x = x';  % Convert X from column to row
end

npts = length(x);

nhops = floor(npts/h);

% Pad x with zeros so that we can extract complete w-length windows
% from it
x = [zeros(1,(w-h)/2),x,zeros(1,(w-h/2))];

a = zeros(nhops, p+1);
g = zeros(nhops, 1);
if ov == 0
  e = zeros(1, npts);
else
  e = zeros(1, (nhops-1)*h+w);
end

pre = [1 -0.9];
x = filter(pre,1,x);

nhops = 3
hop = 3

  % Extract segment of signal
  xx = x((hop - 1)*h + [1:w]);
  % Apply hanning window
  wxx = xx .* hanning(w)';
  % Form autocorrelation (calculates *way* too many points)
  rxx = xcorr(wxx);
  % extract just the points we need (middle p+1 points)
  rxx = rxx(w+[0:p]);% rxx ok
  % Setup the normal equations
  R = toeplitz(rxx(1:p));
inv(R)
//******************************************************************************
and here is the corresponding inverse() Lapack C++ code 

#include <cstdio>
#include "stdafx.h"

extern "C" {
        // LU decomoposition of a general matrix
        void __cdecl dgetrf_(int* M, int *N, double* A, int* lda, int* IPIV,
int* INFO);

        // generate inverse of a matrix given its LU decomposition
        void __cdecl dgetri_(int* N, double* A, int* lda, int* IPIV, double*
WORK, int* lwork, int* INFO);


     
        int __cdecl  dsytrf_(char *uplo, int *n, double *a, int *lda, int
*ipiv, double *work, int *lwork, int *info);
        int __cdecl  dsytri_(char *uplo, int *n, double *a, int *lda, int
*ipiv, double *work, int *info);

    }
    void inverse(double* A, int N)
    {
        int *IPIV = new int[N+1];
        int LWORK = N*N;
        double *WORK = new double[LWORK];
        int INFO;

        dgetrf_(&N,&N,A,&N,IPIV,&INFO);
        dgetri_(&N,A,&N,IPIV,WORK,&LWORK,&INFO);

        delete IPIV;
        delete WORK;
    }
    void inverse2 (double * matrix, int matrix_rank)
    {
        int info = 0;
        int lwork = matrix_rank;
        int * ivpv = new   int [matrix_rank];
        double * work = new   double [lwork];

        dsytrf_("U", & matrix_rank, matrix, & matrix_rank, ivpv, work, &
lwork, & info);
        dsytri_("U", & matrix_rank, matrix, & matrix_rank, ivpv, work, &
info);

        for (int i = 0; i <matrix_rank; i++)
        {
            for (int j = i +1; j &lt;matrix_rank; j++)
                matrix [i * matrix_rank + j] = matrix [j * matrix_rank + i];
        }

        delete [] ivpv;
        delete [] work;
    }


//****************************************************************************************
and here is the C++ code, simply the same Octave's code but in C++


int nsmp = 1024;
    double* d = new double[nsmp+1];
    for (int i = 1; i &lt;= nsmp; i++)
        d[i] = i/double(nsmp);
    float alpha = -0.2;

double** a;
    double* g;
    double* e;
    double* x;

     int p = 20;


    int  h = 128;

    int w = 2*h;

    int ov = 1;

    int sizeax;
    int sizeay;
    int sizee;

    int npts = nsmp;

    int nhops = floor(npts/double(h));

  

    //Pad x with zeros so that we can extract complete w-length windows
    // from it
    double* x2 = new double[nsmp+1];

    for (int i = 0; i &lt;= nsmp; i++)
        x2[i] = xx[i];

    x = new double[(w-h)+nsmp+1];



    //x = [zeros(1,(w-h)/2),x,zeros(1,(w-h/2))];
    for (int i = 1; i &lt;= (w-h)/2; i++)
        x[i] = 0;
    for (int i = (w-h)/2+1; i &lt;= (w-h)/2+nsmp; i++)
    {
        x[i] = x2[i-(w-h)/2];
    }
    for (int i = (w-h)/2+nsmp+1; i &lt;= (w-h)+nsmp; i++)
        x[i] = 0;

   
   
    sizeax = nhops;
              sizeay = p+1;
    //a = zeros(nhops, p+1);
    //g = zeros(nhops, 1);
    a = new double*[nhops+1];
    for (int i = 1; i &lt;= nhops; i++)
    {
        a[i] = new double[(p+1)+1];
        for (int j = 1; j &lt;= p+1; j++)
        {
           a[i][j] = 0;
        }
    }
    g = new double[nhops+1];
    for (int i = 1; i &lt;= nhops; i++)
    {
        g[i] = 0;
    }

    int n1 = (nhops - 1)*h;
    if (ov == 0)
    {
        //e = zeros(1, npts);
        e = new double[npts+1];
        sizee = npts;
        for (int i = 1; i &lt;= npts; i++)
            e[i] = 0;
    }
    else
    {
        //e = zeros(1, (nhops-1)*h+w);
        int n = (nhops-1)*h+w;
        e = new double[n+n1+1]; //on ajoute n1 smp car plus bas on en as
besoin
        sizee = n;
        for (int i = 1; i &lt;= n+n1; i++)
            e[i] = 0;
    }
   
    double * pre = new double[3];
    pre[1] = 1;
    pre[2] = -0.9;

    double prevX = 0;
    double xx2 = 0;
    for (int i = 1; i &lt;= (w-h)+nsmp; i++)
    {
        xx2 = x[i];
        x[i] = pre[1] * xx2 + pre[2]*prevX;
        prevX = xx2;
    }

    free(pre);int n = w;
    double* rxx = new double[n*2+1];
    double* wxx = new double[w+1];double* rs = new double[w+1];

    nhops = 3;

      // Extract segment of signal
      //xx = x((hop - 1)*h + [1,2,3,4,...,w]);
      double* xx = new double[w+1];
      int nn = (hop - 1)*h;
      for (int i = 1; i &lt;= w; i++)
      {
        xx[i] = x[nn + i];
      }

     

      // Apply hanning window
      //wxx = xx .* hanning(w)';
     
      for (int i = 1; i &lt;= w; i++)
      {
        wxx[i] = xx[i];
      }
      fenetrage(wxx, w);//not CPU friendly
     
      double total = 0;

   
      for (int i = 1; i &lt;= n; i++)
      {
          /*n-1 à n-1 * x[0]
          n-2 a n-1 * x[0 à 1]
          ...
          0    à n-1 * x[0 à n-1]*/
         
          total = 0;
          for (int j = n+1-i; j &lt;= n; j++)
          {
               total = total + wxx[j]*wxx[j-(n+1-i)+1];
          }
         
          rxx[i] = total;
      }

     
     
      for (int i = 1; i &lt; n; i++)
      {
     /*   / 0    à n-2 * x[1 à n-1]
          | ...
          \ 0    à 0    * x[n-1]         */
        total = 0;
        for (int j = 1; j &lt;= n - i; j++)
        {
           //i = 0    j = 0 à n-2    1 à n-1
           //i = n-2  j = 0 à 0          n-1
           total = total + wxx[j]* wxx[j + i];
        }
        rxx[i+n] = total;
      }   
           

     
     

      // extract just the points we need (middle p+1 points)
      //rxx = rxx(w+[0:p]);
      for (int i = 0; i &lt;= p; i++)
      {
        rxx[i+1]= rxx[w+i];
      }
     
 
      //Setup the normal equations
      /*R = toeplitz(rxx(1:p));
     
                     a b c d e
                     b a b c d
      T(a,b,c,d,e) = c b a b c
                     d c b a b
                     e d c b a*/
      double** R = new double*[p+1];
      for (int i = 0; i &lt;= p; i++)
      {
        R[i] = new double[p+1];
      }
      for (int i = p; i > 0; i--)
      {
          for (int j = 1; j <= i ; j++)
          {
              R[p-i+1][j+p-i] = rxx[j];
              R[j+p-i][p-i+1] = rxx[j];
          }
      }

      double** R2 = new double*[p];
      for (int i = 0; i < p; i++)
      {
        R2[i] = new double[p];
      }
      for (int i = 0; i < p; i++)
      {
          for (int j = 0; j < p ; j++)
          {
              R2[i][j] = R[i+1][j+1];
          }
      }
      double* R3 = new double[p*p];
      for (int i = 0; i < p; i++)
      {
          for (int j = 0; j < p ; j++)
          {
              R3[i*p+j] = R[i+1][j+1];
          }
         
      }
//****************************************************************************************









--
View this message in context: 
http://octave.1599824.n4.nabble.com/Octave-C-matrix-Inv-comparison-tp4655291p4655510.html
Sent from the Octave - Maintainers mailing list archive at Nabble.com.