Page MenuHomec4science
Files F61461338

cg.c
No OneTemporary
Actions

File Metadata

Created
Mon, May 6, 20:15

cg.c
View Options

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "blas.h"
#include "parameters.h"
const double TOLERANCE = 1.0e-10;
/*
cgsolver solves the linear equation A*x = b where A is
of size m x n
*/
void cgsolver( double *A, double *b, double *x, int m, int n ){
double * R;
double * P;
double * Rold;
double * AP;
double alpha;
double beta;
int incx = 1;
int incy = 1;
int lda = 1;
double al = 1.;
double be = 1.;
int k = 0;
R = (double*) malloc(n* sizeof(double));
P = (double*) malloc(n* sizeof(double));
Rold = (double*) malloc(n* sizeof(double));
AP = (double*) malloc(n* sizeof(double));
// matrixTimesVector( A, X, AP, n );
cblas_dgemv (1, 0, m, n, al, A, lda, x, incx, be, AP, incy);
// subVector(R, B, AP, n);
cblas_dadd(n, R, 1, b, -1, AP);
// copyVector(P,R, n);
cblas_dcopy (n, P, incx, R, incy);
while ( k < n ){
cblas_dcopy (n, Rold, incx, R, incy);
cblas_dgemv (1, 0, m, n, al, A, lda, P, incx, be, AP, incy);
// alpha = dotProduct( R, R, n ) / fmax( dotProduct( P, AP, n ), NEARZERO );
alpha = cblas_ddot (n,R,incx,R,incy) / fmax( cblas_ddot(n, P, incx, AP, incy ), NEARZERO );
cblas_dadd(n,x,1.0,P,alpha,x);
// vectorCombination( 1.0, X, alpha, P, X, n ); // Next estimate of solution
cblas_dadd(n,R,1.0,AP,-alpha,Rold);
// vectorCombination( 1.0, R, -alpha, AP, Rold, n ); // Residual
print_vec("Residuel :\n",Rold,n);
print_vec("X :\n",x,n);
printf("\t[Iteration %d] : vectorNorm(R) = %f and dotProduct = %f\n",k,cblas_dnrm2(n,R,incx), cblas_ddot (n,R,incx,R,incy) );
// if ( vectorNorm( R, n ) < TOLERANCE ) break; // Convergence test
if ( cblas_dnrm2(n,R,incx) < TOLERANCE ) break; // Convergence test
// beta = dotProduct( R, R, n ) / fmax( dotProduct( Rold, Rold, n ), NEARZERO );
beta = cblas_ddot (n,R,incx,R,incy) / fmax( cblas_ddot(n, Rold, incx, Rold, incy ), NEARZERO );
cblas_dadd(n,R,1.0,P,beta,P);
// vectorCombination( 1.0, R, beta, P, P, n ); // Next gradient
printf("\t[Iteration %d] : alpha = %f, beta = %f\n",k,alpha, beta);
k++;
}
free(R);
free(P);
free(Rold);
free(AP);
}

Event Timeline

c4science · Help