function [ctempnew,outputflag] = matsolve(ML,MR,ctempold,R,bcvector);
%%
%% Solve ML*ctempnew = MR*ctempold + (lambda/2) c(x=0,old) + (lambda/2) c(x=0,new).
%% bcvector holds (lambda/2) c(x=0,old) + (lambda/2) c(x=0,new).
%% Use preconditioned conjugate gradient with incomplete Choleski
%% preconditioning to solve the linear system.  initial guess: c_old
RHS = MR*ctempold + bcvector;
numiter = 30;
tol = 1.e-15;
[ctempnew,outputflag]=pcg(ML,RHS,tol,numiter,R',R,ctempold);
return;