#include <stdio.h>
#include <math.h>
#include <lapacke.h>
#define N 3

void function(double F[N],double x[N]){
F[0] = pow(x[0], 0.3e1) + pow(x[0], 0.2e1) * x[1] - 0.1e1 * x[0] * x[2] + 0.6e1;
F[1] = exp(x[0]) + exp(x[1]) - 0.1e1 * x[2];
F[2] = -0.2e1 * x[0] * x[2] + pow(x[1], 0.2e1) - 0.4e1;
}

void jacobian(double J[N][N],double x[N]){
J[0][0] = 0.3e1 * pow(x[0], 0.2e1) + 0.2e1 * x[0] * x[1] - 0.1e1 * x[2];
J[0][1] = pow(x[0], 0.2e1);
J[0][2] = -0.1e1 * x[0];
J[1][0] = exp(x[0]);
J[1][1] = exp(x[1]);
J[1][2] = -0.1e1;
J[2][0] = -0.2e1 * x[2];
J[2][1] = 0.2e1 * x[1];
J[2][2] = -0.2e1 * x[0];
}

void printvec(double x[N]){
    for(int i=0;i<N;i++) printf("%.14e ",x[i]);
    printf("\n");
}

int main(){
    double F[N], J[N][N];
    double x[N] = { -1, -2, 1};
    int ipv[N];
    for(int i=0;i<5;i++){
        function(F,x);
        jacobian(J,x);
        // x = x - J^(-1) F
        LAPACKE_dgetrf(LAPACK_ROW_MAJOR,
            N,N,J[0],N,ipv);
        LAPACKE_dgetrs(LAPACK_ROW_MAJOR,
            'N',N,1,J[0],N,ipv,F,1);
        for(int j=0;j<N;j++)
            x[j]=x[j]-F[j];
        printvec(x);
    }
    function(F,x);
    printf("Error:\n");
    printvec(F);
}