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

Thanks you

I tried both dgetrf/dgetri and dsytrf/dsytri to inverse my matrix but i still have differences between Lapack and Octave. For example in Octave i have : 1.7338e+010 and with dgetrf/dgetri i have : 17351976742.876453 and with dsytrf/dsytri i have : 17352023489.128361 so you see they differ quite a lot...

Do you have an idea what is the Inv() function of Octave in Lapack's formaty ? i tried to look into the code but i couldn't find anything relevant.

Thanks

Jeff

here is the code i used both in Octave and in C++ if you have time to look into it...

#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 <matrix_rank; j++)
                matrix [i * matrix_rank + j] = matrix [j * matrix_rank + i];
        }

        delete [] ivpv;
        delete [] work;
    }
//****************************************************************************************
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 C++ code, simply the same Octave code but in C++

int nsmp = 1024;
    double* d = new double[nsmp+1];
    for (int i = 1; i <= 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 <= 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 <= (w-h)/2; i++)
        x[i] = 0;
    for (int i = (w-h)/2+1; i <= (w-h)/2+nsmp; i++)
    {
        x[i] = x2[i-(w-h)/2];
    }
    for (int i = (w-h)/2+nsmp+1; i <= (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 <= nhops; i++)
    {
        a[i] = new double[(p+1)+1];
        for (int j = 1; j <= p+1; j++)
        {
           a[i][j] = 0;
        }
    }
    g = new double[nhops+1];
    for (int i = 1; i <= 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 <= 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 <= 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 <= (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 <= w; i++)
    {
        xx[i] = x[nn + i];
    }



    // Apply hanning window
    //wxx = xx .* hanning(w)';

    for (int i = 1; i <= w; i++)
    {
        wxx[i] = xx[i];
    }
    fenetrage(wxx, w);//not CPU friendly

    double total = 0;


    for (int i = 1; i <= 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 <= n; j++)
        {
               total = total + wxx[j]*wxx[j-(n+1-i)+1];
        }

        rxx[i] = total;
    }



    for (int i = 1; i < n; i++)
    {
    /*   / 0    à n-2 * x[1 à n-1]
        | ...
        \ 0    à 0    * x[n-1]         */
        total = 0;
        for (int j = 1; j <= 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 <= 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 <= 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];
        }

    }
//****************************************************************************************

Le 08/07/2013 19:45, c.-2 [via Octave] a écrit :

From:	ionone
Subject:	Re: Octave/C++ matrix Inv() comparison
Date:	Tue, 9 Jul 2013 06:31:23 -0700 (PDT)