| [y,nrowy,ncomp,xpts,nxpts,a,nrowa,alpha,nic,b,nrowb,beta,nfc,igofx,re,ae,iflag,work,ndw,iwork,ndiw,neqivp]=bvsup(y,nrowy,ncomp,xpts,nxpts,a,nrowa,alpha,nic,b,nrowb,beta,nfc,igofx,re,ae,iflag,work,ndw,iwork,ndiw,neqivp); |
function [y,nrowy,ncomp,xpts,nxpts,a,nrowa,alpha,nic,b,nrowb,beta,nfc,igofx,re,ae,iflag,work,ndw,iwork,ndiw,neqivp]=bvsup(y,nrowy,ncomp,xpts,nxpts,a,nrowa,alpha,nic,b,nrowb,beta,nfc,igofx,re,ae,iflag,work,ndw,iwork,ndiw,neqivp);
persistent gt200 gt300 is j k kkkcoe kkkcof kkkg kkks kkksto kkksud kkksvc kkku kkkv kkkws kkkyhp kpts lllcof lllip llliws lllsud lllsvc mxnoni mxnonr ndeq nitemp non nrtemp nxptsm xern1 xern2 xern3 xern4 ;
global ml18jr_1; if isempty(ml18jr_1), ml18jr_1=0; end;
global ml8sz_1; if isempty(ml8sz_1), ml8sz_1=0; end;
global ml5mco_3; if isempty(ml5mco_3), ml5mco_3=0; end;
global ml5mco_6; if isempty(ml5mco_6), ml5mco_6=0; end;
global ml15to_2; if isempty(ml15to_2), ml15to_2=0; end;
global ml15to_1; if isempty(ml15to_1), ml15to_1=0; end;
global ml18jr_2; if isempty(ml18jr_2), ml18jr_2=0; end;
global ml5mco_4; if isempty(ml5mco_4), ml5mco_4=0; end;
global ml5mco_2; if isempty(ml5mco_2), ml5mco_2=0; end;
global ml15to_3; if isempty(ml15to_3), ml15to_3=0; end;
global ml18jr_3; if isempty(ml18jr_3), ml18jr_3=0; end;
global ml5mco_5; if isempty(ml5mco_5), ml5mco_5=0; end;
global ml5mco_1; if isempty(ml5mco_1), ml5mco_1=0; end;
global ml15to_4; if isempty(ml15to_4), ml15to_4=0; end;
global ml15to_5; if isempty(ml15to_5), ml15to_5=0; end;
global ml15to_6; if isempty(ml15to_6), ml15to_6=0; end;
global ml15to_8; if isempty(ml15to_8), ml15to_8=0; end;
global ml15to_7; if isempty(ml15to_7), ml15to_7=0; end;
global ml8sz_2; if isempty(ml8sz_2), ml8sz_2=0; end;
global ml18jr_18; if isempty(ml18jr_18), ml18jr_18=0; end;
global ml8sz_3; if isempty(ml8sz_3), ml8sz_3=0; end;
global ml18jr_11; if isempty(ml18jr_11), ml18jr_11=0; end;
global ml15to_9; if isempty(ml15to_9), ml15to_9=zeros(1,15); end;
global ml8sz_4; if isempty(ml8sz_4), ml8sz_4=0; end;
global ml18jr_12; if isempty(ml18jr_12), ml18jr_12=0; end;
if isempty(is), is=0; end;
global ml15to_10; if isempty(ml15to_10), ml15to_10=0; end;
global ml8sz_5; if isempty(ml8sz_5), ml8sz_5=0; end;
if isempty(j), j=0; end;
if isempty(k), k=0; end;
global ml17bw_4; if isempty(ml17bw_4), ml17bw_4=0; end;
global ml17bw_13; if isempty(ml17bw_13), ml17bw_13=0; end;
global ml17bw_14; if isempty(ml17bw_14), ml17bw_14=0; end;
global ml17bw_5; if isempty(ml17bw_5), ml17bw_5=0; end;
global ml17bw_6; if isempty(ml17bw_6), ml17bw_6=0; end;
global ml17bw_7; if isempty(ml17bw_7), ml17bw_7=0; end;
if isempty(gt200), gt200=0; end;
if isempty(gt300), gt300=0; end;
global ml17bw_8; if isempty(ml17bw_8), ml17bw_8=0; end;
global ml17bw_9; if isempty(ml17bw_9), ml17bw_9=0; end;
global ml17bw_10; if isempty(ml17bw_10), ml17bw_10=0; end;
global ml17bw_11; if isempty(ml17bw_11), ml17bw_11=0; end;
global ml17bw_12; if isempty(ml17bw_12), ml17bw_12=0; end;
if isempty(kkkcoe), kkkcoe=0; end;
if isempty(kkkcof), kkkcof=0; end;
if isempty(kkkg), kkkg=0; end;
global ml17bw_17; if isempty(ml17bw_17), ml17bw_17=0; end;
if isempty(kkks), kkks=0; end;
if isempty(kkksto), kkksto=0; end;
if isempty(kkksud), kkksud=0; end;
if isempty(kkksvc), kkksvc=0; end;
if isempty(kkku), kkku=0; end;
if isempty(kkkv), kkkv=0; end;
if isempty(kkkws), kkkws=0; end;
if isempty(kkkyhp), kkkyhp=0; end;
global ml17bw_1; if isempty(ml17bw_1), ml17bw_1=0; end;
global ml15to_11; if isempty(ml15to_11), ml15to_11=0; end;
global ml15to_12; if isempty(ml15to_12), ml15to_12=0; end;
if isempty(kpts), kpts=0; end;
global ml17bw_15; if isempty(ml17bw_15), ml17bw_15=0; end;
global ml17bw_16; if isempty(ml17bw_16), ml17bw_16=0; end;
if isempty(lllcof), lllcof=0; end;
global ml17bw_18; if isempty(ml17bw_18), ml17bw_18=0; end;
if isempty(lllip), lllip=0; end;
if isempty(llliws), llliws=0; end;
if isempty(lllsud), lllsud=0; end;
if isempty(lllsvc), lllsvc=0; end;
global ml15to_13; if isempty(ml15to_13), ml15to_13=0; end;
global ml5mco_7; if isempty(ml5mco_7), ml5mco_7=0; end;
global ml15to_14; if isempty(ml15to_14), ml15to_14=0; end;
global ml18jr_7; if isempty(ml18jr_7), ml18jr_7=0; end;
if isempty(mxnoni), mxnoni=0; end;
if isempty(mxnonr), mxnonr=0; end;
global ml8sz_6; if isempty(ml8sz_6), ml8sz_6=0; end;
if isempty(ndeq), ndeq=0; end;
global ml18jr_8; if isempty(ml18jr_8), ml18jr_8=0; end;
global ml17bw_3; if isempty(ml17bw_3), ml17bw_3=0; end;
global ml17bw_2; if isempty(ml17bw_2), ml17bw_2=0; end;
global ml18jr_10; if isempty(ml18jr_10), ml18jr_10=0; end;
global ml18jr_15; if isempty(ml18jr_15), ml18jr_15=0; end;
global ml18jr_17; if isempty(ml18jr_17), ml18jr_17=0; end;
global ml8sz_7; if isempty(ml8sz_7), ml8sz_7=0; end;
global ml18jr_5; if isempty(ml18jr_5), ml18jr_5=0; end;
if isempty(nitemp), nitemp=0; end;
if isempty(non), non=0; end;
global ml18jr_6; if isempty(ml18jr_6), ml18jr_6=0; end;
global ml18jr_13; if isempty(ml18jr_13), ml18jr_13=0; end;
if isempty(nrtemp), nrtemp=0; end;
global ml15to_15; if isempty(ml15to_15), ml15to_15=0; end;
%***BEGIN PROLOGUE BVSUP
%***PURPOSE Solve a linear two-point boundary value problem using
% superposition coupled with an orthonormalization procedure
% and a variable-step integration scheme.
%***LIBRARY SLATEC
%***CATEGORY I1B1
%***TYPE SINGLE PRECISION (BVSUP-S, DBVSUP-D)
%***KEYWORDS ORTHONORMALIZATION, SHOOTING,
% TWO-POINT BOUNDARY VALUE PROBLEM
%***AUTHOR Scott, M. R., (SNLA)
% Watts, H. A., (SNLA)
%***DESCRIPTION
%
% **********************************************************************
% subroutine BVSUP solves a LINEAR two-point boundary-value problem
% of the form
% dY/dX = MATRIX(X,U)*Y(X) + G(X,U)
% A*Y(Xinitial) = ALPHA , B*Y(Xfinal) = BETA
%
% Coupled with the solution of the initial value problem
%
% dU/dX = F(X,U)
% U(Xinitial) = ETA
%
% **********************************************************************
% Abstract
% The method of solution uses superposition coupled with an
% orthonormalization procedure and a variable-step integration
% scheme. Each time the superposition solutions start to
% lose their numerical linear independence, the vectors are
% reorthonormalized before integration proceeds. The underlying
% principle of the algorithm is then to piece together the
% intermediate (orthogonalized) solutions, defined on the various
% subintervals, to obtain the desired solutions.
%
% **********************************************************************
% INPUT to BVSUP
% **********************************************************************
%
% NROWY = Actual row dimension of Y in calling program.
% NROWY must be .GE. NCOMP
%
% NCOMP = Number of components per solution vector.
% NCOMP is equal to number of original differential
% equations. NCOMP = NIC + NFC.
%
% XPTS = Desired output points for solution. They must be monotonic.
% Xinitial = XPTS(1)
% Xfinal = XPTS(NXPTS)
%
% NXPTS = Number of output points
%
% A(NROWA,NCOMP) = Boundary condition matrix at Xinitial,
% must be contained in (NIC,NCOMP) sub-matrix.
%
% NROWA = Actual row dimension of A in calling program,
% NROWA must be .GE. NIC.
%
% ALPHA(NIC+NEQIVP) = Boundary conditions at Xinitial.
% If NEQIVP .GT. 0 (see below), the boundary
% conditions at Xinitial for the initial value
% equations must be stored starting in
% position (NIC + 1) of ALPHA.
% Thus, ALPHA(NIC+K) = ETA(K).
%
% NIC = Number of boundary conditions at Xinitial.
%
% B(NROWB,NCOMP) = Boundary condition matrix at Xfinal,
% must be contained in (NFC,NCOMP) sub-matrix.
%
% NROWB = Actual row dimension of B in calling program,
% NROWB must be .GE. NFC.
%
% BETA(NFC) = Boundary conditions at Xfinal.
%
% NFC = Number of boundary conditions at Xfinal
%
% IGOFX =0 -- The inhomogeneous term G(X) is identically zero.
% =1 -- The inhomogeneous term G(X) is not identically zero.
% (if IGOFX=1, then subroutine GVEC (or UVEC) must be
% supplied).
%
% RE = Relative error tolerance used by the integrator
% (see one of the integrators)
%
% AE = Absolute error tolerance used by the integrator
% (see one of the integrators)
% **NOTE- RE and AE should not both be zero.
%
% IFLAG = A status parameter used principally for output.
% However, for efficient solution of problems which
% are originally defined as complex valued (but
% converted to real systems to use this code), the
% user must set IFLAG=13 on input. See the comment below
% for more information on solving such problems.
%
% WORK(NDW) = Floating point array used for internal storage.
%
% NDW = Actual dimension of WORK array allocated by user.
% An estimate for NDW can be computed from the following
% NDW = 130 + NCOMP**2 * (6 + NXPTS/2 + expected number of
% orthonormalizations/8)
% For the DISK or TAPE storage mode,
% NDW = 6 * NCOMP**2 + 10 * NCOMP + 130
% However, when the ADAMS integrator is to be used, the estimates are
% NDW = 130 + NCOMP**2 * (13 + NXPTS/2 + expected number of
% orthonormalizations/8)
% and NDW = 13 * NCOMP**2 + 22 * NCOMP + 130 , respectively.
%
% IWORK(NDIW) = Integer array used for internal storage.
%
% NDIW = Actual dimension of IWORK array allocated by user.
% An estimate for NDIW can be computed from the following
% NDIW = 68 + NCOMP * (1 + expected number of
% orthonormalizations)
% **NOTE -- The amount of storage required is problem dependent and may
% be difficult to predict in advance. Experience has shown
% that for most problems 20 or fewer orthonormalizations
% should suffice. If the problem cannot be completed with the
% allotted storage, then a message will be printed which
% estimates the amount of storage necessary. In any case, the
% user can examine the IWORK array for the actual storage
% requirements, as described in the output information below.
%
% NEQIVP = Number of auxiliary initial value equations being added
% to the boundary value problem.
% **NOTE -- Occasionally the coefficients MATRIX and/or G may be
% functions which depend on the independent variable X and
% on U, the solution of an auxiliary initial value problem.
% In order to avoid the difficulties associated with
% interpolation, the auxiliary equations may be solved
% simultaneously with the given boundary value problem.
% This initial value problem may be LINEAR or NONLINEAR.
% See SAND77-1328 for an example.
%
%
% The user must supply subroutines FMAT, GVEC, UIVP and UVEC, when
% needed (they MUST be so named), to evaluate the derivatives
% as follows
%
% A. FMAT must be supplied.
%
% subroutine FMAT(X,Y,YP)
% X = Independent variable (input to FMAT)
% Y = Dependent variable vector (input to FMAT)
% YP = dY/dX = Derivative vector (output from FMAT)
%
% Compute the derivatives for the HOMOGENEOUS problem
% YP(I) = dY(I)/dX = MATRIX(X) * Y(I) , I = 1,...,NCOMP
%
% When (NEQIVP .GT. 0) and MATRIX is dependent on U as
% well as on X, the following common statement must be
% included in FMAT
% COMMON /MLIVP/ NOFST
% For convenience, the U vector is stored at the bottom
% of the Y array. Thus, during any call to FMAT,
% U(I) is referenced by Y(NOFST + I).
%
%
% subroutine BVDER calls FMAT NFC times to evaluate the
% homogeneous equations and, if necessary, it calls FMAT once
% in evaluating the particular solution. Since X remains
% unchanged in this sequence of calls it is possible to
% realize considerable computational savings for complicated
% and expensive evaluations of the MATRIX entries. To do this
% the user merely passes a variable, say XS, via COMMON where
% XS is defined in the main program to be any value except
% the initial X. Then the non-constant elements of MATRIX(X)
% appearing in the differential equations need only be
% computed if X is unequal to XS, whereupon XS is reset to X.
%
%
% B. If NEQIVP .GT. 0 , UIVP must also be supplied.
%
% subroutine UIVP(X,U,UP)
% X = Independent variable (input to UIVP)
% U = Dependent variable vector (input to UIVP)
% UP = dU/dX = Derivative vector (output from UIVP)
%
% Compute the derivatives for the auxiliary initial value eqs
% UP(I) = dU(I)/dX, I = 1,...,NEQIVP.
%
% subroutine BVDER calls UIVP once to evaluate the
% derivatives for the auxiliary initial value equations.
%
%
% C. If NEQIVP = 0 and IGOFX = 1 , GVEC must be supplied.
%
% subroutine GVEC(X,G)
% X = Independent variable (input to GVEC)
% G = Vector of inhomogeneous terms G(X) (output from GVEC)
%
% Compute the inhomogeneous terms G(X)
% G(I) = G(X) values for I = 1,...,NCOMP.
%
% subroutine BVDER calls GVEC in evaluating the particular
% solution provided G(X) is NOT identically zero. Thus, when
% IGOFX=0, the user need NOT write a GVEC subroutine. Also,
% the user does not have to bother with the computational
% savings scheme for GVEC as this is automatically achieved
% via the BVDER subroutine.
%
%
% D. If NEQIVP .GT. 0 and IGOFX = 1 , UVEC must be supplied.
%
% subroutine UVEC(X,U,G)
% X = Independent variable (input to UVEC)
% U = Dependent variable vector from the auxiliary initial
% value problem (input to UVEC)
% G = Array of inhomogeneous terms G(X,U,output from UVEC)
%
% Compute the inhomogeneous terms G(X,U)
% G(I) = G(X,U) values for I = 1,...,NCOMP.
%
% subroutine BVDER calls UVEC in evaluating the particular
% solution provided G(X,U) is NOT identically zero. Thus,
% when IGOFX=0, the user need NOT write a UVEC subroutine.
%
%
%
% The following is optional input to BVSUP to give the user more
% flexibility in use of the code. See SAND75-0198 , SAND77-1328 ,
% SAND77-1690,SAND78-0522, and SAND78-1501 for more information.
%
% ****CAUTION -- The user MUST zero out IWORK(1),...,IWORK(15)
% prior to calling BVSUP. These locations define optional
% input and MUST be zero UNLESS set to special values by
% the user as described below.
%
% IWORK(1) -- Number of orthonormalization points.
% A value need be set only if IWORK(11) = 1
%
% IWORK(9) -- Integrator and orthonormalization parameter
% (default value is 1)
% 1 = RUNGE-KUTTA-FEHLBERG code using GRAM-SCHMIDT test.
% 2 = ADAMS code using GRAM-SCHMIDT TEST.
%
% IWORK(11) -- Orthonormalization points parameter
% (default value is 0)
% 0 - Orthonormalization points not pre-assigned.
% 1 - Orthonormalization points pre-assigned in
% the first IWORK(1) positions of WORK.
%
% IWORK(12) -- Storage parameter
% (default value is 0)
% 0 - All storage IN CORE
% LUN - Homogeneous and inhomogeneous solutions at
% output points and orthonormalization information
% are stored on DISK. The logical unit number to be
% used for DISK I/O (NTAPE) is set to IWORK(12).
%
% WORK(1),... -- Pre-assigned orthonormalization points, stored
% monotonically, corresponding to the direction
% of integration.
%
%
%
% ******************************
% *** COMPLEX VALUED PROBLEM ***
% ******************************
% **NOTE***
% Suppose the original boundary value problem is NC equations
% of the form
% dW/dX = MAT(X,U)*W(X) + H(X,U)
% R*W(Xinitial)=GAMMA , S*W(Xfinal)=DELTA
%
% where all variables are complex valued. The BVSUP code can be
% used by converting to a real system of size 2*NC. To solve the
% larger dimensioned problem efficiently, the user must initialize
% IFLAG=13 on input and order the vector components according to
% Y(1)=real(W(1)),...,Y(NC)=real(W(NC)),Y(NC+1)=imag(W(1)),....,
% Y(2*NC)=imag(W(NC)). Then define
% ...........................
% . real(MAT) -imag(MAT) .
% MATRIX = . .
% . imag(MAT) real(MAT) .
% ...........................
%
% The matrices A,B and vectors G,ALPHA,BETA must be defined
% similarly. Further details can be found in SAND78-1501.
%
%
% **********************************************************************
% OUTPUT from BVSUP
% **********************************************************************
%
% Y(NROWY,NXPTS) = Solution at specified output points.
%
% IFLAG output values
% =-5 Algorithm ,for obtaining starting vectors for the
% special complex problem structure, was unable to obtain
% the initial vectors satisfying the necessary
% independence criteria.
% =-4 Rank of boundary condition matrix A is less than NIC,
% as determined by LSSUDS.
% =-2 Invalid input parameters.
% =-1 Insufficient number of storage locations allocated for
% WORK or IWORK.
%
% =0 Indicates successful solution
%
% =1 A computed solution is returned but UNIQUENESS of the
% solution of the boundary-value problem is questionable.
% For an eigenvalue problem, this should be treated as a
% successful execution since this is the expected mode
% of return.
% =2 A computed solution is returned but the EXISTENCE of the
% solution to the boundary-value problem is questionable.
% =3 A nontrivial solution approximation is returned although
% the boundary condition matrix B*Y(Xfinal) is found to be
% nonsingular (to the desired accuracy level) while the
% right hand side vector is zero. To eliminate this type
% of return, the accuracy of the eigenvalue parameter
% must be improved.
% ***NOTE- We attempt to diagnose the correct problem behavior
% and report possible difficulties by the appropriate
% error flag. However, the user should probably resolve
% the problem using smaller error tolerances and/or
% perturbations in the boundary conditions or other
% parameters. This will often reveal the correct
% interpretation for the problem posed.
%
% =13 Maximum number of orthonormalizations attained before
% reaching Xfinal.
% =20-flag from integrator (DERKF or DEABM) values can range
% from 21 to 25.
% =30 Solution vectors form a dependent set.
%
% WORK(1),...,WORK(IWORK(1)) = Orthonormalization points
% determined by BVPOR.
%
% IWORK(1) = Number of orthonormalizations performed by BVPOR.
%
% IWORK(2) = Maximum number of orthonormalizations allowed as
% calculated from storage allocated by user.
%
% IWORK(3),IWORK(4),IWORK(5),IWORK(6) Give information about
% actual storage requirements for WORK and IWORK
% arrays. In particular,
% required storage for WORK array is
% IWORK(3) + IWORK(4)*(expected number of orthonormalizations)
%
% required storage for IWORK array is
% IWORK(5) + IWORK(6)*(expected number of orthonormalizations)
%
% IWORK(8) = Final value of exponent parameter used in tolerance
% test for orthonormalization.
%
% IWORK(16) = Number of independent vectors returned from MGSBV.
% It is only of interest when IFLAG=30 is obtained.
%
% IWORK(17) = Numerically estimated rank of the boundary
% condition matrix defined from B*Y(Xfinal)
%
% **********************************************************************
%
% Necessary machine constants are defined in the function
% routine R1MACH. The user must make sure that the values
% set in R1MACH are relevant to the computer being used.
%
% **********************************************************************
%
%***REFERENCES M. R. Scott and H. A. Watts, SUPORT - a computer code
% for two-point boundary-value problems via
% orthonormalization, SIAM Journal of Numerical
% Analysis 14, (1977), pp. 40-70.
% B. L. Darlow, M. R. Scott and H. A. Watts, Modifications
% of SUPORT, a linear boundary value problem solver
% Part I - pre-assigning orthonormalization points,
% auxiliary initial value problem, disk or tape storage,
% Report SAND77-1328, Sandia Laboratories, Albuquerque,
% New Mexico, 1977.
% B. L. Darlow, M. R. Scott and H. A. Watts, Modifications
% of SUPORT, a linear boundary value problem solver
% Part II - inclusion of an Adams integrator, Report
% SAND77-1690, Sandia Laboratories, Albuquerque,
% New Mexico, 1977.
% M. E. Lord and H. A. Watts, Modifications of SUPORT,
% a linear boundary value problem solver Part III -
% orthonormalization improvements, Report SAND78-0522,
% Sandia Laboratories, Albuquerque, New Mexico, 1978.
% H. A. Watts, M. R. Scott and M. E. Lord, Computational
% solution of complex*16 valued boundary problems,
% Report SAND78-1501, Sandia Laboratories,
% Albuquerque, New Mexico, 1978.
%***ROUTINES CALLED EXBVP, MACON, XERMSG
%***COMMON BLOCKS ML15TO, ML17BW, ML18JR, ML5MCO, ML8SZ
%***REVISION HISTORY (YYMMDD)
% 750601 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 890831 Modified array declarations. (WRB)
% 890921 Realigned order of variables in certain COMMON blocks.
% (WRB)
% 890921 REVISION DATE from Version 3.2
% 891214 Prologue converted to Version 4.0 format. (BAB)
% 900510 Convert XERRWV calls to XERMSG calls. (RWC)
% 920501 Reformatted the REFERENCES section. (WRB)
%***end PROLOGUE BVSUP
% **********************************************************************
%
%
global ml18jr_9; if isempty(ml18jr_9), ml18jr_9=0; end;
global ml18jr_14; if isempty(ml18jr_14), ml18jr_14=0; end;
global ml18jr_16; if isempty(ml18jr_16), ml18jr_16=0; end;
global ml18jr_4; if isempty(ml18jr_4), ml18jr_4=0; end;
if isempty(nxptsm), nxptsm=0; end;
y_shape=size(y);y=reshape([y(:).',zeros(1,ceil(numel(y)./prod([nrowy])).*prod([nrowy])-numel(y))],nrowy,[]);
a_shape=size(a);a=reshape([a(:).',zeros(1,ceil(numel(a)./prod([nrowa])).*prod([nrowa])-numel(a))],nrowa,[]);
alpha_shape=size(alpha);alpha=reshape(alpha,1,[]);
b_shape=size(b);b=reshape([b(:).',zeros(1,ceil(numel(b)./prod([nrowb])).*prod([nrowb])-numel(b))],nrowb,[]);
beta_shape=size(beta);beta=reshape(beta,1,[]);
work_shape=size(work);work=reshape(work,1,[]);
iwork_shape=size(iwork);iwork=reshape(iwork,1,[]);
xpts_shape=size(xpts);xpts=reshape(xpts,1,[]);
if isempty(xern1), xern1=repmat(' ',1,8); end;
if isempty(xern2), xern2=repmat(' ',1,8); end;
if isempty(xern3), xern3=repmat(' ',1,8); end;
if isempty(xern4), xern4=repmat(' ',1,8); end;
%
% **********************************************************************
% THE COMMON BLOCK BELOW IS USED TO COMMUNICATE WITH SUBROUTINE
% BVDER. THE USER SHOULD NOT ALTER OR USE THIS COMMON BLOCK IN THE
% CALLING PROGRAM.
%
% common :: ;
%% common /ml8sz / c , xsav , igofxd , inhomo , ivp , ncompd , nfcd;
%% common /ml8sz / ml8sz_1 , ml8sz_2 , ml8sz_3 , ml8sz_4 , ml8sz_5 , ml8sz_6 , ml8sz_7;
%
% **********************************************************************
% THESE COMMON BLOCKS AID IN REDUCING THE NUMBER OF SUBROUTINE
% ARGUMENTS PREVALENT IN THIS MODULAR STRUCTURE
%
% common :: ;
%% common /ml18jr/ aed , red , tol , nxptsd , nicd , nopg , mxnon ,ndisk , ntape , neq , indpvt , integ , nps , ntp ,neqivd , numort , nfcc , icoco;
%% common /ml18jr/ ml18jr_1 , ml18jr_2 , ml18jr_3 , ml18jr_4 , ml18jr_5 , ml18jr_6 , ml18jr_7 ,ml18jr_8 , ml18jr_9 , ml18jr_10 , ml18jr_11 , ml18jr_12 , ml18jr_13 , ml18jr_14 ,ml18jr_15 , ml18jr_16 , ml18jr_17 , ml18jr_18;
% common :: ;
%% common /ml17bw/ kkkzpw , needw , neediw , k1 , k2 , k3 , k4 , k5 ,k6 , k7 , k8 , k9 , k10 , k11 , l1 , l2 , kkkint ,lllint;
%% common /ml17bw/ ml17bw_1 , ml17bw_2 , ml17bw_3 , ml17bw_4 , ml17bw_5 , ml17bw_6 , ml17bw_7 , ml17bw_8 ,ml17bw_9 , ml17bw_10 , ml17bw_11 , ml17bw_12 , ml17bw_13 , ml17bw_14 , ml17bw_15 , ml17bw_16 , ml17bw_17 ,ml17bw_18;
%
% **********************************************************************
% THIS COMMON BLOCK IS USED IN SUBROUTINES BVSUP,BVPOR,RKFAB,
% REORT, AND STWAY. IT CONTAINS INFORMATION NECESSARY
% FOR THE ORTHONORMALIZATION TESTING PROCEDURE AND A BACKUP
% RESTARTING CAPABILITY.
%
% common :: ;
%% common /ml15to/ px , pwcnd , tnd , x , xbeg , xend , xot , xop ,info(15) , istkop , knswot , kop , lotjp ,mnswot , nswot;
%% common /ml15to/ ml15to_1 , ml15to_2 , ml15to_3 , ml15to_4 , ml15to_5 , ml15to_6 , ml15to_7 , ml15to_8 ,ml15to_9(15) , ml15to_10 , ml15to_11 , ml15to_12 , ml15to_13 ,ml15to_14 , ml15to_15;
%
% **********************************************************************
% THIS COMMON BLOCK CONTAINS THE MACHINE DEPENDENT PARAMETERS
% USED BY THE CODE
%
% common :: ;
%% common /ml5mco/ uro , sru , eps , sqovfl , twou , fouru , lpar;
%% common /ml5mco/ ml5mco_1 , ml5mco_2 , ml5mco_3 , ml5mco_4 , ml5mco_5 , ml5mco_6 , ml5mco_7;
%
% **********************************************************************
% SET UP MACHINE DEPENDENT CONSTANTS.
%
%***FIRST EXECUTABLE STATEMENT BVSUP
macon;
%
% **********************************************************************
% TEST FOR INVALID INPUT
%
while (1);
gt200=0;
gt300=0;
if( nrowy>=ncomp )
if( ncomp==nic+nfc )
if( nxpts>=2 )
if( nic>0 )
if( nrowa>=nic )
if( nfc>0 )
if( nrowb>=nfc )
if( igofx>=0 && igofx<=1 )
if( re>=0.0 )
if( ae>=0.0 )
if( re~=0.0 || ae~=0.0 )
is = 1;
if( xpts(nxpts)<xpts(1) )
is = 2;
end;
nxptsm = fix(nxpts - 1);
for k = 1 : nxptsm;
if( is==2 )
if( xpts(k)<=xpts(k+1) )
iflag = -2;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
alpha_shape=zeros(alpha_shape);alpha_shape(:)=alpha(1:numel(alpha_shape));alpha=alpha_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
beta_shape=zeros(beta_shape);beta_shape(:)=beta(1:numel(beta_shape));beta=beta_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
xpts_shape=zeros(xpts_shape);xpts_shape(:)=xpts(1:numel(xpts_shape));xpts=xpts_shape;
return;
end;
elseif( xpts(k+1)<=xpts(k) ) ;
iflag = -2;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
alpha_shape=zeros(alpha_shape);alpha_shape(:)=alpha(1:numel(alpha_shape));alpha=alpha_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
beta_shape=zeros(beta_shape);beta_shape(:)=beta(1:numel(beta_shape));beta=beta_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
xpts_shape=zeros(xpts_shape);xpts_shape(:)=xpts(1:numel(xpts_shape));xpts=xpts_shape;
return;
end;
end; k = fix(nxptsm+1);
%
% **********************************************************************
% CHECK FOR DISK STORAGE
%
kpts = fix(nxpts);
ml18jr_8 = 0;
if( iwork(12)~=0 )
ml18jr_9 = fix(iwork(12));
kpts = 1;
ml18jr_8 = 1;
end;
%
% **********************************************************************
% SET INTEG PARAMETER ACCORDING TO CHOICE OF INTEGRATOR.
%
ml18jr_12 = 1;
if( iwork(9)==2 )
ml18jr_12 = 2;
end;
%
% **********************************************************************
% COMPUTE INHOMO
%
while (1);
if( igofx==1 )
ml8sz_4 = 1;
else;
for j = 1 : nic;
if( alpha(j)~=0.0 )
ml8sz_4 = 1;
%2
break;
end;
end;
for j = 1 : nfc;
if( beta(j)~=0.0 )
ml8sz_4 = 2;
%2
break;
end;
end;
ml8sz_4 = 3;
end;
break;
end;
%
% **********************************************************************
% TO TAKE ADVANTAGE OF THE SPECIAL STRUCTURE WHEN SOLVING A
% COMPLEX VALUED PROBLEM,WE INTRODUCE NFCC=NFC WHILE CHANGING
% THE INTERNAL VALUE OF NFC
%
ml18jr_17 = fix(nfc);
if( iflag==13 )
nfc = fix(fix(nfc./2));
end;
%
% **********************************************************************
% DETERMINE NECESSARY STORAGE REQUIREMENTS
%
% FOR BASIC ARRAYS IN BVPOR
kkkyhp = fix(ncomp.*(nfc+1) + neqivp);
kkku = fix(ncomp.*nfc.*kpts);
kkkv = fix(ncomp.*kpts);
kkkcoe = fix(ml18jr_17);
kkks = fix(nfc + 1);
kkksto = fix(ncomp.*(nfc+1) + neqivp + 1);
kkkg = fix(ncomp);
%
% FOR ORTHONORMALIZATION RELATED MATTERS
ml18jr_14 =fix(fix((ml18jr_17.*(ml18jr_17+1))./2));
ml17bw_1 = fix(1 + ml18jr_14 + ml18jr_17);
lllip = fix(ml18jr_17);
%
% FOR ADDITIONAL REQUIRED WORK SPACE
% (LSSUDS)
kkksud = fix(4.*nic +(nrowa+1).*ncomp);
lllsud = fix(nic);
% (SVECS)
kkksvc = fix(1 + 4.*ml18jr_17 + 2.*ml18jr_17.^2);
lllsvc = fix(2.*ml18jr_17);
%
ndeq = fix(ncomp.*nfc + neqivp);
if( ml8sz_4==1 )
ndeq = fix(ndeq + ncomp);
end;
if( ml18jr_12==2 )
% (DEABM)
ml17bw_17 = fix(130 + 21.*ndeq);
ml17bw_18 = 51;
else;
% (DERKF)
ml17bw_17 = fix(33 + 7.*ndeq);
ml17bw_18 = 34;
end;
%
% (COEF)
kkkcof = fix(5.*ml18jr_17 + ml18jr_17.^2);
lllcof = fix(3 + ml18jr_17);
%
kkkws = fix(max([kkksud,kkksvc,ml17bw_17,kkkcof]));
llliws = fix(max([lllsud,lllsvc,ml17bw_18,lllcof]));
%
ml17bw_2 = fix(kkkyhp + kkku + kkkv + kkkcoe + kkks + kkksto +kkkg + ml17bw_1 + kkkws);
ml17bw_3 = fix(17 + lllip + llliws);
% **********************************************************************
% COMPUTE THE NUMBER OF POSSIBLE ORTHONORMALIZATIONS WITH THE
% ALLOTTED STORAGE
%
iwork(3) = fix(ml17bw_2);
iwork(4) = fix(ml17bw_1);
iwork(5) = fix(ml17bw_3);
iwork(6) = fix(lllip);
nrtemp = fix(ndw - ml17bw_2);
nitemp = fix(ndiw - ml17bw_3);
if( nrtemp<0 )
gt200=1;
break;
end;
if( nitemp<0 )
gt200=1;
break;
end;
%
if( ml18jr_8==0 )
%
mxnonr = fix(fix(nrtemp./ml17bw_1));
mxnoni = fix(fix(nitemp./lllip));
ml18jr_7 = fix(min(mxnonr,mxnoni));
non = fix(ml18jr_7);
else;
non = 0;
ml18jr_7 = fix(nrtemp);
end;
%
iwork(2) = fix(ml18jr_7);
%
% **********************************************************************
% CHECK FOR PRE-ASSIGNED ORTHONORMALIZATION POINTS
%
ml18jr_6 = 0;
if( iwork(11)~=1 )
gt300=1;
break;
end;
if( ml18jr_7<iwork(1) )
gt200=1;
break;
end;
ml18jr_6 = 1;
ml18jr_7 = fix(iwork(1));
work(ml18jr_7+1) = 2..*xpts(nxpts) - xpts(1);
gt300=1;
break;
end;
end;
end;
end;
end;
end;
end;
end;
end;
end;
end;
break;
end;
if(gt300==0)
if(gt200==0)
iflag = -2;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
alpha_shape=zeros(alpha_shape);alpha_shape(:)=alpha(1:numel(alpha_shape));alpha=alpha_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
beta_shape=zeros(beta_shape);beta_shape(:)=beta(1:numel(beta_shape));beta=beta_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
xpts_shape=zeros(xpts_shape);xpts_shape(:)=xpts(1:numel(xpts_shape));xpts=xpts_shape;
return;
%
end;
iflag = -1;
if( ml18jr_8~=1 )
xern1=sprintf(['%8i'], ml17bw_2);
xern2=sprintf(['%8i'], ml17bw_1);
xern3=sprintf(['%8i'], ml17bw_3);
xern4=sprintf(['%8i'], lllip);
xermsg('SLATEC','BVSUP',['REQUIRED STORAGE FOR WORK ARRAY IS ',[xern1,[' + ',[xern2,['*(EXPECTED NUMBER OF ORTHONORMALIZATIONS) $$',['REQUIRED STORAGE FOR IWORK ARRAY IS ',[xern3,[' + ',[xern4,'*(EXPECTED NUMBER OF ORTHONORMALIZATIONS)']]]]]]]]],1,0);
else;
xern1=sprintf(['%8i'], ml17bw_2);
xern2=sprintf(['%8i'], ml17bw_3);
xermsg('SLATEC','BVSUP',['REQUIRED STORAGE FOR WORK ARRAY IS ',[xern1,[' + NUMBER OF ORTHONOMALIZATIONS. $$',['REQUIRED STORAGE FOR IWORK ARRAY IS ',xern2]]]],1,0);
end;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
alpha_shape=zeros(alpha_shape);alpha_shape(:)=alpha(1:numel(alpha_shape));alpha=alpha_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
beta_shape=zeros(beta_shape);beta_shape(:)=beta(1:numel(beta_shape));beta=beta_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
xpts_shape=zeros(xpts_shape);xpts_shape(:)=xpts(1:numel(xpts_shape));xpts=xpts_shape;
return;
%
% **********************************************************************
% ALLOCATE STORAGE FROM WORK AND IWORK ARRAYS
%
% (Z)
%gt300
end;
ml17bw_4 = fix(1 +(ml18jr_7+1));
% (P)
ml17bw_5 = fix(ml17bw_4 + ml18jr_14.*(non+1));
% (W)
ml17bw_6 = fix(ml17bw_5 + ml18jr_17.*(non+1));
% (YHP)
ml17bw_7 = fix(ml17bw_6 + kkkyhp);
% (U)
ml17bw_8 = fix(ml17bw_7 + kkku);
% (V)
ml17bw_9 = fix(ml17bw_8 + kkkv);
% (COEF)
ml17bw_10 = fix(ml17bw_9 + kkkcoe);
% (S)
ml17bw_11 = fix(ml17bw_10 + kkks);
% (STOWA)
ml17bw_12 = fix(ml17bw_11 + kkksto);
% (G)
ml17bw_13 = fix(ml17bw_12 + kkkg);
ml17bw_14 = fix(ml17bw_13 + kkkws);
% REQUIRED ADDITIONAL REAL WORK SPACE STARTS AT WORK(K10)
% AND EXTENDS TO WORK(K11-1)
%
% FIRST 17 LOCATIONS OF IWORK ARE USED FOR OPTIONAL
% INPUT AND OUTPUT ITEMS
% (IP)
ml17bw_15 = fix(18 + ml18jr_17.*(non+1));
ml17bw_16 = fix(ml17bw_15 + llliws);
% REQUIRED INTEGER WORK SPACE STARTS AT IWORK(L1)
% AND EXTENDS TO IWORK(L2-1)
%
% **********************************************************************
% SET INDICATOR FOR NORMALIZATION OF PARTICULAR SOLUTION
%
ml18jr_13 = 0;
if( iwork(10)==1 )
ml18jr_13 = 1;
end;
%
% **********************************************************************
% SET PIVOTING PARAMETER
%
ml18jr_11 = 0;
if( iwork(15)==1 )
ml18jr_11 = 1;
end;
%
% **********************************************************************
% SET OTHER COMMON BLOCK PARAMETERS
%
ml8sz_7 = fix(nfc);
ml8sz_6 = fix(ncomp);
ml8sz_3 = fix(igofx);
ml18jr_4 = fix(nxpts);
ml18jr_5 = fix(nic);
ml18jr_2 = re;
ml18jr_1 = ae;
ml18jr_15 = fix(neqivp);
ml15to_14 = 20;
if( iwork(13)==-1 )
ml15to_14 = fix(max(1,iwork(14)));
end;
ml15to_5 = xpts(1);
ml15to_6 = xpts(nxpts);
ml8sz_2 = ml15to_6;
ml18jr_18 = 1;
if( ml8sz_4==3 && ml18jr_6==1 )
work(ml18jr_7+1) = ml15to_6;
end;
%
% **********************************************************************
%
[y,nrowy,xpts,a,nrowa,alpha,b,nrowb,beta,iflag,work,iwork]=exbvp(y,nrowy,xpts,a,nrowa,alpha,b,nrowb,beta,iflag,work,iwork);
nfc = fix(ml18jr_17);
iwork(17) = fix(iwork(ml17bw_15));
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
alpha_shape=zeros(alpha_shape);alpha_shape(:)=alpha(1:numel(alpha_shape));alpha=alpha_shape;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
beta_shape=zeros(beta_shape);beta_shape(:)=beta(1:numel(beta_shape));beta=beta_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_shape;
xpts_shape=zeros(xpts_shape);xpts_shape(:)=xpts(1:numel(xpts_shape));xpts=xpts_shape;
end %subroutine bvsup
%DECK C0LGMC
|
|
