| [y,nrowy,ncomp,xpts,nxpts,a,nrowa,alpha,nic,b,nrowb,beta,nfc,igofx,re,ae,iflag,work,ndw,iwork,ndiw,neqivp]=dbvsup(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]=dbvsup(y,nrowy,ncomp,xpts,nxpts,a,nrowa,alpha,nic,b,nrowb,beta,nfc,igofx,re,ae,iflag,work,ndw,iwork,ndiw,neqivp);
%***BEGIN PROLOGUE DBVSUP
%***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 doubleprecision (BVSUP-S, DBVSUP-D)
%***KEYWORDS ORTHONORMALIZATION, SHOOTING,
% TWO-POINT BOUNDARY VALUE PROBLEM
%***AUTHOR Scott, M. R., (SNLA)
% Watts, H. A., (SNLA)
%***DESCRIPTION
%
% **********************************************************************
%
% subroutine DBVSUP 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 DBVSUP
% **********************************************************************
%
% 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 DGVEC (or DUVEC) 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*16 valued (but
% converted to doubleprecision 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 DFMAT, DGVEC, DUIVP and DUVEC,
% when needed (they must be so named), to evaluate the derivatives
% as follows
%
% A. DFMAT must be supplied.
%
% subroutine DFMAT(X,Y,YP)
% X = independent variable (input to DFMAT)
% Y = dependent variable vector (input to DFMAT)
% YP = DY/DX = derivative vector (output from DFMAT)
%
% 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 DFMAT
% COMMON /DMLIVP/ NOFST
% for convenience, the U vector is stored at the bottom
% of the Y array. Thus, during any call to DFMAT,
% U(I) is referenced by Y(NOFST + I).
%
%
% subroutine DBVDER calls DFMAT NFC times to evaluate the
% homogeneous equations and, if necessary, it calls DFMAT
% 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(dml15t_4)
% 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 , DUIVP must also be supplied.
%
% subroutine DUIVP(X,U,UP)
% X = independent variable (input to DUIVP)
% U = dependent variable vector (input to DUIVP)
% UP = DU/DX = derivative vector (output from DUIVP)
%
% Compute the derivatives for the auxiliary initial value eqs
% UP(I) = DU(I)/DX, I = 1,...,NEQIVP.
%
% subroutine DBVDER calls DUIVP once to evaluate the
% derivatives for the auxiliary initial value equations.
%
%
% C. If NEQIVP = 0 and IGOFX = 1 , DGVEC must be supplied.
%
% subroutine DGVEC(X,G)
% X = independent variable (input to DGVEC)
% G = vector of inhomogeneous terms G(X) (output from
% DGVEC)
%
% Compute the inhomogeneous terms G(X)
% G(I) = G(X) values for I = 1,...,NCOMP.
%
% subroutine DBVDER calls DGVEC in evaluating the particular
% solution provided G(X) is not identically zero. Thus, when
% IGOFX=0, the user need not write a DGVEC subroutine. Also,
% the user does not have to bother with the computational
% savings scheme for DGVEC as this is automatically achieved
% via the DBVDER subroutine.
%
%
% D. If NEQIVP .GT. 0 and IGOFX = 1 , DUVEC must be supplied.
%
% subroutine DUVEC(X,U,G)
% X = independent variable (input to DUVEC)
% U = dependent variable vector from the auxiliary initial
% value problem (input to DUVEC)
% G = array of inhomogeneous terms G(X,U,output from DUVEC)
%
% Compute the inhomogeneous terms G(X,U)
% G(I) = G(X,U) values for I = 1,...,NCOMP.
%
% subroutine DBVDER calls DUVEC in evaluating the particular
% solution provided G(X,U) is not identically zero. Thus,
% when IGOFX=0, the user need not write a DUVEC subroutine.
%
%
%
% The following is optional input to DBVSUP to give user more
% flexibility in use of 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 DBVSUP. 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*16 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*16 valued. The DBVSUP code can be
% used by converting to a doubleprecision 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)=doubleprecision(W(1)),...,Y(NC)=DOUBLE
% PRECISION(W(NC)),Y(NC+1)=IMAG(W(1)),...., Y(2*NC)=IMAG(W(NC)).
% Then define
% ...............................................
% . doubleprecision(MAT) -IMAG(MAT) .
% MATRIX = . .
% . IMAG(MAT) doubleprecision(MAT) .
% ...............................................
%
% The matrices A,B and vectors G,ALPHA,BETA must be defined
% similarly. Further details can be found in SAND78-1501.
%
%
% **********************************************************************
% OUTPUT from DBVSUP
% **********************************************************************
%
% Y(NROWY,NXPTS) = solution at specified output points.
%
% IFLAG Output Values
% =-5 algorithm ,for obtaining starting vectors for the
% special COMPLEX*16 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 DLSSUD.
% =-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 (DDERKF or DDEABM) values can
% range from 21 to 25.
% =30 solution vectors form a dependent set.
%
% WORK(1),...,WORK(IWORK(1)) = orthonormalization points
% determined by DBVPOR.
%
% IWORK(1) = number of orthonormalizations performed by DBVPOR.
%
% 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 DMGSBV.
% 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 D1MACH. The user must make sure that the values
% set in D1MACH 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 DEXBVP, DMACON, XERMSG
%***COMMON BLOCKS DML15T, DML17B, DML18J, DML5MC, DML8SZ
%***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, remove some extraneous
% comments. (RWC)
% 920501 Reformatted the REFERENCES section. (WRB)
%***end PROLOGUE DBVSUP
% **********************************************************************
%
persistent gt 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 dml18j_18; if isempty(dml18j_18), dml18j_18=0; end;
global dml8sz_3; if isempty(dml8sz_3), dml8sz_3=0; end;
global dml18j_11; if isempty(dml18j_11), dml18j_11=0; end;
global dml15t_9; if isempty(dml15t_9), dml15t_9=zeros(1,15); end;
global dml8sz_4; if isempty(dml8sz_4), dml8sz_4=0; end;
global dml18j_12; if isempty(dml18j_12), dml18j_12=0; end;
if isempty(is), is=0; end;
global dml15t_10; if isempty(dml15t_10), dml15t_10=0; end;
global dml8sz_5; if isempty(dml8sz_5), dml8sz_5=0; end;
iwork_shape=size(iwork);iwork=reshape(iwork,1,[]);
if isempty(j), j=0; end;
if isempty(k), k=0; end;
global dml17b_4; if isempty(dml17b_4), dml17b_4=0; end;
global dml17b_13; if isempty(dml17b_13), dml17b_13=0; end;
global dml17b_14; if isempty(dml17b_14), dml17b_14=0; end;
global dml17b_5; if isempty(dml17b_5), dml17b_5=0; end;
global dml17b_6; if isempty(dml17b_6), dml17b_6=0; end;
global dml17b_7; if isempty(dml17b_7), dml17b_7=0; end;
global dml17b_8; if isempty(dml17b_8), dml17b_8=0; end;
global dml17b_9; if isempty(dml17b_9), dml17b_9=0; end;
global dml17b_10; if isempty(dml17b_10), dml17b_10=0; end;
global dml17b_11; if isempty(dml17b_11), dml17b_11=0; end;
global dml17b_12; if isempty(dml17b_12), dml17b_12=0; end;
if isempty(kkkcoe), kkkcoe=0; end;
if isempty(kkkcof), kkkcof=0; end;
if isempty(kkkg), kkkg=0; end;
global dml17b_17; if isempty(dml17b_17), dml17b_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 dml17b_1; if isempty(dml17b_1), dml17b_1=0; end;
global dml15t_11; if isempty(dml15t_11), dml15t_11=0; end;
global dml15t_12; if isempty(dml15t_12), dml15t_12=0; end;
if isempty(kpts), kpts=0; end;
global dml17b_15; if isempty(dml17b_15), dml17b_15=0; end;
global dml17b_16; if isempty(dml17b_16), dml17b_16=0; end;
if isempty(lllcof), lllcof=0; end;
global dml17b_18; if isempty(dml17b_18), dml17b_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 dml15t_13; if isempty(dml15t_13), dml15t_13=0; end;
global dml5mc_7; if isempty(dml5mc_7), dml5mc_7=0; end;
global dml15t_14; if isempty(dml15t_14), dml15t_14=0; end;
global dml18j_7; if isempty(dml18j_7), dml18j_7=0; end;
if isempty(mxnoni), mxnoni=0; end;
if isempty(mxnonr), mxnonr=0; end;
global dml8sz_6; if isempty(dml8sz_6), dml8sz_6=0; end;
if isempty(ndeq), ndeq=0; end;
global dml18j_8; if isempty(dml18j_8), dml18j_8=0; end;
global dml17b_3; if isempty(dml17b_3), dml17b_3=0; end;
global dml17b_2; if isempty(dml17b_2), dml17b_2=0; end;
global dml18j_10; if isempty(dml18j_10), dml18j_10=0; end;
global dml18j_15; if isempty(dml18j_15), dml18j_15=0; end;
global dml18j_17; if isempty(dml18j_17), dml18j_17=0; end;
global dml8sz_7; if isempty(dml8sz_7), dml8sz_7=0; end;
global dml18j_5; if isempty(dml18j_5), dml18j_5=0; end;
if isempty(nitemp), nitemp=0; end;
if isempty(non), non=0; end;
global dml18j_6; if isempty(dml18j_6), dml18j_6=0; end;
global dml18j_13; if isempty(dml18j_13), dml18j_13=0; end;
if isempty(nrtemp), nrtemp=0; end;
global dml15t_15; if isempty(dml15t_15), dml15t_15=0; end;
global dml18j_9; if isempty(dml18j_9), dml18j_9=0; end;
global dml18j_14; if isempty(dml18j_14), dml18j_14=0; end;
global dml18j_16; if isempty(dml18j_16), dml18j_16=0; end;
global dml18j_4; if isempty(dml18j_4), dml18j_4=0; end;
if isempty(nxptsm), nxptsm=0; end;
if isempty(gt), gt=0; end;
a_shape=size(a);a=reshape([a(:).',zeros(1,ceil(numel(a)./prod([nrowa])).*prod([nrowa])-numel(a))],nrowa,[]);
global dml18j_1; if isempty(dml18j_1), dml18j_1=0; end;
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,[]);
global dml8sz_1; if isempty(dml8sz_1), dml8sz_1=0; end;
global dml5mc_3; if isempty(dml5mc_3), dml5mc_3=0; end;
global dml5mc_6; if isempty(dml5mc_6), dml5mc_6=0; end;
global dml15t_2; if isempty(dml15t_2), dml15t_2=0; end;
global dml15t_1; if isempty(dml15t_1), dml15t_1=0; end;
global dml18j_2; if isempty(dml18j_2), dml18j_2=0; end;
global dml5mc_4; if isempty(dml5mc_4), dml5mc_4=0; end;
global dml5mc_2; if isempty(dml5mc_2), dml5mc_2=0; end;
global dml15t_3; if isempty(dml15t_3), dml15t_3=0; end;
global dml18j_3; if isempty(dml18j_3), dml18j_3=0; end;
global dml5mc_5; if isempty(dml5mc_5), dml5mc_5=0; end;
global dml5mc_1; if isempty(dml5mc_1), dml5mc_1=0; end;
global dml15t_4; if isempty(dml15t_4), dml15t_4=0; end;
global dml15t_5; if isempty(dml15t_5), dml15t_5=0; end;
global dml15t_6; if isempty(dml15t_6), dml15t_6=0; end;
global dml15t_8; if isempty(dml15t_8), dml15t_8=0; end;
global dml15t_7; if isempty(dml15t_7), dml15t_7=0; end;
xpts_shape=size(xpts);xpts=reshape(xpts,1,[]);
global dml8sz_2; if isempty(dml8sz_2), dml8sz_2=0; end;
y_shape=size(y);y=reshape([y(:).',zeros(1,ceil(numel(y)./prod([nrowy])).*prod([nrowy])-numel(y))],nrowy,[]);
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
% DBVDER. THE USER SHOULD NOT ALTER OR USE THIS COMMON BLOCK IN
% THE CALLING PROGRAM.
%
% common :: ;
%% common /dml8sz/ c , xsav , igofxd , inhomo , ivp , ncompd , nfcd;
%% common /dml8sz/ dml8sz_1 , dml8sz_2 , dml8sz_3 , dml8sz_4 , dml8sz_5 , dml8sz_6 , dml8sz_7;
%
% ******************************************************************
% THESE COMMON BLOCKS AID IN REDUCING THE NUMBER OF SUBROUTINE
% ARGUMENTS PREVALENT IN THIS MODULAR STRUCTURE
%
% common :: ;
%% common /dml18j/ aed , red , tol , nxptsd , nicd , nopg , mxnon ,ndisk , ntape , neq , indpvt , integ , nps , ntp ,neqivd , numort , nfcc , icoco;
%% common /dml18j/ dml18j_1 , dml18j_2 , dml18j_3 , dml18j_4 , dml18j_5 , dml18j_6 , dml18j_7 ,dml18j_8 , dml18j_9 , dml18j_10 , dml18j_11 , dml18j_12 , dml18j_13 , dml18j_14 ,dml18j_15 , dml18j_16 , dml18j_17 , dml18j_18;
% common :: ;
%% common /dml17b/ kkkzpw , needw , neediw , k1 , k2 , k3 , k4 , k5 ,k6 , k7 , k8 , k9 , k10 , k11 , l1 , l2 , kkkint ,lllint;
%% common /dml17b/ dml17b_1 , dml17b_2 , dml17b_3 , dml17b_4 , dml17b_5 , dml17b_6 , dml17b_7 , dml17b_8 ,dml17b_9 , dml17b_10 , dml17b_11 , dml17b_12 , dml17b_13 , dml17b_14 , dml17b_15 , dml17b_16 , dml17b_17 ,dml17b_18;
%
% ******************************************************************
% THIS COMMON BLOCK IS USED IN SUBROUTINES DBVSUP,DBVPOR,DRKFAB,
% DREORT, AND DSTWAY. IT CONTAINS INFORMATION NECESSARY
% FOR THE ORTHONORMALIZATION TESTING PROCEDURE AND A BACKUP
% RESTARTING CAPABILITY.
%
% common :: ;
%% common /dml15t/ px , pwcnd , tnd , x , xbeg , xend , xot , xop ,info(15) , istkop , knswot , kop , lotjp ,mnswot , nswot;
%% common /dml15t/ dml15t_1 , dml15t_2 , dml15t_3 , dml15t_4 , dml15t_5 , dml15t_6 , dml15t_7 , dml15t_8 ,dml15t_9(15) , dml15t_10 , dml15t_11 , dml15t_12 , dml15t_13 ,dml15t_14 , dml15t_15;
%
% ******************************************************************
% THIS COMMON BLOCK CONTAINS THE MACHINE DEPENDENT PARAMETERS
% USED BY THE CODE
%
% common :: ;
%% common /dml5mc/ uro , sru , eps , sqovfl , twou , fouru , lpar;
%% common /dml5mc/ dml5mc_1 , dml5mc_2 , dml5mc_3 , dml5mc_4 , dml5mc_5 , dml5mc_6 , dml5mc_7;
%
% *****************************************************************
% SET UP MACHINE DEPENDENT CONSTANTS.
%
%***FIRST EXECUTABLE STATEMENT DBVSUP
dmacon;
%
% ************************************************
% TEST FOR INVALID INPUT
%
if( nrowy<ncomp )
iflag = -2;
elseif( ncomp~=nic+nfc ) ;
iflag = -2;
elseif( nxpts<2 ) ;
iflag = -2;
elseif( nic<=0 ) ;
iflag = -2;
elseif( nrowa<nic ) ;
iflag = -2;
elseif( nfc<=0 ) ;
iflag = -2;
elseif( nrowb<nfc ) ;
iflag = -2;
elseif( igofx<0 || igofx>1 ) ;
iflag = -2;
elseif( re<0.0d0 ) ;
iflag = -2;
elseif( ae<0.0d0 ) ;
iflag = -2;
else;
while (1);
if( re==0.0d0 && ae==0.0d0 )
iflag = -2;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_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;
xpts_shape=zeros(xpts_shape);xpts_shape(:)=xpts(1:numel(xpts_shape));xpts=xpts_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
else;
% BEGIN BLOCK PERMITTING ...EXITS TO 70
is = 1;
if( xpts(nxpts)<xpts(1) )
is = 2;
end;
nxptsm = fix(nxpts - 1);
for k = 1 : nxptsm;
if( is==2 )
% .........EXIT
if( xpts(k)<=xpts(k+1) )
iflag = -2;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_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;
xpts_shape=zeros(xpts_shape);xpts_shape(:)=xpts(1:numel(xpts_shape));xpts=xpts_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
% .........EXIT
elseif( xpts(k+1)<=xpts(k) ) ;
iflag = -2;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_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;
xpts_shape=zeros(xpts_shape);xpts_shape(:)=xpts(1:numel(xpts_shape));xpts=xpts_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
end; k = fix(nxptsm+1);
%
% ******************************************
% CHECK FOR DISK STORAGE
%
kpts = fix(nxpts);
dml18j_8 = 0;
if( iwork(12)~=0 )
dml18j_9 = fix(iwork(12));
kpts = 1;
dml18j_8 = 1;
end;
%
% ******************************************
% SET INTEG PARAMETER ACCORDING TO
% CHOICE OF INTEGRATOR.
%
dml18j_12 = 1;
if( iwork(9)==2 )
dml18j_12 = 2;
end;
%
% ******************************************
% COMPUTE INHOMO
%
% ............EXIT
gt=0;
while (1);
if( igofx==1 )
% ..................EXIT
dml8sz_4 = 1;
else;
for j = 1 : nic;
% ...............EXIT
if( alpha(j)~=0.0d0 )
dml8sz_4 = 1;
gt=1;
break;
end;
end;
if(gt==1)
break;
end;
for j = 1 : nfc;
% ............EXIT
if( beta(j)~=0.0d0 )
dml8sz_4 = 2;
% ......EXIT
gt=1;
break;
end;
end;
if(gt==1)
break;
end;
% ...............EXIT
dml8sz_4 = 3;
end;
break;
end;
%
% *********************************************************
% TO TAKE ADVANTAGE OF THE SPECIAL STRUCTURE WHEN
% SOLVING A COMPLEX*16 VALUED PROBLEM,WE INTRODUCE
% NFCC=NFC WHILE CHANGING THE INTERNAL VALUE OF NFC
%
dml18j_17 = fix(nfc);
if( iflag==13 )
nfc = fix(fix(nfc./2));
end;
%
% *********************************************************
% DETERMINE NECESSARY STORAGE REQUIREMENTS
%
% FOR BASIC ARRAYS IN DBVPOR
kkkyhp = fix(ncomp.*(nfc+1) + neqivp);
kkku = fix(ncomp.*nfc.*kpts);
kkkv = fix(ncomp.*kpts);
kkkcoe = fix(dml18j_17);
kkks = fix(nfc + 1);
kkksto = fix(ncomp.*(nfc+1) + neqivp + 1);
kkkg = fix(ncomp);
%
% FOR ORTHONORMALIZATION RELATED MATTERS
dml18j_14 =fix(fix((dml18j_17.*(dml18j_17+1))./2));
dml17b_1 = fix(1 + dml18j_14 + dml18j_17);
lllip = fix(dml18j_17);
%
% FOR ADDITIONAL REQUIRED WORK SPACE
% (DLSSUD)
kkksud = fix(4.*nic +(nrowa+1).*ncomp);
lllsud = fix(nic);
% (DVECS)
kkksvc = fix(1 + 4.*dml18j_17 + 2.*dml18j_17.^2);
lllsvc = fix(2.*dml18j_17);
%
ndeq = fix(ncomp.*nfc + neqivp);
if( dml8sz_4==1 )
ndeq = fix(ndeq + ncomp);
end;
if( dml18j_12==2 )
% (DDEABM)
dml17b_17 = fix(130 + 21.*ndeq);
dml17b_18 = 51;
else;
% (DDERKF)
dml17b_17 = fix(33 + 7.*ndeq);
dml17b_18 = 34;
end;
%
% (COEF)
kkkcof = fix(5.*dml18j_17 + dml18j_17.^2);
lllcof = fix(3 + dml18j_17);
%
kkkws = fix(max([kkksud,kkksvc,dml17b_17,kkkcof]));
llliws = fix(max([lllsud,lllsvc,dml17b_18,lllcof]));
%
dml17b_2 = fix(kkkyhp + kkku + kkkv + kkkcoe + kkks + kkksto + kkkg +dml17b_1 + kkkws);
dml17b_3 = fix(17 + lllip + llliws);
% *********************************************************
% COMPUTE THE NUMBER OF POSSIBLE ORTHONORMALIZATIONS
% WITH THE ALLOTTED STORAGE
%
iwork(3) = fix(dml17b_2);
iwork(4) = fix(dml17b_1);
iwork(5) = fix(dml17b_3);
iwork(6) = fix(lllip);
nrtemp = fix(ndw - dml17b_2);
nitemp = fix(ndiw - dml17b_3);
% ...EXIT
if( nrtemp>=0 )
% ...EXIT
if( nitemp>=0 )
%
if( dml18j_8==0 )
%
mxnonr = fix(fix(nrtemp./dml17b_1));
mxnoni = fix(fix(nitemp./lllip));
dml18j_7 = fix(min(mxnonr,mxnoni));
non = fix(dml18j_7);
else;
non = 0;
dml18j_7 = fix(nrtemp);
end;
%
iwork(2) = fix(dml18j_7);
%
% *********************************************************
% CHECK FOR PRE-ASSIGNED ORTHONORMALIZATION POINTS
%
dml18j_6 = 0;
% ......EXIT
if( iwork(11)~=1 )
break;
end;
if( dml18j_7>=iwork(1) )
dml18j_6 = 1;
dml18j_7 = fix(iwork(1));
work(dml18j_7+1) = 2.0d0.*xpts(nxpts) - xpts(1);
% .........EXIT
break;
end;
end;
end;
%
iflag = -1;
if( dml18j_8~=1 )
xern1=sprintf(['%8i'], dml17b_2);
xern2=sprintf(['%8i'], dml17b_1);
xern3=sprintf(['%8i'], dml17b_3);
xern4=sprintf(['%8i'], lllip);
xermsg('SLATEC','DBVSUP',['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'], dml17b_2);
xern2=sprintf(['%8i'], dml17b_3);
xermsg('SLATEC','DBVSUP',['REQUIRED STORAGE FOR WORK ARRAY IS ',[xern1,[' + NUMBER OF ORTHONOMALIZATIONS. $$',['REQUIRED STORAGE FOR IWORK ARRAY IS ',xern2]]]],1,0);
end;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_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;
xpts_shape=zeros(xpts_shape);xpts_shape(:)=xpts(1:numel(xpts_shape));xpts=xpts_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
return;
end;
break;
end;
%
% ***************************************************************
% ALLOCATE STORAGE FROM WORK AND IWORK ARRAYS
%
% (Z)
dml17b_4 = fix(1 +(dml18j_7+1));
% (P)
dml17b_5 = fix(dml17b_4 + dml18j_14.*(non+1));
% (W)
dml17b_6 = fix(dml17b_5 + dml18j_17.*(non+1));
% (YHP)
dml17b_7 = fix(dml17b_6 + kkkyhp);
% (U)
dml17b_8 = fix(dml17b_7 + kkku);
% (V)
dml17b_9 = fix(dml17b_8 + kkkv);
% (COEF)
dml17b_10 = fix(dml17b_9 + kkkcoe);
% (S)
dml17b_11 = fix(dml17b_10 + kkks);
% (STOWA)
dml17b_12 = fix(dml17b_11 + kkksto);
% (G)
dml17b_13 = fix(dml17b_12 + kkkg);
dml17b_14 = fix(dml17b_13 + kkkws);
% REQUIRED ADDITIONAL doubleprecision 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)
dml17b_15 = fix(18 + dml18j_17.*(non+1));
dml17b_16 = fix(dml17b_15 + llliws);
% REQUIRED INTEGER WORK SPACE STARTS AT IWORK(L1)
% AND EXTENDS TO IWORK(L2-1)
%
% ***************************************************************
% SET INDICATOR FOR NORMALIZATION OF PARTICULAR SOLUTION
%
dml18j_13 = 0;
if( iwork(10)==1 )
dml18j_13 = 1;
end;
%
% ***************************************************************
% SET PIVOTING PARAMETER
%
dml18j_11 = 0;
if( iwork(15)==1 )
dml18j_11 = 1;
end;
%
% ***************************************************************
% SET OTHER COMMON BLOCK PARAMETERS
%
dml8sz_7 = fix(nfc);
dml8sz_6 = fix(ncomp);
dml8sz_3 = fix(igofx);
dml18j_4 = fix(nxpts);
dml18j_5 = fix(nic);
dml18j_2 = re;
dml18j_1 = ae;
dml18j_15 = fix(neqivp);
dml15t_14 = 20;
if( iwork(13)==-1 )
dml15t_14 = fix(max(1,iwork(14)));
end;
dml15t_5 = xpts(1);
dml15t_6 = xpts(nxpts);
dml8sz_2 = dml15t_6;
dml18j_18 = 1;
if( dml8sz_4==3 && dml18j_6==1 )
work(dml18j_7+1) = dml15t_6;
end;
%
% ***************************************************************
%
[y,nrowy,xpts,a,nrowa,alpha,b,nrowb,beta,iflag,work,iwork]=dexbvp(y,nrowy,xpts,a,nrowa,alpha,b,nrowb,beta,iflag,work,iwork);
nfc = fix(dml18j_17);
iwork(17) = fix(iwork(dml17b_15));
end;
iwork_shape=zeros(iwork_shape);iwork_shape(:)=iwork(1:numel(iwork_shape));iwork=iwork_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;
xpts_shape=zeros(xpts_shape);xpts_shape(:)=xpts(1:numel(xpts_shape));xpts=xpts_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
end %subroutine dbvsup
%DECK DCBRT
|
|
