/* solving the matrix equation A*x=b using LAPACK */
#include <vecLib/clapack.h>
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include "nrutil.h"
#include <time.h>

#define NR_END 1


float *vector(long nl, long nh)
/* allocate a float vector with subscript range v[nl..nh] */
{
    float *v;

    v=(float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(float)));

    return v-nl+NR_END;
}

float **matrix(long nrl, long nrh, long ncl, long nch)
/* allocate a float matrix with cubscript range m[nrl..nrh][ncl..nch] */
{
    long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
    float **m;

    /* allocate pointers to rows */
    m=(float **) malloc((size_t)((nrow+NR_END)*sizeof(float*)));

    m += NR_END;
    m -= nrl;

    /* allocate rows and set pointers to them */
    m[nrl]=(float *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float)));

    m[nrl] += NR_END;
    m[nrl] -= nrl;

    for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;

    /* return pointer to array of pointers to rows */
    return m;
}



int solve_lin_eq_LAPACK ( float **A, float *B, int n ){

    double Afortran[n*n], Btmp[n];
    int i,j, nb=1, check, Bpivot[n];
    int fail=0;

    for (i=0; i<n; i++){		/* to call a Fortran routine from C we */
        /* have to transform the matrix to an array*/
        Btmp[i]=(double)B[i+1];
        for(j=0; j<n; j++){
            Afortran[j+n*i]=(double)A[j+1][i+1];
        }
    }

    /*here for general matrix, but if A is symmetric (no-force case) one could
     actually use a special subroutine for those like ssysv_(...) (LAPACK)*/
    dgesv_(&n, &nb, Afortran, &n, Bpivot, Btmp, &n, &check);
    printf("check value=%i\n",check);
    if (check!=0) {
        printf("dgesv failed!\n");
        printf("check value=%i\n",check);
        fail=1;
    }

    if (fail==0){
        for (i=0; i<n; i++)
            B[i+1]=(float)Btmp[i];
    }
    return(fail);

}

int main(){

    float **A, *B;
    int fail;

    A=matrix(1,3,1,3);
    B=vector(1,3);

    A[1][1]=1674.456299;
    A[1][2]=-874.579834;
    A[1][3]=-799.876465;
    A[2][1]=-874.579834;
    A[2][2]=1799.875854;
    A[2][3]=-925.296021;
    A[3][1]=-799.876465;
    A[3][2]=-925.296021;
    A[3][3]=1725.172485;
    B[1]=74.703453;
    B[2]=-1799.875854;
    B[3]=1725.172241;



    fail=solve_lin_eq_LAPACK ( A, B, 3 );
    printf("%f %f %f\n",B[1],B[2],B[3]);
    return(0);
}