segmentation fault when calling dgetri_() from a c++ program

[enid]

Hi there,

I am trying to call dgetri_() from a c++ program on my MAC OS X. It is a very simple program (I copied it below), and I feel I allocated the memory correctly.. No idea why it fails. It compiles and links alright. Just when I run it, it has a "Segmentation fault" occured when calling dgetri_(). Does someone have some clue? Thanks a million!

cheers,
Enid

/********** main.cpp ************/
#include <stdio.h>
#include <stdlib.h>
#include <iostream.h>

using namespace std;
#include <vector.h>

#include <vecLib/vecLib.h>
#include <vecLib/clapack.h>

#define size 3

int main(int agrc, char** argv)
{
double *work = (double*)malloc(size*size*sizeof(double));
double *tmpmatrix = (double *)malloc(size*size*sizeof(double));
int *index = (int*)malloc(size * sizeof(int));
int info;

tmpmatrix[0] = 1.0; tmpmatrix[1] = 0.0; tmpmatrix[2] = 1.0;
tmpmatrix[3] = 0.0; tmpmatrix[4] = 1.0; tmpmatrix[5] = 0.0;
tmpmatrix[6] = 1.0; tmpmatrix[7] = 0.0; tmpmatrix[8] = 1.0;

int s = size;
int LDA = size;
int lwork = size;
dgetri_((__CLPK_integer *)(&s), (__CLPK_doublereal *)(tmpmatrix), \ (__CLPK_integer *)(&LDA), (__CLPK_integer *)(index), \ (__CLPK_doublereal *)(work), (__CLPK_integer *)(&lwork), \ (__CLPK_integer *)(&info));

if(info != 0) {// failed.
cout << "It failed! (info = " << info << " ) " << endl;
}

free(tmpmatrix);
free(work);
free(index);
}

[Stan Tomov]

Hello,
dgetri inverts a matrix using its LU factorization computed by DGETRF,
so if you wanted to invert tmpmatrix the following will work on linux

Code: Select all

#include <stdio.h>
#include <stdlib.h>
#include <iostream.h>

using namespace std;
#include <vector.h>

#define size 3

extern "C" void dgetri_(int &, double *, int &, int *, double *, int &, int &);
extern "C" void dgetrf_(int &, int &, double *, int &, int *, int &);

int main(int agrc, char** argv)
{
  double *work = new double[size*size];
  double *tmpmatrix = new double[size*size];
  int *index = new int[size];
  int info;

  tmpmatrix[0] = 1.0; tmpmatrix[1] = 0.0; tmpmatrix[2] = 1.0;
  tmpmatrix[3] = 0.0; tmpmatrix[4] = 1.0; tmpmatrix[5] = 0.0;
  tmpmatrix[6] = 0.0; tmpmatrix[7] = 0.0; tmpmatrix[8] = 1.0;

  int s = size;
  int LDA = size;
  int lwork = size*size;
  dgetrf_(s,LDA,tmpmatrix,LDA,index,info);
  dgetri_(s,tmpmatrix, LDA, index, work, lwork, info);

  if(info != 0) {// failed.
    cout << "It failed! (info = " << info << " ) " << endl;
  }
 
  delete [] tmpmatrix;
  delete [] work;
  delete [] index;
}


I only changed tmpmatrix[6] so that the matrix is invertible.
Only dgetri should also work if you initialize the pivoting indices,
for example

Code: Select all

 index[0]=1;
 index[1]=2;
 index[2]=3;

Note that the enumeration starts from 1.

Hope that this helps.
Stan

[enid]

Aha, it works! That was a stupid mistake. I should have read the description of the function more carefully!

Thanks a lot!
Enid.