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                                        */
/*                                                        */
/* Copyright (c) 2016 Tomonori Kouya, All rights reserved */
/* 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 ADD_PREC 1
#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[])
{
	unsigned long prec, add_prec;
	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
	add_prec = prec + diff_prec; // additional precision 

	printf("Original Prec in bits: %ld\n", prec);
	printf("Additional Prec      : %ld\n", add_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);
	mpfr_inits2(add_prec, alpha[ADD_PREC], beta[ADD_PREC], tmp_val[ADD_PREC], init_r_norm2[ADD_PREC], r_new_norm2[ADD_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));
	r[ADD_PREC] = (mpfr_t *)calloc(dim, sizeof(mpfr_t));
	p[ADD_PREC] = (mpfr_t *)calloc(dim, sizeof(mpfr_t));
	w[ADD_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(r[ADD_PREC], dim, add_prec);
	mpfr_init2_array(p[ADD_PREC], dim, add_prec);
	mpfr_init2_array(w[ADD_PREC], dim, add_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);
	mpfr_set0(x[ADD_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);
	mpfr_mynorm2(init_r_norm2[ADD_PREC], r[ADD_PREC], dim);

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

	// main loop
	for(itimes = 0; itimes < maxitimes; itimes++)
	{
		// w := A * p
		mpfr_mymv(w[ORG_PREC], A, p[ORG_PREC], dim);
		mpfr_mymv(w[ADD_PREC], A, p[ADD_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);
		mpfr_mydotp(alpha[ADD_PREC], r[ADD_PREC], p[ADD_PREC], dim);
		mpfr_mydotp(tmp_val[ADD_PREC], p[ADD_PREC], w[ADD_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);
		mpfr_div(alpha[ADD_PREC], alpha[ADD_PREC], tmp_val[ADD_PREC], _tk_default_rmode);

		// x := x + alpha * p
		mpfr_myaxpy(x[ORG_PREC], alpha[ORG_PREC], p[ORG_PREC], x[ORG_PREC], dim);
		mpfr_myaxpy(x[ADD_PREC], alpha[ADD_PREC], p[ADD_PREC], x[ADD_PREC], dim);
		mpfr_diff_array(diff_x, x[ORG_PREC], x[ADD_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);
		//get_coef_ev(d_coef, d_tmp, x[ADD_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);
		mpfr_mydotp(beta[ADD_PREC], r[ADD_PREC], r[ADD_PREC], dim);
		mpfr_neg(alpha[ADD_PREC], alpha[ADD_PREC], _tk_default_rmode); // alpha := -alpha
		mpfr_myaxpy(r[ADD_PREC], alpha[ADD_PREC], w[ADD_PREC], r[ADD_PREC], dim);
		mpfr_diff_array(diff_r, r[ORG_PREC], r[ADD_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);
		//get_coef_ev(d_coef, d_tmp, r[ADD_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);
		mpfr_mynorm2(r_new_norm2[ADD_PREC], r[ADD_PREC], dim);
		mpfr_sqr(tmp_val[ADD_PREC], r_new_norm2[ADD_PREC], _tk_default_rmode);
		mpfr_div(beta[ADD_PREC], tmp_val[ADD_PREC], beta[ADD_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);
		mpfr_add(tmp_val[ORG_PREC], tmp_val[ORG_PREC], abs_tol, _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_div(tmp_val[ADD_PREC], r_new_norm2[ADD_PREC], init_r_norm2[ADD_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);
		mpfr_myaxpy(p[ADD_PREC], beta[ADD_PREC], p[ADD_PREC], r[ADD_PREC], dim);
		mpfr_diff_array(diff_p, p[ORG_PREC], p[ADD_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);
		get_coef_ev(d_coef, d_tmp, p[ADD_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(r[ADD_PREC], dim);
	mpfr_clear_array(p[ADD_PREC], dim);
	mpfr_clear_array(w[ADD_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(r[ADD_PREC]);
	free(p[ADD_PREC]);
	free(w[ADD_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);
	mpfr_clears(alpha[ADD_PREC], beta[ADD_PREC], tmp_val[ADD_PREC], init_r_norm2[ADD_PREC], r_new_norm2[ADD_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");

	}

	// run conjugate-gradient routine
	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.

Delete this entry (admin only).