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

typedef int32_t integer;
typedef double doublereal;

int dgetrf_(integer *m, integer *n, doublereal *a, integer *lda,
	    integer *ipiv, integer *info);
 
int dgetri_(integer *n, doublereal *a, integer *lda, integer 
	    *ipiv, doublereal *work, integer *lwork, integer *info);

int invert_general_matrix (double *val, int dim)
{
    integer n = dim;
    integer info;
    integer lwork;
    integer *ipiv;
    double *work;
    int err = 0;

    ipiv = malloc(dim * sizeof *ipiv);
    work = malloc(sizeof *work);
    
    if (ipiv == NULL || work == NULL) {
	return 1;
    }

    fprintf(stderr, "calling dgetrf: val=%p, ipiv=%p, work=%p\n",
            (void *) val, (void *) ipiv, (void *) work);
    dgetrf_(&n, &n, val, &n, ipiv, &info);   

    if (info != 0) {
	fprintf(stderr, "dgetrf: matrix is singular\n");
	return 1;
    }
    
    fprintf(stderr, "after dgetrf: val=%p, ipiv=%p, work=%p\n",
	    (void *) val, (void *) ipiv, (void *) work);

    lwork = -1;
    dgetri_(&n, val, &n, ipiv, work, &lwork, &info);

    if (info != 0 || work[0] <= 0.0) {
	return 1;
    }

    fprintf(stderr, "after dgetri workspace query:\n"
	    " val=%p, ipiv=%p, work=%p\n",
	    (void *) val, (void *) ipiv, (void *) work);

    lwork = (integer) work[0];

    printf("dgetri: workspace = %d\n", (int) lwork);

    work = realloc(work, lwork * sizeof *work);
    if (work == NULL) {
	return 1;
    }

    fprintf(stderr, "before 'real' dgetri call:\n"
	    " val=%p, ipiv=%p, work=%p\n",
	    (void *) val, (void *) ipiv, (void *) work);    

    dgetri_(&n, val, &n, ipiv, work, &lwork, &info);
    printf("dgetri: info = %d\n", (int) info);

    fprintf(stderr, "after 'real' dgetri call:\n"
	    " val=%p, ipiv=%p, work=%p\n",
	    (void *) val, (void *) ipiv, (void *) work);    

    free(work);
    free(ipiv);

    if (info != 0) {
	fprintf(stderr, "dgetri: matrix is singular\n");
	err = 1;
    }

    return err;
}

int main (void)
{
    double *A;
    int dim = 10;
    int i, j;
    int err = 0;

    A = malloc(dim * dim * sizeof *A);

    for (i=0; i<dim; i++) {
	for (j=0; j<dim; j++) {
	    A[j*dim+i] = 0.0/0.0;
	}
    }

    invert_general_matrix(A, dim);
    free(A);

    return 0;
}