pdepe/pde1dm's vectorized option

In pdepe, when vectorized is turned on, the elapsed time is the same as that when vectorized is off. Why no improvement?

The doc of pde1dm says that for vectorized mode, we return coefficients at multiple x locations, so in the pdefunc function

nx=length(x); c = ones(1,nx);

need to replace the original

c=1

While pde1dm requires the above change in the pdefunc function code, when c and s in the pdefunc function aren't changed to their vector form, pdepe still allows vectorized to be turned on. Why such a difference?

[feynman feynman]

any expert of pde1dm/pdepe please?

[Torsten]

Where did you see that you can use the Vectorized option for "pdepe" ? I only see the following options available according to the documentation:

Option structure, specified as a structure array. Use the odeset function to create or modify the option structure. pdepe supports these options:

• 'RelTol'
• 'AbsTol'
• 'NormControl'
• 'InitialStep'
• 'MaxStep'
• 'Events'

In most cases, default values for these options provide satisfactory solutions.

[feynman feynman]

I didn't see the vectorized option for pdepe anywhere either. I just tried turning it on and no error was reported.

I think actually pdefun, icfun etc accept tensors as arguments so maybe vectorization is already adopted?

[Torsten]

I think actually pdefun, icfun etc accept tensors as arguments so maybe vectorization is already adopted?

The inputs to "pdefun" and "icfun" are pointwise in x, not vectorwise. The functions are called in a loop from "pdepe" for every single x-value, and the same is expected for the returned output.

type pdepe

function varargout = pdepe(m,pde,ic,bc,xmesh,t,options,varargin)
%PDEPE  Solve initial-boundary value problems for parabolic-elliptic PDEs in 1-D.
%   SOL = PDEPE(M,PDEFUN,ICFUN,BCFUN,XMESH,TSPAN) solves initial-boundary
%   value problems for small systems of parabolic and elliptic PDEs in one
%   space variable x and time t to modest accuracy. There are npde unknown
%   solution components that satisfy a system of npde equations of the form
%
%   c(x,t,u,Du/Dx) * Du/Dt = x^(-m) * D(x^m * f(x,t,u,Du/Dx))/Dx + s(x,t,u,Du/Dx)
%
%   Here f(x,t,u,Du/Dx) is a flux and s(x,t,u,Du/Dx) is a source term. m must
%   be 0, 1, or 2, corresponding to slab, cylindrical, or spherical symmetry,
%   respectively. The coupling of the partial derivatives with respect to
%   time is restricted to multiplication by a diagonal matrix c(x,t,u,Du/Dx).
%   The diagonal elements of c are either identically zero or positive.
%   An entry that is identically zero corresponds to an elliptic equation and
%   otherwise to a parabolic equation. There must be at least one parabolic
%   equation. An entry of c corresponding to a parabolic equation is permitted
%   to vanish at isolated values of x provided they are included in the mesh
%   XMESH, and in particular, is always allowed to vanish at the ends of the
%   interval. The PDEs hold for t0 <= t <= tf and a <= x <= b. The interval
%   [a,b] must be finite. If m > 0, it is required that 0 <= a. The solution
%   components are to have known values at the initial time t = t0, the
%   initial conditions. The solution components are to satisfy boundary
%   conditions at x=a and x=b for all t of the form
%
%       p(x,t,u) + q(x,t) * f(x,t,u,Du/Dx) = 0
%
%   q(x,t) is a diagonal matrix. The diagonal elements of q must be either
%   identically zero or never zero. Note that the boundary conditions are
%   expressed in terms of the flux rather than Du/Dx. Also, of the two
%   coefficients, only p can depend on u.
%
%   The input argument M defines the symmetry of the problem. PDEFUN, ICFUN,
%   and BCFUN are function handles.
%
%   [C,F,S] = PDEFUN(X,T,U,DUDX) evaluates the quantities defining the
%   differential equation. The input arguments are scalars X and T and
%   vectors U and DUDX that approximate the solution and its partial
%   derivative with respect to x, respectively. PDEFUN returns column
%   vectors: C (containing the diagonal of the matrix c(x,t,u,Dx/Du)),
%   F, and S (representing the flux and source term, respectively).
%
%   U = ICFUN(X) evaluates the initial conditions. For a scalar X, ICFUN
%   must return a column vector, corresponding to the initial values of
%   the solution components at X.
%
%   [PL,QL,PR,QR] = BCFUN(XL,UL,XR,UR,T) evaluates the components of the
%   boundary conditions at time T. XL and XR are scalars representing the
%   left and right boundary points. UL and UR are column vectors with the
%   solution at the left and right boundary, respectively. PL and QL are
%   column vectors corresponding to p and the diagonal of q, evaluated at
%   the left boundary, similarly PR and QR correspond to the right boundary.
%   When m > 0 and a = 0, boundedness of the solution near x = 0 requires
%   that the flux f vanish at a = 0. PDEPE imposes this boundary condition
%   automatically.
%
%   PDEPE returns values of the solution on a mesh provided as the input
%   array XMESH. The entries of XMESH must satisfy
%       a = XMESH(1) < XMESH(2) < ... < XMESH(NX) = b
%   for some NX >= 3. Discontinuities in c and/or s due to material
%   interfaces are permitted if the problem requires the flux f to be
%   continuous at the interfaces and a mesh point is placed at each
%   interface. The ODEs resulting from discretization in space are integrated
%   to obtain approximate solutions at times specified in the input array
%   TSPAN. The entries of TSPAN must satisfy
%       t0 = TSPAN(1) < TSPAN(2) < ... < TSPAN(NT) = tf
%   for some NT >= 3. The arrays XMESH and TSPAN do not play the same roles
%   in PDEPE: The time integration is done with an ODE solver that selects
%   both the time step and formula dynamically. The cost depends weakly on
%   the length of TSPAN. Second order approximations to the solution are made
%   on the mesh specified in XMESH. Generally it is best to use closely
%   spaced points where the solution changes rapidly. PDEPE does not select
%   the mesh in x automatically like it does in t; you must choose an
%   appropriate fixed mesh yourself. The discretization takes into account
%   the coordinate singularity at x = 0 when m > 0, so it is not necessary to
%   use a fine mesh near x = 0 for this reason. The cost depends strongly on
%   the length of XMESH.
%
%   The solution is returned as a multidimensional array SOL. UI = SOL(:,:,i)
%   is an approximation to component i of the solution vector u for
%   i = 1:npde. The entry UI(j,k) = SOL(j,k,i) approximates UI at
%   (t,x) = (TSPAN(j),XMESH(k)).
%
%   SOL = PDEPE(M,PDEFUN,ICFUN,BCFUN,XMESH,TSPAN,OPTIONS) solves as above
%   with default integration parameters replaced by values in OPTIONS, an
%   argument created with the ODESET function. Only some of the options of
%   the underlying ODE solver are available in PDEPE - RelTol, AbsTol,
%   NormControl, InitialStep, and MaxStep. See ODESET for details.
%
%   [SOL,TSOL,SOLE,TE,IE] = PDEPE(M,PDEFUN,ICFUN,BCFUN,XMESH,TSPAN,OPTIONS)
%   with the 'Events' property in OPTIONS set to a function handle EVENTS,
%   solves as above while also finding where event functions g(t,u(x,t))
%   are zero. For each function you specify whether the integration is to
%   terminate at a zero and whether the direction of the zero crossing
%   matters. These are the three column vectors returned by EVENTS:
%   [VALUE,ISTERMINAL,DIRECTION] = EVENTS(M,T,XMESH,UMESH).
%   XMESH contains the spatial mesh and UMESH is the solution at the mesh
%   points. Use PDEVAL to evaluate the solution between mesh points.
%   For the I-th event function: VALUE(I) is the value of the function,
%   ISTERMINAL(I) = 1 if the integration is to terminate at a zero of this
%   event function and 0 otherwise. DIRECTION(I) = 0 if all zeros are to be
%   computed (the default), +1 if only zeros where the event function is
%   increasing, and -1 if only zeros where the event function is decreasing.
%   Output TSOL is a column vector of times specified in TSPAN, prior to
%   first terminal event. SOL(j,:,:) is the solution at T(j). TE is a vector
%   of times at which events occur. SOLE(j,:,:) is the solution at TE(j) and
%   indices in vector IE specify which event occurred.
%
%   If UI = SOL(j,:,i) approximates component i of the solution at time TSPAN(j)
%   and mesh points XMESH, PDEVAL evaluates the approximation and its partial
%   derivative Dui/Dx at the array of points XOUT and returns them in UOUT
%   and DUOUTDX:  [UOUT,DUOUTDX] = PDEVAL(M,XMESH,UI,XOUT)
%   NOTE that the partial derivative Dui/Dx is evaluated here rather than the
%   flux. The flux is continuous, but at a material interface the partial
%   derivative may have a jump.
%
%   Example
%         x = linspace(0,1,20);
%         t = [0 0.5 1 1.5 2];
%         sol = pdepe(0,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
%     solve the problem defined by the function pdex1pde with the initial and
%     boundary conditions provided by the functions pdex1ic and pdex1bc,
%     respectively. The solution is obtained on a spatial mesh of 20 equally
%     spaced points in [0 1] and it is returned at times t = [0 0.5 1 1.5 2].
%     Often a good way to study a solution is plot it as a surface and use
%     Rotate 3D. The first unknown, u1, is extracted from sol and plotted with
%         u1 = sol(:,:,1);
%         surf(x,t,u1);
%     PDEX1 shows how this problem can be coded using subfunctions. For more
%     examples see PDEX2, PDEX3, PDEX4, PDEX5.  The examples can be read
%     separately, but read in order they form a mini-tutorial on using PDEPE.
%
%   See also PDEVAL, ODE15S, ODESET, ODEGET, FUNCTION_HANDLE.

%   The spatial discretization used is that of R.D. Skeel and M. Berzins, A
%   method for the spatial discretization of parabolic equations in one space
%   variable, SIAM J. Sci. Stat. Comput., 11 (1990) 1-32. The time
%   integration is done with ODE15S. PDEPE exploits the capabilities this
%   code has for solving the differential-algebraic equations that arise when
%   there are elliptic equations and for handling Jacobians with specified
%   sparsity pattern.

%   Elliptic equations give rise to algebraic equations after discretization.
%   Initial values taken from the given initial conditions may not be
%   "consistent" with the discretization, so PDEPE may have to adjust these
%   values slightly before beginning the time integration. For this reason
%   the values output for these components at t0 may have a discretization
%   error comparable to that at any other time. No adjustment is necessary
%   for solution components corresponding to parabolic equations. If the mesh
%   is sufficiently fine, the program will be able to find consistent initial
%   values close to the given ones, so if it has difficulty initializing,
%   try refining the mesh.

%   Lawrence F. Shampine and Jacek Kierzenka
%   Copyright 1984-2022 The MathWorks, Inc.
%   $Revision: 1.8.4.9.12.1 $  $Date: 2013/09/27 03:10:20 $

if nargin < 7
    options = [];
end

switch m
    case {0, 1, 2}
    otherwise
        error(message('MATLAB:pdepe:InvalidM'))
end

nt = length(t);
if nt < 3
    error(message('MATLAB:pdepe:TSPANnotEnoughPts'))
end
if any(diff(t) <= 0)
    error(message('MATLAB:pdepe:TSPANnotIncreasing'))
end

xmesh = xmesh(:);
if m > 0 && xmesh(1) < 0
    error(message('MATLAB:pdepe:NegXMESHwithPosM'))
end
nx = length(xmesh);
if nx < 3
    error(message('MATLAB:pdepe:XMESHnotEnoughPts'))
end
if any(diff(xmesh) <= 0)
    error(message('MATLAB:pdepe:XMESHnotIncreasing'))
end

% Initialize the nx-1 points xi where functions will be evaluated
% and as many coefficients in the difference formulas as possible.
% For problems with coordinate singularity, use 'singular' formula
% on all subintervals.
singular = (xmesh(1) == 0 && m > 0);
xL = xmesh(1:end-1);
xR = xmesh(2:end);
xM = xL + 0.5*(xR-xL);
switch m
    case 0
        xi = xM;
        zeta = xi;
    case 1
        if singular
            xi = (2/3)*(xL.^2 + xL.*xR + xR.^2) ./ (xL+xR);
        else
            xi = (xR-xL) ./ log(xR./xL);
        end
        zeta = (xi .* xM).^(1/2);
    case 2
        if singular
            xi = (2/3)*(xL.^2 + xL.*xR + xR.^2) ./ (xL+xR);
        else
            xi = xL .* xR .* log(xR./xL) ./ (xR-xL);
        end
        zeta = (xL .* xR .* xM).^(1/3);
end
if singular
    xim = (zeta .^ (m+1))./xi;
else
    xim = xi .^ m;
end
zxmp1 = (zeta.^(m+1) - xL.^(m+1)) / (m+1);
xzmp1 = (xR.^(m+1) - zeta.^(m+1)) / (m+1);

% Form the initial values with a column of unknowns at each
% mesh point and then reshape into a column vector for dae.
temp = feval(ic,xmesh(1),varargin{:});
if size(temp,1) ~= numel(temp)
    error(message('MATLAB:pdepe:InvalidOutputICFUN'))
end

npde = length(temp);
y0 = zeros(npde,nx);
y0(:,1) = temp;
for j = 2:nx
    y0(:,j) = feval(ic,xmesh(j),varargin{:});
end

% Classify the equations so that a constant, diagonal mass matrix
% can be formed. The entries of c are to be identically zero or
% vanish only at the entries of xmesh, so need be checked only at one
% point not in the mesh.
[U,Ux] = pdentrp(singular,m,xmesh(1),y0(:,1),xmesh(2),y0(:,2),xi(1));
[c,f,s] = feval(pde,xi(1),t(1),U,Ux,varargin{:});
if any([size(c,1),size(f,1),size(s,1)]~=npde) || ...
        any([numel(c),numel(f),numel(s)]~=npde)
    error(message('MATLAB:pdepe:UnexpectedOutputPDEFUN',sprintf('%d',npde)))
end
[pL,qL,pR,qR] = feval(bc,xmesh(1),y0(:,1),xmesh(nx),y0(:,nx),t(1),varargin{:});
if any([size(pL,1),size(qL,1),size(pR,1),size(qR,1)]~=npde) || ...
        any([numel(pL),numel(qL),numel(pR),numel(qR)]~=npde)
    error(message('MATLAB:pdepe:UnexpectedOutputBCFUN',sprintf('%d',npde)))
end

D = ones(npde,nx);
D( c == 0, 2:nx-1) = 0;
if ~singular
    D( qL == 0, 1) = 0;
end
D( qR == 0, nx) = 0;
M = spdiags( D(:), 0, npde*nx, npde*nx);

% Construct block-diagonal pattern of Jacobian matrix
S = kron( spdiags(ones(nx,3),-1:1,nx,nx), ones(npde,npde));

% Extract relevant options and augment with new ones
reltol = odeget(options,'RelTol',[]);
abstol = odeget(options,'AbsTol',[]);
normcontrol = odeget(options,'NormControl',[]);
initialstep = odeget(options,'InitialStep',[]);
maxstep = odeget(options,'MaxStep',[]);
events = odeget(options,'Events',[]);  % events(m,t,xmesh,umesh)
hasEvents = ~isempty(events);
if hasEvents
    eventfcn = @(t,y) events(m,t,xmesh,y,varargin{:});
else
    eventfcn = [];
end
opts = odeset('RelTol',reltol,'AbsTol',abstol,'NormControl',normcontrol,...
    'InitialStep',initialstep,'MaxStep',maxstep,'Events',eventfcn,...
    'Mass',M,'JPattern',S);

% Call DAE solver
tfinal = t(end);
try
    if hasEvents
        [t,y,te,ye,ie] = ode15s(@pdeodes,t,y0(:),opts);
    else
        [t,y] = ode15s(@pdeodes,t,y0(:),opts);
    end
catch ME
    if strcmp(ME.identifier,'MATLAB:daeic12:IndexGTOne')
        error(message('MATLAB:pdepe:SpatialDiscretizationFailed'))
    else
        rethrow(ME);
    end
end

% Verify and process the solution
if t(end) ~= tfinal
    nt = length(t);
    if ~hasEvents || (te(end) ~= t(end))  % did not stop on a terminal event
        warning(message('MATLAB:pdepe:TimeIntegrationFailed',sprintf('%e',t(end))));
    end
end

sol = zeros(nt,nx,npde);
for i = 1:npde
    sol(:,:,i) = y(:,i:npde:end);
end
varargout{1} = sol;
if hasEvents % [sol,t,sole,te,ie] = pdepe(...)
    varargout{2} = t(:);
    sole = zeros(length(te),nx,npde);
    for i = 1:npde
        sole(:,:,i) = ye(:,i:npde:end);
    end
    varargout{3} = sole;
    varargout{4} = te(:);
    varargout{5} = ie(:);
end

%---------------------------------------------------------------------------
% Nested functions
%---------------------------------------------------------------------------

    function dudt = pdeodes(tnow,y)
        %PDEODES  Assemble the difference equations and evaluate the time derivative
        %   for the ODE system.

        u = reshape(y,npde,nx);
        up = zeros(npde,nx);
        [U,Ux] = pdentrp(singular,m,xmesh(1),u(:,1),xmesh(2),u(:,2),xi(1));
        [cL,fL,sL] = feval(pde,xi(1),tnow,U,Ux,varargin{:});

        %  Evaluate the boundary conditions
        [pL,qL,pR,qR] = feval(bc,xmesh(1),u(:,1),xmesh(nx),u(:,nx),tnow,varargin{:});

        %  Left boundary
        if singular
            denom = cL;
            denom(denom == 0) = 1;
            up(:,1) = (sL + (m+1) * fL / xi(1)) ./ denom;
        else
            up(:,1) = pL;
            idx = (qL ~= 0);
            denom = (qL(idx)/xmesh(1)^m) .* (zxmp1(1)*cL(idx));
            denom(denom == 0) = 1;
            up(idx,1) = ( pL(idx) + (qL(idx)/xmesh(1)^m) .* ...
                (xim(1)*fL(idx) + zxmp1(1)*sL(idx))) ./ denom;
        end
        %  Interior points
        for ii = 2:nx-1
            [U,Ux] = pdentrp(singular,m,xmesh(ii),u(:,ii),xmesh(ii+1),u(:,ii+1),xi(ii));
            [cR,fR,sR] = feval(pde,xi(ii),tnow,U,Ux,varargin{:});

            denom = zxmp1(ii) * cR + xzmp1(ii-1) * cL;
            denom(denom == 0) = 1;
            up(:,ii) = ((xim(ii) * fR - xim(ii-1) * fL) + ...
                (zxmp1(ii) * sR + xzmp1(ii-1) * sL)) ./ denom;

            cL = cR;
            fL = fR;
            sL = sR;
        end
        %  Right boundary
        up(:,nx) = pR;
        idx = (qR ~= 0);
        denom = -(qR(idx)/xmesh(nx)^m) .* (xzmp1(nx-1)*cL(idx));
        denom(denom == 0) = 1;
        up(idx,nx) = ( pR(idx) + (qR(idx)/xmesh(nx)^m) .* ...
            (xim(nx-1)*fL(idx) - xzmp1(nx-1)*sL(idx))) ./ denom;

        dudt = up(:);
    end  % pdeodes

% --------------------------------------------------------------------------

end  % pdepe

[feynman feynman]

Oh yes. Then why is there no error reported when vectorized is turned on as in?

options = odeset('Vectorized','on');
sol = pdepe(0, @pde, @ic, @bc, x, t, options);

[Torsten]

Because this option is simply ignored and it is assumed that you read the documentation. That's why no error message is reported. Only the documented options are tested for in the options structure supplied by the user.

But just look through the code - that's why I included it:

% Extract relevant options and augment with new ones 
reltol = odeget(options,'RelTol',[]); 
abstol = odeget(options,'AbsTol',[]); 
normcontrol = odeget(options,'NormControl',[]); 
initialstep = odeget(options,'InitialStep',[]); 
maxstep = odeget(options,'MaxStep',[]); 
events = odeget(options,'Events',[]);  %events(m,t,xmesh,umesh) 
hasEvents = ~isempty(events); 
if hasEvents 
    eventfcn = @(t,y) events(m,t,xmesh,y,varargin{:}); 
else 
    eventfcn = []; 
end 
opts = odeset('RelTol',reltol,'AbsTol',abstol,'NormControl',normcontrol,... 
              'InitialStep',initialstep,'MaxStep',maxstep,'Events',eventfcn,... 
              'Mass',M,'JPattern',S);
% Call DAE solver
tfinal = t(end);
try
    if hasEvents
        [t,y,te,ye,ie] = ode15s(@pdeodes,t,y0(:),opts);
    else
        [t,y] = ode15s(@pdeodes,t,y0(:),opts);
    end
catch ME
    if strcmp(ME.identifier,'MATLAB:daeic12:IndexGTOne')
        error(message('MATLAB:pdepe:SpatialDiscretizationFailed'))
    else
        rethrow(ME);
    end
end

[feynman feynman]

thank you so much for finding it out! Then I guess there won't be too much difference in efficiency between looping and vectorization unless the computation is really huge. It's interesting that pde1dm provides vectorization as an additional option.

[Torsten]

And you think it's not possible to copy the "pdepe" code and change the two calls to ode15s by two calls to ode23t ?

[feynman feynman]

Why do you mention ode23t? When should it replace ode15s?

[Torsten]

I just remembered one of your thousand questions when you asked whether it's possible to replace the integrator within "pdepe" ... :-)

[feynman feynman]

:) I don't think anything in a copied pdepe code can be changed, for otherwise it reports errors

[Torsten]

If you rename the "pdepe" into something else and also copy the function "pdentrp" called from the "pdepe" code, it doesn't work ?

[feynman feynman]

right, since pdentrp isn't recognized by matlab.

[Torsten]

type pdentrp

function [U,Ux] = pdentrp(singular,m,xL,uL,xR,uR,xout)
%PDENTRP  Interpolation helper function for PDEPE.
%   [U,UX] = PDENTRP(M,XL,UL,XR,UR,XOUT) uses solution values UL at XL and UR at XR
%   for successive mesh points XL < XR to interpolate the solution values U and
%   the partial derivative with respect to x, UX, at arguments XOUT(i) with
%   XL <= XOUT(i) <= XR.  UL and UR are column vectors. Column i of the output
%   arrays U, UX correspond to XOUT(i).
%
%   See also PDEPE, PDEVAL, PDEODES.

%   Lawrence F. Shampine and Jacek Kierzenka
%   Copyright 1984-2023 The MathWorks, Inc.

xout = xout(:).';
nout = length(xout);

uRL = uR - uL;

% Use singular interpolant on all subintervals.
if singular
    U  = uL + uRL*((xout .^ 2 - xL^2) / (xR^2 - xL^2));
    Ux =      uRL*(2*xout / (xR^2 - xL^2));
else
    if m == 0
        U  = uL + uRL*( (xout - xL) / (xR - xL));
        Ux =      uRL*(ones(1,nout) / (xR - xL));
    elseif m == 1
        U  = uL + uRL*(log(xout/xL) / log(xR/xL));
        Ux =      uRL*( (1 ./ xout) / log(xR/xL));
    elseif m == 2
        U  = uL + uRL*((xR ./ xout) .* ((xout - xL)/(xR - xL)));
        Ux =      uRL*((xR ./ xout) .* (xL ./ xout)/(xR - xL));
    end
end

[feynman feynman]

oh thanks so much!! I didn't know type worked when open and help didn't return anything. But for other functions, all these commands show results.

[feynman feynman]

@Torsten"The inputs to "pdefun" and "icfun" are pointwise in x, not vectorwise. The functions are called in a loop from "pdepe" for every single x-value, and the same is expected for the returned output."

If so is there no point vectorising pdefun and bcfun etc?

[Torsten]

Sure, if you can vectorize this loop in the "pdepe" code and rewrite "pdentrp" such that it accepts array inputs, you are done:

%  Interior points 
for ii = 2:nx-1
    [U,Ux] = pdentrp(singular,m,xmesh(ii),u(:,ii),xmesh(ii+1),u(:,ii+1),xi(ii)); 
    [cR,fR,sR] = feval(pde,xi(ii),tnow,U,Ux,varargin{:}); 
    denom = zxmp1(ii) * cR + xzmp1(ii-1) * cL; 
    denom(denom == 0) = 1; 
    up(:,ii) = ((xim(ii) * fR - xim(ii-1) * fL) + ... 
                (zxmp1(ii) * sR + xzmp1(ii-1) * sL)) ./ denom; 
    cL = cR; 
    fL = fR; 
    sL = sR; 
end

[feynman feynman]

thank you very much for the hint!