| [iflg,np,n,an,bn,cn,mp,m,am,bm,cm,idimy,y,ierror,w]=cblktr(iflg,np,n,an,bn,cn,mp,m,am,bm,cm,idimy,y,ierror,w); |
function [iflg,np,n,an,bn,cn,mp,m,am,bm,cm,idimy,y,ierror,w]=cblktr(iflg,np,n,an,bn,cn,mp,m,am,bm,cm,idimy,y,ierror,w);
persistent iw1 iw2 iw3 iwah iwbh iwd iwu iww m2 nh nl ;
global ccblk_4; if isempty(ccblk_4), ccblk_4=0; end;
global ccblk_3; if isempty(ccblk_3), ccblk_3=0; end;
global ccblk_7; if isempty(ccblk_7), ccblk_7=0; end;
if isempty(iw1), iw1=0; end;
if isempty(iw2), iw2=0; end;
if isempty(iw3), iw3=0; end;
if isempty(iwah), iwah=0; end;
if isempty(iwbh), iwbh=0; end;
if isempty(iwd), iwd=0; end;
if isempty(iwu), iwu=0; end;
if isempty(iww), iww=0; end;
global ccblk_2; if isempty(ccblk_2), ccblk_2=0; end;
if isempty(m2), m2=0; end;
global ccblk_6; if isempty(ccblk_6), ccblk_6=0; end;
if isempty(nh), nh=0; end;
if isempty(nl), nl=0; end;
global ccblk_5; if isempty(ccblk_5), ccblk_5=0; end;
global ccblk_1; if isempty(ccblk_1), ccblk_1=0; end;
%***BEGIN PROLOGUE CBLKTR
%***PURPOSE Solve a block tridiagonal system of linear equations
% (usually resulting from the discretization of separable
% two-dimensional elliptic equations).
%***LIBRARY SLATEC (FISHPACK)
%***CATEGORY I2B4B
%***TYPE COMPLEX (BLKTRI-S, CBLKTR-C)
%***KEYWORDS ELLIPTIC PDE, FISHPACK, TRIDIAGONAL LINEAR SYSTEM
%***AUTHOR Adams, J., (NCAR)
% Swarztrauber, P. N., (NCAR)
% Sweet, R., (NCAR)
%***DESCRIPTION
%
% subroutine CBLKTR is a complex version of subroutine BLKTRI.
% Both subroutines solve a system of linear equations of the form
%
% AN(J)*X(I,J-1) + AM(I)*X(I-1,J) + (BN(J)+BM(I))*X(I,J)
%
% + CN(J)*X(I,J+1) + CM(I)*X(I+1,J) = Y(I,J)
%
% For I = 1,2,...,M and J = 1,2,...,N.
%
% I+1 and I-1 are evaluated modulo M and J+1 and J-1 modulo N, i.e.,
%
% X(I,0) = X(I,N), X(I,N+1) = X(I,1),
% X(0,J) = X(M,J), X(M+1,J) = X(1,J).
%
% These equations usually result from the discretization of
% separable elliptic equations. Boundary conditions may be
% Dirichlet, Neumann, or periodic.
%
%
% * * * * * * * * * * On INPUT * * * * * * * * * *
%
% IFLG
% = 0 Initialization only. Certain quantities that depend on NP,
% N, AN, BN, and CN are computed and stored in the work
% array W.
% = 1 The quantities that were computed in the initialization are
% used to obtain the solution X(I,J).
%
% NOTE A call with IFLG=0 takes approximately one half the time
% time as a call with IFLG = 1. However, the
% initialization does not have to be repeated unless NP, N,
% AN, BN, or CN change.
%
% NP
% = 0 If AN(1) and CN(N) are not zero, which corresponds to
% periodic boundary conditions.
% = 1 If AN(1) and CN(N) are zero.
%
% N
% The number of unknowns in the J-direction. N must be greater
% than 4. The operation count is proportional to MNlog2(N), hence
% N should be selected less than or equal to M.
%
% AN,BN,CN
% Real one-dimensional arrays of length N that specify the
% coefficients in the linear equations given above.
%
% MP
% = 0 If AM(1) and CM(M) are not zero, which corresponds to
% periodic boundary conditions.
% = 1 If AM(1) = CM(M) = 0 .
%
% M
% The number of unknowns in the I-direction. M must be greater
% than 4.
%
% AM,BM,CM
% Complex one-dimensional arrays of length M that specify the
% coefficients in the linear equations given above.
%
% IDIMY
% The row (or first) dimension of the two-dimensional array Y as
% it appears in the program calling BLKTRI. This parameter is
% used to specify the variable dimension of Y. IDIMY must be at
% least M.
%
% Y
% A complex two-dimensional array that specifies the values of
% the right side of the linear system of equations given above.
% Y must be dimensioned Y(IDIMY,N) with IDIMY .GE. M.
%
% W
% A one-dimensional array that must be provided by the user for
% work space.
% If NP=1 define K=INT(log2(N))+1 and set L=2**(K+1) then
% W must have dimension (K-2)*L+K+5+MAX(2N,12M)
%
% If NP=0 define K=INT(log2(N-1))+1 and set L=2**(K+1) then
% W must have dimension (K-2)*L+K+5+2N+MAX(2N,12M)
%
% **IMPORTANT** For purposes of checking, the required dimension
% of W is computed by BLKTRI and stored in W(1)
% in floating point format.
%
% * * * * * * * * * * On Output * * * * * * * * * *
%
% Y
% contains the solution X.
%
% IERROR
% An error flag that indicates invalid input parameters. Except
% for number zero, a solution is not attempted.
%
% = 0 No error.
% = 1 M is less than 5.
% = 2 N is less than 5.
% = 3 IDIMY is less than M.
% = 4 BLKTRI failed while computing results that depend on the
% coefficient arrays AN, BN, CN. Check these arrays.
% = 5 AN(J)*CN(J-1) is less than 0 for some J. Possible reasons
% for this condition are
% 1. The arrays AN and CN are not correct.
% 2. Too large a grid spacing was used in the discretization
% of the elliptic equation.
% 3. The linear equations resulted from a partial
% differential equation which was not elliptic.
%
% W
% contains intermediate values that must not be destroyed if
% CBLKTR will be called again with IFLG=1. W(1) contains the
% number of locations required by W in floating point format.
%
% *Long Description:
%
% * * * * * * * program Specifications * * * * * * * * * * * *
%
% Dimension of AN(N),BN(N),CN(N),AM(M),BM(M),CM(M),Y(IDIMY,N)
% Arguments W(see argument list)
%
% Latest June 1979
% Revision
%
% Required CBLKTR,CBLKT1,PROC,PROCP,CPROC,CPROCP,CCMPB,INXCA,
% Subprograms INXCB,INXCC,CPADD,PGSF,PPGSF,PPPSF,BCRH,TEVLC,
% R1MACH
%
% Special The algorithm may fail if ABS(BM(I)+BN(J)) is less
% Conditions than ABS(AM(I))+ABS(AN(J))+ABS(CM(I))+ABS(CN(J))
% for some I and J. The algorithm will also fail if
% AN(J)*CN(J-1) is less than zero for some J.
% See the description of the output parameter IERROR.
%
% Common CCBLK
% Blocks
%
% I/O NONE
%
% Precision Single
%
% Specialist Paul Swarztrauber
%
% Language FORTRAN
%
% History CBLKTR is a complex version of BLKTRI (version 3)
%
% Algorithm Generalized Cyclic Reduction (see reference below)
%
% Space
% Required CONTROL DATA 7600
%
% Portability American National Standards Institute FORTRAN.
% The machine accuracy is set using function R1MACH.
%
% Required NONE
% Resident
% Routines
%
% References Swarztrauber,P. and R. SWEET, 'Efficient Fortran
% Subprograms for the solution of elliptic equations'
% NCAR TN/IA-109, July, 1975, 138 PP.
%
% SWARZTRAUBER P. ,'A Direct Method for The Discrete
% Solution of Separable Elliptic Equations', SIAM
% J. Numer. Anal.,11(1974) PP. 1136-1150.
%
% * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
%
%***REFERENCES P. N. Swarztrauber and R. Sweet, Efficient Fortran
% subprograms for the solution of elliptic equations,
% NCAR TN/IA-109, July 1975, 138 pp.
% P. N. Swarztrauber, A direct method for the discrete
% solution of separable elliptic equations, SIAM Journal
% on Numerical Analysis 11, (1974), pp. 1136-1150.
%***ROUTINES CALLED CBLKT1, CCMPB, CPROC, CPROCP, PROC, PROCP
%***COMMON BLOCKS CCBLK
%***REVISION HISTORY (YYMMDD)
% 801001 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 890531 REVISION DATE from Version 3.2
% 891214 Prologue converted to Version 4.0 format. (BAB)
% 920501 Reformatted the REFERENCES section. (WRB)
%***end PROLOGUE CBLKTR
%
an_shape=size(an);an=reshape(an,1,[]);
bn_shape=size(bn);bn=reshape(bn,1,[]);
cn_shape=size(cn);cn=reshape(cn,1,[]);
am_shape=size(am);am=reshape(am,1,[]);
bm_shape=size(bm);bm=reshape(bm,1,[]);
cm_shape=size(cm);cm=reshape(cm,1,[]);
y_shape=size(y);y=reshape([y(:).',zeros(1,ceil(numel(y)./prod([idimy])).*prod([idimy])-numel(y))],idimy,[]);
w_shape=size(w);w=reshape(w,1,[]);
% common :: ;
%% common /ccblk / npp , k , eps , cnv , nm , ncmplx , ik;
%% common /ccblk / ccblk_1 , ccblk_2 , ccblk_3 , ccblk_4 , ccblk_5 , ccblk_6 , ccblk_7;
%***FIRST EXECUTABLE STATEMENT CBLKTR
ccblk_5 = fix(n);
m2 = fix(m + m);
ierror = 0;
if( m<5 )
ierror = 1;
elseif( ccblk_5<3 ) ;
ierror = 2;
elseif( idimy<m ) ;
ierror = 3;
else;
nh = fix(n);
ccblk_1 = fix(np);
if( ccblk_1~=0 )
nh = fix(nh + 1);
end;
ccblk_7 = 2;
ccblk_2 = 1;
while (1);
ccblk_7 = fix(ccblk_7 + ccblk_7);
ccblk_2 = fix(ccblk_2 + 1);
if( nh<=ccblk_7 )
break;
end;
end;
nl = fix(ccblk_7);
ccblk_7 = fix(ccblk_7 + ccblk_7);
nl = fix(nl - 1);
iwah =fix((ccblk_2-2).*ccblk_7 + ccblk_2 + 6);
if( ccblk_1~=0 )
%
% DIVIDE W INTO WORKING SUB ARRAYS
%
iw1 = fix(iwah);
iwbh = fix(iw1 + ccblk_5);
w(1) = iw1 - 1 + max(2.*ccblk_5,12.*m);
else;
iwbh = fix(iwah + ccblk_5 + ccblk_5);
iw1 = fix(iwbh);
w(1) = iw1 - 1 + max(2.*ccblk_5,12.*m);
ccblk_5 = fix(ccblk_5 - 1);
end;
%
% subroutine CCMPB COMPUTES THE ROOTS OF THE B POLYNOMIALS
%
if( ierror==0 )
iw2 = fix(iw1 + m2);
iw3 = fix(iw2 + m2);
iwd = fix(iw3 + m2);
iww = fix(iwd + m2);
iwu = fix(iww + m2);
if( iflg==0 )
[nl,ierror,an,bn,cn,dumvar6,dumvar7,dumvar8]=ccmpb(nl,ierror,an,bn,cn,w(sub2ind(size(w),max(2,1)):end),w(sub2ind(size(w),max(iwah,1)):end),w(sub2ind(size(w),max(iwbh,1)):end)); dumvar6i=find((w(sub2ind(size(w),max(2,1)):end))~=(dumvar6));dumvar7i=find((w(sub2ind(size(w),max(iwah,1)):end))~=(dumvar7));dumvar8i=find((w(sub2ind(size(w),max(iwbh,1)):end))~=(dumvar8)); w(2-1+dumvar6i)=dumvar6(dumvar6i); w(iwah-1+dumvar7i)=dumvar7(dumvar7i); w(iwbh-1+dumvar8i)=dumvar8(dumvar8i);
elseif( mp~=0 ) ;
%
% subroutine CBLKT1 SOLVES THE LINEAR SYSTEM
%
[nl,an,bn,cn,m,am,bm,cm,idimy,y,dumvar11,dumvar12,dumvar13,dumvar14,dumvar15,dumvar16,dumvar17]=cblkt1(nl,an,bn,cn,m,am,bm,cm,idimy,y,w(sub2ind(size(w),max(2,1)):end),w(sub2ind(size(w),max(iw1,1)):end),w(sub2ind(size(w),max(iw2,1)):end),w(sub2ind(size(w),max(iw3,1)):end),w(sub2ind(size(w),max(iwd,1)):end),w(sub2ind(size(w),max(iww,1)):end),w(sub2ind(size(w),max(iwu,1)):end),@proc,@cproc); dumvar11i=find((w(sub2ind(size(w),max(2,1)):end))~=(dumvar11));dumvar12i=find((w(sub2ind(size(w),max(iw1,1)):end))~=(dumvar12));dumvar13i=find((w(sub2ind(size(w),max(iw2,1)):end))~=(dumvar13));dumvar14i=find((w(sub2ind(size(w),max(iw3,1)):end))~=(dumvar14));dumvar15i=find((w(sub2ind(size(w),max(iwd,1)):end))~=(dumvar15));dumvar16i=find((w(sub2ind(size(w),max(iww,1)):end))~=(dumvar16));dumvar17i=find((w(sub2ind(size(w),max(iwu,1)):end))~=(dumvar17)); w(2-1+dumvar11i)=dumvar11(dumvar11i); w(iw1-1+dumvar12i)=dumvar12(dumvar12i); w(iw2-1+dumvar13i)=dumvar13(dumvar13i); w(iw3-1+dumvar14i)=dumvar14(dumvar14i); w(iwd-1+dumvar15i)=dumvar15(dumvar15i); w(iww-1+dumvar16i)=dumvar16(dumvar16i); w(iwu-1+dumvar17i)=dumvar17(dumvar17i);
else;
[nl,an,bn,cn,m,am,bm,cm,idimy,y,dumvar11,dumvar12,dumvar13,dumvar14,dumvar15,dumvar16,dumvar17]=cblkt1(nl,an,bn,cn,m,am,bm,cm,idimy,y,w(sub2ind(size(w),max(2,1)):end),w(sub2ind(size(w),max(iw1,1)):end),w(sub2ind(size(w),max(iw2,1)):end),w(sub2ind(size(w),max(iw3,1)):end),w(sub2ind(size(w),max(iwd,1)):end),w(sub2ind(size(w),max(iww,1)):end),w(sub2ind(size(w),max(iwu,1)):end),@procp,@cprocp); dumvar11i=find((w(sub2ind(size(w),max(2,1)):end))~=(dumvar11));dumvar12i=find((w(sub2ind(size(w),max(iw1,1)):end))~=(dumvar12));dumvar13i=find((w(sub2ind(size(w),max(iw2,1)):end))~=(dumvar13));dumvar14i=find((w(sub2ind(size(w),max(iw3,1)):end))~=(dumvar14));dumvar15i=find((w(sub2ind(size(w),max(iwd,1)):end))~=(dumvar15));dumvar16i=find((w(sub2ind(size(w),max(iww,1)):end))~=(dumvar16));dumvar17i=find((w(sub2ind(size(w),max(iwu,1)):end))~=(dumvar17)); w(2-1+dumvar11i)=dumvar11(dumvar11i); w(iw1-1+dumvar12i)=dumvar12(dumvar12i); w(iw2-1+dumvar13i)=dumvar13(dumvar13i); w(iw3-1+dumvar14i)=dumvar14(dumvar14i); w(iwd-1+dumvar15i)=dumvar15(dumvar15i); w(iww-1+dumvar16i)=dumvar16(dumvar16i); w(iwu-1+dumvar17i)=dumvar17(dumvar17i);
end;
end;
end;
an_shape=zeros(an_shape);an_shape(:)=an(1:numel(an_shape));an=an_shape;
bn_shape=zeros(bn_shape);bn_shape(:)=bn(1:numel(bn_shape));bn=bn_shape;
cn_shape=zeros(cn_shape);cn_shape(:)=cn(1:numel(cn_shape));cn=cn_shape;
am_shape=zeros(am_shape);am_shape(:)=am(1:numel(am_shape));am=am_shape;
bm_shape=zeros(bm_shape);bm_shape(:)=bm(1:numel(bm_shape));bm=bm_shape;
cm_shape=zeros(cm_shape);cm_shape(:)=cm(1:numel(cm_shape));cm=cm_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
end %subroutine cblktr
%DECK CBRT
|
|
