-
Notifications
You must be signed in to change notification settings - Fork 4
/
linear_eq_dsgesv.c
executable file
·102 lines (86 loc) · 2.47 KB
/
linear_eq_dsgesv.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
/*************************************************/
/* LAPACK/BLAS Tutorial */
/* Solver for Linear equation with DSGESV */
/* Last Update: 2016-11-30 (Wed) T.Kouya */
/*************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "lapacke.h"
#include "cblas.h"
int main()
{
lapack_int i, j, dim, itimes;
lapack_int inc_vec_x, inc_vec_b;
lapack_int *pivot, info;
double *mat_a, *vec_b, *vec_x, *vec_x_approx;
double alpha, beta;
double running_time;
// input dimension of a linear equation to be solved
printf("Dim = "); scanf("%d", &dim);
if(dim <= 0)
{
printf("Illegal dimension! (dim = %d)\n", dim);
return EXIT_FAILURE;
}
// inisialize a matrix and vectors
mat_a = (double *)calloc(dim * dim, sizeof(double));
vec_x = (double *)calloc(dim, sizeof(double));
vec_x_approx = (double *)calloc(dim, sizeof(double));
vec_b = (double *)calloc(dim, sizeof(double));
// input mat_a and vec_x
for(i = 0; i < dim; i++)
{
for(j = 0; j < dim; j++)
{
//mat_a[i * dim + j] = (double)rand() / (double)RAND_MAX;
mat_a[i * dim + j] = 1.0 / (double)(i + j + 1);
if((i + j + 1) % 2 != 0)
mat_a[i * dim + j] = -mat_a[i * dim + j];
}
mat_a[i * dim + i] += 2.0;
vec_x[i] = 1.0 / (double)(i + 1);
}
// size(vec_x) == size(vec_b)
inc_vec_x = inc_vec_b = 1;
// vec_b := 1.0 * mat_a * vec_x + 0.0 * vec_b
alpha = 1.0;
beta = 0.0;
cblas_dgemv(CblasRowMajor, CblasNoTrans, dim, dim, alpha, mat_a, dim, vec_x, inc_vec_x, beta, vec_b, inc_vec_b);
// print
for(i = 0; i < dim; i++)
{
printf("[");
for(j = 0; j < dim; j++)
printf("%10.3f ", mat_a[i * dim + j]);
printf("] %10.3f = %10.3f\n", vec_x[i], vec_b[i]);
}
// initialize pivot
pivot = (lapack_int *)calloc(dim, sizeof(lapack_int));
// single-double mixed presicion iterative refinement method
info = LAPACKE_dsgesv(LAPACK_ROW_MAJOR, dim, 1, mat_a, dim, pivot, vec_b, 1, vec_x_approx, 1, &itimes);
printf("DSGESV info = %d, Iterative Times = %d\n", info, itimes);
// print
printf("calculated x = \n");
for(i = 0; i < dim; i++)
{
printf("%3d -> %3d: ", i, pivot[i]);
printf("%25.17e ", vec_x_approx[i]);
printf("\n");
}
// diff
printf("x - calculated x = \n");
for(i = 0; i < dim; i++)
{
printf("%3d: ", i);
printf("%10.2e ", fabs((vec_x[i] - vec_x_approx[i]) / vec_x[i]));
printf("\n");
}
// free
free(mat_a);
free(vec_x);
free(vec_x_approx);
free(vec_b);
free(pivot);
return EXIT_SUCCESS;
}