#include <stdio.h> 
#include <math.h> 
#include <strings.h>
#include"mkl.h"
#define N 5
#define NB 36
#define PI acos(-1.0)
void main()
{

inline double delta(double r)
{
	if (r <= 2)
		return (1 + cos(PI*r / 2)) / 4;
	else
		return 0;
}

	int i, j, k, b, kb[NB], jb[NB], ib[NB];
	int ipiv[NB],info;
	double  a[NB][NB],  A1[NB][64];
	double body_coord[NB][3];
	char t[20];
	FILE *fpB = fopen("body.log", "r+");
	for (i = 0; i < NB; i++)
	{
		fscanf(fpB, "%lf %lf %lf", &body_coord[i][0], &body_coord[i][1], &body_coord[i][2]);
		body_coord[i][0] += 20; body_coord[i][1] += 20; body_coord[i][2] +=20;
	}
	fclose(fpB);

	for (b = 0; b < NB; b++)
	{
		ib[b] = (int)(floor(body_coord[b][0]));
		jb[b] = (int)(floor(body_coord[b][1]));
		kb[b] = (int)(floor(body_coord[b][2]));
		for (k = 0; k < 4; k++) for (j = 0; j < 4; j++) for (i = 0; i < 4; i++)
		{
			A1[b][i + 4 * j + 16 * k] = delta(ib[b] - 1 + i - body_coord[b][0])*delta(jb[b] - 1 + j - body_coord[b][1])*delta(kb[b] - 1 + k - body_coord[b][2]);
		}
	}

	for (j = 0; j < NB; j++)for (i = 0; i <= j; i++)
	{
		a[j][i] = 0.0;
		for (k = 0; k < 64; k++)a[j][i] +=(double) A1[j][k] * A1[i][k];//A is a symmetric matrix
		//hmli begin
		sprintf(t,"%lf",a[j][i]);
		printf("%s",t);
		a[j][i]=atof(t);
		//hmli end
		a[i][j] = (double)a[j][i];
	}

	info=LAPACKE_dgetrf(LAPACK_ROW_MAJOR, NB, NB, &a[0][0], NB, &ipiv[0]);
	printf("%i\n",info);
	info=LAPACKE_dgetri(LAPACK_ROW_MAJOR, NB, &a[0][0], NB, &ipiv[0]);
	printf("%i\n",info);

	FILE *fp = fopen("ainv.dat", "w");
	for (j = 0; j < NB; j++)
	{
		for (i = 0; i < NB; i++)fprintf(fp, "%30.6f ", a[j][i]);
		fprintf(fp, "\n");
	}
	fclose(fp);
}