# Demo entry 6498721

mpfr_cg.c

Submitted by T.Kouya on May 30, 2017 at 09:02
Language: C. Code size: 12.1 kB.

/**********************************************************/
/* Conjugate-Gradient Linear Solver with MPFR/GMP         */
/* ANSI C standard                                        */
/*                                                        */
/* Version 0.0: 2016-11-17 (Thu) First published          */
/**********************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "cblas.h"
#include "lapacke.h"

//#define MIN(i, j) ((i) < (j) ? (i) : (j))
//#define MAX(i, j) ((i) > (j) ? (i) : (j))

// MPFR/GMP is used
#include "gmp.h"
#include "mpfr.h"

// C linear compucation with double, QD, MPFR/GMP
#include "linear_c.h"

#define ORG_PREC 0
#define DIF_PREC 2

// d_matrix := (double)matrix
void mpfr_set_d_array(double d_array[], mpfr_t array[], int dim)
{
int i;

for(i = 0; i < dim; i++)
d_array[i] = mpfr_get_d(array[i], _tk_default_rmode);
}

// diff := x - y
void mpfr_diff_array(mpfr_t diff[], mpfr_t x[], mpfr_t y[], int dim)
{
unsigned long prec;
mpfr_t minus_1;

prec = mpfr_get_max_prec_array(y, dim);
mpfr_init2(minus_1, prec); mpfr_set_ui(minus_1, 1UL, _tk_default_rmode); mpfr_neg(minus_1, minus_1, _tk_default_rmode);

mpfr_myaxpy(diff, minus_1, y, x, dim);

mpfr_clear(minus_1);
}

// coef := V^T * x -> x = sum_i coef[i] * v_i
void get_coef_ev(double d_coef[], double d_x[], mpfr_t x[], double d_ev[], int dim)
{
int i;

// d_x := (double)x
mpfr_set_d_array(d_x, x, dim);

// d_coef := V^T * x
cblas_dgemv(CblasRowMajor, CblasTrans, dim, dim, 1.0, d_ev, dim, d_x, 1, 0.0, d_coef, 1);

for(i = 0; i < dim - 1; i++)
printf("%10.3e, ", fabs(d_coef[i]));
printf("%10.3e\n", fabs(d_coef[dim - 1]));
}

// Solve A * x = b
int mpfr_conjugate_gradient_error(mpfr_t *x[3], mpfr_t A[], mpfr_t b[], int dim, mpfr_t rel_tol, mpfr_t abs_tol, int maxitimes, unsigned long diff_prec, double d_ev[])
{
int itimes, i;

// vectors
mpfr_t *r[2], *p[2], *w[2], *diff_p, *diff_r, *diff_x;
double *d_coef, *d_tmp;

// constant
mpfr_t alpha[2], beta[2], tmp_val[2], init_r_norm2[2], r_new_norm2[2];

// Initialize
prec = mpfr_get_max_prec_array(x[ORG_PREC], dim); // x's precision is locally used

printf("Original Prec in bits: %ld\n", prec);

mpfr_inits2(prec, alpha[ORG_PREC], beta[ORG_PREC], tmp_val[ORG_PREC], init_r_norm2[ORG_PREC], r_new_norm2[ORG_PREC], (mpfr_ptr)0);

r[ORG_PREC] = (mpfr_t *)calloc(dim, sizeof(mpfr_t));
p[ORG_PREC] = (mpfr_t *)calloc(dim, sizeof(mpfr_t));
w[ORG_PREC] = (mpfr_t *)calloc(dim, sizeof(mpfr_t));
diff_p = (mpfr_t *)calloc(dim, sizeof(mpfr_t));
diff_r = (mpfr_t *)calloc(dim, sizeof(mpfr_t));
diff_x = (mpfr_t *)calloc(dim, sizeof(mpfr_t));
d_coef = (double *)calloc(dim, sizeof(double));
d_tmp = (double *)calloc(dim, sizeof(double));

mpfr_init2_array(r[ORG_PREC], dim, prec);
mpfr_init2_array(p[ORG_PREC], dim, prec);
mpfr_init2_array(w[ORG_PREC], dim, prec);
mpfr_init2_array(diff_p, dim, diff_prec);
mpfr_init2_array(diff_r, dim, diff_prec);
mpfr_init2_array(diff_x, dim, diff_prec);

// x_0 := 0
mpfr_set0(x[ORG_PREC], dim);

// r := b - A * x_0;
mpfr_mycopy(r[ORG_PREC], b, dim); // x_0 = 0
mpfr_mycopy(r[ADD_PREC], b, dim); // x_0 = 0

// init_r_norm2 = ||r||_2
mpfr_mynorm2(init_r_norm2[ORG_PREC], r[ORG_PREC], dim);

// p := r
mpfr_mycopy(p[ORG_PREC], r[ORG_PREC], dim);

// main loop
for(itimes = 0; itimes < maxitimes; itimes++)
{
// w := A * p
mpfr_mymv(w[ORG_PREC], A, p[ORG_PREC], dim);

// alpha := (r, p) / (p, A * p)
mpfr_mydotp(alpha[ORG_PREC], r[ORG_PREC], p[ORG_PREC], dim);
mpfr_mydotp(tmp_val[ORG_PREC], p[ORG_PREC], w[ORG_PREC], dim);

if(mpfr_cmp_ui(tmp_val[ORG_PREC], 0UL) == 0)
{
itimes = -1;
break;
}
//alpha /= tmp_val;
mpfr_div(alpha[ORG_PREC], alpha[ORG_PREC], tmp_val[ORG_PREC], _tk_default_rmode);

// x := x + alpha * p
mpfr_myaxpy(x[ORG_PREC], alpha[ORG_PREC], p[ORG_PREC], x[ORG_PREC], dim);
mpfr_mycopy(x[DIF_PREC], diff_x, dim);
//get_coef_ev(d_coef, d_tmp, x[DIF_PREC], d_ev, dim);
//get_coef_ev(d_coef, d_tmp, x[ORG_PREC], d_ev, dim);

// r_new := r - alpha * A * p
mpfr_mydotp(beta[ORG_PREC], r[ORG_PREC], r[ORG_PREC], dim);
mpfr_neg(alpha[ORG_PREC], alpha[ORG_PREC], _tk_default_rmode); // alpha := -alpha
mpfr_myaxpy(r[ORG_PREC], alpha[ORG_PREC], w[ORG_PREC], r[ORG_PREC], dim);
//get_coef_ev(d_coef, d_tmp, diff_r, d_ev, dim);
//get_coef_ev(d_coef, d_tmp, r[ORG_PREC], d_ev, dim);

// beta := ||r_new||_2^2 / ||r||_2^2
//beta = r_new_norm2 * r_new_norm2 / beta;
mpfr_mynorm2(r_new_norm2[ORG_PREC], r[ORG_PREC], dim);
mpfr_sqr(tmp_val[ORG_PREC], r_new_norm2[ORG_PREC], _tk_default_rmode);
mpfr_div(beta[ORG_PREC], tmp_val[ORG_PREC], beta[ORG_PREC], _tk_default_rmode);

// check residual
//if(r_new_norm2 <= rel_tol * init_r_norm2 + abs_tol)
mpfr_mul(tmp_val[ORG_PREC], rel_tol, init_r_norm2[ORG_PREC], _tk_default_rmode);
if(mpfr_cmp(r_new_norm2[ORG_PREC], tmp_val[ORG_PREC]) <= 0)
break;

// tmp_val := r_new_norm2 / init_r_norm2
mpfr_div(tmp_val[ORG_PREC], r_new_norm2[ORG_PREC], init_r_norm2[ORG_PREC], _tk_default_rmode);
//mpfr_printf("%3d %10.3RNe %10.3RNe\n", itimes, tmp_val[ORG_PREC], tmp_val[ADD_PREC]);

// p := r + beta * p
mpfr_myaxpy(p[ORG_PREC], beta[ORG_PREC], p[ORG_PREC], r[ORG_PREC], dim);
//get_coef_ev(d_coef, d_tmp, diff_p, d_ev, dim);
//get_coef_ev(d_coef, d_tmp, p[ORG_PREC], d_ev, dim);

}

// clean
mpfr_clear_array(r[ORG_PREC], dim);
mpfr_clear_array(p[ORG_PREC], dim);
mpfr_clear_array(w[ORG_PREC], dim);
mpfr_clear_array(diff_p, dim);
mpfr_clear_array(diff_r, dim);
mpfr_clear_array(diff_x, dim);
free(r[ORG_PREC]);
free(p[ORG_PREC]);
free(w[ORG_PREC]);
free(diff_p);
free(diff_r);
free(diff_x);
free(d_coef);
free(d_tmp);

mpfr_clears(alpha[ORG_PREC], beta[ORG_PREC], tmp_val[ORG_PREC], init_r_norm2[ORG_PREC], r_new_norm2[ORG_PREC], (mpfr_ptr)0);

return itimes;
}

int main(int argc, char *argv[])
{
unsigned long prec;
int i, j, dimension, cg_itimes;
mpfr_t *matrix, *true_x, *b, *x[3];
mpfr_t rel_tol, abs_tol, tmp_val;
// lapack
int info;
double *d_matrix, *d_eig;

if(argc <= 2)
{
fprintf(stderr, "USAGE: %s [dimension] [prec]\n", argv[0]);
return EXIT_SUCCESS;
}

// get precitsion in bits
prec = (unsigned long)atoi(argv[2]);

dimension = atoi(argv[1]);

if(dimension <= 1)
{
fprintf(stderr, "ERROR: dimension = %d is illegal!", dimension);
return EXIT_FAILURE;
}

// initialize
mpfr_inits2(prec, rel_tol, abs_tol, tmp_val, (mpfr_ptr)NULL);

matrix = (mpfr_t *)calloc(dimension * dimension, sizeof(mpfr_t));
true_x = (mpfr_t *)calloc(dimension, sizeof(mpfr_t));
x[0]   = (mpfr_t *)calloc(dimension, sizeof(mpfr_t));
x[1]   = (mpfr_t *)calloc(dimension, sizeof(mpfr_t));
x[2]   = (mpfr_t *)calloc(dimension, sizeof(mpfr_t));
b      = (mpfr_t *)calloc(dimension, sizeof(mpfr_t));
d_matrix = (double *)calloc(dimension * dimension, sizeof(double));
d_eig = (double *)calloc(dimension, sizeof(double));

printf("callocs, dim = %d, prec = %ld\n", dimension, prec);

mpfr_init2_array(matrix, dimension * dimension, prec);
mpfr_init2_array(true_x, dimension, prec);
mpfr_init2_array(x[0]  , dimension, prec);
mpfr_init2_array(x[1]  , dimension, prec);
mpfr_init2_array(x[2]  , dimension, prec);
mpfr_init2_array(b     , dimension, prec);

// set test problem
set_test_mpfr_linear_eq(matrix, true_x, b, dimension);
printf("set_test_mpfr_linear_eq, dim = %d, prec = %ld\n", dimension, prec);

// d_matrix := (double)matrix
mpfr_set_d_array(d_matrix, matrix, dimension * dimension);

// solve A * V = \lambda * V
// 'V' ... get eigenvectors
info = LAPACKE_dsyev(LAPACK_ROW_MAJOR, 'V', 'U', dimension, d_matrix, dimension, d_eig);

// error occurs if info > 0
if(info > 0)
{
printf("QR decomposition failed! (%d) \n", info);
return EXIT_FAILURE;
}
else if(info < 0)
{
printf("%d-th argument of DSYEV is illegal!\n", info);
return EXIT_FAILURE;
}

printf("Envenvectors: \n");
for(i = 0; i < dimension; i++)
{
printf("%3d: ", i);
printf("%25.17e", d_eig[i]);
/*		for(j = 0; j < dimension; j++)
printf("%5g ", d_matrix[i * dim + j]);
*/
printf("\n");

}

mpfr_set_str(rel_tol, "1.0e-30", 10, _tk_default_rmode);
mpfr_set_str(abs_tol, "1.0e-100", 10, _tk_default_rmode);
//cg_itimes = mpfr_conjugate_gradient_error(x, matrix, b, dimension, rel_tol, abs_tol, dimension * 5, prec * 4);
cg_itimes = mpfr_conjugate_gradient_error(x, matrix, b, dimension, rel_tol, abs_tol, dimension * 5, prec * 4, d_matrix);

// print solution
printf("cg iterative times: %d\n", cg_itimes);
for(i = 0; i < dimension; i++)
{
get_mpfr_relerr(tmp_val, x[0][i], true_x[i]);
//mpfr_printf("%3d %25.17RNe %10.3RNe\n", i, x[0][i], tmp_val);
mpfr_printf("%3d %25.17RNe %25.17RNe %10.3RNe\n", i, x[0][i], x[1][i], x[2][i]);
}

// free
mpfr_clears(rel_tol, abs_tol, tmp_val, (mpfr_ptr)NULL);

mpfr_clear_array(matrix, dimension * dimension);
mpfr_clear_array(true_x, dimension);
mpfr_clear_array(x[0], dimension);
mpfr_clear_array(x[1], dimension);
mpfr_clear_array(x[2], dimension);
mpfr_clear_array(b, dimension);

free(matrix);
free(true_x);
free(x[0]);
free(x[1]);
free(x[2]);
free(b);
free(d_matrix);
free(d_eig);

return EXIT_SUCCESS;
}

This snippet took 0.02 seconds to highlight.

Back to the Entry List or Home.