| [n,nelt,ia,ja,a,isym,nl,il,jl,l,dinv,nu,iu,ju,u,nrow,ncol]=dsilus(n,nelt,ia,ja,a,isym,nl,il,jl,l,dinv,nu,iu,ju,u,nrow,ncol); |
function [n,nelt,ia,ja,a,isym,nl,il,jl,l,dinv,nu,iu,ju,u,nrow,ncol]=dsilus(n,nelt,ia,ja,a,isym,nl,il,jl,l,dinv,nu,iu,ju,u,nrow,ncol);
%***BEGIN PROLOGUE DSILUS
%***PURPOSE Incomplete LU Decomposition Preconditioner SLAP Set Up.
% Routine to generate the incomplete LDU decomposition of a
% matrix. The unit lower triangular factor L is stored by
% rows and the unit upper triangular factor U is stored by
% columns. The inverse of the diagonal matrix D is stored.
% No fill in is allowed.
%***LIBRARY SLATEC (SLAP)
%***CATEGORY D2E
%***TYPE doubleprecision (SSILUS-S, DSILUS-D)
%***KEYWORDS INCOMPLETE LU FACTORIZATION, ITERATIVE PRECONDITION,
% NON-SYMMETRIC LINEAR SYSTEM, SLAP, SPARSE
%***AUTHOR Greenbaum, Anne, (Courant Institute)
% Seager, Mark K., (LLNL)
% Lawrence Livermore National Laboratory
% PO BOX 808, L-60
% Livermore, CA 94550 (510) 423-3141
% seager@llnl.gov
%***DESCRIPTION
%
% *Usage:
% INTEGER N, NELT, IA(NELT), JA(NELT), ISYM
% INTEGER NL, IL(NL), JL(NL), NU, IU(NU), JU(NU)
% INTEGER NROW(N), NCOL(N)
% doubleprecision A(NELT), L(NL), DINV(N), U(NU)
%
% CALL DSILUS( N, NELT, IA, JA, A, ISYM, NL, IL, JL, L,
% $ DINV, NU, IU, JU, U, NROW, NCOL )
%
% *Arguments:
% N :IN Integer
% Order of the Matrix.
% NELT :IN Integer.
% Number of elements in arrays IA, JA, and A.
% IA :IN Integer IA(NELT).
% JA :IN Integer JA(NELT).
% A :IN doubleprecision A(NELT).
% These arrays should hold the matrix A in the SLAP Column
% format. See 'Description', below.
% ISYM :IN Integer.
% Flag to indicate symmetric storage format.
% If ISYM=0, all non-zero entries of the matrix are stored.
% If ISYM=1, the matrix is symmetric, and only the lower
% triangle of the matrix is stored.
% NL :OUT Integer.
% Number of non-zeros in the L array.
% IL :OUT Integer IL(NL).
% JL :OUT Integer JL(NL).
% L :OUT doubleprecision L(NL).
% IL, JL, L contain the unit lower triangular factor of the
% incomplete decomposition of some matrix stored in SLAP
% Row format. The Diagonal of ones *IS* stored. See
% 'DESCRIPTION', below for more details about the SLAP format.
% NU :OUT Integer.
% Number of non-zeros in the U array.
% IU :OUT Integer IU(NU).
% JU :OUT Integer JU(NU).
% U :OUT doubleprecision U(NU).
% IU, JU, U contain the unit upper triangular factor of the
% incomplete decomposition of some matrix stored in SLAP
% Column format. The Diagonal of ones *IS* stored. See
% 'Description', below for more details about the SLAP
% format.
% NROW :WORK Integer NROW(N).
% NROW(I) is the number of non-zero elements in the I-th row
% of L.
% NCOL :WORK Integer NCOL(N).
% NCOL(I) is the number of non-zero elements in the I-th
% column of U.
%
% *Description
% IL, JL, L should contain the unit lower triangular factor of
% the incomplete decomposition of the A matrix stored in SLAP
% Row format. IU, JU, U should contain the unit upper factor
% of the incomplete decomposition of the A matrix stored in
% SLAP Column format This ILU factorization can be computed by
% the DSILUS routine. The diagonals (which are all one's) are
% stored.
%
% =================== S L A P Column format ==================
%
% This routine requires that the matrix A be stored in the
% SLAP Column format. In this format the non-zeros are stored
% counting down columns (except for the diagonal entry, which
% must appear first in each 'column') and are stored in the
% doubleprecision array A. In other words, for each column
% in the matrix put the diagonal entry in A. Then put in the
% other non-zero elements going down the column (except the
% diagonal) in order. The IA array holds the row index for
% each non-zero. The JA array holds the offsets into the IA,
% A arrays for the beginning of each column. That is,
% IA(JA(ICOL)), A(JA(ICOL)) points to the beginning of the
% ICOL-th column in IA and A. IA(JA(ICOL+1)-1),
% A(JA(ICOL+1)-1) points to the end of the ICOL-th column.
% Note that we always have JA(N+1) = NELT+1, where N is the
% number of columns in the matrix and NELT is the number of
% non-zeros in the matrix.
%
% Here is an example of the SLAP Column storage format for a
% 5x5 Matrix (in the A and IA arrays '|' denotes the end of a
% column):
%
% 5x5 Matrix SLAP Column format for 5x5 matrix on left.
% 1 2 3 4 5 6 7 8 9 10 11
% |11 12 0 0 15| A: 11 21 51 | 22 12 | 33 53 | 44 | 55 15 35
% |21 22 0 0 0| IA: 1 2 5 | 2 1 | 3 5 | 4 | 5 1 3
% | 0 0 33 0 35| JA: 1 4 6 8 9 12
% | 0 0 0 44 0|
% |51 0 53 0 55|
%
% ==================== S L A P Row format ====================
%
% This routine requires that the matrix A be stored in the
% SLAP Row format. In this format the non-zeros are stored
% counting across rows (except for the diagonal entry, which
% must appear first in each 'row') and are stored in the
% doubleprecision array A. In other words, for each row in
% the matrix put the diagonal entry in A. Then put in the
% other non-zero elements going across the row (except the
% diagonal) in order. The JA array holds the column index for
% each non-zero. The IA array holds the offsets into the JA,
% A arrays for the beginning of each row. That is,
% JA(IA(IROW)),A(IA(IROW)) are the first elements of the IROW-
% th row in JA and A, and JA(IA(IROW+1)-1), A(IA(IROW+1)-1)
% are the last elements of the IROW-th row. Note that we
% always have IA(N+1) = NELT+1, where N is the number of rows
% in the matrix and NELT is the number of non-zeros in the
% matrix.
%
% Here is an example of the SLAP Row storage format for a 5x5
% Matrix (in the A and JA arrays '|' denotes the end of a row):
%
% 5x5 Matrix SLAP Row format for 5x5 matrix on left.
% 1 2 3 4 5 6 7 8 9 10 11
% |11 12 0 0 15| A: 11 12 15 | 22 21 | 33 35 | 44 | 55 51 53
% |21 22 0 0 0| JA: 1 2 5 | 2 1 | 3 5 | 4 | 5 1 3
% | 0 0 33 0 35| IA: 1 4 6 8 9 12
% | 0 0 0 44 0|
% |51 0 53 0 55|
%
%***SEE ALSO SILUR
%***REFERENCES 1. Gene Golub and Charles Van Loan, Matrix Computations,
% Johns Hopkins University Press, Baltimore, Maryland,
% 1983.
%***ROUTINES CALLED (NONE)
%***REVISION HISTORY (YYMMDD)
% 890404 DATE WRITTEN
% 890404 Previous REVISION DATE
% 890915 Made changes requested at July 1989 CML Meeting. (MKS)
% 890922 Numerous changes to prologue to make closer to SLATEC
% standard. (FNF)
% 890929 Numerous changes to reduce SP/DP differences. (FNF)
% 910411 Prologue converted to Version 4.0 format. (BAB)
% 920511 Added complete declaration section. (WRB)
% 920929 Corrected format of reference. (FNF)
% 930701 Updated CATEGORY section. (FNF, WRB)
%***end PROLOGUE DSILUS
% .. Scalar Arguments ..
% .. Array Arguments ..
% .. Local Scalars ..
persistent i ibgn icol iend indx indx1 indx2 indxc1 indxc2 indxr1 indxr2 irow itemp j jbgn jend jtemp k kc kr temp ;
if isempty(temp), temp=0; end;
if isempty(i), i=0; end;
if isempty(ibgn), ibgn=0; end;
if isempty(icol), icol=0; end;
if isempty(iend), iend=0; end;
if isempty(indx), indx=0; end;
if isempty(indx1), indx1=0; end;
if isempty(indx2), indx2=0; end;
if isempty(indxc1), indxc1=0; end;
if isempty(indxc2), indxc2=0; end;
if isempty(indxr1), indxr1=0; end;
if isempty(indxr2), indxr2=0; end;
if isempty(irow), irow=0; end;
if isempty(itemp), itemp=0; end;
if isempty(j), j=0; end;
if isempty(jbgn), jbgn=0; end;
if isempty(jend), jend=0; end;
if isempty(jtemp), jtemp=0; end;
if isempty(k), k=0; end;
if isempty(kc), kc=0; end;
if isempty(kr), kr=0; end;
%***FIRST EXECUTABLE STATEMENT DSILUS
%
% Count number of elements in each row of the lower triangle.
%
for i = 1 : n;
nrow(i) = 0;
ncol(i) = 0;
end; i = fix(n+1);
%VD$R NOCONCUR
%VD$R NOVECTOR
for icol = 1 : n;
jbgn = fix(ja(icol) + 1);
jend = fix(ja(icol+1) - 1);
if( jbgn<=jend )
for j = jbgn : jend;
if( ia(j)<icol )
ncol(icol) = fix(ncol(icol) + 1);
else;
nrow(ia(j)) = fix(nrow(ia(j)) + 1);
if( isym~=0 )
ncol(ia(j)) = fix(ncol(ia(j)) + 1);
end;
end;
end; j = fix(jend+1);
end;
end; icol = fix(n+1);
ju(1) = 1;
il(1) = 1;
for icol = 1 : n;
il(icol+1) = fix(il(icol) + nrow(icol));
ju(icol+1) = fix(ju(icol) + ncol(icol));
nrow(icol) = fix(il(icol));
ncol(icol) = fix(ju(icol));
end; icol = fix(n+1);
%
% Copy the matrix A into the L and U structures.
for icol = 1 : n;
dinv(icol) = a(ja(icol));
jbgn = fix(ja(icol) + 1);
jend = fix(ja(icol+1) - 1);
if( jbgn<=jend )
for j = jbgn : jend;
irow = fix(ia(j));
if( irow<icol )
% Part of the upper triangle.
iu(ncol(icol)) = fix(irow);
u(ncol(icol)) = a(j);
ncol(icol) = fix(ncol(icol) + 1);
else;
% Part of the lower triangle (stored by row).
jl(nrow(irow)) = fix(icol);
l(nrow(irow)) = a(j);
nrow(irow) = fix(nrow(irow) + 1);
if( isym~=0 )
% Symmetric...Copy lower triangle into upper triangle as well.
iu(ncol(irow)) = fix(icol);
u(ncol(irow)) = a(j);
ncol(irow) = fix(ncol(irow) + 1);
end;
end;
end; j = fix(jend+1);
end;
end; icol = fix(n+1);
%
% Sort the rows of L and the columns of U.
for k = 2 : n;
jbgn = fix(ju(k));
jend = fix(ju(k+1) - 1);
if( jbgn<jend )
for j = jbgn : jend - 1;
for i = j + 1 : jend;
if( iu(j)>iu(i) )
itemp = fix(iu(j));
iu(j) = fix(iu(i));
iu(i) = fix(itemp);
temp = u(j);
u(j) = u(i);
u(i) = temp;
end;
end; i = fix(jend+1);
end; j = fix(jend - 1+1);
end;
ibgn = fix(il(k));
iend = fix(il(k+1) - 1);
if( ibgn<iend )
for i = ibgn : iend - 1;
for j = i + 1 : iend;
if( jl(i)>jl(j) )
jtemp = fix(ju(i));
ju(i) = fix(ju(j));
ju(j) = fix(jtemp);
temp = l(i);
l(i) = l(j);
l(j) = temp;
end;
end; j = fix(iend+1);
end; i = fix(iend - 1+1);
end;
end; k = fix(n+1);
%
% Perform the incomplete LDU decomposition.
for i = 2 : n;
%
% I-th row of L
indx1 = fix(il(i));
indx2 = fix(il(i+1) - 1);
if( indx1<=indx2 )
for indx = indx1 : indx2;
if( indx~=indx1 )
indxr1 = fix(indx1);
indxr2 = fix(indx - 1);
indxc1 = fix(ju(jl(indx)));
indxc2 = fix(ju(jl(indx)+1) - 1);
if( indxc1<=indxc2 )
kr = fix(jl(indxr1));
while( true );
kc = fix(iu(indxc1));
if( kr>kc )
indxc1 = fix(indxc1 + 1);
if( indxc1>indxc2 )
break;
end;
elseif( kr<kc ) ;
indxr1 = fix(indxr1 + 1);
if( indxr1>indxr2 )
break;
end;
kr = fix(jl(indxr1));
elseif( kr==kc ) ;
l(indx) = l(indx) - l(indxr1).*dinv(kc).*u(indxc1);
indxr1 = fix(indxr1 + 1);
indxc1 = fix(indxc1 + 1);
if( indxr1>indxr2 || indxc1>indxc2 )
break;
end;
kr = fix(jl(indxr1));
else;
break;
end;
end;
end;
end;
l(indx) = l(indx)./dinv(jl(indx));
end;
end;
%
% I-th column of U
indx1 = fix(ju(i));
indx2 = fix(ju(i+1) - 1);
if( indx1<=indx2 )
for indx = indx1 : indx2;
if( indx~=indx1 )
indxc1 = fix(indx1);
indxc2 = fix(indx - 1);
indxr1 = fix(il(iu(indx)));
indxr2 = fix(il(iu(indx)+1) - 1);
if( indxr1<=indxr2 )
kr = fix(jl(indxr1));
while( true );
kc = fix(iu(indxc1));
if( kr>kc )
indxc1 = fix(indxc1 + 1);
if( indxc1>indxc2 )
break;
end;
elseif( kr<kc ) ;
indxr1 = fix(indxr1 + 1);
if( indxr1>indxr2 )
break;
end;
kr = fix(jl(indxr1));
elseif( kr==kc ) ;
u(indx) = u(indx) - l(indxr1).*dinv(kc).*u(indxc1);
indxr1 = fix(indxr1 + 1);
indxc1 = fix(indxc1 + 1);
if( indxr1>indxr2 || indxc1>indxc2 )
break;
end;
kr = fix(jl(indxr1));
else;
break;
end;
end;
end;
end;
u(indx) = u(indx)./dinv(iu(indx));
end;
end;
%
% I-th diagonal element
indxr1 = fix(il(i));
indxr2 = fix(il(i+1) - 1);
if( indxr1<=indxr2 )
%
indxc1 = fix(ju(i));
indxc2 = fix(ju(i+1) - 1);
if( indxc1<=indxc2 )
kr = fix(jl(indxr1));
while( true );
kc = fix(iu(indxc1));
if( kr>kc )
indxc1 = fix(indxc1 + 1);
if( indxc1>indxc2 )
break;
end;
elseif( kr<kc ) ;
indxr1 = fix(indxr1 + 1);
if( indxr1>indxr2 )
break;
end;
kr = fix(jl(indxr1));
elseif( kr==kc ) ;
dinv(i) = dinv(i) - l(indxr1).*dinv(kc).*u(indxc1);
indxr1 = fix(indxr1 + 1);
indxc1 = fix(indxc1 + 1);
if( indxr1>indxr2 || indxc1>indxc2 )
break;
end;
kr = fix(jl(indxr1));
else;
break;
end;
end;
end;
end;
end;
%
% Replace diagonal elements by their inverses.
%VD$ VECTOR
for i = 1 : n;
dinv(i) = 1.0d0./dinv(i);
end; i = fix(n+1);
%
%------------- LAST LINE OF DSILUS FOLLOWS ----------------------------
end %subroutine dsilus
%DECK DSINDG
|
|
