| [n,lgmr,x,xl,zl,hes,maxlp1,q,v,r0nrm,wk,sz,jscal,jpre,msolve,nmsl,rpar,ipar,nelt,ia,ja,a,isym]=sxlcal(n,lgmr,x,xl,zl,hes,maxlp1,q,v,r0nrm,wk,sz,jscal,jpre,msolve,nmsl,rpar,ipar,nelt,ia,ja,a,isym); |
function [n,lgmr,x,xl,zl,hes,maxlp1,q,v,r0nrm,wk,sz,jscal,jpre,msolve,nmsl,rpar,ipar,nelt,ia,ja,a,isym]=sxlcal(n,lgmr,x,xl,zl,hes,maxlp1,q,v,r0nrm,wk,sz,jscal,jpre,msolve,nmsl,rpar,ipar,nelt,ia,ja,a,isym);
%***BEGIN PROLOGUE SXLCAL
%***SUBSIDIARY
%***PURPOSE Internal routine for SGMRES.
%***LIBRARY SLATEC (SLAP)
%***CATEGORY D2A4, D2B4
%***TYPE SINGLE PRECISION (SXLCAL-S, DXLCAL-D)
%***KEYWORDS GENERALIZED MINIMUM RESIDUAL, ITERATIVE PRECONDITION,
% NON-SYMMETRIC LINEAR SYSTEM, SLAP, SPARSE
%***AUTHOR Brown, Peter, (LLNL), pnbrown@llnl.gov
% Hindmarsh, Alan, (LLNL), alanh@llnl.gov
% Seager, Mark K., (LLNL), seager@llnl.gov
% Lawrence Livermore National Laboratory
% PO Box 808, L-60
% Livermore, CA 94550 (510) 423-3141
%***DESCRIPTION
% This routine computes the solution XL, the current SGMRES
% iterate, given the V(I)'s and the QR factorization of the
% Hessenberg matrix HES. This routine is only called when
% ITOL=11.
%
% *Usage:
% INTEGER N, LGMR, MAXLP1, JSCAL, JPRE, NMSL, IPAR(USER DEFINED)
% INTEGER NELT, IA(NELT), JA(NELT), ISYM
% REAL X(N), XL(N), ZL(N), HES(MAXLP1,MAXL), Q(2*MAXL),
% $ V(N,MAXLP1), R0NRM, WK(N), SZ(N), RPAR(USER DEFINED),
% $ A(NELT)
% EXTERNAL MSOLVE
%
% CALL SXLCAL(N, LGMR, X, XL, ZL, HES, MAXLP1, Q, V, R0NRM,
% $ WK, SZ, JSCAL, JPRE, MSOLVE, NMSL, RPAR, IPAR,
% $ NELT, IA, JA, A, ISYM)
%
% *Arguments:
% N :IN Integer
% The order of the matrix A, and the lengths
% of the vectors SR, SZ, R0 and Z.
% LGMR :IN Integer
% The number of iterations performed and
% the current order of the upper Hessenberg
% matrix HES.
% X :IN Real X(N)
% The current approximate solution as of the last restart.
% XL :OUT Real XL(N)
% An array of length N used to hold the approximate
% solution X(L).
% Warning: XL and ZL are the same array in the calling routine.
% ZL :IN Real ZL(N)
% An array of length N used to hold the approximate
% solution Z(L).
% HES :IN Real HES(MAXLP1,MAXL)
% The upper triangular factor of the QR decomposition
% of the (LGMR+1) by LGMR upper Hessenberg matrix whose
% entries are the scaled inner-products of A*V(*,i) and V(*,k).
% MAXLP1 :IN Integer
% MAXLP1 = MAXL + 1, used for dynamic dimensioning of HES.
% MAXL is the maximum allowable order of the matrix HES.
% Q :IN Real Q(2*MAXL)
% A real array of length 2*MAXL containing the components
% of the Givens rotations used in the QR decomposition
% of HES. It is loaded in SHEQR.
% V :IN Real V(N,MAXLP1)
% The N by(LGMR+1) array containing the LGMR
% orthogonal vectors V(*,1) to V(*,LGMR).
% R0NRM :IN Real
% The scaled norm of the initial residual for the
% current call to SPIGMR.
% WK :IN Real WK(N)
% A real work array of length N.
% SZ :IN Real SZ(N)
% A vector of length N containing the non-zero
% elements of the diagonal scaling matrix for Z.
% JSCAL :IN Integer
% A flag indicating whether arrays SR and SZ are used.
% JSCAL=0 means SR and SZ are not used and the
% algorithm will perform as if all
% SR(i) = 1 and SZ(i) = 1.
% JSCAL=1 means only SZ is used, and the algorithm
% performs as if all SR(i) = 1.
% JSCAL=2 means only SR is used, and the algorithm
% performs as if all SZ(i) = 1.
% JSCAL=3 means both SR and SZ are used.
% JPRE :IN Integer
% The preconditioner type flag.
% MSOLVE :EXT External.
% Name of the routine which solves a linear system Mz = r for
% z given r with the preconditioning matrix M (M is supplied via
% RPAR and IPAR arrays. The name of the MSOLVE routine must
% be declared external in the calling program. The calling
% sequence to MSOLVE is:
% CALL MSOLVE(N, R, Z, NELT, IA, JA, A, ISYM, RPAR, IPAR)
% Where N is the number of unknowns, R is the right-hand side
% vector and Z is the solution upon return. NELT, IA, JA, A and
% ISYM are defined as below. RPAR is a real array that can be
% used to pass necessary preconditioning information and/or
% workspace to MSOLVE. IPAR is an integer work array for the
% same purpose as RPAR.
% NMSL :IN Integer
% The number of calls to MSOLVE.
% RPAR :IN Real RPAR(USER DEFINED)
% Real workspace passed directly to the MSOLVE routine.
% IPAR :IN Integer IPAR(USER DEFINED)
% Integer workspace passed directly to the MSOLVE routine.
% NELT :IN Integer
% The length of arrays IA, JA and A.
% IA :IN Integer IA(NELT)
% An integer array of length NELT containing matrix data.
% It is passed directly to the MATVEC and MSOLVE routines.
% JA :IN Integer JA(NELT)
% An integer array of length NELT containing matrix data.
% It is passed directly to the MATVEC and MSOLVE routines.
% A :IN Real A(NELT)
% A real array of length NELT containing matrix data.
% It is passed directly to the MATVEC and MSOLVE routines.
% ISYM :IN Integer
% A flag to indicate symmetric matrix storage.
% If ISYM=0, all non-zero entries of the matrix are
% stored. If ISYM=1, the matrix is symmetric and
% only the upper or lower triangular part is stored.
%
%***SEE ALSO SGMRES
%***ROUTINES CALLED SAXPY, SCOPY, SHELS
%***REVISION HISTORY (YYMMDD)
% 871001 DATE WRITTEN
% 881213 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)
% 910502 Removed MSOLVE from ROUTINES CALLED list. (FNF)
% 910506 Made subsidiary to SGMRES. (FNF)
% 920511 Added complete declaration section. (WRB)
%***end PROLOGUE SXLCAL
% The following is for optimized compilation on LLNL/LTSS Crays.
%LLL. OPTIMIZE
% .. Scalar Arguments ..
% .. Array Arguments ..
persistent i k ll llp1 ;
hes_shape=size(hes);hes=reshape([hes(:).',zeros(1,ceil(numel(hes)./prod([maxlp1])).*prod([maxlp1])-numel(hes))],maxlp1,[]);
q_shape=size(q);q=reshape(q,1,[]);
rpar_shape=size(rpar);rpar=reshape(rpar,1,[]);
sz_shape=size(sz);sz=reshape(sz,1,[]);
v_shape=size(v);v=reshape([v(:).',zeros(1,ceil(numel(v)./prod([n])).*prod([n])-numel(v))],n,[]);
ipar_shape=size(ipar);ipar=reshape(ipar,1,[]);
% .. subroutine Arguments ..
% .. Local Scalars ..
if isempty(i), i=0; end;
if isempty(k), k=0; end;
if isempty(ll), ll=0; end;
if isempty(llp1), llp1=0; end;
% .. External Subroutines ..
%***FIRST EXECUTABLE STATEMENT SXLCAL
ll = fix(lgmr);
llp1 = fix(ll + 1);
for k = 1 : llp1;
wk(k) = 0;
end; k = fix(llp1+1);
wk(1) = r0nrm;
[hes,maxlp1,ll,q,wk]=shels(hes,maxlp1,ll,q,wk);
for k = 1 : n;
zl(k) = 0;
end; k = fix(n+1);
for i = 1 : ll;
[n,wk(i),v(sub2ind(size(v),1,i):end),dumvar4,zl]=saxpy(n,wk(i),v(sub2ind(size(v),1,i):end),1,zl,1);
end; i = fix(ll+1);
if((jscal==1) ||(jscal==3) )
for k = 1 : n;
zl(k) = zl(k)./sz(k);
end; k = fix(n+1);
end;
if( jpre>0 )
[n,zl,dumvar3,wk]=scopy(n,zl,1,wk,1);
[n,wk,zl,nelt,ia,ja,a,isym,rpar,ipar]=msolve(n,wk,zl,nelt,ia,ja,a,isym,rpar,ipar);
nmsl = fix(nmsl + 1);
end;
% calculate XL from X and ZL.
for k = 1 : n;
xl(k) = x(k) + zl(k);
end; k = fix(n+1);
%------------- LAST LINE OF SXLCAL FOLLOWS ----------------------------
hes_shape=zeros(hes_shape);hes_shape(:)=hes(1:numel(hes_shape));hes=hes_shape;
q_shape=zeros(q_shape);q_shape(:)=q(1:numel(q_shape));q=q_shape;
rpar_shape=zeros(rpar_shape);rpar_shape(:)=rpar(1:numel(rpar_shape));rpar=rpar_shape;
sz_shape=zeros(sz_shape);sz_shape(:)=sz(1:numel(sz_shape));sz=sz_shape;
v_shape=zeros(v_shape);v_shape(:)=v(1:numel(v_shape));v=v_shape;
ipar_shape=zeros(ipar_shape);ipar_shape(:)=ipar(1:numel(ipar_shape));ipar=ipar_shape;
end
%DECK TEVLC
|
|
