/* This function is designed for the Poisson problem discretized by a mesh with N points in each spatial dimension and a total of n mesh points in the mesh, that is, n = N^2 in two and n = N^3 in three dimensions, for instance. The problem is assumed to be split in the last spatial dimension and l_N = N / np and l_n = n / np. The inputs to this function include l_n, l_N, N, and np, which are assumed to be supplied consistently. The vectors l_x, l_r, l_p, l_q, gl, gr should be defined and allocated in the calling routine. The vectors l_x, l_r, l_p, l_q should be of length l_n and gl and gr should be of length l_n / l_N. This function solves Ax = b for x using the (unpreconditioned) conjugatre gradient (CG) method. The matrix-vector product y = A * x for vectors x and y stored in l_x and l_y on each MPI process must be supplied in a function Ax() in a file Ax.c with header file Ax.h and have the prototype void Ax(double *l_y, double *l_x, int l_n, int l_N, int N, int id, int idleft, int idright, int np, MPI_Comm comm, double *gl, double *gr); The user must also supply a parallel dot product with prototype double parallel_dot(double *l_x, double *l_y, int l_n, MPI_Comm comm); We solve for l_x. The initial guess is supplied to the function as the initial content of l_x. flag, relres, and iter are output variables that describe the result of the call. l_r is the right-hand side l_b when the function is called and contains output data of l_r when the call is complete . This implementation of the CG algorithm uses only 4 vectors: l_x = initial guess on input, solution on output, l_r = right-hand side l_b on input, residual on output, l_p = auxiliary variable = search direction, l_q = auxiliary variable = local portion of q = A * p These vectors are of length l_n. Two auxiliary vectors gl and gr of length l_n / l_N must be supplied. */ #include "cg.h" void cg(double *l_x, int *flag, double *relres, int *iter, /* output */ double *l_r, double tol, int maxit, /* input */ double *l_p, double *l_q, int l_n, int l_N, int N, int id, int idleft, int idright, int np, MPI_Comm comm, double *gl, double *gr) { int it, l_i; double n2b, tolb, normr, alpha, pq, beta, rho, rho1; n2b = sqrt(parallel_dot(l_r,l_r,l_n,comm)); /* Norm of the rhs vector b */ if(n2b <= 1.0e-14) { /* if the rhs vector is all zeros... */ for (l_i=0; l_i tolb) && (it < maxit)) { it++; if(it == 1) { for (l_i=0; l_i