| [nperod,n,mperod,m,a,b,c,idimy,y,ierror,w]=cmgnbn(nperod,n,mperod,m,a,b,c,idimy,y,ierror,w); |
function [nperod,n,mperod,m,a,b,c,idimy,y,ierror,w]=cmgnbn(nperod,n,mperod,m,a,b,c,idimy,y,ierror,w);
persistent a1 gt200 gt300 gt400 gt500 i ipstor irev iwb2 iwb3 iwba iwbb iwbc iwd iwp iwtcos iww1 iww2 iww3 j k mh mhm1 mhmi mhpi modd mp mskip nby2 np ;
if isempty(i), i=0; end;
if isempty(ipstor), ipstor=0; end;
if isempty(irev), irev=0; end;
if isempty(iwb2), iwb2=0; end;
if isempty(iwb3), iwb3=0; end;
if isempty(iwba), iwba=0; end;
if isempty(iwbb), iwbb=0; end;
if isempty(iwbc), iwbc=0; end;
if isempty(iwd), iwd=0; end;
if isempty(iwp), iwp=0; end;
if isempty(iwtcos), iwtcos=0; end;
if isempty(iww1), iww1=0; end;
if isempty(iww2), iww2=0; end;
if isempty(iww3), iww3=0; end;
if isempty(j), j=0; end;
if isempty(k), k=0; end;
if isempty(mh), mh=0; end;
if isempty(gt200), gt200=0; end;
if isempty(gt300), gt300=0; end;
if isempty(gt400), gt400=0; end;
if isempty(gt500), gt500=0; end;
%***BEGIN PROLOGUE CMGNBN
%***PURPOSE Solve a complex block tridiagonal linear system of
% equations by a cyclic reduction algorithm.
%***LIBRARY SLATEC (FISHPACK)
%***CATEGORY I2B4B
%***TYPE COMPLEX (GENBUN-S, CMGNBN-C)
%***KEYWORDS CYCLIC REDUCTION, ELLIPTIC PDE, FISHPACK,
% TRIDIAGONAL LINEAR SYSTEM
%***AUTHOR Adams, J., (NCAR)
% Swarztrauber, P. N., (NCAR)
% Sweet, R., (NCAR)
%***DESCRIPTION
%
% subroutine CMGNBN solves the complex linear system of equations
%
% A(I)*X(I-1,J) + B(I)*X(I,J) + C(I)*X(I+1,J)
%
% + X(I,J-1) - 2.*X(I,J) + X(I,J+1) = Y(I,J)
%
% For I = 1,2,...,M and J = 1,2,...,N.
%
% The indices I+1 and I-1 are evaluated modulo M, i.e.,
% X(0,J) = X(M,J) and X(M+1,J) = X(1,J), and X(I,0) may be equal to
% 0, X(I,2), or X(I,N) and X(I,N+1) may be equal to 0, X(I,N-1), or
% X(I,1) depending on an input parameter.
%
%
% * * * * * * * * Parameter Description * * * * * * * * * *
%
% * * * * * * On Input * * * * * *
%
% NPEROD
% Indicates the values that X(I,0) and X(I,N+1) are assumed to
% have.
%
% = 0 If X(I,0) = X(I,N) and X(I,N+1) = X(I,1).
% = 1 If X(I,0) = X(I,N+1) = 0 .
% = 2 If X(I,0) = 0 and X(I,N+1) = X(I,N-1).
% = 3 If X(I,0) = X(I,2) and X(I,N+1) = X(I,N-1).
% = 4 If X(I,0) = X(I,2) and X(I,N+1) = 0.
%
% N
% The number of unknowns in the J-direction. N must be greater
% than 2.
%
% MPEROD
% = 0 If A(1) and C(M) are not zero
% = 1 If A(1) = C(M) = 0
%
% M
% The number of unknowns in the I-direction. N must be greater
% than 2.
%
% A,B,C
% One-dimensional complex arrays of length M that specify the
% coefficients in the linear equations given above. If MPEROD = 0
% the array elements must not depend upon the index I, but must be
% constant. Specifically, the subroutine checks the following
% condition
%
% A(I) = C(1)
% C(I) = C(1)
% B(I) = B(1)
%
% For I=1,2,...,M.
%
% IDIMY
% The row (or first) dimension of the two-dimensional array Y as
% it appears in the program calling CMGNBN. This parameter is
% used to specify the variable dimension of Y. IDIMY must be at
% least M.
%
% Y
% A two-dimensional complex array that specifies the values of the
% right side of the linear system of equations given above. Y
% must be dimensioned at least M*N.
%
% W
% A one-dimensional complex array that must be provided by the
% user for work space. W may require up to 4*N +
% (10 + INT(log2(N)))*M LOCATIONS. The actual number of locations
% used is computed by CMGNBN and is returned in location W(1).
%
%
% * * * * * * On Output * * * * * *
%
% Y
% contains the solution X.
%
% IERROR
% An error flag which indicates invalid input parameters. Except
% for number zero, a solution is not attempted.
%
% = 0 No error.
% = 1 M .LE. 2
% = 2 N .LE. 2
% = 3 IDIMY .LT. M
% = 4 NPEROD .LT. 0 or NPEROD .GT. 4
% = 5 MPEROD .LT. 0 or MPEROD .GT. 1
% = 6 A(I) .NE. C(1) or C(I) .NE. C(1) or B(I) .NE. B(1) for
% some I=1,2,...,M.
% = 7 A(1) .NE. 0 or C(M) .NE. 0 and MPEROD = 1
%
% W
% W(1) contains the required length of W.
%
% *Long Description:
%
% * * * * * * * program Specifications * * * * * * * * * * * *
%
% Dimension of A(M),B(M),C(M),Y(IDIMY,N),W(see parameter list)
% Arguments
%
% Latest June 1979
% Revision
%
% Subprograms CMGNBN,CMPOSD,CMPOSN,CMPOSP,CMPCSG,CMPMRG,
% Required CMPTRX,CMPTR3,PIMACH
%
% Special None
% Conditions
%
% Common None
% Blocks
%
% I/O None
%
% Precision Single
%
% Specialist Roland Sweet
%
% Language FORTRAN
%
% History Written by Roland Sweet at NCAR in June, 1977
%
% Algorithm The linear system is solved by a cyclic reduction
% algorithm described in the reference.
%
% Space 4944(DECIMAL) = 11520(octal) locations on the NCAR
% Required Control Data 7600
%
% Timing and The execution time T on the NCAR Control Data
% Accuracy 7600 for subroutine CMGNBN is roughly proportional
% to M*N*log2(N), but also depends on the input
% parameter NPEROD. Some typical values are listed
% in the table below.
% To measure the accuracy of the algorithm a
% uniform random number generator was used to create
% a solution array X for the system given in the
% 'PURPOSE' with
%
% A(I) = C(I) = -0.5*B(I) = 1, I=1,2,...,M
%
% and, when MPEROD = 1
%
% A(1) = C(M) = 0
% A(M) = C(1) = 2.
%
% The solution X was substituted into the given sys-
% tem and a right side Y was computed. Using this
% array Y subroutine CMGNBN was called to produce an
% approximate solution Z. Then the relative error,
% defined as
%
% E = MAX(ABS(Z(I,J)-X(I,J)))/MAX(ABS(X(I,J)))
%
% where the two maxima are taken over all I=1,2,...,M
% and J=1,2,...,N, was computed. The value of E is
% given in the table below for some typical values of
% M and N.
%
%
% M (=N) MPEROD NPEROD T(MSECS) E
% ------ ------ ------ -------- ------
%
% 31 0 0 77 1.0E-12
% 31 1 1 45 4.0E-13
% 31 1 3 91 2.0E-12
% 32 0 0 59 7.0E-14
% 32 1 1 65 5.0E-13
% 32 1 3 97 2.0E-13
% 33 0 0 80 6.0E-13
% 33 1 1 67 5.0E-13
% 33 1 3 76 3.0E-12
% 63 0 0 350 5.0E-12
% 63 1 1 215 6.0E-13
% 63 1 3 412 1.0E-11
% 64 0 0 264 1.0E-13
% 64 1 1 287 3.0E-12
% 64 1 3 421 3.0E-13
% 65 0 0 338 2.0E-12
% 65 1 1 292 5.0E-13
% 65 1 3 329 1.0E-11
%
% Portability American National Standards Institute Fortran.
% The machine dependent constant PI is defined in
% function PIMACH.
%
% Required COS
% Resident
if isempty(mhm1), mhm1=0; end;
if isempty(mhmi), mhmi=0; end;
if isempty(mhpi), mhpi=0; end;
if isempty(modd), modd=0; end;
if isempty(mp), mp=0; end;
if isempty(mskip), mskip=0; end;
if isempty(nby2), nby2=0; end;
if isempty(np), np=0; end;
% Routines
%
% Reference Sweet, R., 'A Cyclic Reduction Algorithm for
% Solving Block Tridiagonal Systems Of Arbitrary
% Dimensions,' SIAM J. on Numer. Anal.,
% 14(SEPT., 1977), PP. 706-720.
%
% * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
%
%***REFERENCES R. Sweet, A cyclic reduction algorithm for solving
% block tridiagonal systems of arbitrary dimensions,
% SIAM Journal on Numerical Analysis 14, (September
% 1977), pp. 706-720.
%***ROUTINES CALLED CMPOSD, CMPOSN, CMPOSP
%***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 CMGNBN
%
%
if isempty(a1), a1=0; end;
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,[]);
b_shape=size(b);b=reshape(b,1,[]);
a_shape=size(a);a=reshape(a,1,[]);
c_shape=size(c);c=reshape(c,1,[]);
%***FIRST EXECUTABLE STATEMENT CMGNBN
ierror = 0;
if( m<=2 )
ierror = 1;
end;
if( n<=2 )
ierror = 2;
end;
if( idimy<m )
ierror = 3;
end;
if( nperod<0 || nperod>4 )
ierror = 4;
end;
if( mperod<0 || mperod>1 )
ierror = 5;
end;
if( mperod==1 )
if( abs(a(1))~=0. && abs(c(m))~=0. )
ierror = 7;
end;
else;
for i = 2 : m;
if( abs(a(i)-c(1))~=0. )
ierror = 6;
break;
elseif( abs(c(i)-c(1))~=0. ) ;
ierror = 6;
break;
elseif( abs(b(i)-b(1))~=0. ) ;
ierror = 6;
break;
end;
end;
end;
if( ierror~=0 )
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;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
return;
end;
iwba = fix(m + 1);
iwbb = fix(iwba + m);
iwbc = fix(iwbb + m);
iwb2 = fix(iwbc + m);
iwb3 = fix(iwb2 + m);
iww1 = fix(iwb3 + m);
iww2 = fix(iww1 + m);
iww3 = fix(iww2 + m);
iwd = fix(iww3 + m);
iwtcos = fix(iwd + m);
iwp = fix(iwtcos + 4.*n);
for i = 1 : m;
k = fix(iwba + i - 1);
w(k) = -a(i);
k = fix(iwbc + i - 1);
w(k) = -c(i);
k = fix(iwbb + i - 1);
w(k) = 2. - b(i);
for j = 1 : n;
y(i,j) = -y(i,j);
end; j = fix(n+1);
end; i = fix(m+1);
mp = fix(mperod + 1);
np = fix(nperod + 1);
gt200=0;
gt300=0;
gt400=0;
gt500=0;
if( mp~=1 )
gt400=1;
end;
while (1);
if(gt500==0)
if(gt400==0)
if(gt300==0)
if(gt200==0)
if( np==2 )
[m,n,dumvar3,dumvar4,dumvar5,dumvar6,y,idimy,w,dumvar10,dumvar11,dumvar12,dumvar13]=cmposd(m,n,1,w(sub2ind(size(w),max(iwba,1)):end),w(sub2ind(size(w),max(iwbb,1)):end),w(sub2ind(size(w),max(iwbc,1)):end),y,idimy,w,w(sub2ind(size(w),max(iww1,1)):end),w(sub2ind(size(w),max(iwd,1)):end),w(sub2ind(size(w),max(iwtcos,1)):end),w(sub2ind(size(w),max(iwp,1)):end)); dumvar4i=find((w(sub2ind(size(w),max(iwba,1)):end))~=(dumvar4));dumvar5i=find((w(sub2ind(size(w),max(iwbb,1)):end))~=(dumvar5));dumvar6i=find((w(sub2ind(size(w),max(iwbc,1)):end))~=(dumvar6));dumvar10i=find((w(sub2ind(size(w),max(iww1,1)):end))~=(dumvar10));dumvar11i=find((w(sub2ind(size(w),max(iwd,1)):end))~=(dumvar11));dumvar12i=find((w(sub2ind(size(w),max(iwtcos,1)):end))~=(dumvar12));dumvar13i=find((w(sub2ind(size(w),max(iwp,1)):end))~=(dumvar13)); w(iwba-1+dumvar4i)=dumvar4(dumvar4i); w(iwbb-1+dumvar5i)=dumvar5(dumvar5i); w(iwbc-1+dumvar6i)=dumvar6(dumvar6i); w(iww1-1+dumvar10i)=dumvar10(dumvar10i); w(iwd-1+dumvar11i)=dumvar11(dumvar11i); w(iwtcos-1+dumvar12i)=dumvar12(dumvar12i); w(iwp-1+dumvar13i)=dumvar13(dumvar13i);
elseif( np==3 ) ;
[m,n,dumvar3,dumvar4,dumvar5,dumvar6,dumvar7,y,idimy,w,dumvar11,dumvar12,dumvar13,dumvar14,dumvar15,dumvar16,dumvar17,dumvar18]=cmposn(m,n,1,2,w(sub2ind(size(w),max(iwba,1)):end),w(sub2ind(size(w),max(iwbb,1)):end),w(sub2ind(size(w),max(iwbc,1)):end),y,idimy,w,w(sub2ind(size(w),max(iwb2,1)):end),w(sub2ind(size(w),max(iwb3,1)):end),w(sub2ind(size(w),max(iww1,1)):end),w(sub2ind(size(w),max(iww2,1)):end),w(sub2ind(size(w),max(iww3,1)):end),w(sub2ind(size(w),max(iwd,1)):end),w(sub2ind(size(w),max(iwtcos,1)):end),w(sub2ind(size(w),max(iwp,1)):end)); dumvar5i=find((w(sub2ind(size(w),max(iwba,1)):end))~=(dumvar5));dumvar6i=find((w(sub2ind(size(w),max(iwbb,1)):end))~=(dumvar6));dumvar7i=find((w(sub2ind(size(w),max(iwbc,1)):end))~=(dumvar7));dumvar11i=find((w(sub2ind(size(w),max(iwb2,1)):end))~=(dumvar11));dumvar12i=find((w(sub2ind(size(w),max(iwb3,1)):end))~=(dumvar12));dumvar13i=find((w(sub2ind(size(w),max(iww1,1)):end))~=(dumvar13));dumvar14i=find((w(sub2ind(size(w),max(iww2,1)):end))~=(dumvar14));dumvar15i=find((w(sub2ind(size(w),max(iww3,1)):end))~=(dumvar15));dumvar16i=find((w(sub2ind(size(w),max(iwd,1)):end))~=(dumvar16));dumvar17i=find((w(sub2ind(size(w),max(iwtcos,1)):end))~=(dumvar17));dumvar18i=find((w(sub2ind(size(w),max(iwp,1)):end))~=(dumvar18)); w(iwba-1+dumvar5i)=dumvar5(dumvar5i); w(iwbb-1+dumvar6i)=dumvar6(dumvar6i); w(iwbc-1+dumvar7i)=dumvar7(dumvar7i); w(iwb2-1+dumvar11i)=dumvar11(dumvar11i); w(iwb3-1+dumvar12i)=dumvar12(dumvar12i); w(iww1-1+dumvar13i)=dumvar13(dumvar13i); w(iww2-1+dumvar14i)=dumvar14(dumvar14i); w(iww3-1+dumvar15i)=dumvar15(dumvar15i); w(iwd-1+dumvar16i)=dumvar16(dumvar16i); w(iwtcos-1+dumvar17i)=dumvar17(dumvar17i); w(iwp-1+dumvar18i)=dumvar18(dumvar18i);
elseif( np==4 ) ;
[m,n,dumvar3,dumvar4,dumvar5,dumvar6,dumvar7,y,idimy,w,dumvar11,dumvar12,dumvar13,dumvar14,dumvar15,dumvar16,dumvar17,dumvar18]=cmposn(m,n,1,1,w(sub2ind(size(w),max(iwba,1)):end),w(sub2ind(size(w),max(iwbb,1)):end),w(sub2ind(size(w),max(iwbc,1)):end),y,idimy,w,w(sub2ind(size(w),max(iwb2,1)):end),w(sub2ind(size(w),max(iwb3,1)):end),w(sub2ind(size(w),max(iww1,1)):end),w(sub2ind(size(w),max(iww2,1)):end),w(sub2ind(size(w),max(iww3,1)):end),w(sub2ind(size(w),max(iwd,1)):end),w(sub2ind(size(w),max(iwtcos,1)):end),w(sub2ind(size(w),max(iwp,1)):end)); dumvar5i=find((w(sub2ind(size(w),max(iwba,1)):end))~=(dumvar5));dumvar6i=find((w(sub2ind(size(w),max(iwbb,1)):end))~=(dumvar6));dumvar7i=find((w(sub2ind(size(w),max(iwbc,1)):end))~=(dumvar7));dumvar11i=find((w(sub2ind(size(w),max(iwb2,1)):end))~=(dumvar11));dumvar12i=find((w(sub2ind(size(w),max(iwb3,1)):end))~=(dumvar12));dumvar13i=find((w(sub2ind(size(w),max(iww1,1)):end))~=(dumvar13));dumvar14i=find((w(sub2ind(size(w),max(iww2,1)):end))~=(dumvar14));dumvar15i=find((w(sub2ind(size(w),max(iww3,1)):end))~=(dumvar15));dumvar16i=find((w(sub2ind(size(w),max(iwd,1)):end))~=(dumvar16));dumvar17i=find((w(sub2ind(size(w),max(iwtcos,1)):end))~=(dumvar17));dumvar18i=find((w(sub2ind(size(w),max(iwp,1)):end))~=(dumvar18)); w(iwba-1+dumvar5i)=dumvar5(dumvar5i); w(iwbb-1+dumvar6i)=dumvar6(dumvar6i); w(iwbc-1+dumvar7i)=dumvar7(dumvar7i); w(iwb2-1+dumvar11i)=dumvar11(dumvar11i); w(iwb3-1+dumvar12i)=dumvar12(dumvar12i); w(iww1-1+dumvar13i)=dumvar13(dumvar13i); w(iww2-1+dumvar14i)=dumvar14(dumvar14i); w(iww3-1+dumvar15i)=dumvar15(dumvar15i); w(iwd-1+dumvar16i)=dumvar16(dumvar16i); w(iwtcos-1+dumvar17i)=dumvar17(dumvar17i); w(iwp-1+dumvar18i)=dumvar18(dumvar18i);
elseif( np==5 ) ;
%
% REVERSE COLUMNS WHEN NPEROD = 4
%
irev = 1;
nby2 = fix(fix(n./2));
gt500=1;
%to 500
continue;
else;
[m,n,dumvar3,dumvar4,dumvar5,y,idimy,w,dumvar9,dumvar10,dumvar11,dumvar12,dumvar13,dumvar14,dumvar15,dumvar16]=cmposp(m,n,w(sub2ind(size(w),max(iwba,1)):end),w(sub2ind(size(w),max(iwbb,1)):end),w(sub2ind(size(w),max(iwbc,1)):end),y,idimy,w,w(sub2ind(size(w),max(iwb2,1)):end),w(sub2ind(size(w),max(iwb3,1)):end),w(sub2ind(size(w),max(iww1,1)):end),w(sub2ind(size(w),max(iww2,1)):end),w(sub2ind(size(w),max(iww3,1)):end),w(sub2ind(size(w),max(iwd,1)):end),w(sub2ind(size(w),max(iwtcos,1)):end),w(sub2ind(size(w),max(iwp,1)):end)); dumvar3i=find((w(sub2ind(size(w),max(iwba,1)):end))~=(dumvar3));dumvar4i=find((w(sub2ind(size(w),max(iwbb,1)):end))~=(dumvar4));dumvar5i=find((w(sub2ind(size(w),max(iwbc,1)):end))~=(dumvar5));dumvar9i=find((w(sub2ind(size(w),max(iwb2,1)):end))~=(dumvar9));dumvar10i=find((w(sub2ind(size(w),max(iwb3,1)):end))~=(dumvar10));dumvar11i=find((w(sub2ind(size(w),max(iww1,1)):end))~=(dumvar11));dumvar12i=find((w(sub2ind(size(w),max(iww2,1)):end))~=(dumvar12));dumvar13i=find((w(sub2ind(size(w),max(iww3,1)):end))~=(dumvar13));dumvar14i=find((w(sub2ind(size(w),max(iwd,1)):end))~=(dumvar14));dumvar15i=find((w(sub2ind(size(w),max(iwtcos,1)):end))~=(dumvar15));dumvar16i=find((w(sub2ind(size(w),max(iwp,1)):end))~=(dumvar16)); w(iwba-1+dumvar3i)=dumvar3(dumvar3i); w(iwbb-1+dumvar4i)=dumvar4(dumvar4i); w(iwbc-1+dumvar5i)=dumvar5(dumvar5i); w(iwb2-1+dumvar9i)=dumvar9(dumvar9i); w(iwb3-1+dumvar10i)=dumvar10(dumvar10i); w(iww1-1+dumvar11i)=dumvar11(dumvar11i); w(iww2-1+dumvar12i)=dumvar12(dumvar12i); w(iww3-1+dumvar13i)=dumvar13(dumvar13i); w(iwd-1+dumvar14i)=dumvar14(dumvar14i); w(iwtcos-1+dumvar15i)=dumvar15(dumvar15i); w(iwp-1+dumvar16i)=dumvar16(dumvar16i);
end;
%gt200
end;
gt200=0;
ipstor = fix(real(w(iww1)));
irev = 2;
if( nperod==4 )
gt500=1;
continue;
end;
%gt300
end;
gt300=0;
%to 600
if( mp==1 )
break;
end;
if( mp==2 )
w(1) = complex(real(ipstor+iwp-1),0.);
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;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
return;
end;
%gt400
end;
gt400=0;
%
% REORDER UNKNOWNS WHEN MP =0
%
mh =fix(fix((m+1)./2));
mhm1 = fix(mh - 1);
modd = 1;
if( mh.*2==m )
modd = 2;
end;
for j = 1 : n;
for i = 1 : mhm1;
mhpi = fix(mh + i);
mhmi = fix(mh - i);
w(i) = y(mhmi,j) - y(mhpi,j);
w(mhpi) = y(mhmi,j) + y(mhpi,j);
end; i = fix(mhm1+1);
w(mh) = 2..*y(mh,j);
if( modd~=1 )
w(m) = 2..*y(m,j);
end;
for i = 1 : m;
y(i,j) = w(i);
end; i = fix(m+1);
end; j = fix(n+1);
k = fix(iwbc + mhm1 - 1);
i = fix(iwba + mhm1);
w(k) =complex(0.,0.);
w(i) =complex(0.,0.);
w(k+1) = 2..*w(k+1);
if( modd==2 )
w(iwbb-1) = w(k+1);
else;
k = fix(iwbb + mhm1 - 1);
w(k) = w(k) - w(i-1);
w(iwbc-1) = w(iwbc-1) + w(iwbb-1);
end;
%to 100
continue;
%gt500
end;
gt500=0;
for j = 1 : nby2;
mskip = fix(n + 1 - j);
for i = 1 : m;
a1 = y(i,j);
y(i,j) = y(i,mskip);
y(i,mskip) = a1;
end; i = fix(m+1);
end; j = fix(nby2+1);
if( irev==1 )
[m,n,dumvar3,dumvar4,dumvar5,dumvar6,dumvar7,y,idimy,w,dumvar11,dumvar12,dumvar13,dumvar14,dumvar15,dumvar16,dumvar17,dumvar18]=cmposn(m,n,1,2,w(sub2ind(size(w),max(iwba,1)):end),w(sub2ind(size(w),max(iwbb,1)):end),w(sub2ind(size(w),max(iwbc,1)):end),y,idimy,w,w(sub2ind(size(w),max(iwb2,1)):end),w(sub2ind(size(w),max(iwb3,1)):end),w(sub2ind(size(w),max(iww1,1)):end),w(sub2ind(size(w),max(iww2,1)):end),w(sub2ind(size(w),max(iww3,1)):end),w(sub2ind(size(w),max(iwd,1)):end),w(sub2ind(size(w),max(iwtcos,1)):end),w(sub2ind(size(w),max(iwp,1)):end)); dumvar5i=find((w(sub2ind(size(w),max(iwba,1)):end))~=(dumvar5));dumvar6i=find((w(sub2ind(size(w),max(iwbb,1)):end))~=(dumvar6));dumvar7i=find((w(sub2ind(size(w),max(iwbc,1)):end))~=(dumvar7));dumvar11i=find((w(sub2ind(size(w),max(iwb2,1)):end))~=(dumvar11));dumvar12i=find((w(sub2ind(size(w),max(iwb3,1)):end))~=(dumvar12));dumvar13i=find((w(sub2ind(size(w),max(iww1,1)):end))~=(dumvar13));dumvar14i=find((w(sub2ind(size(w),max(iww2,1)):end))~=(dumvar14));dumvar15i=find((w(sub2ind(size(w),max(iww3,1)):end))~=(dumvar15));dumvar16i=find((w(sub2ind(size(w),max(iwd,1)):end))~=(dumvar16));dumvar17i=find((w(sub2ind(size(w),max(iwtcos,1)):end))~=(dumvar17));dumvar18i=find((w(sub2ind(size(w),max(iwp,1)):end))~=(dumvar18)); w(iwba-1+dumvar5i)=dumvar5(dumvar5i); w(iwbb-1+dumvar6i)=dumvar6(dumvar6i); w(iwbc-1+dumvar7i)=dumvar7(dumvar7i); w(iwb2-1+dumvar11i)=dumvar11(dumvar11i); w(iwb3-1+dumvar12i)=dumvar12(dumvar12i); w(iww1-1+dumvar13i)=dumvar13(dumvar13i); w(iww2-1+dumvar14i)=dumvar14(dumvar14i); w(iww3-1+dumvar15i)=dumvar15(dumvar15i); w(iwd-1+dumvar16i)=dumvar16(dumvar16i); w(iwtcos-1+dumvar17i)=dumvar17(dumvar17i); w(iwp-1+dumvar18i)=dumvar18(dumvar18i);
gt200=1;
continue;
elseif( irev==2 ) ;
gt300=1;
continue;
end;
break;
end;
for j = 1 : n;
for i = 1 : mhm1;
mhmi = fix(mh - i);
mhpi = fix(mh + i);
w(mhmi) = .5.*(y(mhpi,j)+y(i,j));
w(mhpi) = .5.*(y(mhpi,j)-y(i,j));
end; i = fix(mhm1+1);
w(mh) = .5.*y(mh,j);
if( modd~=1 )
w(m) = .5.*y(m,j);
end;
for i = 1 : m;
y(i,j) = w(i);
end; i = fix(m+1);
end; j = fix(n+1);
%
% RETURN STORAGE REQUIREMENTS FOR W ARRAY.
%
w(1) = complex(real(ipstor+iwp-1),0.);
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;
b_shape=zeros(b_shape);b_shape(:)=b(1:numel(b_shape));b=b_shape;
a_shape=zeros(a_shape);a_shape(:)=a(1:numel(a_shape));a=a_shape;
c_shape=zeros(c_shape);c_shape(:)=c(1:numel(c_shape));c=c_shape;
end %subroutine cmgnbn
%DECK CMLRI
|
|
