| [nperod,n,mperod,m,a,b,c,idimy,y,ierror,w]=genbun(nperod,n,mperod,m,a,b,c,idimy,y,ierror,w); |
function [nperod,n,mperod,m,a,b,c,idimy,y,ierror,w]=genbun(nperod,n,mperod,m,a,b,c,idimy,y,ierror,w);
persistent a1 gt i ipstor irev iwb2 iwb3 iwba iwbb iwbc iwd iwp iwtcos iww1 iww2 iww3 j k mh mhm1 mhmi mhpi modd mp mp1 mskip nby2 np ;
if isempty(a1), a1=0; end;
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(gt), gt=zeros(1,6); end;
%***BEGIN PROLOGUE GENBUN
%***PURPOSE Solve by a cyclic reduction algorithm the linear system
% of equations that results from a finite difference
% approximation to certain 2-d elliptic PDE's on a centered
% grid .
%***LIBRARY SLATEC (FISHPACK)
%***CATEGORY I2B4B
%***TYPE SINGLE PRECISION (GENBUN-S, CMGNBN-C)
%***KEYWORDS ELLIPTIC, FISHPACK, PDE, TRIDIAGONAL
%***AUTHOR Adams, J., (NCAR)
% Swarztrauber, P. N., (NCAR)
% Sweet, R., (NCAR)
%***DESCRIPTION
%
% subroutine GENBUN solves the 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. M must be greater
% than 2.
%
% A,B,C
% One-dimensional 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 GENBUN. This parameter is
% used to specify the variable dimension of Y. IDIMY must be at
% least M.
%
% Y
% A two-dimensional 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 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
% GENBUN and is returned in location W(1).
%
%
% * * * * * * 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 .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 1, 1976
% Revision
%
% Subprograms GENBUN,POISD2,POISN2,POISP2,COSGEN,MERGE,TRIX,TRI3,
% Required PIMACH
%
% Special NONE
% Conditions
%
% Common NONE
% Blocks
%
% I/O NONE
%
% Precision Single
%
% Specialist Roland Sweet
%
% Language FORTRAN
%
% History Standardized April 1, 1973
% Revised August 20,1973
% Revised January 1, 1976
%
% 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 GENBUN 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. More comprehensive timing
% charts may be found in the reference.
% 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, using doubleprecision, a right side Y was
% computed. Using this array Y subroutine GENBUN 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 36 6.0E-14
% 31 1 1 21 4.0E-13
% 31 1 3 41 3.0E-13
% 32 0 0 29 9.0E-14
% 32 1 1 32 3.0E-13
% 32 1 3 48 1.0E-13
% 33 0 0 36 9.0E-14
% 33 1 1 30 4.0E-13
% 33 1 3 34 1.0E-13
% 63 0 0 150 1.0E-13
% 63 1 1 91 1.0E-12
% 63 1 3 173 2.0E-13
% 64 0 0 122 1.0E-13
% 64 1 1 128 1.0E-12
% 64 1 3 199 6.0E-13
% 65 0 0 143 2.0E-13
% 65 1 1 120 1.0E-12
% 65 1 3 138 4.0E-13
%
% Portability American National Standards Institute Fortran.
% The machine dependent constant PI is defined in
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(mp1), mp1=0; end;
if isempty(mskip), mskip=0; end;
if isempty(nby2), nby2=0; end;
if isempty(np), np=0; end;
% function PIMACH.
%
% Required COS
% Resident
% 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 POISD2, POISN2, POISP2
%***REVISION HISTORY (YYMMDD)
% 801001 DATE WRITTEN
% 861211 REVISION DATE from Version 3.2
% 891214 Prologue converted to Version 4.0 format. (BAB)
% 920501 Reformatted the REFERENCES section. (WRB)
%***end PROLOGUE GENBUN
%
%
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 GENBUN
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( a(1)~=0. || c(m)~=0. )
ierror = 7;
end;
else;
for i = 2 : m;
if( a(i)~=c(1) )
ierror = 6;
break;
elseif( c(i)~=c(1) ) ;
ierror = 6;
break;
elseif( b(i)~=b(1) ) ;
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;
mp1 = fix(m + 1);
iwba = fix(mp1);
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);
gt(:)=0;
if( mp==1 )
gt(4)=1;
end;
while (1);
if(gt(6)==0)
if(gt(5)==0)
if(gt(4)==0)
if(gt(3)==0)
if(gt(2)==0)
if( np==2 )
[m,n,dumvar3,dumvar4,dumvar5,dumvar6,y,idimy,w,dumvar10,dumvar11,dumvar12,dumvar13]=poisd2(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]=poisn2(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]=poisn2(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));
gt(5)=1;
continue;
else;
[m,n,dumvar3,dumvar4,dumvar5,y,idimy,w,dumvar9,dumvar10,dumvar11,dumvar12,dumvar13,dumvar14,dumvar15,dumvar16]=poisp2(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;
end;
gt(2)=0;
ipstor = fix(w(iww1));
irev = 2;
if( nperod==4 )
gt(5)=1;
continue;
end;
end;
gt(3)=0;
if( mp==1 )
gt(6)=1;
continue;
end;
if( mp==2 )
break;
end;
%
% REORDER UNKNOWNS WHEN MP =0
%
end;
gt(4)=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) = 0.;
w(i) = 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;
continue;
end;
gt(5)=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]=poisn2(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);
gt(2)=1;
continue;
elseif( irev==2 ) ;
gt(3)=1;
continue;
end;
end;
gt(6)=0;
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);
break;
end;
%
% RETURN STORAGE REQUIREMENTS FOR W ARRAY.
%
w(1) = ipstor + iwp - 1;
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 genbun
%DECK GVEC
|
|
