/*************************************************/
/* LAPACK Tutorial */
/* Sample Program of DGESV */
/* with Intel Math Kernel */
/* Last Update: 2013-01-21 (Mon) T.Kouya */
// icc -I/opt/intel/include -I/opt/intel/mkl/include -L/opt/intel/mkl/lib/intel64 -L/opt/intel/lib/intel64 blas1_mkl.c -Wl,--start-group /opt/intel/mkl/lib/intel64/libmkl_intel_lp64.a /opt/intel/mkl/lib/intel64/libmkl_sequential.a /opt/intel/mkl/lib/intel64/libmkl_core.a -Wl,--end-group -lpthread
/*************************************************/
#include <stdio.h>
#include <stdlib.h>
#include "mkl.h" // for Intel Math Kernel Library
#include "mkl_cblas.h"
#include "mkl_lapacke.h"
#define min(i, j) (((i) < (j)) ? (i) : (j))
#define max(i, j) (((i) > (j)) ? (i) : (j))
int main()
{
MKL_INT i, j, dim; // dimension of vectors
MKL_INT *pivot, info;
double *ma, *vx, *vb;
double alpha, beta;
char str_mkl_version[1024];
double running_time;
// print MKL version
MKL_Get_Version_String(str_mkl_version, 1024);
printf("%s\n", str_mkl_version);
// input dimension
printf("Dim = "); scanf("%d", &dim);
if(dim <= 0)
{
printf("Illigal dimension! (dim = %d)\n", dim);
return EXIT_FAILURE;
}
// Initialize
pivot = (MKL_INT *)calloc(sizeof(MKL_INT), dim);
ma = (double *)calloc(sizeof(double), dim * dim);
vx = (double *)calloc(sizeof(double), dim);
vb = (double *)calloc(sizeof(double), dim);
// input va and vb
for(i = 0; i < dim; i++)
{
for(j = 0; j < dim; j++)
{
// A = hirbert matrix
//*(ma + i * dim + j) = 1.0 / (double)(i + j + 1);
// A = frank matrix
*(ma + i * dim + j) = (double)(dim - max(i, j));
}
// x = answer
*(vx + i) = (double) (i + 1);
}
// b := A * x
alpha = 1.0;
beta = 0.0;
cblas_dgemv(CblasRowMajor, CblasNoTrans, dim, dim, alpha, ma, dim, vx, 1, beta, vb, 1);
// print
/* printf("A = \n");
for(i = 0; i < dim; i++)
{
printf("%3d: ", i);
for(j = 0; j < dim; j++)
printf("%10f ", *(ma + i * dim + j));
printf("\n");
}
// print
printf("x = \n");
for(i = 0; i < dim; i++)
{
printf("%3d: ", i);
printf("%10f ", *(vx + i));
printf("\n");
}
// print
printf("b = \n");
for(i = 0; i < dim; i++)
{
printf("%3d: ", i);
printf("%10f ", *(vb + i));
printf("\n");
}
*/
// solve A * X = C -> C := X
running_time = dsecnd(); // start
info = LAPACKE_dgesv(LAPACK_ROW_MAJOR, dim, 1, ma, dim, pivot, vb, 1);
running_time = dsecnd() - running_time; // end
// error occurs if info > 0
if(info > 0)
{
printf("pivot ma[%d][%d] is zero!\n", info, info);
return EXIT_FAILURE;
}
printf("Dimension: %d, Running Time(s): %g\n", dim, running_time);
// print
/* printf("calculated x = \n");
for(i = 0; i < dim; i++)
{
printf("%3d: ", i);
printf("%10f ", *(vb + i));
printf("\n");
}
// diff
printf("x - calculated x = \n");
for(i = 0; i < dim; i++)
{
printf("%3d: ", i);
printf("%10f ", fabs(*(vx + i) - *(vb + i)));
printf("\n");
}
*/
// free
free(pivot);
free(ma);
free(vb);
free(vx);
return EXIT_SUCCESS;
}