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 :
On 8 Jul 2013, at 18:35, Ed Meyer <
[hidden email]>
wrote:
>
> >> the numbers agree in the first two places so there
may not be programming errors
> >> but the lapack functions used by octave do
equilibration and pivoting to improve
> >> the solution - why not just use lapack (dgesvx)?
>
>
>
> > thank you for your answer
>
> > it "deblocks" me a little ^^. So what you said is that
the C++ routine i use are not as precise > as in Octave ? it
makes sense now that you say it...
>
> > So I looked at the subroutine dgesvx as you told me
>
> > I see that there are a lot of parameters
>
> > Could you help me with those parameters ? I don't
understand a clue
>
> > All i need to do i take the inverse of a <double**
matrix> with <int n> the size of the matrix
>
> > but i don't understand what are all those parameters,
even with the explanation of the routine
>
> > mucho thanks for the help )))
>
> > Jeff
>
> If you haven't called fortran routines from C++ before
there are a few issues to
> be aware of; see e.g.
>
http://www.yolinux.com/TUTORIALS/LinuxTutorialMixingFortranAndC.html
> the main things you have to be aware of is how fortran
stores matrices, all args
> are pointers, and the naming convention (e.g. on linux
dgesv is dgesv_).
>
> It would probably be worthwhile trying the simpler routine
dgesv which doesn't do
> equilibration but has a simpler interface:
>
> dgesv (n, nrhs, a, lda, ipiv, b, ldb, info)
>
> where "b" is an (n,n) identity matrix (I like using stl
vectors for temp arrays):
>
> vector<double> b(n*n, 0.0);
> for (i=0; i<n; i++)
> b[i*(n+1)] = 1.0;
>
> ipiv is an output array of ints:
>
> vector<int> ipiv(n);
> and you call with pointers for ints:
>
> dgesv_ (&n, &n, &a[0], &n, &ipiv[0],
&b[0], &ldb, &info);
>
> Calling fortran from C++ is a pain but it does open up a
huge amount
> of code once you figure out how to do it.
>
> --
> Ed Meyer
Otherwise,
If you want a linear algebra library with a more typical C++
interface you may want to have a look at eigen [1]
c.
[1]
http://eigen.tuxfamily.org
To unsubscribe from Octave/C++ matrix Inv() comparison,
click
here.
NAML
View this message in context: Re: Octave/C++ matrix Inv() comparison