#include <cstdio>
#include <cstring>
#include "lapacke.h"
#include "cblas.h"

int main()
{
	const int sz = 10;
	double M[sz*sz] = {		// main matrix
		5, 1.00251875, 5.03E-03, 1.99248125, 1.989975, 1.25E-02, 1.51E-02, 1.98245625, 1.97995, 0.02256875, 
		1.00251875, 5, 2.99248125, 0.0100375, 0.01254375, 1.9849625, 1.98245625, 0.0200625, 0.02256875, 1.9749375, 
		0.00503, 2.99248125, 5, 1.01254375, 0.01505, 1.98245625, 1.97995, 0.02256875, 0.025075, 1.97243125, 
		1.99248125, 0.0100375, 1.01254375, 5, 2.98245625, 0.0200625, 0.02256875, 1.9749375, 1.97243125, 0.0300875, 
		1.989975, 0.01254375, 0.01505, 2.98245625, 5, 1.02256875, 0.025075, 1.97243125, 1.969925, 0.03259375, 
		0.0125, 1.9849625, 1.98245625, 0.0200625, 1.02256875, 5, 2.97243125, 0.0300875, 0.03259375, 1.9649125, 
		0.0151, 1.98245625, 1.97995, 0.02256875, 0.025075, 2.97243125, 5, 1.03259375, 0.0351, 1.96240625, 
		1.98245625, 0.0200625, 0.02256875, 1.9749375, 1.97243125, 0.0300875, 1.03259375, 5, 2.96240625, 0.0401125,
		1.97995, 0.02256875, 0.025075, 1.97243125, 1.969925, 0.03259375, 0.0351, 2.96240625, 5, 1.04261875, 
		0.02256875, 1.9749375, 1.97243125, 0.0300875, 0.03259375, 1.9649125, 1.96240625, 0.0401125, 1.04261875, 5 };

	double Minv[sz*sz] = {	// inverse matrix
  		0.310457947, -0.121891883, 0.068132496, -0.071131317, -0.04345154, 0.012819451, 0.016009522, -0.050265929, -0.052250891, 0.020555943, 
  		-0.121891883, 0.394775269, -0.203727041, 0.079337243, -0.006988023, -0.029946088, -0.04093194, 0.015678336, 0.021025629, -0.052122264, 
  		0.068132496, -0.203727041, 0.407473554, -0.159204457, 0.07269563, -0.060348885, -0.026047026, 0.003527077, 0.01482109, -0.049276304, 
  		-0.071131317, 0.079337243, -0.159204457, 0.407769592, -0.203767941, 0.074264539, -0.013522603, -0.02609559, -0.040100594, 0.01535636,
  		-0.04345154, -0.006988023, 0.07269563, -0.203767941, 0.404897082, -0.155631253, 0.073398233, -0.059451286, -0.028991606, 0.01174095, 
  		0.012819451, -0.029946088, -0.060348885, 0.074264539, -0.155631253, 0.403495193, -0.201745512, 0.071213714, -0.007202764, -0.042309677, 
  		0.016009522, -0.04093194, -0.026047026, -0.013522603, 0.073398233, -0.201745512, 0.404584084, -0.15630474, 0.075667198, -0.068060222, 
  		-0.050265929, 0.015678336, 0.003527077, -0.02609559, -0.059451286, 0.071213714, -0.15630474, 0.401900334, -0.197398696, 0.064486344, 
  		-0.052250891, 0.021025629, 0.01482109, -0.040100594, -0.028991606, -0.007202764, 0.075667198, -0.197398696, 0.389243217, -0.119935656, 
  		0.020555943, -0.052122264, -0.049276304, 0.01535636, 0.01174095, -0.042309677, -0.068060222, 0.064486344, -0.119935656, 0.307596154 };		
	
	double C[sz*sz];		// matrix for checking
	
	// calculate the two norms
	double norm1 = LAPACKE_dlange(LAPACK_COL_MAJOR, '1', sz, sz, M, sz);
	double norm1_inv = LAPACKE_dlange(LAPACK_COL_MAJOR, '1', sz, sz, Minv, sz);
	printf("1-norm(M) = %g\n", norm1);	
	printf("1-norm(Minv) = %g\n", norm1_inv);	
	
	cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, sz, sz, sz, 1, M, sz, Minv, sz, 0, C, sz);	// C = M*Minv
	for (int i = 0; i < sz; i++)
		C[i + sz*i] -= 1;				// C = M*Minv - I
	printf("1-norm(M*Minv-I) = %g\t*** just in case check Minv = inverse(M) with reasonable precision\n", LAPACKE_dlange(LAPACK_COL_MAJOR, '1', sz, sz, C, sz));
	
	// decompose the two matrices
	int info = LAPACKE_dpotrf(LAPACK_COL_MAJOR, 'U', sz, M, sz);
	if (info != 0)
		printf("DPOTRF info = %d\n", info);
	
	info = LAPACKE_dpotrf(LAPACK_COL_MAJOR, 'U', sz, Minv, sz);
	if (info != 0)
		printf("DPOTRF info = %d\n", info);	
	
	// find reciprocal condition numbers
	double rcond, rcond_inv;
	info = LAPACKE_dpocon(LAPACK_COL_MAJOR, 'U', sz, M, sz, norm1, &rcond);
	if (info != 0)
		printf("DPOCON info = %d\n", info);
	
	info = LAPACKE_dpocon(LAPACK_COL_MAJOR, 'U', sz, Minv, sz, norm1_inv, &rcond_inv);
	if (info != 0)
		printf("DPOCON info = %d\n", info);	
	
	printf("1-norm cond(M) = %g\n", 1/rcond);
	printf("1-norm cond(Minv) = %g\n", 1/rcond_inv);
	printf("1-norm(M) * 1-norm(Minv) = %g\n", norm1*norm1_inv);
	
	return 0;
}