solve.m

五点差分型多重网格方法：各种插值算子的比较）

%SOLVE  Solve a linear system.
%
%       [X,RESIDS,ITS] = SOLVE(A,B,X0,RTOL,PRTOL,MAX_IT,MAX_TIME,MAX_MFLOP,...
%       RESTART) applies a solver defined by "solver_flag", with the given
%       tolerances and limits, to a linear system AX=B.  The solution X, 
%       residual history RESIDS, and iterations ITS are returned.
%
%       Accesses global variables in "include_flags"

% James Bordner and Faisal Saied
% Department of Computer Science
% University of Illinois at Urbana-Champaign
% 10 April 1995

% Modified for Matlab Version 6 Compatability
% Ryan McKenzie
% University of Kentucky Center for Computational Sciences
% April 2004
%
% For some reason, locally generated variables cannot be seen outside the scope of a
% particular function in version 6 unless they have a global reference. This seems to 
% only occur when the newly generated variable is passed as a parameter. I have taken
% locally generated variables throughout MGLab and "bridged" them to their destination
% routines using global references. It's an ugly fix, so maybe someone should come up
% with a more centralized solution.
%
% Removed stopping criteria and other solver information from the parameter list since
% they are all stored in global constants anyway.
%
% The if-elseif-end statement that was below tried to evaluate flags that shared their
% names with subroutines. This caused each condition in the switch to attempt to execute 
% each of those solver routines on any run regardless of the selected parameters. I just
% did a global replacement on the flags whose names I changed. Also changed it to a switch
% statement to make it more readable.

function [x,resids,its] = solve(A,b,x0)

include_globals
include_bridge_globals
include_flags

% gobally referencing variables from the run routine "bridge"
b = b_in_solve;
A = A_in_solve;
x0 = x0_in_solve;

% "bridging" variables to solution method routines
b_in_sol_method = b;
A_in_sol_method = A;
x_in_sol_method = x0;

x=x0;

disp(sprintf('Running...\n'));

switch solver_flag
    case VMG_SOLVE
        disp(sprintf('VMG Solver\n'));
        [x,resids,its]=vmg(A,b,x0);
    case FMG_SOLVE
        disp(sprintf('FMG Solver\n'));
        [x,resids,its]=fmg(A,b);
    case PCG_SOLVE
        disp(sprintf('PCG Solver\n'));
        [x,resids,its]=pcg(A,b,x0);
    case BICG_STAB_SOLVE
        disp(sprintf('BICG_STAB Solver\n'));
        [x,resids,its]=pbicgstab(A,b,x0);
    case CGS_SOLVE
        disp(sprintf('CGS Solver\n'));
        [x,resids,its]=pcgs(A,b,x0);
    case GMRES_SOLVE
        disp(sprintf('GMRES Solver\n'));
        [x,resids,its]=pgmres(A,b,x0);
    case SOR_SOLVE
        disp(sprintf('SOR Solver\n'));
        [x,resids,its]=sor(A,b,x0);
end

fprintf('Initial Residual = %g \n', norm(b-A*x0));
fprintf('Final residual = %g \n', norm(b-A*x));
disp (sprintf('Finished.\n'));