| [intl,iorder,a,b,m,mbdcnd,bda,alpha,bdb,beta,c,d,n,nbdcnd,bdc,gama,bdd,xnu,cofx,cofy,an,bn,cn,dn,un,zn,am,bm,cm,dm,um,zm,grhs,usol,idmn,w,pertrb,ierror]=spelip(intl,iorder,a,b,m,mbdcnd,bda,alpha,bdb,beta,c,d,n,nbdcnd,bdc,gama,bdd,xnu,cofx,cofy,an,bn,cn,dn |
function [intl,iorder,a,b,m,mbdcnd,bda,alpha,bdb,beta,c,d,n,nbdcnd,bdc,gama,bdd,xnu,cofx,cofy,an,bn,cn,dn,un,zn,am,bm,cm,dm,um,zm,grhs,usol,idmn,w,pertrb,ierror]=spelip(intl,iorder,a,b,m,mbdcnd,bda,alpha,bdb,beta,c,d,n,nbdcnd,bdc,gama,bdd,xnu,cofx,cofy,an,bn,cn,dn,un,zn,am,bm,cm,dm,um,zm,grhs,usol,idmn,w,pertrb,ierror);
persistent ai ax1 axi bi bxi ci cxi cxm dj dy1 dyj ej eyj fj fyj fyn i i1 iord j mp np prtrb singlr xi yj ;
if isempty(ai), ai=0; end;
global splpcm_5; if isempty(splpcm_5), splpcm_5=0; end;
if isempty(ax1), ax1=0; end;
if isempty(axi), axi=0; end;
if isempty(bi), bi=0; end;
global splpcm_6; if isempty(splpcm_6), splpcm_6=0; end;
if isempty(bxi), bxi=0; end;
if isempty(ci), ci=0; end;
global splpcm_7; if isempty(splpcm_7), splpcm_7=0; end;
if isempty(cxi), cxi=0; end;
if isempty(cxm), cxm=0; end;
global splpcm_8; if isempty(splpcm_8), splpcm_8=0; end;
if isempty(dj), dj=0; end;
global splpcm_15; if isempty(splpcm_15), splpcm_15=0; end;
global splpcm_19; if isempty(splpcm_19), splpcm_19=0; end;
global splpcm_16; if isempty(splpcm_16), splpcm_16=0; end;
global splpcm_20; if isempty(splpcm_20), splpcm_20=0; end;
if isempty(dy1), dy1=0; end;
if isempty(dyj), dyj=0; end;
if isempty(ej), ej=0; end;
if isempty(eyj), eyj=0; end;
if isempty(fj), fj=0; end;
if isempty(fyj), fyj=0; end;
if isempty(fyn), fyn=0; end;
if isempty(prtrb), prtrb=0; end;
global splpcm_17; if isempty(splpcm_17), splpcm_17=0; end;
global splpcm_18; if isempty(splpcm_18), splpcm_18=0; end;
if isempty(xi), xi=0; end;
if isempty(yj), yj=0; end;
if isempty(i), i=0; end;
if isempty(i1), i1=0; end;
if isempty(iord), iord=0; end;
global splpcm_11; if isempty(splpcm_11), splpcm_11=0; end;
if isempty(j), j=0; end;
global splpcm_13; if isempty(splpcm_13), splpcm_13=0; end;
global splpcm_3; if isempty(splpcm_3), splpcm_3=0; end;
global splpcm_1; if isempty(splpcm_1), splpcm_1=0; end;
global splpcm_2; if isempty(splpcm_2), splpcm_2=0; end;
global splpcm_4; if isempty(splpcm_4), splpcm_4=0; end;
global splpcm_9; if isempty(splpcm_9), splpcm_9=0; end;
if isempty(mp), mp=0; end;
global splpcm_12; if isempty(splpcm_12), splpcm_12=0; end;
global splpcm_10; if isempty(splpcm_10), splpcm_10=0; end;
if isempty(np), np=0; end;
global splpcm_14; if isempty(splpcm_14), splpcm_14=0; end;
%***BEGIN PROLOGUE SPELIP
%***SUBSIDIARY
%***PURPOSE Subsidiary to SEPELI
%***LIBRARY SLATEC
%***TYPE SINGLE PRECISION (SPELIP-S)
%***AUTHOR (UNKNOWN)
%***DESCRIPTION
%
% SPELIP sets up vectors and arrays for input to BLKTRI
% and computes a second order solution in USOL. A return jump to
% SEPELI occurs if IORDER=2. If IORDER=4 a fourth order
% solution splpcm_11 generated in USOL.
%
%***SEE ALSO SEPELI
%***ROUTINES CALLED BLKTRI, CHKSNG, DEFER, MINSOL, ORTHOG, TRISP
%***COMMON BLOCKS SPLPCM
%***REVISION HISTORY (YYMMDD)
% 801001 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 891214 Prologue converted to Version 4.0 format. (BAB)
% 900402 Added TYPE section. (WRB)
%***end PROLOGUE SPELIP
%
bda_shape=size(bda);bda=reshape(bda,1,[]);
bdb_shape=size(bdb);bdb=reshape(bdb,1,[]);
bdc_shape=size(bdc);bdc=reshape(bdc,1,[]);
bdd_shape=size(bdd);bdd=reshape(bdd,1,[]);
w_shape=size(w);w=reshape(w,1,[]);
grhs_shape=size(grhs);grhs=reshape([grhs(:).',zeros(1,ceil(numel(grhs)./prod([idmn])).*prod([idmn])-numel(grhs))],idmn,[]);
usol_shape=size(usol);usol=reshape([usol(:).',zeros(1,ceil(numel(usol)./prod([idmn])).*prod([idmn])-numel(usol))],idmn,[]);
an_shape=size(an);an=reshape(an,1,[]);
bn_shape=size(bn);bn=reshape(bn,1,[]);
cn_shape=size(cn);cn=reshape(cn,1,[]);
dn_shape=size(dn);dn=reshape(dn,1,[]);
un_shape=size(un);un=reshape(un,1,[]);
zn_shape=size(zn);zn=reshape(zn,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,[]);
dm_shape=size(dm);dm=reshape(dm,1,[]);
um_shape=size(um);um=reshape(um,1,[]);
zm_shape=size(zm);zm=reshape(zm,1,[]);
% common :: ;
%% common /splpcm/ kswx , kswy , k , l , ait , bit , cit , dit ,mit , nit , is , ms , js , ns , dlx , dly ,tdlx3 , tdly3 , dlx4 , dly4;
%% common /splpcm/ splpcm_1 , splpcm_2 , splpcm_3 , splpcm_4 , splpcm_5 , splpcm_6 , splpcm_7 , splpcm_8 ,splpcm_9 , splpcm_10 , splpcm_11 , splpcm_12 , splpcm_13 , splpcm_14 , splpcm_15 , splpcm_16 ,splpcm_17 , splpcm_18 , splpcm_19 , splpcm_20;
if isempty(singlr), singlr=false; end;
%***FIRST EXECUTABLE STATEMENT SPELIP
splpcm_1 = fix(mbdcnd + 1);
splpcm_2 = fix(nbdcnd + 1);
splpcm_3 = fix(m + 1);
splpcm_4 = fix(n + 1);
splpcm_5 = a;
splpcm_6 = b;
splpcm_7 = c;
splpcm_8 = d;
%
% SET RIGHT HAND SIDE VALUES FROM GRHS IN USOL ON THE INTERIOR
% AND NON-SPECIFIED BOUNDARIES.
%
for i = 2 : m;
for j = 2 : n;
usol(i,j) = grhs(i,j);
end; j = fix(n+1);
end; i = fix(m+1);
if( splpcm_1~=2 && splpcm_1~=3 )
for j = 2 : n;
usol(1,j) = grhs(1,j);
end; j = fix(n+1);
end;
if( splpcm_1~=2 && splpcm_1~=5 )
for j = 2 : n;
usol(splpcm_3,j) = grhs(splpcm_3,j);
end; j = fix(n+1);
end;
if( splpcm_2~=2 && splpcm_2~=3 )
for i = 2 : m;
usol(i,1) = grhs(i,1);
end; i = fix(m+1);
end;
if( splpcm_2~=2 && splpcm_2~=5 )
for i = 2 : m;
usol(i,splpcm_4) = grhs(i,splpcm_4);
end; i = fix(m+1);
end;
if( splpcm_1~=2 && splpcm_1~=3 && splpcm_2~=2 && splpcm_2~=3 )
usol(1,1)= grhs(1,1);
end;
if( splpcm_1~=2 && splpcm_1~=5 && splpcm_2~=2 && splpcm_2~=3 )
usol(splpcm_3,1)= grhs(splpcm_3,1);
end;
if( splpcm_1~=2 && splpcm_1~=3 && splpcm_2~=2 && splpcm_2~=5 )
usol(1,splpcm_4)= grhs(1,splpcm_4);
end;
if( splpcm_1~=2 && splpcm_1~=5 && splpcm_2~=2 && splpcm_2~=5 )
usol(splpcm_3,splpcm_4)= grhs(splpcm_3,splpcm_4);
end;
i1 = 1;
%
% SET SWITCHES FOR PERIODIC OR NON-PERIODIC BOUNDARIES
%
mp = 1;
np = 1;
if( splpcm_1==1 )
mp = 0;
end;
if( splpcm_2==1 )
np = 0;
end;
%
% SET DLX,DLY AND SIZE OF BLOCK TRI-DIAGONAL SYSTEM GENERATED
% IN NINT,MINT
%
splpcm_15 =(splpcm_6-splpcm_5)./m;
splpcm_9 = fix(splpcm_3 - 1);
if( splpcm_1==2 )
splpcm_9 = fix(splpcm_3 - 2);
end;
if( splpcm_1==4 )
splpcm_9 = fix(splpcm_3);
end;
splpcm_16 =(splpcm_8-splpcm_7)./n;
splpcm_10 = fix(splpcm_4 - 1);
if( splpcm_2==2 )
splpcm_10 = fix(splpcm_4 - 2);
end;
if( splpcm_2==4 )
splpcm_10 = fix(splpcm_4);
end;
splpcm_17 = 2.0.*splpcm_15.^3;
splpcm_19 = splpcm_15.^4;
splpcm_18 = 2.0.*splpcm_16.^3;
splpcm_20 = splpcm_16.^4;
%
% SET SUBSCRIPT LIMITS FOR PORTION OF ARRAY TO INPUT TO BLKTRI
%
splpcm_11 = 1;
splpcm_13 = 1;
if( splpcm_1==2 || splpcm_1==3 )
splpcm_11 = 2;
end;
if( splpcm_2==2 || splpcm_2==3 )
splpcm_13 = 2;
end;
splpcm_14 = fix(splpcm_10 + splpcm_13 - 1);
splpcm_12 = fix(splpcm_9 + splpcm_11 - 1);
%
% SET X - DIRECTION
%
for i = 1 : splpcm_9;
xi = splpcm_5 +(splpcm_11+i-2).*splpcm_15;
[xi,ai,bi,ci]=cofx(xi,ai,bi,ci);
axi =(ai./splpcm_15-0.5.*bi)./splpcm_15;
bxi = -2..*ai./splpcm_15.^2 + ci;
cxi =(ai./splpcm_15+0.5.*bi)./splpcm_15;
am(i) = axi;
bm(i) = bxi;
cm(i) = cxi;
end; i = fix(splpcm_9+1);
%
% SET Y DIRECTION
%
for j = 1 : splpcm_10;
yj = splpcm_7 +(splpcm_13+j-2).*splpcm_16;
[yj,dj,ej,fj]=cofy(yj,dj,ej,fj);
dyj =(dj./splpcm_16-0.5.*ej)./splpcm_16;
eyj =(-2..*dj./splpcm_16.^2+fj);
fyj =(dj./splpcm_16+0.5.*ej)./splpcm_16;
an(j) = dyj;
bn(j) = eyj;
cn(j) = fyj;
end; j = fix(splpcm_10+1);
%
% ADJUST EDGES IN X DIRECTION UNLESS PERIODIC
%
ax1 = am(1);
cxm = cm(splpcm_9);
if( splpcm_1~=1 )
if( splpcm_1==3 )
%
% DIRICHLET-MIXED IN X DIRECTION
%
am(1) = 0.0;
am(splpcm_9) = am(splpcm_9) + cxm;
bm(splpcm_9) = bm(splpcm_9) - 2..*beta.*splpcm_15.*cxm;
cm(splpcm_9) = 0.0;
elseif( splpcm_1==4 ) ;
%
% MIXED - MIXED IN X DIRECTION
%
am(1) = 0.0;
bm(1) = bm(1) + 2..*splpcm_15.*alpha.*ax1;
cm(1) = cm(1) + ax1;
am(splpcm_9) = am(splpcm_9) + cxm;
bm(splpcm_9) = bm(splpcm_9) - 2..*splpcm_15.*beta.*cxm;
cm(splpcm_9) = 0.0;
elseif( splpcm_1==5 ) ;
%
% MIXED-DIRICHLET IN X DIRECTION
%
am(1) = 0.0;
bm(1) = bm(1) + 2..*alpha.*splpcm_15.*ax1;
cm(1) = cm(1) + ax1;
cm(splpcm_9) = 0.0;
else;
%
% DIRICHLET-DIRICHLET IN X DIRECTION
%
am(1) = 0.0;
cm(splpcm_9) = 0.0;
end;
end;
%
% ADJUST IN Y DIRECTION UNLESS PERIODIC
%
dy1 = an(1);
fyn = cn(splpcm_10);
if( splpcm_2~=1 )
if( splpcm_2==3 )
%
% DIRICHLET-MIXED IN Y DIRECTION
%
an(1) = 0.0;
an(splpcm_10) = an(splpcm_10) + fyn;
bn(splpcm_10) = bn(splpcm_10) - 2..*splpcm_16.*xnu.*fyn;
cn(splpcm_10) = 0.0;
elseif( splpcm_2==4 ) ;
%
% MIXED - MIXED DIRECTION IN Y DIRECTION
%
an(1) = 0.0;
bn(1) = bn(1) + 2..*splpcm_16.*gama.*dy1;
cn(1) = cn(1) + dy1;
an(splpcm_10) = an(splpcm_10) + fyn;
bn(splpcm_10) = bn(splpcm_10) - 2.0.*splpcm_16.*xnu.*fyn;
cn(splpcm_10) = 0.0;
elseif( splpcm_2==5 ) ;
%
% MIXED-DIRICHLET IN Y DIRECTION
%
an(1) = 0.0;
bn(1) = bn(1) + 2..*splpcm_16.*gama.*dy1;
cn(1) = cn(1) + dy1;
cn(splpcm_10) = 0.0;
else;
%
% DIRICHLET-DIRICHLET IN Y DIRECTION
%
an(1) = 0.0;
cn(splpcm_10) = 0.0;
end;
end;
if( splpcm_1~=1 )
%
% ADJUST USOL ALONG X EDGE
%
for j = splpcm_13 : splpcm_14;
if( splpcm_1~=2 && splpcm_1~=3 )
usol(splpcm_11,j) = usol(splpcm_11,j) + 2.0.*splpcm_15.*ax1.*bda(j);
else;
usol(splpcm_11,j) = usol(splpcm_11,j) - ax1.*usol(1,j);
end;
if( splpcm_1~=2 && splpcm_1~=5 )
usol(splpcm_12,j) = usol(splpcm_12,j) - 2.0.*splpcm_15.*cxm.*bdb(j);
else;
usol(splpcm_12,j) = usol(splpcm_12,j) - cxm.*usol(splpcm_3,j);
end;
end; j = fix(splpcm_14+1);
end;
if( splpcm_2~=1 )
%
% ADJUST USOL ALONG Y EDGE
%
for i = splpcm_11 : splpcm_12;
if( splpcm_2~=2 && splpcm_2~=3 )
usol(i,splpcm_13) = usol(i,splpcm_13) + 2.0.*splpcm_16.*dy1.*bdc(i);
else;
usol(i,splpcm_13) = usol(i,splpcm_13) - dy1.*usol(i,1);
end;
if( splpcm_2~=2 && splpcm_2~=5 )
usol(i,splpcm_14) = usol(i,splpcm_14) - 2.0.*splpcm_16.*fyn.*bdd(i);
else;
usol(i,splpcm_14) = usol(i,splpcm_14) - fyn.*usol(i,splpcm_4);
end;
end; i = fix(splpcm_12+1);
end;
%
% SAVE ADJUSTED EDGES IN GRHS IF IORDER=4
%
if( iorder==4 )
for j = splpcm_13 : splpcm_14;
grhs(splpcm_11,j) = usol(splpcm_11,j);
grhs(splpcm_12,j) = usol(splpcm_12,j);
end; j = fix(splpcm_14+1);
for i = splpcm_11 : splpcm_12;
grhs(i,splpcm_13) = usol(i,splpcm_13);
grhs(i,splpcm_14) = usol(i,splpcm_14);
end; i = fix(splpcm_12+1);
end;
iord = fix(iorder);
pertrb = 0.0;
%
% CHECK IF OPERATOR IS SINGULAR
%
[mbdcnd,nbdcnd,alpha,beta,gama,xnu,cofx,cofy,singlr]=chksng(mbdcnd,nbdcnd,alpha,beta,gama,xnu,cofx,cofy,singlr);
%
% COMPUTE NON-ZERO EIGENVECTOR IN NULL SPACE OF TRANSPOSE
% IF SINGULAR
%
if( singlr )
[splpcm_9,am,bm,cm,dm,um,zm]=trisp(splpcm_9,am,bm,cm,dm,um,zm);
end;
if( singlr )
[splpcm_10,an,bn,cn,dn,un,zn]=trisp(splpcm_10,an,bn,cn,dn,un,zn);
end;
%
% MAKE INITIALIZATION CALL TO BLKTRI
%
if( intl==0 )
[intl,np,splpcm_10,an,bn,cn,mp,splpcm_9,am,bm,cm,idmn,usol(sub2ind(size(usol),splpcm_11,splpcm_13):end),ierror,w]=blktri(intl,np,splpcm_10,an,bn,cn,mp,splpcm_9,am,bm,cm,idmn,usol(sub2ind(size(usol),splpcm_11,splpcm_13):end),ierror,w);
end;
if( ierror==0 )
%
% ADJUST RIGHT HAND SIDE IF NECESSARY
%
while( true );
if( singlr )
[usol,idmn,zn,zm,pertrb]=orthog(usol,idmn,zn,zm,pertrb);
end;
%
% COMPUTE SOLUTION
%
[i1,np,splpcm_10,an,bn,cn,mp,splpcm_9,am,bm,cm,idmn,usol(sub2ind(size(usol),splpcm_11,splpcm_13):end),ierror,w]=blktri(i1,np,splpcm_10,an,bn,cn,mp,splpcm_9,am,bm,cm,idmn,usol(sub2ind(size(usol),splpcm_11,splpcm_13):end),ierror,w);
if( ierror~=0 )
break;
end;
%
% SET PERIODIC BOUNDARIES IF NECESSARY
%
if( splpcm_1==1 )
for j = 1 : splpcm_4;
usol(splpcm_3,j) = usol(1,j);
end; j = fix(splpcm_4+1);
end;
if( splpcm_2==1 )
for i = 1 : splpcm_3;
usol(i,splpcm_4) = usol(i,1);
end; i = fix(splpcm_3+1);
end;
%
% MINIMIZE SOLUTION WITH RESPECT TO WEIGHTED LEAST SQUARES
% NORM IF OPERATOR IS SINGULAR
%
if( singlr )
[usol,idmn,zn,zm,prtrb]=minsol(usol,idmn,zn,zm,prtrb);
end;
%
% RETURN IF DEFERRED CORRECTIONS AND A FOURTH ORDER SOLUTION ARE
% NOT FLAGGED
%
if( iord==2 )
break;
end;
iord = 2;
%
% COMPUTE NEW RIGHT HAND SIDE FOR FOURTH ORDER SOLUTION
%
[cofx,cofy,idmn,usol,grhs]=defer(cofx,cofy,idmn,usol,grhs);
end;
end;
bda_shape=zeros(bda_shape);bda_shape(:)=bda(1:numel(bda_shape));bda=bda_shape;
bdb_shape=zeros(bdb_shape);bdb_shape(:)=bdb(1:numel(bdb_shape));bdb=bdb_shape;
bdc_shape=zeros(bdc_shape);bdc_shape(:)=bdc(1:numel(bdc_shape));bdc=bdc_shape;
bdd_shape=zeros(bdd_shape);bdd_shape(:)=bdd(1:numel(bdd_shape));bdd=bdd_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
grhs_shape=zeros(grhs_shape);grhs_shape(:)=grhs(1:numel(grhs_shape));grhs=grhs_shape;
usol_shape=zeros(usol_shape);usol_shape(:)=usol(1:numel(usol_shape));usol=usol_shape;
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;
dn_shape=zeros(dn_shape);dn_shape(:)=dn(1:numel(dn_shape));dn=dn_shape;
un_shape=zeros(un_shape);un_shape(:)=un(1:numel(un_shape));un=un_shape;
zn_shape=zeros(zn_shape);zn_shape(:)=zn(1:numel(zn_shape));zn=zn_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;
dm_shape=zeros(dm_shape);dm_shape(:)=dm(1:numel(dm_shape));dm=dm_shape;
um_shape=zeros(um_shape);um_shape(:)=um(1:numel(um_shape));um=um_shape;
zm_shape=zeros(zm_shape);zm_shape(:)=zm(1:numel(zm_shape));zm=zm_shape;
end %subroutine spelip
%DECK SPENC
|
|
