# Date: Fri, 16 Sep 1994 12:16:27 +0200 # From: "Michela Redivo Zaglia -tel.+39 49 8287625" <ELEN@ELETT1.DEI.UNIPD.IT> # Subject: NA5 package # This is a shell archive. # Remove everything above and including the cut line. # Then run the rest of the file through sh. #----cut here-----cut here-----cut here-----cut here----# #!/bin/sh # shar: Shell Archiver # Run the following text with /bin/sh to create: # disk$site:[elen.cgs]aaaread.me # disk$site:[elen.cgs]makefile # disk$site:[elen.cgs]bsmrzs.f # disk$site:[elen.cgs]chksig.f # disk$site:[elen.cgs]compo1.f # disk$site:[elen.cgs]compo2.f # disk$site:[elen.cgs]compol.f # disk$site:[elen.cgs]compp1.f # disk$site:[elen.cgs]cresol.f # disk$site:[elen.cgs]exa1.f # disk$site:[elen.cgs]exa2.f # disk$site:[elen.cgs]exa3.f # disk$site:[elen.cgs]exa4.f # disk$site:[elen.cgs]exa5.f # disk$site:[elen.cgs]fh.f # disk$site:[elen.cgs]gameta.f # disk$site:[elen.cgs]gsolvd.f # disk$site:[elen.cgs]matvec.f # disk$site:[elen.cgs]pinner.f # disk$site:[elen.cgs]polmul.f # disk$site:[elen.cgs]rxsumm.f # disk$site:[elen.cgs]solshf.f # disk$site:[elen.cgs]ssumm.f # disk$site:[elen.cgs]stomhf.f # disk$site:[elen.cgs]storhf.f # disk$site:[elen.cgs]sunorm.f # disk$site:[elen.cgs]zsumm.f # This archive created: Fri Sep 16 11:19:31 1994 echo shar: extracting aaaread.me sed 's/^X//' << \SHAR_EOF > aaaread.me X *************************************************************************** X * All the software contained in this library is protected by copyright. * X * Permission to use, copy, modify, and distribute this software for any * X * purpose without fee is hereby granted, provided that this entire notice * X * is included in all copies of any software which is or includes a copy * X * or modification of this software and in all copies of the supporting * X * documentation for such software. * X *************************************************************************** X * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED * X * WARRANTY. IN NO EVENT, NEITHER THE AUTHORS NOR THE PUBLISHER, BE LIABLE * X * FOR ANY ERROR IN THE SOFTWARE, ANY MISUSE OF IT OR ANY DAMAGE ARISING * X * OUT OF ITS USE. THE ENTIRE RISK OF USING THE SOFTWARE LIES WITH THE * X * PARTY DOING SO. * X *************************************************************************** X * ANY USE OF THE SOFTWARE CONSTITUTES ACCEPTANCE OF THE TERMS OF THE * X * ABOVE STATEMENT. * X *************************************************************************** X X * AUTHORS: X X C. BREZINSKI X UNIVERSITY OF LILLE - FRANCE X X M. REDIVO ZAGLIA X UNIVERSITY OF PADOVA - ITALY X X * REFERENCE: X X TREATMENT OF NEAR-BREAKDOWN IN THE CGS ALGORITHM X NUMERICAL ALGORITHMS 7 (1994) PP. 33-73. X X * SOFTWARE REVISION DATE: X X FEBRUARY, 1993 X X * MODULES: X X - MAIN PROGRAMS EXA1, EXA2, EXA3, EXA4, EXA5 X X - SUBROUTINES BSMRZS, CHKSIG, COMPO1, COMPO2, COMPOL, X COMPP1, CRESOL, GSOLVD, POLMUL, SOLSHF, X STOMHF, STORHF. X X - DOUBLE PRECISION FUNCTIONS FH, GAMETA, RXSUMM, SSUMM, SUNORM, ZSUMM. X X X * MODULES GIVEN IN NETLIB (NA1 FROM NUMERALGO): X X - SUBROUTINE MATVEC. X X - DOUBLE PRECISION FUNCTION PINNER. X X * REMARK: A makefile IS ALSO PROVIDED X SHAR_EOF if test 2189 -ne "`wc -c aaaread.me`" then echo shar: error transmitting aaaread.me '(should have been 2189 characters)' fi echo shar: extracting makefile sed 's/^X//' << \SHAR_EOF > makefile X.KEEP_STATE: XFC=f77 XOBJS1 = bsmrzs.o chksig.o compo1.o compo2.o compol.o compp1.o \ X cresol.o gsolvd.o polmul.o solshf.o stomhf.o storhf.o \ X fh.o gameta.o rxsumm.o ssumm.o sunorm.o zsumm.o matvec.o pinner.o X XFLIBS = -lm X Xall : exa1 exa2 exa3 exa4 exa5 X Xexa1 : $(OBJS1) exa1.o X $(FC) $(OBJS1) exa1.o -o exa1 $(FLIBS) X Xexa2 : $(OBJS1) exa2.o X $(FC) $(OBJS1) exa2.o -o exa2 $(FLIBS) X Xexa3 : $(OBJS1) exa3.o X $(FC) $(OBJS1) exa3.o -o exa3 $(FLIBS) X Xexa4 : $(OBJS1) exa4.o X $(FC) $(OBJS1) exa4.o -o exa4 $(FLIBS) X Xexa5 : $(OBJS1) exa5.o X $(FC) $(OBJS1) exa5.o -o exa5 $(FLIBS) X X X SHAR_EOF if test 597 -ne "`wc -c makefile`" then echo shar: error transmitting makefile '(should have been 597 characters)' fi echo shar: extracting bsmrzs.f sed 's/^X//' << \SHAR_EOF > bsmrzs.f XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC SUBROUTINE NAME - BSMRZS + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC DESCRIPTION: + XC + XC COMPUTES THE SOLUTION OF A SYSTEM OF LINEAR EQUATIONS BY THE + XC BSMRZS ALGORITHM WHICH AVOIDS BREAKDOWN AND NEAR-BREAKDOWN IN THE + XC CGS ALGORITHM. + XC + XC USAGE: + XC CALL BSMRZS (N, A, AB, Y, X, R, NORM, MAXBCK, MAXN, V, W, U, + XC LMAT, RMAT, ETA, ETAP, GAMMA, GAMMAP, D, C, P, P1, + XC PWK, IR, NK, EPS, EPS1, EPS2, INIT, IFLAG, IER) + XC + XC ARGUMENTS: + XC + XC N - INPUT INTEGER VALUE, DIMENSION OF THE SYSTEM. + XC + XC A - INPUT REAL MATRIX OF DIMENSION (N,N) CONTAINING THE + XC MATRIX OF THE SYSTEM. + XC + XC AB - INPUT/WORKING REAL VECTOR OF DIMENSION N CONTAINING, + XC BEFORE THE FIRST CALL, THE RIGHT HAND SIDE OF THE SYSTEM. + XC THE VECTOR IS USED AS A WORKING AREA IN THE NEXT CALLS + XC AND THEN THE INPUT VALUES ARE ALWAYS DESTROYED BY THE + XC SUBROUTINE. + XC + XC Y - INPUT/OUTPUT REAL VECTOR OF DIMENSION N CONTAINING, + XC BEFORE THE FIRST CALL, THE AUXILIARY VECTOR y. + XC IF IFLAG .EQ. 0 THE VECTOR Y IS SET TO r = b - A x + XC 0 0 + XC BY THE SUBROUTINE DURING THE FIRST CALL. + XC + XC X - INPUT/OUTPUT REAL VECTOR OF DIMENSION N CONTAINING + XC AFTER THE (k+1)-TH CALL, THE SOLUTION x . + XC k+1 + XC BEFORE THE FIRST CALL, X MUST CONTAIN x . + XC 0 + XC + XC R - OUTPUT REAL VECTOR OF DIMENSION N CONTAINING, AFTER THE + XC (k+1)-TH CALL, THE RESIDUAL VECTOR r . + XC k+1 + XC IF, AT THE FIRST CALL, THE NORM OF THE VECTOR R = r IS + XC 0 + XC LESS THAN EPS2 THEN THE VECTOR R CONTAINS r IN OUTPUT. + XC 0 + XC + XC NORM - OUTPUT REAL VALUE CONTAINING, AFTER THE (k+1)-TH CALL, + XC THE NORM OF THE RESIDUAL VECTOR r . + XC k+1 + XC + XC MAXBCK - INPUT INTEGER VALUE. IT REPRESENTS THE MAXIMUM VALUE + XC ALLOWED FOR THE JUMP m . IT MUST BE GREATER OR EQUAL TO + XC k + XC 2 AND LESS OR EQUAL TO MAXN. IT IS USED TO CONTROL THE + XC DIMENSION OF THE WORKING ARRAYS. + XC + XC MAXN - INPUT INTEGER VALUE. IT REPRESENTS THE MAXIMUM VALUE + XC ALLOWED, GREATER OR EQUAL TO N, FOR n + m . + XC k k + XC + XC V - WORKING REAL MATRIX OF DIMENSION (N,0:2*MAXBCK). + XC V(I), I=0,2*MAXBCK, CONTAINS A**I * z . + XC k + XC + XC W - WORKING REAL MATRIX OF DIMENSION (N,0:2*MAXBCK-2). + XC W(I), I=0,2*MAXBCK-2, CONTAINS A**I * r . + XC k + XC + XC U - WORKING REAL MATRIX OF DIMENSION (N,0:2*MAXBCK-1). + XC W(I), I=0,2*MAXBCK-1, CONTAINS A**I * s . + XC k + XC + XC LMAT - WORKING REAL VECTOR OF DIMENSION: + XC (8*MAXBCK) IF MAXBCK <= 2 + XC ((2*MAXBCK-1)*(2*MAXBCK-1)) IF MAXBCK > 2 + XC + XC RMAT - WORKING REAL MATRIX OF DIMENSION (2,0:2*MAXBCK-1). + XC + XC ETA - WORKING REAL VECTOR CONTAINING THE ETA . + XC i + XC + XC ETAP - WORKING REAL VECTOR CONTAINING THE ETA'. + XC i + XC + XC GAMMA - WORKING REAL VECTOR CONTAINING THE GAMMA . + XC i + XC + XC GAMMAP - WORKING REAL VECTOR CONTAINING THE GAMMA'. + XC i + XC + XC D - WORKING REAL VECTOR OF DIMENSION (0:2*MAXBCK-1). + XC + XC C - WORKING REAL VECTOR OF DIMENSION (0:2*MAXBCK-1). + XC + XC P - INPUT/OUTPUT REAL VECTOR OF DIMENSION (0:MAXN). + XC IT CONTAINS, BEFORE THE (k+1)-TH CALL, THE COEFFICIENTS + XC OF THE POLYNOMIAL P OF DEGREE NK, AND AFTER THE (k+1)-TH + XC k + XC CALL, THE COEFFICIENTS OF THE POLYNOMIAL P OF DEGREE + XC k+1 + XC NK+MK. THE POLYNOMIAL IS INITIALIZED BY THE SUBROUTINE + XC DURING THE FIRST CALL. + XC + XC P1 - INPUT/OUTPUT REAL VECTOR OF DIMENSION (0:MAXN). + XC IT CONTAINS, BEFORE THE (k+1)-TH CALL, THE COEFFICIENTS + XC (1) + XC OF THE POLYNOMIAL P OF DEGREE NK, AND AFTER THE (k+1)-TH + XC k (1) + XC CALL, THE COEFFICIENTS OF THE POLYNOMIAL P OF DEGREE + XC k+1 + XC NK+MK. THE POLYNOMIAL IS INITIALIZED BY THE SUBROUTINE + XC DURING THE FIRST CALL. + XC + XC PWK - WORKING REAL VECTOR OF DIMENSION (0:2*MAXN+1). + XC + XC IR - INPUT ROWS DIMENSION EXACTLY AS SPECIFIED IN THE + XC DIMENSION STATEMENTS IN THE CALLING PROGRAM FOR + XC THE MATRICES A, V, W AND U. + XC + XC NK - OUTPUT INTEGER VALUE, DIMENSION OF THE INTERMEDIATE + XC KRYLOV SUBSPACE. + XC + XC EPS - INPUT REAL VALUE USED FOR TESTING THE NEAR BREAKDOWN. + XC THE SUBROUTINE DOES NOT CONTROL WHETHER OR NOT EPS IS + XC A NEGATIVE OR ZERO REAL NUMBER. + XC + XC EPS1 - INPUT REAL VALUE USED FOR TESTING THE PIVOTS IN GAUSSIAN + XC ELIMINATION FOR SOLVING THE AUXILIARY SYSTEMS. + XC IF DABS(X) .LE. EPS1, THEN X IS CONSIDERED TO BE ZERO. + XC THE SUBROUTINE DOES NOT CONTROL WHETHER OR NOT EPS1 IS + XC A NEGATIVE OR ZERO REAL NUMBER. + XC + XC EPS2 - INPUT REAL VALUE USED FOR TESTING THE FINAL SOLUTION. + XC THE FINAL SOLUTION IS OBTAINED WHEN THE NORM OF THE + XC RESIDUAL VECTOR IS SMALLER THAN EPS2. + XC THE SUBROUTINE DOES NOT CONTROL WHETHER OR NOT EPS2 IS + XC A NEGATIVE OR ZERO REAL NUMBER. + XC + XC INIT - INPUT/OUTPUT INTEGER TO BE SET TO ZERO BEFORE THE + XC FIRST CALL OF THE SUBROUTINE. ITS VALUE IS CHANGED + XC TO 1 BY THE SUBROUTINE DURING THE FIRST CALL. + XC FOR A NEW APPLICATION OF THE METHOD, INIT MUST SET + XC AGAIN TO ZERO. + XC + XC IFLAG - INPUT INTEGER. + XC IF IFLAG .EQ. 0 THEN THE AUXILIARY VECTOR y IS CHOSEN + XC TO COINCIDE WITH THE VECTOR z = r . + XC 0 0 + XC IF IFLAG .NE. 0 THEN THE USER MUST DEFINE y IN THE MAIN + XC PROGRAM BEFORE THE FIRST CALL. + XC + XC IER - OUTPUT INDEX WARNING/ERROR. + XC IER = 100 CALL OF THE SUBROUTINE WITH A NON ZERO IER + XC VALUE. + XC IER = 200 THE NORM OF THE RESIDUAL VECTOR IS LESS OR + XC EQUAL TO EPS2. THE EXACT SOLUTION HAS BEEN + XC OBTAINED. + XC IER = 300 SOLUTION NOT OBTAINED AFTER REACHING THE + XC VALUE NK+MK=MAXN. + XC IER = 400 SOLUTION NOT OBTAINED AFTER REACHING THE + XC DIMENSION N OF THE SYSTEM, DUE TO THE + XC PRECISION OF THE COMPUTER. + XC THE COMPUTATION CONTINUES UNTIL NK+MK=MAXN + XC AT MOST. + XC IER = 500 THE VALUE OF MAXBCK IS GREATER THAN MAXN OR + XC LESS THAN 2. + XC IER = 600 THE VALUE OF m EXCEEDS THE VALUE OF MAXBCK.+ XC k + XC IER = 700 JUMP FOR BREAKDOWN C(0) = 0 + XC (mk) + XC IER = 800 JUMP FOR BREAKDOWN sigma = 0 + XC k+1 (mk) + XC IER = 900 JUMP FOR NEAR-BREAKDOWN |sigma | <= EPS + XC k+1 + XC IER = 1000 JUMP FOR NEAR-BREAKDOWN |D(0)/C(0)|>= 1/EPS + XC IER = 1100 JUMP FOR NEAR-BREAKDOWN + XC NK=0, |C(1)/C(0)| >= 1/EPS + XC IER = 1200 JUMP FOR NEAR-BREAKDOWN + XC NK<>0, |D(0)/C(0)| <= EPS + XC + XC + XC EXTERNAL MODULES: + XC + XC - DOUBLE PRECISION FUNCTIONS SUNORM, FH, GAMETA, PINNER, RXSUMM, + XC SSUMM, ZSUMM + XC - SUBROUTINES CHKSIG, COMPO1, COMPO2, COMPOL, COMPP1, + XC CRESOL, GSOLVD, MATVEC, POLMUL, SOLSHF, + XC STOMHF, STORHF + XC + XC REMARKS: + XC + XC - ALL THE REAL ARGUMENTS PASSED FROM THE CALLING PROGRAM MUST BE + XC IN DOUBLE PRECISION. + XC - THE VARIABLE NK, THE VECTORS Y, X, R, AB, D, C, P, P1, PWK, LMAT + XC AND THE ARRAYS V, W, U, AND RMAT MUST NOT BE MODIFIED BY THE USER + XC BETWEEN TWO CONSECUTIVE CALLS OF THE SUBROUTINE. + XC - IN THE CALLING PROCEDURE, THE ROWS DIMENSION FOR THE ARRAYS A, V, W,+ XC AND U MUST BE THE SAME. + XC - IN THE CALLING PROCEDURE THE MINIMAL DIMENSIONS FOR THE MATRIX AND + XC VECTOR ARGUMENTS (IR=N) ARE: + XC A(N,N) + XC AB(N) + XC Y(N) + XC X(N) + XC R(N) + XC V(N,0:2*MAXBCK) + XC W(N,0:2*MAXBCK-2) + XC U(N,0:2*MAXBCK-1) + XC LMAT(8*MAXBCK) IF MAXBCK <= 2 + XC OR + XC LMAT((2*MAXBCK-1)*(2*MAXBCK-1)) IF MAXBCK > 2 + XC RMAT(2,0:2*MAXBCK-1) + XC ETA(0:MAXBCK) + XC ETAP(0:MAXBCK-1) + XC GAMMA(0:MAXBCK) + XC GAMMAP(0:MAXBCK-1) + XC D(0:2*MAXBCK-1) + XC C(0:2*MAXBCK-1) + XC P(0:MAXN) + XC P1(0:MAXN) + XC PWK(0:2*MAXN+1) + XC + XC AUTHORS: + XC + XC C. BREZINSKI + XC UNIVERSITY OF LILLE - FRANCE + XC M. REDIVO ZAGLIA + XC UNIVERSITY OF PADOVA - ITALY + XC + XC REFERENCE: + XC + XC - TREATMENT OF NEAR-BREAKDOWN IN THE CGS ALGORITHM + XC NUMERICAL ALGORITHMS 7 (1994) 33-73 + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC X SUBROUTINE BSMRZS (N, A, AB, Y, X, R, NORM, MAXBCK, MAXN, X * V, W, U, LMAT, RMAT, ETA, ETAP, GAMMA, GAMMAP, D, C, X * P, P1, PWK, IR, NK, EPS, EPS1, EPS2, INIT, IFLAG, IER) XC XC ... SPECIFICATIONS FOR ARGUMENTS XC X INTEGER IER, IFLAG, INIT, IR, MAXBCK, MAXN, N, NK X DOUBLE PRECISION EPS, EPS1, EPS2, NORM X DOUBLE PRECISION A(IR,1), W(IR,0:1), V(IR,0:1), U(IR,0:1), X * D(0:1), C(0:1), Y(1), AB(1), X(1), R(1), LMAT(1), X * RMAT(2,0:1), ETA(0:1), ETAP(0:1), GAMMA(0:1), X * GAMMAP(0:1), P(0:1), P1(0:1), PWK(0:1) XC XC ... SPECIFICATIONS FOR EXTERNAL MODULES XC X DOUBLE PRECISION SUNORM, PINNER, RXSUMM, SSUMM, ZSUMM XC XC ... SPECIFICATIONS FOR LOCAL VARIABLES XC X INTEGER I, INDS, J, MK, NFIND X DOUBLE PRECISION SIGMA X LOGICAL JUMP X SAVE NFIND, SIGMA XC XC ... EXECUTABLE STATEMENTS XC XC XC ... FIRST CALL OF THE SUBROUTINE XC X IF (INIT .EQ. 0) THEN X IER = 0 X INIT = 1 X NK = 0 X NFIND = 0 XC XC ... STORE THE VALUES FOR P1(0) AND P(0) XC X P1(0) = 1.0D0 X P(0) = 1.0D0 XC XC ... CONTROL THE VALUE OF MAXN AND OF MAXBCK XC X IF (MAXBCK .GT. MAXN .OR. MAXBCK .LT. 2) THEN X IER = 500 X RETURN X END IF XC XC ... COMPUTATION OF THE VECTORS r AND z XC 0 0 XC X CALL MATVEC (N, A, X, 1, R, 1, IR, 0) X DO 10 I = 1, N X R(I) = AB(I) - R(I) XC XC ... STORE THE VALUES FOR V(0), W(0) AND U(0) XC X V(I,0) = R(I) X W(I,0) = R(I) X U(I,0) = R(I) XC XC ... TEST FOR y = z XC 0 X IF (IFLAG .EQ. 0) Y(I) = R(I) X 10 CONTINUE XC XC ... COMPUTE THE NORM OF THE RESIDUAL VECTOR r XC 0 XC X NORM = SUNORM(N,R) XC XC (0) XC ... INITIALIZE THE VALUE OF sigma XC 0 XC X SIGMA = PINNER(N, Y, R, 1, IR) XC XC ... TEST FOR THE SOLUTION XC X IF (NORM .LE. EPS2) THEN X IER = 200 X NK = NK + 1 X RETURN X END IF X END IF XC XC ... FIRST AND NEXT CALLS OF THE SUBROUTINE XC X IF (IER .NE. 0) THEN X IER = 100 X RETURN X END IF XC XC ... INITIALIZATION OF m XC k X MK = 1 XC XC ... COMPUTATION OF V(1), V(2) AND U(1) XC X CALL MATVEC (N, A, V, 1, V, 2, IR, 0) X CALL MATVEC (N, A, V, 2, V, 3, IR, 0) X CALL MATVEC (N, A, U, 1, U, 2, IR, 0) XC XC (m ) XC k-1 XC ... SET D(0)= sigma XC k XC X D(0) = SIGMA XC XC ... COMPUTATION OF C(0)=(y,A*z ) XC k X C(0) = PINNER(N, Y, V, 2, IR) XC XC ... COMPUTATION OF D(1) AND OF C(1) XC X D(1) = PINNER(N, Y, U, 2, IR) X C(1) = PINNER(N, Y, V, 3, IR) XC XC ... INITIALIZATION OF THE FLAG JUMP XC X JUMP = .FALSE. XC XC ... CASE MK = 1 XC XC ... CONTROL FOR BREAKDOWN XC (IMPOSSIBLE TO COMPUTE BETA(0)) XC X IF (C(0) .EQ. 0.0D0) THEN X IER = 700 X JUMP = .TRUE. X END IF XC XC ... IF NO JUMP ... X IF (.NOT. JUMP) THEN XC XC ... COMPUTE BETA(0), ALPHA(0) XC AND (IF NK<>0) ALPHA'(0) XC X CALL SOLSHF (LMAT, RMAT, C, D, P1, NK, MK, EPS1, INDS) XC XC ... COMPUTE GAMMA(0), ETA(0) XC AND (IF NK<>0) ETA'(0) XC X CALL CRESOL (RMAT, GAMMA, GAMMAP, ETA, ETAP, P1, NK, MK) XC XC ... COMPUTE THE VECTOR s XC k+1 XC X IF (NK .EQ. 0) THEN X DO 20 J = 1, N X AB(J) = ETA(0)*U(J,0) + U(J,1) - ETA(0)*GAMMA(0)*V(J,1) X * - GAMMA(0)*V(J,2) X 20 CONTINUE X ELSE X DO 30 J = 1, N X AB(J) = ETA(0)*U(J,0) + U(J,1) - ETA(0)*GAMMA(0)*V(J,1) X * - GAMMA(0)*V(J,2) + ETAP(0)*W(J,0) + U(J,1) X 30 CONTINUE X END IF XC XC (m ) XC k XC ... COMPUTATION OF sigma = (y, s ) XC k+1 k+1 XC X SIGMA = PINNER(N, Y, AB, 1, IR) XC XC ... COMPUTE THE RESIDUAL VECTOR r XC k+1 XC X DO 40 J = 1, N X R(J) = W(J,0)-2.0D0*GAMMA(0)*U(J,1)+GAMMA(0)**2*V(J,2) X 40 CONTINUE XC XC ... COMPUTE THE NORM OF THE RESIDUAL VECTOR r XC k+1 XC X NORM = SUNORM(N,R) XC XC ... CONTROL FOR BREAKDOWN OR NEAR-BREAKDOWN XC OR POSSIBLE SOLUTION (CASE MK = 1) XC X CALL CHKSIG (C, D, NK, MK, EPS, SIGMA, JUMP, IER) XC XC ... TEST FOR THE JUMP XC X IF ((NORM .LE. EPS2) .OR. (.NOT. JUMP)) THEN XC XC ... COMPUTE THE SOLUTION VECTOR x XC k+1 XC X DO 50 J = 1, N X X(J) = X(J)+2.0D0*GAMMA(0)*U(J,0)-GAMMA(0)**2*V(J,1) X 50 CONTINUE XC XC ... TEST FOR THE SOLUTION XC X IF (NORM .LE. EPS2) THEN X IER = 200 XC XC ... COMPUTE THE VALUE OF n XC k+1 X NK = NK + 1 X RETURN X END IF X END IF X END IF XC XC ... START OF THE LOOP TO COMPUTE BETA , BETA' AND XC i i XC ALPHA , ALPHA' UNTIL THE SYSTEM IS NON SINGULAR XC i i XC OR NO BREAKDOWN AND NO NEAR-BREAKDOWN AT THE XC NEXT STEP XC X 60 CONTINUE X IF (JUMP .OR. (INDS .NE. 0)) THEN XC XC ... SET THE FLAG JUMP XC X JUMP = .FALSE. XC XC ... INCREASE THE VALUE OF m XC k X MK = MK + 1 XC XC ... CONTROL THE VALUE OF MAXN XC X IF (MK .EQ. MAXN-NK+1) THEN X NK = MAXN X IER = 300 X RETURN X END IF XC XC ... CONTROL THE VALUE OF m XC k X IF (MK .GT. MAXBCK) THEN X NK = NK + MK - 1 X IER = 600 X RETURN X END IF XC XC ... FOR THE NEW m XC k XC X DO 70 I = -1, 0 XC XC ... COMPUTATION OF V(2*MK-1) AND OF V(2*MK) XC X CALL MATVEC (N, A, V, 2*MK+I, V, 2*MK+I+1, IR, 0) XC XC ... COMPUTATION OF W(2*MK-3), W(2*MK-2) XC U(2*MK-2), U(2*MK-1) XC ONLY IF m <= n + 1 XC k k XC ELSE COMPUTATION OF U(NK+MK) XC X IF (MK .LE. NK+1) THEN X CALL MATVEC (N, A, W, 2*MK+I-2, W, 2*MK+I-1, IR, 0) X CALL MATVEC (N, A, U, 2*MK+I-1, U, 2*MK+I, IR, 0) X ELSE IF (I .EQ. -1) THEN X CALL MATVEC (N, A, U, NK+MK, U, NK+MK+1, IR, 0) X END IF X 70 CONTINUE XC X DO 80 I = 2, 1, -1 XC XC ... COMPUTATION OF C(2m -2) AND C(2m -1), XC k k XC ... COMPUTATION OF D(2m -2) AND D(2m -1) XC k k XC IF m <= n XC k k XC ... COMPUTATION OF D(n + m - 1) XC k k XC IF m > n XC k k XC X C(2*MK-I) = PINNER(N, Y, V, 2*MK-I+2, IR) X IF (MK .LE. NK) THEN X D(2*MK-I) = PINNER(N, Y, U, 2*MK-I+1, IR) X ELSE IF (I .EQ. 2) THEN X D(MK+NK-1) = PINNER(N, Y, U, MK+NK, IR) X END IF XC X 80 CONTINUE XC XC ... SOLVE THE SYSTEM XC XC ... COMPUTATION OF BETA , BETA' AND OF XC i i XC ALPHA , ALPHA' XC i i XC X CALL SOLSHF (LMAT, RMAT, C, D, P1, NK, MK, EPS1, INDS) XC XC ... IF NON SINGULAR SYSTEM THEN ... XC X IF (INDS .EQ. 0) THEN XC XC ... CONSTRUCT THE VECTORS OF THE SOLUTIONS XC GAMMA , GAMMA' , ETA AND ETA' XC i i i i XC X CALL CRESOL (RMAT, GAMMA, GAMMAP, ETA, ETAP, P1, NK, MK) XC XC ... COMPUTATION OF THE VECTOR s XC k+1 XC X CALL COMPO2 (GAMMA, GAMMAP, ETA, ETAP, NK, MK, X * LMAT(1), LMAT(2*MAXBCK), LMAT(4*MAXBCK), X * LMAT(6*MAXBCK+1)) X DO 90 J = 1, N X AB(J) = SSUMM(LMAT(1), LMAT(2*MAXBCK), LMAT(4*MAXBCK), X * LMAT(6*MAXBCK+1), V, W, U, IR, J, NK, MK) X 90 CONTINUE XC XC (m ) XC k XC ... COMPUTATION OF sigma = (y, s ) XC k+1 k+1 XC X SIGMA = PINNER(N, Y, AB, 1, IR) XC XC XC ... COMPUTE THE RESIDUAL VECTOR r XC k+1 XC X CALL COMPOL (GAMMA, GAMMAP, NK, MK, LMAT(1), LMAT(2*MAXBCK), X * LMAT(4*MAXBCK), LMAT(6*MAXBCK+1)) X DO 100 J = 1, N X R(J) = RXSUMM(LMAT(2*MAXBCK), LMAT(4*MAXBCK), X * LMAT(6*MAXBCK+1), V, W, U, IR, J, NK, MK, 1) X 100 CONTINUE XC XC ... COMPUTE THE NORM OF THE RESIDUAL VECTOR r XC k+1 XC X NORM = SUNORM(N,R) XC XC ... CONTROL FOR BREAKDOWN OR NEAR-BREAKDOWN XC OR POSSIBLE SOLUTION (CASE MK > 1) XC X CALL CHKSIG (C, D, NK, MK, EPS, SIGMA, JUMP, IER) XC XC ... TEST FOR THE JUMP XC X IF ((NORM .LE. EPS2) .OR. (.NOT. JUMP)) THEN XC XC ... COMPUTE THE SOLUTION VECTOR x XC k+1 XC X DO 110 J = 1, N X X(J) = X(J) - RXSUMM(LMAT(1), LMAT(4*MAXBCK), X * LMAT(6*MAXBCK+1), V, W, U, IR, J, NK, MK, 0) X 110 CONTINUE XC XC ... TEST FOR THE SOLUTION XC X IF (NORM .LE. EPS2) THEN X IER = 200 XC XC ... COMPUTE THE VALUE OF n XC k+1 X NK = NK + MK X RETURN X END IF X END IF X END IF XC XC ... END OF THE LOOP TO COMPUTE BETA , BETA' AND XC i i XC ALPHA , ALPHA' XC i i X GO TO 60 X END IF XC XC ... COMPUTE THE VALUE OF n XC k+1 X NK = NK + MK XC XC ... CONTROL THE VALUE OF n XC k XC X IF (NK .GE. N .AND. NFIND .EQ. 0) THEN X IER = 400 X NFIND = 1 X END IF XC XC ... CONTROL THE VALUE OF MAXN XC X IF (NK .EQ. MAXN) THEN X IER = 300 X RETURN X END IF XC XC ... COMPUTE THE VECTOR z XC k+1 XC X CALL COMPO1 (ETA, ETAP, NK-MK, MK, LMAT(1), LMAT(2*MAXBCK), X * LMAT(4*MAXBCK)) X DO 120 J = 1, N X V(J,0) = ZSUMM(LMAT(1), LMAT(2*MAXBCK), LMAT(4*MAXBCK), X * V, W, U, IR, J, NK-MK, MK) X 120 CONTINUE XC XC ... SAVE THE VALUES OF s AND OF r XC k+1 k+1 X DO 130 J = 1, N X U(J,0) = AB(J) X W(J,0) = R(J) X 130 CONTINUE XC XC ... COMPUTE THE COEFFICIENTS OF XC (1) XC THE POLYNOMIALS P AND P XC k+1 k+1 XC X CALL COMPP1 (ETA, ETAP, GAMMA, GAMMAP, P, P1, NK-MK, MK, PWK) XC X RETURN X END SHAR_EOF if test 33381 -ne "`wc -c bsmrzs.f`" then echo shar: error transmitting bsmrzs.f '(should have been 33381 characters)' fi echo shar: extracting chksig.f sed 's/^X//' << \SHAR_EOF > chksig.f XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC SUBROUTINE NAME - CHKSIG + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC DESCRIPTION: + XC + XC CONTROLS THE CONDITIONS THAT COULD PRODUCE A BREAKDOWN, + XC A NEAR-BREAKDOWN OR IF THE SOLUTION IS ABOUT TO BE OBTAINED. + XC + XC USAGE: + XC + XC CALL CHKSIG (C, D, NK, MK, EPS, SIGMA, JUMP, IER) + XC + XC ARGUMENTS: + XC + XC C - INPUT REAL VECTOR CONTAINING THE c 'S. + XC i + XC + XC D - INPUT REAL VECTOR CONTAINING THE d 'S. + XC i + XC + XC NK - INPUT INTEGER EQUAL TO THE DIMENSION OF THE PREVIOUS + XC MATRIX OF THE SYSTEM. + XC + XC MK - INPUT INTEGER EQUAL TO THE LENGTH OF THE NEW JUMP. + XC + XC EPS - INPUT REAL VALUE USED FOR TESTING THE NEAR BREAKDOWN. + XC THE SUBROUTINE DOES NOT CONTROL WHETHER OR NOT EPS IS + XC A NEGATIVE OR ZERO REAL NUMBER. + XC + XC (m ) + XC k + XC SIGMA - INPUT REAL VALUE EQUAL TO THE sigma . + XC k+1 + XC + XC JUMP - OUTPUT LOGICAL VALUE. + XC IF JUMP = .FALSE. THEN THERE IS NO BREAKDOWN, NO + XC NEAR-BREAKDOWN NOR POSSIBLE SOLUTION. + XC IF JUMP = .TRUE. THEN THERE IS A BREAKDOWN OR A + XC NEAR-BREAKDOWN OR A POSSIBLE SOLUTION. + XC + XC IER - OUTPUT INDEX WARNING/ERROR + XC (mk) + XC IER = 800 BREAKDOWN OR POSSIBLE SOLUTION. sigma = 0 + XC k+1 + XC IER = 900 NEAR-BREAKDOWN OR POSSIBLE SOLUTION. + XC (mk) + XC |sigma | <= EPS + XC k+1 + XC IER = 1000 JUMP FOR NEAR-BREAKDOWN. |D(0)/C(0)| >= 1/EPS. + XC IER = 1100 JUMP FOR NEAR-BREAKDOWN. + XC NK=0, |C(1)/C(0)| >= 1/EPS. + XC IER = 1200 JUMP FOR NEAR-BREAKDOWN. + XC NK<>0, |D(0)/C(0)| <= EPS. + XC + XC REMARKS: + XC + XC - ALL THE REAL ARGUMENTS PASSED FROM THE CALLING PROGRAM MUST BE + XC IN DOUBLE PRECISION. + XC - IN THE CALLING PROCEDURE THE DIMENSIONS FOR THE ARGUMENTS ARE: + XC C(0:*) + XC D(0:*) + XC + XC AUTHORS: + XC + XC C. BREZINSKI + XC UNIVERSITY OF LILLE - FRANCE + XC M. REDIVO ZAGLIA + XC UNIVERSITY OF PADOVA - ITALY + XC + XC REFERENCES: + XC + XC - TREATMENT OF NEAR-BREAKDOWN IN THE CGS ALGORITHMS + XC NUMERICAL ALGORITHMS 7 (1994) 33-73 + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC X SUBROUTINE CHKSIG (C, D, NK, MK, EPS, SIGMA, JUMP, IER) XC XC ... SPECIFICATIONS FOR ARGUMENTS XC X INTEGER IER, NK X LOGICAL JUMP X DOUBLE PRECISION EPS, SIGMA X DOUBLE PRECISION C(0:1), D(0:1) XC XC ... EXECUTABLE STATEMENTS XC XC ... CONTROL FOR BREAKDOWN OR NEAR BREAKDOWN XC OR IF THE SOLUTION IS ABOUT TO BE OBTAINED XC X IF (DABS(SIGMA) .LE. EPS) THEN X JUMP = .TRUE. X IF (SIGMA .EQ. 0.0D0) THEN X IER = 800 X ELSE X IER = 900 X END IF X RETURN X END IF XC XC ... ONLY FOR MK = 1 XC X IF (MK .EQ. 1) THEN XC XC ... CONTROL FOR NEAR BREAKDOWN ... XC XC ... IN COMPUTING BETA(0) XC X IF (DABS(D(0)/C(0)) .GE. 1.0D0/EPS) THEN X JUMP = .TRUE. X IER = 1000 XC XC OR ... XC ... IN COMPUTING ALPHA(0) XC FOR THE CASE MK > NK = 0 XC X ELSE IF (NK.EQ.0 .AND. DABS(C(1)/C(0)) .GE. 1.0D0/EPS) THEN X JUMP = .TRUE. X IER = 1100 XC XC OR ... XC ... IN COMPUTING ALPHA'(0) XC FOR THE CASE MK <= NK, NK<>0 XC X ELSE IF (NK.NE.0 .AND. DABS(D(0)/C(0)) .LE. EPS) THEN X JUMP = .TRUE. X IER = 1200 X END IF X END IF XC X RETURN X END SHAR_EOF if test 7834 -ne "`wc -c chksig.f`" then echo shar: error transmitting chksig.f '(should have been 7834 characters)' fi echo shar: extracting compo1.f sed 's/^X//' << \SHAR_EOF > compo1.f XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC SUBROUTINE NAME - COMPO1 + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC DESCRIPTION: + XC + XC COMPUTES THE COEFFICIENTS OF THE POLYNOMIALS APPEARING IN THE + XC FORMULA FOR THE COMPUTATION OF z . + XC k+1 + XC + XC USAGE: + XC CALL COMPO1 (ETA, ETAP, NK, MK, POL1, POL2, POL3) + XC + XC ARGUMENTS: + XC + XC ETA - INPUT REAL VECTOR CONTAINING THE ETA . + XC i + XC + XC ETAP - INPUT REAL VECTOR CONTAINING THE ETA'. + XC i + XC + XC NK - INPUT INTEGER EQUAL TO THE DIMENSION OF THE PREVIOUS + XC MATRIX OF THE SYSTEM. + XC + XC MK - INPUT INTEGER EQUAL TO THE DIMENSION OF THE NEW JUMP. + XC + XC POL1 - OUTPUT REAL VECTOR CONTAINING THE COEFFICIENTS OF THE + XC POLYNOMIAL t(xi)**2. + XC + XC POL2 - OUTPUT REAL VECTOR CONTAINING THE COEFFICIENTS OF THE + XC POLYNOMIAL t(xi) * q(xi). + XC + XC POL3 - OUTPUT REAL VECTOR CONTAINING THE COEFFICIENTS OF THE + XC POLYNOMIAL q(xi)**2. + XC + XC EXTERNAL MODULES: + XC + XC - SUBROUTINE POLMUL + XC + XC REMARKS: + XC + XC - ALL THE REAL ARGUMENTS PASSED FROM THE CALLING PROGRAM MUST BE + XC IN DOUBLE PRECISION. + XC - IN THE CALLING PROCEDURE THE DIMENSIONS FOR THE MATRIX AND VECTOR + XC ARGUMENTS ARE: + XC ETA(0:*) + XC ETAP(0:*) + XC POL1(0:*) + XC POL2(0:*) + XC POL3(0:*) + XC + XC AUTHORS: + XC + XC C. BREZINSKI + XC UNIVERSITY OF LILLE - FRANCE + XC M. REDIVO ZAGLIA + XC UNIVERSITY OF PADOVA - ITALY + XC + XC REFERENCE: + XC + XC - TREATMENT OF NEAR-BREAKDOWN IN THE CGS ALGORITHM + XC NUMERICAL ALGORITHMS 7 (1994) 33-73 + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC X SUBROUTINE COMPO1 (ETA, ETAP, NK, MK, POL1, POL2, POL3) XC XC ... SPECIFICATIONS FOR ARGUMENTS XC X INTEGER NK, MK X DOUBLE PRECISION ETA(0:1), ETAP(0:1), POL1(0:1), POL2(0:1), X * POL3(0:2) XC XC ... SPECIFICATIONS FOR LOCAL VARIABLES XC X INTEGER IER, NV XC XC ... EXECUTABLE STATEMENTS XC XC XC ... CASE MK = 1 XC X IF (MK .EQ. 1) THEN X POL3(0) = ETA(0)**2 X POL3(1) = 2.0D0 * ETA(0) X POL3(2) = 1.0D0 X IF (MK .LE. NK) THEN X POL1(0) = ETAP(0)**2 X POL2(0) = ETA(0) * ETAP(0) X POL2(1) = ETAP(0) X END IF X RETURN X END IF XC XC ... SET THE DEGREE OF THE POLYNOMIAL t(xi) XC (COEFFICIENTS STORED IN ETAP) XC XC ... CASE MK <= NK XC X IF (MK .LE. NK) THEN X NV = MK - 1 X ELSE XC XC ... CASE MK > NK XC X NV = NK - 1 X END IF XC X IF (NK .NE. 0) THEN XC XC ... COMPUTES THE COEFFICIENTS OF THE XC POLYNOMIAL t(xi)**2 XC X CALL POLMUL(ETAP, NV, ETAP, NV, POL1, IER) XC XC ... COMPUTES THE COEFFICIENTS OF THE XC POLYNOMIAL t(xi) * q(xi) XC X CALL POLMUL(ETAP, NV, ETA, MK, POL2, IER) X END IF XC XC ... COMPUTES THE COEFFICIENTS OF THE XC POLYNOMIAL q(xi)**2 XC X CALL POLMUL(ETA, MK, ETA, MK, POL3, IER) XC X RETURN X END SHAR_EOF if test 6954 -ne "`wc -c compo1.f`" then echo shar: error transmitting compo1.f '(should have been 6954 characters)' fi echo shar: extracting compo2.f sed 's/^X//' << \SHAR_EOF > compo2.f XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC SUBROUTINE NAME - COMPO2 + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC DESCRIPTION: + XC + XC COMPUTES THE COEFFICIENTS OF THE POLYNOMIALS APPEARING IN THE + XC FORMULA FOR THE COMPUTATION OF s . + XC k+1 + XC + XC USAGE: + XC CALL COMPO2 (GAMMA, GAMMAP, ETA, ETAP, NK, MK, POL1, POL2, + XC POL3, POL4) + XC + XC ARGUMENTS: + XC + XC GAMMA - INPUT REAL VECTOR CONTAINING THE GAMMA . + XC i + XC + XC GAMMAP - INPUT REAL VECTOR CONTAINING THE GAMMA' . + XC i + XC + XC ETA - INPUT REAL VECTOR CONTAINING THE ETA . + XC i + XC + XC ETAP - INPUT REAL VECTOR CONTAINING THE ETA'. + XC i + XC + XC NK - INPUT INTEGER EQUAL TO THE DIMENSION OF THE PREVIOUS + XC MATRIX OF THE SYSTEM. + XC + XC MK - INPUT INTEGER EQUAL TO THE DIMENSION OF THE NEW JUMP. + XC + XC POL1 - OUTPUT REAL VECTOR CONTAINING THE COEFFICIENTS OF THE + XC POLYNOMIAL t(xi) * (1 - xi*v(xi)). + XC + XC POL2 - OUTPUT REAL VECTOR CONTAINING THE COEFFICIENTS OF THE + XC POLYNOMIAL (1 - xi*v(xi)) * q(xi). + XC + XC POL3 - OUTPUT REAL VECTOR CONTAINING THE COEFFICIENTS OF THE + XC POLYNOMIAL t(xi) * w(xi). + XC + XC POL4 - OUTPUT REAL VECTOR CONTAINING THE COEFFICIENTS OF THE + XC POLYNOMIAL w(xi) * q(xi). + XC + XC EXTERNAL MODULES: + XC + XC - SUBROUTINE POLMUL + XC + XC REMARKS: + XC + XC - ALL THE REAL ARGUMENTS PASSED FROM THE CALLING PROGRAM MUST BE + XC IN DOUBLE PRECISION. + XC - IN THE CALLING PROCEDURE THE DIMENSIONS FOR THE MATRIX AND VECTOR + XC ARGUMENTS ARE: + XC GAMMA(0:*) + XC GAMMAP(0:*) + XC ETA(0:*) + XC ETAP(0:*) + XC POL1(0:*) + XC POL2(0:*) + XC POL3(0:*) + XC POL4(0:*) + XC + XC AUTHORS: + XC + XC C. BREZINSKI + XC UNIVERSITY OF LILLE - FRANCE + XC M. REDIVO ZAGLIA + XC UNIVERSITY OF PADOVA - ITALY + XC + XC REFERENCE: + XC + XC - TREATMENT OF NEAR-BREAKDOWN IN THE CGS ALGORITHM + XC NUMERICAL ALGORITHMS 7 (1994) 33-73 + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC X SUBROUTINE COMPO2 (GAMMA, GAMMAP, ETA, ETAP, NK, MK, POL1, X * POL2, POL3, POL4) XC XC ... SPECIFICATIONS FOR ARGUMENTS XC X INTEGER NK, MK X DOUBLE PRECISION GAMMA(0:1), GAMMAP(0:1), ETA(0:1), ETAP(0:1), X * POL1(0:1), POL2(0:1), POL3(0:1), POL4(0:1) XC XC ... SPECIFICATIONS FOR LOCAL VARIABLES XC X INTEGER I, IER, NV, NV1 XC XC ... EXECUTABLE STATEMENTS XC XC XC ... CASE MK = 1 XC X IF (MK .EQ. 1) THEN X POL2(0) = ETA(0) X POL4(0) = GAMMA(0) * ETA(0) X POL4(1) = GAMMA(0) X IF (MK .LE. NK) THEN X POL1(0) = ETAP(0) X POL3(0) = ETAP(0) * GAMMA(0) X END IF X RETURN X END IF XC XC ... SET THE DEGREE OF THE POLYNOMIAL v(xi) XC (COEFFICIENTS STORED IN GAMMAP) XC AND OF THE POLYNOMIAL t(xi) XC (COEFFICIENTS STORED IN ETAP) XC XC ... CASE MK <= NK XC X IF (MK .LE. NK) THEN X NV = MK - 2 X NV1 = MK - 1 X ELSE XC XC ... CASE MK > NK XC X NV = NK - 1 X NV1 = NK - 1 X END IF XC X POL4(0) = 1.0D0 X IF (NK .NE. 0) THEN XC XC ... STORES THE COEFFICIENTS OF THE XC POLYNOMIAL (1 - xi*v(xi)) XC X DO 10 I = 1, NV+1 X POL4(I) = -GAMMAP(I-1) X 10 CONTINUE XC XC ... COMPUTES THE COEFFICIENTS OF THE XC POLYNOMIAL t(xi) * (1 - xi*v(xi)) XC X CALL POLMUL(ETAP, NV1, POL4, NV+1, POL1, IER) XC XC ... COMPUTES THE COEFFICIENTS OF THE XC POLYNOMIAL t(xi) * w(xi) XC X CALL POLMUL(ETAP, NV1, GAMMA, MK-1, POL3, IER) X END IF XC XC ... COMPUTES THE COEFFICIENTS OF THE XC POLYNOMIAL (1 - xi*v(xi)) * q(xi) XC X CALL POLMUL(POL4, NV+1, ETA, MK, POL2, IER) XC XC ... COMPUTES THE COEFFICIENTS OF THE XC POLYNOMIAL w(xi) * q(xi) XC X CALL POLMUL(GAMMA, MK-1, ETA, MK, POL4, IER) XC X RETURN X END SHAR_EOF if test 8655 -ne "`wc -c compo2.f`" then echo shar: error transmitting compo2.f '(should have been 8655 characters)' fi echo shar: extracting compol.f sed 's/^X//' << \SHAR_EOF > compol.f XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC SUBROUTINE NAME - COMPOL + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC DESCRIPTION: + XC + XC COMPUTES THE COEFFICIENTS OF THE POLYNOMIALS APPEARING IN THE + XC FORMULAE FOR THE COMPUTATION OF x (THE NEW SOLUTION) AND OF r + XC k+1 k+1 + XC (THE NEW RESIDUAL). + XC + XC USAGE: + XC CALL COMPOL (GAMMA, GAMMAP, NK, MK, POL1, POL2, POL3, POL4) + XC + XC ARGUMENTS: + XC + XC GAMMA - INPUT REAL VECTOR CONTAINING THE GAMMA . + XC i + XC + XC GAMMAP - INPUT REAL VECTOR CONTAINING THE GAMMA' . + XC i + XC + XC NK - INPUT INTEGER EQUAL TO THE DIMENSION OF THE PREVIOUS + XC MATRIX OF THE SYSTEM. + XC + XC MK - INPUT INTEGER EQUAL TO THE DIMENSION OF THE NEW JUMP. + XC + XC POL1 - OUTPUT REAL VECTOR CONTAINING THE COEFFICIENTS OF THE + XC POLYNOMIAL v(xi) * (xi*v(xi) - 2). + XC + XC POL2 - OUTPUT REAL VECTOR CONTAINING THE COEFFICIENTS OF THE + XC POLYNOMIAL (1 - xi*v(xi))**2. + XC + XC POL3 - OUTPUT REAL VECTOR CONTAINING THE COEFFICIENTS OF THE + XC POLYNOMIAL (1 - xi*v(xi)) * w(xi). + XC + XC POL4 - OUTPUT REAL VECTOR CONTAINING THE COEFFICIENTS OF THE + XC POLYNOMIAL w(xi)**2. + XC + XC EXTERNAL MODULES: + XC + XC - SUBROUTINE POLMUL + XC + XC REMARKS: + XC + XC - ALL THE REAL ARGUMENTS PASSED FROM THE CALLING PROGRAM MUST BE + XC IN DOUBLE PRECISION. + XC - IN THE CALLING PROCEDURE THE DIMENSIONS FOR THE MATRIX AND VECTOR + XC ARGUMENTS ARE: + XC GAMMA(0:*) + XC GAMMAP(0:*) + XC POL1(0:*) + XC POL2(0:*) + XC POL3(0:*) + XC POL4(0:*) + XC + XC AUTHORS: + XC + XC C. BREZINSKI + XC UNIVERSITY OF LILLE - FRANCE + XC M. REDIVO ZAGLIA + XC UNIVERSITY OF PADOVA - ITALY + XC + XC REFERENCE: + XC + XC - TREATMENT OF NEAR-BREAKDOWN IN THE CGS ALGORITHM + XC NUMERICAL ALGORITHMS 7 (1994) 33-73 + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC X SUBROUTINE COMPOL (GAMMA, GAMMAP, NK, MK, POL1, POL2, POL3, POL4) XC XC ... SPECIFICATIONS FOR ARGUMENTS XC X INTEGER NK, MK X DOUBLE PRECISION GAMMA(0:1), GAMMAP(0:1), POL1(0:1), POL2(0:1), X * POL3(0:1), POL4(0:1) XC XC ... SPECIFICATIONS FOR LOCAL VARIABLES XC X INTEGER I, IER, NV XC XC ... EXECUTABLE STATEMENTS XC XC XC ... CASE MK = 1 XC X IF (MK .EQ. 1) THEN X POL3(0) = GAMMA(0) X POL4(0) = GAMMA(0)**2 X RETURN X END IF XC XC ... SET THE DEGREE OF THE POLYNOMIAL v(xi) XC (COEFFICIENTS STORED IN GAMMAP) XC XC ... CASE MK <= NK XC X IF (MK .LE. NK) THEN X NV = MK - 2 X ELSE XC XC ... CASE MK > NK XC X NV = NK - 1 X END IF X IF (NK .NE. 0) THEN XC XC ... STORES THE COEFFICIENTS OF THE XC POLYNOMIAL (xi*v(xi) - 2) XC X POL4(0) = -2.0D0 X DO 10 I = 1, NV+1 X POL4(I) = GAMMAP(I-1) X 10 CONTINUE XC XC ... COMPUTES THE COEFFICIENTS OF THE XC POLYNOMIAL v(xi) * (xi*v(xi) - 2) XC X CALL POLMUL(GAMMAP, NV, POL4, NV+1, POL1, IER) XC XC ... STORES THE COEFFICIENTS OF THE XC POLYNOMIAL (1 - xi*v(xi)) XC X POL4(0) = 1.0D0 X DO 20 I = 1, NV+1 X POL4(I) = -GAMMAP(I-1) X 20 CONTINUE XC XC ... COMPUTES THE COEFFICIENTS OF THE XC POLYNOMIAL (1 - xi*v(xi))**2 XC X CALL POLMUL(POL4, NV+1, POL4, NV+1, POL2, IER) XC X ELSE X POL4(0) = 1.0D0 X POL2(0) = 1.0D0 X END IF XC XC ... COMPUTES THE COEFFICIENTS OF THE XC POLYNOMIAL (1 - xi*v(xi)) * w(xi) XC X CALL POLMUL(POL4, NV+1, GAMMA, MK-1, POL3, IER) XC XC ... COMPUTES THE COEFFICIENTS OF THE XC POLYNOMIAL w(xi)**2 XC X CALL POLMUL(GAMMA, MK-1, GAMMA, MK-1, POL4, IER) XC X RETURN X END SHAR_EOF if test 7943 -ne "`wc -c compol.f`" then echo shar: error transmitting compol.f '(should have been 7943 characters)' fi echo shar: extracting compp1.f sed 's/^X//' << \SHAR_EOF > compp1.f XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC SUBROUTINE NAME - COMPP1 + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC DESCRIPTION: + XC + XC (1) + XC COMPUTES THE COEFFICIENTS OF THE POLYNOMIALS P AND P . + XC k+1 k+1 + XC + XC USAGE: + XC CALL COMPP1 (ETA, ETAP, GAMMA, GAMMAP, P, P1, NK, MK, POL) + XC + XC ARGUMENTS: + XC + XC ETA - INPUT REAL VECTOR CONTAINING THE ETA . + XC i + XC + XC ETAP - INPUT REAL VECTOR CONTAINING THE ETA'. + XC i + XC + XC GAMMA - INPUT/WORKING REAL VECTOR CONTAINING, IN INPUT, THE + XC GAMMA . + XC i + XC + XC GAMMAP - INPUT/WORKING REAL VECTOR CONTAINING, IN INPUT, THE + XC GAMMA' . + XC i + XC + XC P - INPUT/OUTPUT REAL VECTOR CONTAINING, IN INPUT, THE + XC COEFFICIENTS OF THE POLYNOMIAL P OF DEGREE NK AND, IN + XC k + XC OUTPUT, THE COEFFICIENTS OF THE POLYNOMIAL P . + XC k+1 + XC + XC P1 - INPUT/OUTPUT REAL VECTOR CONTAINING, IN INPUT, THE + XC (1) + XC COEFFICIENTS OF THE POLYNOMIAL P OF DEGREE NK AND, + XC k (1) + XC IN OUTPUT, THE COEFFICIENTS OF THE POLYNOMIAL P . + XC k+1 + XC + XC NK - INPUT INTEGER EQUAL TO THE DIMENSION OF THE PREVIOUS + XC MATRIX OF THE SYSTEM. + XC + XC MK - INPUT INTEGER EQUAL TO THE DIMENSION OF THE NEW JUMP. + XC + XC POL - WORKING REAL VECTOR OF DIMENSION AT LEAST 2*(NK+MK+1). + XC + XC EXTERNAL MODULES: + XC + XC - SUBROUTINE POLMUL + XC + XC REMARKS: + XC + XC - ALL THE REAL ARGUMENTS PASSED FROM THE CALLING PROGRAM MUST BE + XC IN DOUBLE PRECISION. + XC - IN THE CALLING PROCEDURE THE DIMENSIONS FOR THE MATRIX AND VECTOR + XC ARGUMENTS ARE: + XC ETA(0:*) + XC ETAP(0:*) + XC GAMMA(0:*) + XC GAMMAP(0:*) + XC P(0:*) + XC P1(0:*) + XC POL(0:*) + XC + XC AUTHORS: + XC + XC C. BREZINSKI + XC UNIVERSITY OF LILLE - FRANCE + XC M. REDIVO ZAGLIA + XC UNIVERSITY OF PADOVA - ITALY + XC + XC REFERENCE: + XC + XC - TREATMENT OF NEAR-BREAKDOWN IN THE CGS ALGORITHM + XC NUMERICAL ALGORITHMS 7 (1994) 33-73 + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC X SUBROUTINE COMPP1 (ETA, ETAP, GAMMA, GAMMAP, P, P1, NK, MK, POL) XC XC ... SPECIFICATIONS FOR ARGUMENTS XC X INTEGER NK, MK X DOUBLE PRECISION ETA(0:1), ETAP(0:1), GAMMA(0:1), GAMMAP(0:1), X * P(0:1), P1(0:1), POL(0:1) XC XC ... SPECIFICATIONS FOR LOCAL VARIABLES XC X INTEGER I, IER, NM, NV, NV1 XC XC ... EXECUTABLE STATEMENTS XC XC XC ... CASE MK = 1 XC X IF (MK .EQ. 1) THEN X POL(0) = ETA(0)*P1(0) X POL(NK+1) = 1.0D0 X P(NK+1) = -GAMMA(0) X IF (MK .LE. NK) THEN X POL(0) = POL(0) + ETAP(0) X DO 10 I = 1, NK X POL(I) = ETA(0)*P1(I) + P1(I-1) + ETAP(0)*P(I) X P(I) = P(I) - GAMMA(0)*P1(I-1) X 10 CONTINUE X END IF X DO 20 I = 0, NK+1 X P1(I) = POL(I) X 20 CONTINUE X RETURN X END IF XC X NM = NK + MK + 1 XC XC ... SET THE DEGREE OF THE POLYNOMIAL t(xi) XC (COEFFICIENTS STORED IN ETAP) AND XC OF THE POLYNOMIAL v(xi) (COEFFICIENTS XC STORED IN GAMMAP) XC XC ... CASE MK <= NK XC X IF (MK .LE. NK) THEN X NV = MK - 1 X NV1 = MK - 2 X ELSE XC XC ... CASE MK > NK XC X NV = NK - 1 X NV1 = NK - 1 X END IF XC XC ... COMPUTES THE COEFFICIENTS OF THE XC POLYNOMIAL t(xi) * P (xi) (NK .NE. 0) XC k XC X IF (NK .NE. 0) CALL POLMUL(ETAP, NV, P, NK, POL, IER) XC XC ... COMPUTES THE COEFFICIENTS OF THE XC (1) XC POLYNOMIAL q(xi) * P (xi) AND XC k XC SAVE IN POL THIS TERM OF THE XC (1) XC POLYNOMIAL P (xi) XC k+1 XC X IF (MK .LE. NK) THEN X CALL POLMUL(ETA, MK, P1, NK, POL(NM), IER) X ELSE X CALL POLMUL(P1, NK, ETA, MK, POL(NM), IER) X END IF XC XC ... CASE NK .NE. 0 XC X IF (NK .NE. 0) THEN XC XC ... ADDS IN POL THE SECOND TERM TO COMPLETE XC THE COMPUTATION OF THE COEFFICIENTS OF THE XC (1) XC POLYNOMIAL P (xi) XC k+1 XC X DO 30 I = 0, NK + NV X POL(NM+I) = POL(NM+I) + POL(I) X 30 CONTINUE XC XC ... STORES THE COEFFICIENTS OF THE XC POLYNOMIAL (1 - xi*v(xi)) XC X DO 40 I = NV1 + 1, 1, -1 X GAMMAP(I) = -GAMMAP(I-1) X 40 CONTINUE X GAMMAP(0) = 1.0D0 XC XC ... COMPUTES THE COEFFICIENTS OF THE XC POLYNOMIAL (1 - xi*v(xi))* P (xi) XC k XC X CALL POLMUL(GAMMAP, NV1+1, P, NK, POL, IER) XC X ELSE X POL(0) = P(0) X END IF XC XC ... STORES THE COEFFICIENTS OF THE XC POLYNOMIAL -xi*w(xi) XC X DO 50 I = MK, 1, -1 X GAMMA(I) = -GAMMA(I-1) X 50 CONTINUE X GAMMA(0) = 0.0D0 XC XC ... COMPUTES THE COEFFICIENTS OF THE XC (1) XC POLYNOMIAL -xi*w(xi) * P (xi) XC k XC X IF (MK .LE. NK) THEN X CALL POLMUL(GAMMA, MK, P1, NK, P, IER) X ELSE X CALL POLMUL(P1, NK, GAMMA, MK, P, IER) X END IF XC XC ... COMPUTES THE COEFFICIENTS OF THE XC POLYNOMIAL P (xi) XC k+1 XC X DO 60 I = 0, NK + NV1 + 1 X P(I) = P(I) + POL(I) X 60 CONTINUE XC XC ... STORES THE COEFFICIENTS OF THE XC (1) XC POLYNOMIAL P (xi) XC k+1 XC X DO 70 I = 0, NK + MK X P1(I) = POL(NM+I) X 70 CONTINUE XC X RETURN X END SHAR_EOF if test 11012 -ne "`wc -c compp1.f`" then echo shar: error transmitting compp1.f '(should have been 11012 characters)' fi echo shar: extracting cresol.f sed 's/^X//' << \SHAR_EOF > cresol.f XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC SUBROUTINE NAME - CRESOL + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC DESCRIPTION: + XC + XC STORES IN SEPARATE VECTORS THE VALUES OBTAINED FOR THE + XC ALPHA , ALPHA', BETA , BETA' BY SOLVING THE AUXILIARY SYSTEM, + XC i i i i + XC AND COMPUTES THE ETA , ETA' , GAMMA AND GAMMA' . + XC i i i i + XC + XC USAGE: + XC CALL CRESOL (RMAT, GAMMA, GAMMAP, ETA, ETAP, P1, NK, MK) + XC + XC ARGUMENTS: + XC + XC RMAT - INPUT REAL MATRIX OF DIMENSION (2,0:*) WHICH CONTAINS + XC THE SOLUTIONS OF THE AUXILIARY SYSTEM. + XC + XC GAMMA - OUTPUT REAL VECTOR CONTAINING THE GAMMA . + XC i + XC + XC GAMMAP - OUTPUT REAL VECTOR CONTAINING THE GAMMA' . + XC i + XC + XC ETA - OUTPUT REAL VECTOR CONTAINING THE ETA . + XC i + XC + XC ETAP - OUTPUT REAL VECTOR CONTAINING THE ETA' . + XC i + XC + XC P1 - INPUT REAL VECTOR CONTAINING THE COEFFICIENTS OF THE + XC (1) + XC POLYNOMIAL P OF DEGREE NK. + XC k + XC + XC NK - INPUT INTEGER EQUAL TO THE DIMENSION OF THE PREVIOUS + XC MATRIX OF THE SYSTEM. + XC + XC MK - INPUT INTEGER EQUAL TO THE DIMENSION OF THE NEW JUMP. + XC + XC EXTERNAL MODULES: + XC + XC DOUBLE PRECISION FUNCTION GAMETA + XC + XC REMARKS: + XC + XC - ALL THE REAL ARGUMENTS PASSED FROM THE CALLING PROGRAM MUST BE + XC IN DOUBLE PRECISION. + XC - IN THE CALLING PROCEDURE THE DIMENSIONS FOR THE MATRIX AND VECTOR + XC ARGUMENTS ARE: + XC RMAT(2,0:*) + XC GAMMA(0:*) + XC GAMMAP(0:*) + XC ETA(0:*) + XC ETAP(0:*) + XC P1(0:*) + XC + XC AUTHORS: + XC + XC C. BREZINSKI + XC UNIVERSITY OF LILLE - FRANCE + XC M. REDIVO ZAGLIA + XC UNIVERSITY OF PADOVA - ITALY + XC + XC REFERENCE: + XC + XC - TREATMENT OF NEAR-BREAKDOWN IN THE CGS ALGORITHM + XC NUMERICAL ALGORITHMS 7 (1994) 33-73 + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC X SUBROUTINE CRESOL (RMAT, GAMMA, GAMMAP, ETA, ETAP, P1, NK, MK) XC XC ... SPECIFICATIONS FOR ARGUMENTS XC X INTEGER NK, MK X DOUBLE PRECISION RMAT(2,0:1), GAMMA(0:1), GAMMAP(0:1), ETA(0:1), X * ETAP(0:1), P1(0:1) XC XC ... SPECIFICATIONS FOR EXTERNAL FUNCTIONS XC X DOUBLE PRECISION GAMETA XC XC ... SPECIFICATIONS FOR LOCAL VARIABLES XC X INTEGER I, I2, NC, NC1, NCC XC XC ... EXECUTABLE STATEMENTS XC XC XC ... CASE MK = 1 XC X IF (MK .EQ. 1) THEN XC XC ... STORES GAMMA = BETA XC 0 0 XC X GAMMA(0) = RMAT(1,0) XC XC ... STORES ETA = ALPHA + P1(NK-1) XC 0 0 XC X ETA(0) = RMAT(2,0) XC XC ... STORES ETA = ALPHA = 1 XC 1 1 XC X ETA(1) = 1.0D0 XC XC ... STORES ETA' = ALPHA' (IF MK <= NK) XC 0 0 XC X IF (MK .LE. NK) ETAP(0) = RMAT(2,1) XC X RETURN X END IF XC XC ... CASE MK <= NK XC X IF (MK .LE. NK) THEN X NC = MK - 2 X NC1 = MK - 1 X ELSE XC XC ... CASE MK > NK XC X NC = NK - 1 X NC1 = NK - 1 X END IF XC XC ... STORES THE BETA AND BETA' XC i i XC AND THE ALPHA AND ALPHA' XC i i XC XC ... BOTH CASES: MK <= NK AND MK > NK XC X IF (NK .NE. 0) THEN X DO 10 I = 0, NC X I2 = I*2 X GAMMA(I) = RMAT(1,I2) X GAMMAP(I) = RMAT(1,I2+1) X ETA(I) = RMAT(2,I2) X ETAP(I) = RMAT(2,I2+1) X 10 CONTINUE X END IF XC X ETA(MK) = 1.0D0 XC XC ... CASE MK <= NK XC X IF (MK .LE. NK) THEN X GAMMA(MK-1) = RMAT(1,2*MK-2) X ETA(MK-1) = RMAT(2,2*MK-2) X ETAP(MK-1) = RMAT(2,2*MK-1) XC XC ... CASE MK > NK XC X ELSE X DO 20 I = NK, MK-1 X NCC = NK + I X GAMMA(I) = RMAT(1,NCC) X ETA(I) = RMAT(2,NCC) X 20 CONTINUE X END IF XC XC ... COMPUTES THE GAMMA AND GAMMA' XC i i XC AND THE ETA AND ETA' XC i i XC XC ... BOTH CASES: MK <= NK AND MK > NK XC X IF (NK .NE. 0) THEN X DO 30 I = 0, NC X I2 = I*2 X RMAT(1,I2) = GAMETA(MK-1,I,GAMMA,P1,NK) X RMAT(1,I2+1) = GAMETA(NC,I,GAMMAP,P1,NK) X RMAT(2,I2) = GAMETA(MK-1,I,ETA,P1,NK) X IF (NK .GE. MK-I) RMAT(2,I2) = RMAT(2,I2) + P1(NK-MK+I) X RMAT(2,I2+1) = GAMETA(NC1,I,ETAP,P1,NK) X 30 CONTINUE X END IF XC XC ... CASE MK <= NK XC X IF (MK .LE. NK) THEN X IF (NK .GE. 1) THEN X RMAT(2,2*MK-2) = ETA(MK-1) + P1(NK-1) X ELSE X RMAT(2,2*MK-2) = ETA(MK-1) X END IF XC XC ... CASE MK > NK XC X ELSE X DO 40 I = NK, MK-1 X NCC = NK + I X RMAT(1,NCC) = GAMETA(MK-1,I,GAMMA,P1,NK) X RMAT(2,NCC) = GAMETA(MK-1,I,ETA,P1,NK) X IF (NK .GE. MK-I) RMAT(2,NCC) = RMAT(2,NCC)+P1(NK-MK+I) X 40 CONTINUE X END IF XC XC ... STORES THE GAMMA AND GAMMA' XC i i XC AND THE ETA AND ETA' XC i i XC XC ... BOTH CASES: MK <= NK AND MK > NK XC X IF (NK .NE. 0) THEN X DO 50 I = 0, NC X I2 = I*2 X GAMMA(I) = RMAT(1,I2) X GAMMAP(I) = RMAT(1,I2+1) X ETA(I) = RMAT(2,I2) X ETAP(I) = RMAT(2,I2+1) X 50 CONTINUE X END IF XC XC ... CASE MK <= NK XC X IF (MK .LE. NK) THEN X ETA(MK-1) = RMAT(2,2*MK-2) XC XC ... CASE MK > NK XC X ELSE X DO 60 I = NK, MK-1 X NCC = NK + I X GAMMA(I) = RMAT(1,NCC) X ETA(I) = RMAT(2,NCC) X 60 CONTINUE X END IF XC X RETURN X END SHAR_EOF if test 10510 -ne "`wc -c cresol.f`" then echo shar: error transmitting cresol.f '(should have been 10510 characters)' fi echo shar: extracting exa1.f sed 's/^X//' << \SHAR_EOF > exa1.f XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC PROGRAM NAME - EXA1 + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC DESCRIPTION: + XC + XC MAIN PROGRAM TO TEST THE SUBROUTINE BSMRZS. + XC + XC MODULES REQUIRED: + XC + XC - SUBROUTINES BSMRZS, CHKSIG, COMPO1, COMPO2, COMPOL, + XC COMPP1, CRESOL, GSOLVD, POLMUL, SOLSHF, + XC STOMHF, STORHF. + XC - DOUBLE PRECISION FUNCTIONS FH, GAMETA, RXSUMM, SSUMM, SUNORM, + XC ZSUMM. + XC + XC MODULES REQUIRED GIVEN IN NETLIB (NA1 FROM NUMERALGO): + XC + XC - SUBROUTINES MATVEC. + XC - DOUBLE PRECISION FUNCTIONS PINNER. + XC + XC REMARK: + XC + XC - IF MAXBCK IS SET EQUAL TO A VALUE LESS THAN 3, THEN THE DIMENSION + XC OF THE WORKING VECTOR LMAT MUST BE SET EQUAL TO (8*MAXBCK). + XC + XC AUTHORS: + XC + XC C. BREZINSKI + XC UNIVERSITY OF LILLE - FRANCE + XC M. REDIVO ZAGLIA + XC UNIVERSITY OF PADOVA - ITALY + XC + XC REFERENCE: + XC + XC - TREATMENT OF NEAR-BREAKDOWN IN THE CGS ALGORITHM + XC NUMERICAL ALGORITHMS 7 (1994) 33-73 + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC X PROGRAM EXA1 XC XC ... SPECIFICATIONS FOR VARIABLES AND CONSTANTS XC X INTEGER I, IER, IFLAG, INIT, IR, J, K, KSAV, MAXBCK, MAXN, X * N, NBC, NK, NKSAV X DOUBLE PRECISION EPS, EPS1, EPS2, RBEST, TMP, ACONST, NORM XC XC ... SPECIFICATIONS FOR PARAMETER CONSTANTS XC X PARAMETER (N=40, IR=40, MAXBCK=40, MAXN=40, NBC=40, X * EPS=1.0D-8, EPS1=1.0D-8, EPS2=1.D-9, ACONST=0.95) XC XC ... SPECIFICATIONS FOR VECTORS AND MATRICES XC X DOUBLE PRECISION A(IR,IR), AB(IR), Y(IR), X(IR), R(IR), X * C(0:2*MAXBCK-1), D(0:2*MAXBCK-1), V(IR,0:2*MAXBCK), X * W(IR,0:2*MAXBCK-2), U(IR,0:2*MAXBCK-1), X * RMAT(2,0:2*MAXBCK-1), LMAT((2*MAXBCK-1)*(2*MAXBCK-1)), X * ETA(0:MAXBCK), ETAP(0:MAXBCK-1), GAMMA(0:MAXBCK), X * GAMMAP(0:MAXBCK-1), P(0:MAXN), P1(0:MAXN), PWK(0:2*MAXN+1), X * ADIG(IR), XBEST(IR), ACT(IR), CB(IR) XC XC ... SPECIFICATIONS FOR EXTERNAL FUNCTIONS XC X DOUBLE PRECISION SUNORM XC XC ... EXECUTABLE STATEMENTS XC XC XC ... CHOICE OF IFLAG XC XC IFLAG .EQ. 0 Y IS SET TO r BY THE SUBROUTINE XC 0 XC IFLAG .NE. 0 THE USER MUST DEFINE Y IN XC THE MAIN PROGRAM XC X IFLAG = 0 XC X WRITE (*,*) ' OUTPUT RESULTS FROM EXA1' X WRITE (*,*) X WRITE (*,*) ' DIMENSION OF THE SYSTEM = ',N X WRITE (*,*) ' NUMBER OF CALLS = ',NBC X WRITE (*,*) ' MAXBCK (MAX JUMP MK) = ',MAXBCK X WRITE (*,*) ' MAXN (MAX NK+MK) = ',MAXN X WRITE (*,*) ' IFLAG = ',IFLAG X WRITE (*,*) ' EPS (NEAR-BREAKDOWN) = ',EPS X WRITE (*,*) ' EPS1 (GAUSS) = ',EPS1 X WRITE (*,*) ' EPS2 (SOLUTION) = ',EPS2 X WRITE (*,*) ' ACONST = ',ACONST X WRITE (*,*) X WRITE (*,*) ' ITER. NK NORM' X DO 20 I = 1, N X X(I) = 0.0D0 X Y(I) = 1.0D0 X DO 10 J = 1, N X A(I,J) = 0.0D0 X 10 CONTINUE X 20 CONTINUE XC X A(N,1) = ACONST X AB(N) = ACONST X CB(N) = AB(N) X DO 30 I = 2, N X A(I-1,I) = 1.0D0 X AB(I-1) = 1.0D0 X CB(I-1) = AB(I-1) X 30 CONTINUE XC X INIT = 0 XC X DO 50 K = 0, NBC X CALL BSMRZS (N, A, AB, Y, X, R, NORM, MAXBCK, MAXN, V, W, X * U, LMAT, RMAT, ETA, ETAP, GAMMA, GAMMAP, D, C, P, P1, X * PWK, IR, NK, EPS, EPS1, EPS2, INIT, IFLAG, IER) X IF (K .EQ. 0 .OR. NORM .LT. RBEST) THEN X RBEST = NORM X KSAV = K X NKSAV = NK X DO 40 I = 1, N X XBEST(I) = X(I) X 40 CONTINUE X END IF X WRITE (*,100) K, NK, NORM X IF (IER .EQ. 400) THEN X WRITE (*,*) X WRITE (*,*) ' SOLUTION NOT OBTAINED AFTER REACHING' X WRITE (*,*) ' THE DIMENSION N OF THE SYSTEM. THE' X WRITE (*,*) ' COMPUTATION CONTINUES UNTIL MAXN AT MOST' X WRITE (*,*) X IER = 0 X ELSE IF (IER .GE. 700) THEN X WRITE (*,*) X WRITE (*,*) ' JUMP FOR BREAKDOWN OR NEAR-BREAKDOWN ',IER X WRITE (*,*) X IER = 0 X END IF X IF (IER .NE. 0) THEN X WRITE (*,*) X IF (IER .EQ. 200) THEN X WRITE (*,*) ' THE SOLUTION HAS BEEN OBTAINED' X ELSE IF (IER .EQ. 300) THEN X WRITE (*,*) ' SOLUTION NOT OBTAINED AFTER' X WRITE (*,*) ' REACHING THE VALUE OF MAXN' X ELSE IF (IER .EQ. 600) THEN X WRITE (*,*) ' THE VALUE OF THE JUMP MK' X WRITE (*,*) ' IS GREATER THAN MAXBCK' X ELSE X WRITE (*,*) ' STOP DUE TO ERROR ', IER X STOP X END IF X GO TO 60 X END IF X 50 CONTINUE X WRITE (*,*) X WRITE (*,300) ' SOLUTION NOT OBTAINED AFTER', NBC X WRITE (*,*) ' ITERATIONS OF THE SUBROUTINE' X 60 CONTINUE X WRITE (*,*) X DO 70 I= 1, N X TMP = DABS(X(I)-1.0D0) X IF (TMP .EQ. 0.0D0) THEN X ADIG(I) = 999.00 X ELSE X ADIG(I) = -DLOG10(TMP) X END IF X 70 CONTINUE XC XC ... COMPUTE THE ACTUAL RESIDUAL XC X CALL MATVEC (N, A, X, 1, ACT, 1, IR, 0) X DO 80 I = 1, N X ACT(I) = ACT(I)-CB(I) X 80 CONTINUE X WRITE(*,500) ' ACTUAL RESIDUAL',SUNORM(N,ACT) X WRITE (*,*) XC X IF (IER .EQ. 200) THEN X WRITE (*,*) X * ' FINAL SOLUTION N. OF SIGN. DIG. ' X WRITE (*,200) (X(I),ADIG(I),I=1,N) X ELSE X WRITE (*,*) X * ' LAST SOLUTION N. OF SIGN. DIG. ' X WRITE (*,200) (X(I),ADIG(I),I=1,N) X WRITE (*,*) X WRITE (*,*) ' SOLUTION OF MINIMAL RESIDUAL OBTAINED AT' X WRITE (*,*) ' ITER. NK NORM' X WRITE (*,100) KSAV, NKSAV, RBEST X WRITE (*,*) X WRITE (*,400) (XBEST(I),I=1,N) X END IF X STOP X 100 FORMAT (I5,I8,D28.15) X 200 FORMAT (D25.15,F15.2) X 300 FORMAT (T2,A,I4) X 400 FORMAT (D25.15) X 500 FORMAT (T2,A,D24.15) X END SHAR_EOF if test 8688 -ne "`wc -c exa1.f`" then echo shar: error transmitting exa1.f '(should have been 8688 characters)' fi echo shar: extracting exa2.f sed 's/^X//' << \SHAR_EOF > exa2.f XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC PROGRAM NAME - EXA2 + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC DESCRIPTION: + XC + XC MAIN PROGRAM TO TEST THE SUBROUTINE BSMRZS. + XC + XC MODULES REQUIRED: + XC + XC - SUBROUTINES BSMRZS, CHKSIG, COMPO1, COMPO2, COMPOL, + XC COMPP1, CRESOL, GSOLVD, POLMUL, SOLSHF, + XC STOMHF, STORHF. + XC - DOUBLE PRECISION FUNCTIONS FH, GAMETA, RXSUMM, SSUMM, SUNORM, + XC ZSUMM. + XC + XC MODULES REQUIRED GIVEN IN NETLIB (NA1 FROM NUMERALGO): + XC + XC - SUBROUTINES MATVEC. + XC - DOUBLE PRECISION FUNCTIONS PINNER. + XC + XC REMARK: + XC + XC - IF MAXBCK IS SET EQUAL TO A VALUE LESS THAN 3, THEN THE DIMENSION + XC OF THE WORKING VECTOR LMAT MUST BE SET EQUAL TO (8*MAXBCK). + XC + XC AUTHORS: + XC + XC C. BREZINSKI + XC UNIVERSITY OF LILLE - FRANCE + XC M. REDIVO ZAGLIA + XC UNIVERSITY OF PADOVA - ITALY + XC + XC REFERENCE: + XC + XC - TREATMENT OF NEAR-BREAKDOWN IN THE CGS ALGORITHM + XC NUMERICAL ALGORITHMS 7 (1994) 33-73 + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC X PROGRAM EXA2 X XC XC ... SPECIFICATIONS FOR VARIABLES AND CONSTANTS XC X INTEGER I, IER, IFLAG, INIT, IR, J, K, KSAV, MAXBCK, MAXN, X * N, NBC, NK, NKSAV X DOUBLE PRECISION EPS, EPS1, EPS2, RBEST, TMP, ACONST, SOL, NORM XC XC ... SPECIFICATIONS FOR PARAMETER CONSTANTS XC X PARAMETER (N=40, IR=40, MAXBCK=10, MAXN=40, NBC=40, X * EPS=1.D-4, EPS1=1.D-12, EPS2=1.D-14, ACONST=1.D-4) XC XC ... SPECIFICATIONS FOR VECTORS AND MATRICES XC X DOUBLE PRECISION A(IR,IR), AB(IR), Y(IR), X(IR), R(IR), X * C(0:2*MAXBCK-1), D(0:2*MAXBCK-1), V(IR,0:2*MAXBCK), X * W(IR,0:2*MAXBCK-2), U(IR,0:2*MAXBCK-1), X * RMAT(2,0:2*MAXBCK-1), LMAT((2*MAXBCK-1)*(2*MAXBCK-1)), X * ETA(0:MAXBCK), ETAP(0:MAXBCK-1), GAMMA(0:MAXBCK), X * GAMMAP(0:MAXBCK-1), P(0:MAXN), P1(0:MAXN), PWK(0:2*MAXN+1), X * ADIG(IR), XBEST(IR), SOLU(IR), ACT(IR), CB(IR) XC XC ... SPECIFICATIONS FOR EXTERNAL FUNCTIONS XC X DOUBLE PRECISION SUNORM XC XC ... EXECUTABLE STATEMENTS XC XC XC ... CHOICE OF IFLAG XC XC IFLAG .EQ. 0 Y IS SET TO r BY THE SUBROUTINE XC 0 XC IFLAG .NE. 0 THE USER MUST DEFINE Y IN XC THE MAIN PROGRAM XC X IFLAG = 0 XC X WRITE (*,*) ' OUTPUT RESULTS FROM EXA2' X WRITE (*,*) X WRITE (*,*) ' DIMENSION OF THE SYSTEM = ',N X WRITE (*,*) ' NUMBER OF CALLS = ',NBC X WRITE (*,*) ' MAXBCK (MAX JUMP MK) = ',MAXBCK X WRITE (*,*) ' MAXN (MAX NK+MK) = ',MAXN X WRITE (*,*) ' IFLAG = ',IFLAG X WRITE (*,*) ' EPS (NEAR-BREAKDOWN) = ',EPS X WRITE (*,*) ' EPS1 (GAUSS) = ',EPS1 X WRITE (*,*) ' EPS2 (SOLUTION) = ',EPS2 X WRITE (*,*) ' ACONST = ',ACONST X WRITE (*,*) X WRITE (*,*) ' ITER. NK NORM' X DO 20 I = 1, N X X(I) = 0.0D0 X Y(I) = 1.0D0 X AB(I) = 0.0D0 X DO 10 J = 1, N X A(I,J) = 0.0D0 X 10 CONTINUE X 20 CONTINUE XC X AB(1) = 5.0D0 X AB(2) = -3.0D0 X AB(3) = 4.0D0 X AB(4) = -4.0D0 X DO 30 I = 1, N X CB(I) = AB(I) X A(I,I) = -(-1.0D0)**I X 30 CONTINUE X DO 40 I = 1, N/2 X A(2*I-1, 2*I) = I - 1.0D0 + ACONST X 40 CONTINUE XC X INIT = 0 XC X DO 60 K = 0, NBC X CALL BSMRZS (N, A, AB, Y, X, R, NORM, MAXBCK, MAXN, V, W, X * U, LMAT, RMAT, ETA, ETAP, GAMMA, GAMMAP, D, C, P, P1, X * PWK, IR, NK, EPS, EPS1, EPS2, INIT, IFLAG, IER) X IF (K .EQ. 0 .OR. NORM .LT. RBEST) THEN X RBEST = NORM X KSAV = K X NKSAV = NK X DO 50 I = 1, N X XBEST(I) = X(I) X 50 CONTINUE X END IF X WRITE (*,100) K, NK, NORM X IF (IER .EQ. 400) THEN X WRITE (*,*) X WRITE (*,*) ' SOLUTION NOT OBTAINED AFTER REACHING' X WRITE (*,*) ' THE DIMENSION N OF THE SYSTEM. THE' X WRITE (*,*) ' COMPUTATION CONTINUES UNTIL MAXN AT MOST' X WRITE (*,*) X IER = 0 X ELSE IF (IER .GE. 700) THEN X WRITE (*,*) X WRITE (*,*) ' JUMP FOR BREAKDOWN OR NEAR-BREAKDOWN ',IER X WRITE (*,*) X IER = 0 X END IF X IF (IER .NE. 0) THEN X WRITE (*,*) X IF (IER .EQ. 200) THEN X WRITE (*,*) ' THE SOLUTION HAS BEEN OBTAINED' X ELSE IF (IER .EQ. 300) THEN X WRITE (*,*) ' SOLUTION NOT OBTAINED AFTER' X WRITE (*,*) ' REACHING THE VALUE OF MAXN' X ELSE IF (IER .EQ. 600) THEN X WRITE (*,*) ' THE VALUE OF THE JUMP MK' X WRITE (*,*) ' IS GREATER THAN MAXBCK' X ELSE X WRITE (*,*) ' STOP DUE TO ERROR ', IER X STOP X END IF X GO TO 70 X END IF X 60 CONTINUE X WRITE (*,*) X WRITE (*,300) ' SOLUTION NOT OBTAINED AFTER', NBC X WRITE (*,*) ' ITERATIONS OF THE SUBROUTINE' XC X 70 CONTINUE X WRITE (*,*) XC XC ... COMPUTE THE EXACT SOLUTION XC X DO 80 I = 1, N/2 X SOLU(2*I) = -CB(2*I) X SOLU(2*I-1) = CB(2*I-1) + (I-1.0D0+ACONST) * CB(2*I) X 80 CONTINUE X DO 90 I= 1, N X TMP = DABS(X(I)-SOLU(I)) X IF (TMP .EQ. 0.0D0) THEN X ADIG(I) = 999.00 X ELSE X SOL = SOLU(I) X IF (SOLU(I) .EQ. 0.0D0) SOL=1.0D0 X ADIG(I) = -DLOG10(TMP/DABS(SOL)) X END IF X 90 CONTINUE XC XC ... COMPUTE THE ACTUAL RESIDUAL XC X CALL MATVEC (N, A, X, 1, ACT, 1, IR, 0) X DO 95 I = 1, N X ACT(I) = ACT(I)-CB(I) X 95 CONTINUE X WRITE(*,500) ' ACTUAL RESIDUAL',SUNORM(N,ACT) X WRITE (*,*) XC X IF (IER .EQ. 200) THEN X WRITE (*,*) X * ' FINAL SOLUTION N. OF SIGN. DIG. ' X WRITE (*,200) (X(I),ADIG(I),I=1,N) X ELSE X WRITE (*,*) X * ' LAST SOLUTION N. OF SIGN. DIG. ' X WRITE (*,200) (X(I),ADIG(I),I=1,N) X WRITE (*,*) X WRITE (*,*) ' SOLUTION OF MINIMAL RESIDUAL OBTAINED AT' X WRITE (*,*) ' ITER. NK NORM' X WRITE (*,100) KSAV, NKSAV, RBEST X WRITE (*,*) X WRITE (*,400) (XBEST(I),I=1,N) X END IF X STOP X 100 FORMAT (I5,I8,D28.15) X 200 FORMAT (D25.15,F15.2) X 300 FORMAT (T2,A,I4) X 400 FORMAT (D25.15) X 500 FORMAT (T2,A,D24.15) X END SHAR_EOF if test 9082 -ne "`wc -c exa2.f`" then echo shar: error transmitting exa2.f '(should have been 9082 characters)' fi echo shar: extracting exa3.f sed 's/^X//' << \SHAR_EOF > exa3.f XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC PROGRAM NAME - EXA3 + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC DESCRIPTION: + XC + XC MAIN PROGRAM TO TEST THE SUBROUTINE BSMRZS. + XC + XC MODULES REQUIRED: + XC + XC - SUBROUTINES BSMRZS, CHKSIG, COMPO1, COMPO2, COMPOL, + XC COMPP1, CRESOL, GSOLVD, POLMUL, SOLSHF, + XC STOMHF, STORHF. + XC - DOUBLE PRECISION FUNCTIONS FH, GAMETA, RXSUMM, SSUMM, SUNORM, + XC ZSUMM. + XC + XC MODULES REQUIRED GIVEN IN NETLIB (NA1 FROM NUMERALGO): + XC + XC - SUBROUTINES MATVEC. + XC - DOUBLE PRECISION FUNCTIONS PINNER. + XC + XC REMARK: + XC + XC - IF MAXBCK IS SET EQUAL TO A VALUE LESS THAN 3, THEN THE DIMENSION + XC OF THE WORKING VECTOR LMAT MUST BE SET EQUAL TO (8*MAXBCK). + XC + XC AUTHORS: + XC + XC C. BREZINSKI + XC UNIVERSITY OF LILLE - FRANCE + XC M. REDIVO ZAGLIA + XC UNIVERSITY OF PADOVA - ITALY + XC + XC REFERENCE: + XC + XC - TREATMENT OF NEAR-BREAKDOWN IN THE CGS ALGORITHM + XC NUMERICAL ALGORITHMS 7 (1994) 33-73 + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + X PROGRAM EXA3 X XC XC ... SPECIFICATIONS FOR VARIABLES AND CONSTANTS XC X INTEGER I, IER, IFLAG, INIT, IR, J, K, KSAV, MAXBCK, MAXN, X * N, NBC, NK, NKSAV X DOUBLE PRECISION EPS, EPS1, EPS2, RBEST, TMP, ACONST, SOL, NORM XC XC ... SPECIFICATIONS FOR PARAMETER CONSTANTS XC X PARAMETER (N=40, IR=40, MAXBCK=10, MAXN=40, NBC=40, X * EPS=1.0D-8, EPS1=1.0D-13, EPS2=1.D-13, ACONST=1.D-9) XC XC XC ... SPECIFICATIONS FOR VECTORS AND MATRICES XC X DOUBLE PRECISION A(IR,IR), AB(IR), Y(IR), X(IR), R(IR), X * C(0:2*MAXBCK-1), D(0:2*MAXBCK-1), V(IR,0:2*MAXBCK), X * W(IR,0:2*MAXBCK-2), U(IR,0:2*MAXBCK-1), X * RMAT(2,0:2*MAXBCK-1), LMAT((2*MAXBCK-1)*(2*MAXBCK-1)), X * ETA(0:MAXBCK), ETAP(0:MAXBCK-1), GAMMA(0:MAXBCK), X * GAMMAP(0:MAXBCK-1), P(0:MAXN), P1(0:MAXN), PWK(0:2*MAXN+1), X * ADIG(IR), XBEST(IR), SOLU(IR), ACT(IR), CB(IR) XC XC ... SPECIFICATIONS FOR EXTERNAL FUNCTIONS XC X DOUBLE PRECISION SUNORM XC XC ... EXECUTABLE STATEMENTS XC XC XC ... CHOICE OF IFLAG XC XC IFLAG .EQ. 0 Y IS SET TO r BY THE SUBROUTINE XC 0 XC IFLAG .NE. 0 THE USER MUST DEFINE Y IN XC THE MAIN PROGRAM XC X IFLAG = 0 XC X WRITE (*,*) ' OUTPUT RESULTS FROM EXA3' X WRITE (*,*) X WRITE (*,*) ' DIMENSION OF THE SYSTEM = ',N X WRITE (*,*) ' NUMBER OF CALLS = ',NBC X WRITE (*,*) ' MAXBCK (MAX JUMP MK) = ',MAXBCK X WRITE (*,*) ' MAXN (MAX NK+MK) = ',MAXN X WRITE (*,*) ' IFLAG = ',IFLAG X WRITE (*,*) ' EPS (NEAR-BREAKDOWN) = ',EPS X WRITE (*,*) ' EPS1 (GAUSS) = ',EPS1 X WRITE (*,*) ' EPS2 (SOLUTION) = ',EPS2 X WRITE (*,*) ' ACONST = ',ACONST X WRITE (*,*) X WRITE (*,*) ' ITER. NK NORM' X DO 20 I = 1, N X X(I) = 0.0D0 X Y(I) = 1.0D0 X AB(I) = DBLE(I) X CB(I) = AB(I) X DO 10 J = 1, N X A(I,J) = 0.0D0 X 10 CONTINUE X 20 CONTINUE XC X DO 30 I = 1, N/2 X A(2*I-1, 2*I) = 1.0D0 + I*ACONST X A(2*I, 2*I-1) = -1.0D0 X 30 CONTINUE XC X INIT = 0 XC X DO 50 K = 0, NBC X CALL BSMRZS (N, A, AB, Y, X, R, NORM, MAXBCK, MAXN, V, W, X * U, LMAT, RMAT, ETA, ETAP, GAMMA, GAMMAP, D, C, P, P1, X * PWK, IR, NK, EPS, EPS1, EPS2, INIT, IFLAG, IER) X IF (K .EQ. 0 .OR. NORM .LT. RBEST) THEN X RBEST = NORM X KSAV = K X NKSAV = NK X DO 40 I = 1, N X XBEST(I) = X(I) X 40 CONTINUE X END IF X WRITE (*,100) K, NK, NORM X IF (IER .EQ. 400) THEN X WRITE (*,*) X WRITE (*,*) ' SOLUTION NOT OBTAINED AFTER REACHING' X WRITE (*,*) ' THE DIMENSION N OF THE SYSTEM. THE' X WRITE (*,*) ' COMPUTATION CONTINUES UNTIL MAXN AT MOST' X WRITE (*,*) X IER = 0 X ELSE IF (IER .GE. 700) THEN X WRITE (*,*) X WRITE (*,*) ' JUMP FOR BREAKDOWN OR NEAR-BREAKDOWN ',IER X WRITE (*,*) X IER = 0 X END IF X IF (IER .NE. 0) THEN X WRITE (*,*) X IF (IER .EQ. 200) THEN X WRITE (*,*) ' THE SOLUTION HAS BEEN OBTAINED' X ELSE IF (IER .EQ. 300) THEN X WRITE (*,*) ' SOLUTION NOT OBTAINED AFTER' X WRITE (*,*) ' REACHING THE VALUE OF MAXN' X ELSE IF (IER .EQ. 600) THEN X WRITE (*,*) ' THE VALUE OF THE JUMP MK' X WRITE (*,*) ' IS GREATER THAN MAXBCK' X ELSE X WRITE (*,*) ' STOP DUE TO ERROR ', IER X STOP X END IF X GO TO 60 X END IF X 50 CONTINUE X WRITE (*,*) X WRITE (*,300) ' SOLUTION NOT OBTAINED AFTER', NBC X WRITE (*,*) ' ITERATIONS OF THE SUBROUTINE' XC X 60 CONTINUE X WRITE (*,*) XC XC ... COMPUTE THE EXACT SOLUTION XC X DO 80 I = 1, N/2 X SOLU(2*I-1) = -CB(2*I) X SOLU(2*I) = CB(2*I-1) /(1.0D0 + I*ACONST) X 80 CONTINUE X DO 90 I= 1, N X TMP = DABS(X(I)-SOLU(I)) X IF (TMP .EQ. 0.0D0) THEN X ADIG(I) = 999.00 X ELSE X SOL = SOLU(I) X IF (SOLU(I) .EQ. 0.0D0) SOL=1.0D0 X ADIG(I) = -DLOG10(TMP/DABS(SOL)) X END IF X 90 CONTINUE XC XC ... COMPUTE THE ACTUAL RESIDUAL XC X CALL MATVEC (N, A, X, 1, ACT, 1, IR, 0) X DO 95 I = 1, N X ACT(I) = ACT(I)-CB(I) X 95 CONTINUE X WRITE(*,500) ' ACTUAL RESIDUAL',SUNORM(N,ACT) X WRITE (*,*) XC X IF (IER .EQ. 200) THEN X WRITE (*,*) X * ' FINAL SOLUTION N. OF SIGN. DIG. ' X WRITE (*,200) (X(I),ADIG(I),I=1,N) X ELSE X WRITE (*,*) X * ' LAST SOLUTION N. OF SIGN. DIG. ' X WRITE (*,200) (X(I),ADIG(I),I=1,N) X WRITE (*,*) X WRITE (*,*) ' SOLUTION OF MINIMAL RESIDUAL OBTAINED AT' X WRITE (*,*) ' ITER. NK NORM' X WRITE (*,100) KSAV, NKSAV, RBEST X WRITE (*,*) X WRITE (*,400) (XBEST(I),I=1,N) X END IF X STOP X 100 FORMAT (I5,I8,D28.15) X 200 FORMAT (D25.15,F15.2) X 300 FORMAT (T2,A,I4) X 400 FORMAT (D25.15) X 500 FORMAT (T2,A,D24.15) X END X SHAR_EOF if test 8964 -ne "`wc -c exa3.f`" then echo shar: error transmitting exa3.f '(should have been 8964 characters)' fi echo shar: extracting exa4.f sed 's/^X//' << \SHAR_EOF > exa4.f XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC PROGRAM NAME - EXA4 + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC DESCRIPTION: + XC + XC MAIN PROGRAM TO TEST THE SUBROUTINE BSMRZS. + XC + XC MODULES REQUIRED: + XC + XC - SUBROUTINES BSMRZS, CHKSIG, COMPO1, COMPO2, COMPOL, + XC COMPP1, CRESOL, GSOLVD, POLMUL, SOLSHF, + XC STOMHF, STORHF. + XC - DOUBLE PRECISION FUNCTIONS FH, GAMETA, RXSUMM, SSUMM, SUNORM, + XC ZSUMM. + XC + XC MODULES REQUIRED GIVEN IN NETLIB (NA1 FROM NUMERALGO): + XC + XC - SUBROUTINES MATVEC. + XC - DOUBLE PRECISION FUNCTIONS PINNER. + XC + XC REMARK: + XC + XC - IF MAXBCK IS SET EQUAL TO A VALUE LESS THAN 3, THEN THE DIMENSION + XC OF THE WORKING VECTOR LMAT MUST BE SET EQUAL TO (8*MAXBCK). + XC + XC AUTHORS: + XC + XC C. BREZINSKI + XC UNIVERSITY OF LILLE - FRANCE + XC M. REDIVO ZAGLIA + XC UNIVERSITY OF PADOVA - ITALY + XC + XC REFERENCE: + XC + XC - TREATMENT OF NEAR-BREAKDOWN IN THE CGS ALGORITHM + XC NUMERICAL ALGORITHMS 7 (1994) 33-73 + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC X PROGRAM EXA4 XC XC ... SPECIFICATIONS FOR VARIABLES AND CONSTANTS XC X INTEGER I, IER, IFLAG, INIT, IR, J, K, KSAV, MAXBCK, MAXN, X * N, NBC, NK, NKSAV X DOUBLE PRECISION EPS, EPS1, EPS2, RBEST, TMP, NORM, ACONST XC XC ... SPECIFICATIONS FOR PARAMETER CONSTANTS XC X PARAMETER (N=200, IR=200, MAXBCK=10, MAXN=200, NBC=200, X * EPS=1.0D-5, EPS1=1.0D-14, EPS2=1.0D-9, ACONST=1.D-8) XC XC ... SPECIFICATIONS FOR VECTORS AND MATRICES XC X DOUBLE PRECISION A(IR,IR), AB(IR), Y(IR), X(IR), R(IR), X * C(0:2*MAXBCK-1), D(0:2*MAXBCK-1), V(IR,0:2*MAXBCK), X * W(IR,0:2*MAXBCK-2), U(IR,0:2*MAXBCK-1), X * RMAT(2,0:2*MAXBCK-1), LMAT((2*MAXBCK-1)*(2*MAXBCK-1)), X * ETA(0:MAXBCK), ETAP(0:MAXBCK-1), GAMMA(0:MAXBCK), X * GAMMAP(0:MAXBCK-1), P(0:MAXN), P1(0:MAXN), PWK(0:2*MAXN+1), X * ADIG(IR), XBEST(IR), ACT(IR), CB(IR) XC XC ... SPECIFICATIONS FOR EXTERNAL FUNCTIONS XC X DOUBLE PRECISION SUNORM XC XC ... EXECUTABLE STATEMENTS XC XC XC ... CHOICE OF IFLAG XC XC IFLAG .EQ. 0 Y IS SET TO r BY THE SUBROUTINE XC 0 XC IFLAG .NE. 0 THE USER MUST DEFINE Y IN XC THE MAIN PROGRAM XC X IFLAG = 0 XC X WRITE (*,*) ' OUTPUT RESULTS FROM EXA4' X WRITE (*,*) X WRITE (*,*) ' DIMENSION OF THE SYSTEM = ',N X WRITE (*,*) ' NUMBER OF CALLS = ',NBC X WRITE (*,*) ' MAXBCK (MAX JUMP MK) = ',MAXBCK X WRITE (*,*) ' MAXN (MAX NK+MK) = ',MAXN X WRITE (*,*) ' IFLAG = ',IFLAG X WRITE (*,*) ' EPS (NEAR-BREAKDOWN) = ',EPS X WRITE (*,*) ' EPS1 (GAUSS) = ',EPS1 X WRITE (*,*) ' EPS2 (SOLUTION) = ',EPS2 X WRITE (*,*) ' A CONSTANT = ',ACONST X WRITE (*,*) X WRITE (*,*) ' ITER. NK NORM' X DO 20 I = 1, N X X(I) = 0.0D0 X Y(I) = 1.0D0 X DO 10 J = 1, N X A(I,J) = 0.0D0 X 10 CONTINUE X 20 CONTINUE XC X A(1,1) = ACONST X AB(1) = 1.0D0 + ACONST X CB(1) = AB(1) X AB(N) = -1.0D0 + ACONST X CB(N) = AB(N) X DO 30 I = 2, N X A(I-1,I) = 1.0D0 X A(I,I-1) = -1.0D0 X A(I,I) = ACONST X IF (I .NE. N) THEN X AB(I) = ACONST X CB(I) = AB(I) X END IF X 30 CONTINUE XC X INIT = 0 XC X DO 50 K = 0, NBC X CALL BSMRZS (N, A, AB, Y, X, R, NORM, MAXBCK, MAXN, V, W, X * U, LMAT, RMAT, ETA, ETAP, GAMMA, GAMMAP, D, C, P, P1, X * PWK, IR, NK, EPS, EPS1, EPS2, INIT, IFLAG, IER) X IF (K .EQ. 0 .OR. NORM .LT. RBEST) THEN X RBEST = NORM X KSAV = K X NKSAV = NK X DO 40 I = 1, N X XBEST(I) = X(I) X 40 CONTINUE X END IF X WRITE (*,100) K, NK, NORM X IF (IER .EQ. 400) THEN X WRITE (*,*) X WRITE (*,*) ' SOLUTION NOT OBTAINED AFTER REACHING' X WRITE (*,*) ' THE DIMENSION N OF THE SYSTEM. THE' X WRITE (*,*) ' COMPUTATION CONTINUES UNTIL MAXN AT MOST' X WRITE (*,*) X IER = 0 X ELSE IF (IER .GE. 700) THEN X WRITE (*,*) X WRITE (*,*) ' JUMP FOR BREAKDOWN OR NEAR-BREAKDOWN ',IER X WRITE (*,*) X IER = 0 X END IF X IF (IER .NE. 0) THEN X WRITE (*,*) X IF (IER .EQ. 200) THEN X WRITE (*,*) ' THE SOLUTION HAS BEEN OBTAINED' X ELSE IF (IER .EQ. 300) THEN X WRITE (*,*) ' SOLUTION NOT OBTAINED AFTER' X WRITE (*,*) ' REACHING THE VALUE OF MAXN' X ELSE IF (IER .EQ. 600) THEN X WRITE (*,*) ' THE VALUE OF THE JUMP MK' X WRITE (*,*) ' IS GREATER THAN MAXBCK' X ELSE X WRITE (*,*) ' STOP DUE TO ERROR ', IER X STOP X END IF X GO TO 60 X END IF X 50 CONTINUE X WRITE (*,*) X WRITE (*,300) ' SOLUTION NOT OBTAINED AFTER', NBC X WRITE (*,*) ' ITERATIONS OF THE SUBROUTINE' XC X 60 CONTINUE X WRITE (*,*) X DO 70 I= 1, N X TMP = DABS(X(I)-1.0D0) X IF (TMP .EQ. 0.0D0) THEN X ADIG(I) = 999.00 X ELSE X ADIG(I) = -DLOG10(TMP/1.0D0) X END IF X 70 CONTINUE XC XC ... COMPUTE THE ACTUAL RESIDUAL XC X CALL MATVEC (N, A, X, 1, ACT, 1, IR, 0) X DO 80 I = 1, N X ACT(I) = ACT(I)-CB(I) X 80 CONTINUE X WRITE(*,500) ' ACTUAL RESIDUAL',SUNORM(N,ACT) X WRITE (*,*) XC X IF (IER .EQ. 200) THEN X WRITE (*,*) X * ' FINAL SOLUTION N. OF SIGN. DIG. ' X WRITE (*,200) (X(I),ADIG(I),I=1,N) X ELSE X WRITE (*,*) X * ' LAST SOLUTION N. OF SIGN. DIG. ' X WRITE (*,200) (X(I),ADIG(I),I=1,N) X WRITE (*,*) X WRITE (*,*) ' SOLUTION OF MINIMAL RESIDUAL OBTAINED AT' X WRITE (*,*) ' ITER. NK NORM' X WRITE (*,100) KSAV, NKSAV, RBEST X WRITE (*,*) X WRITE (*,400) (XBEST(I),I=1,N) X END IF X STOP X 100 FORMAT (I5,I8,D28.15) X 200 FORMAT (D25.15,F15.2) X 300 FORMAT (T2,A,I4) X 400 FORMAT (D25.15) X 500 FORMAT (T2,A,D24.15) X END SHAR_EOF if test 8860 -ne "`wc -c exa4.f`" then echo shar: error transmitting exa4.f '(should have been 8860 characters)' fi echo shar: extracting exa5.f sed 's/^X//' << \SHAR_EOF > exa5.f XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC PROGRAM NAME - EXA5 + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC DESCRIPTION: + XC + XC MAIN PROGRAM TO TEST THE SUBROUTINE BSMRZS. + XC + XC MODULES REQUIRED: + XC + XC - SUBROUTINES BSMRZS, CHKSIG, COMPO1, COMPO2, COMPOL, + XC COMPP1, CRESOL, GSOLVD, POLMUL, SOLSHF, + XC STOMHF, STORHF. + XC - DOUBLE PRECISION FUNCTIONS FH, GAMETA, RXSUMM, SSUMM, SUNORM, + XC ZSUMM. + XC + XC MODULES REQUIRED GIVEN IN NETLIB (NA1 FROM NUMERALGO): + XC + XC - SUBROUTINES MATVEC. + XC - DOUBLE PRECISION FUNCTIONS PINNER. + XC + XC REMARK: + XC + XC - IF MAXBCK IS SET EQUAL TO A VALUE LESS THAN 3, THEN THE DIMENSION + XC OF THE WORKING VECTOR LMAT MUST BE SET EQUAL TO (8*MAXBCK). + XC + XC AUTHORS: + XC + XC C. BREZINSKI + XC UNIVERSITY OF LILLE - FRANCE + XC M. REDIVO ZAGLIA + XC UNIVERSITY OF PADOVA - ITALY + XC + XC REFERENCE: + XC + XC - TREATMENT OF NEAR-BREAKDOWN IN THE CGS ALGORITHM + XC NUMERICAL ALGORITHMS 7 (1994) 33-73 + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC X PROGRAM EXA5 XC XC ... SPECIFICATIONS FOR VARIABLES AND CONSTANTS XC X INTEGER I, IER, IFLAG, INIT, IR, J, K, KSAV, MAXBCK, MAXN, X * N, NBC, NK, NKSAV X DOUBLE PRECISION EPS, EPS1, EPS2, RBEST, TMP, NORM XC XC ... SPECIFICATIONS FOR PARAMETER CONSTANTS XC X PARAMETER (N=401, IR=401, MAXBCK=10, MAXN=401, NBC=100, X * EPS=1.D-7, EPS1=1.D-7, EPS2=1.D-14) XC XC ... SPECIFICATIONS FOR VECTORS AND MATRICES XC X DOUBLE PRECISION A(IR,IR), AB(IR), Y(IR), X(IR), R(IR), X * C(0:2*MAXBCK-1), D(0:2*MAXBCK-1), V(IR,0:2*MAXBCK), X * W(IR,0:2*MAXBCK-2), U(IR,0:2*MAXBCK-1), X * RMAT(2,0:2*MAXBCK-1), LMAT((2*MAXBCK-1)*(2*MAXBCK-1)), X * ETA(0:MAXBCK), ETAP(0:MAXBCK-1), GAMMA(0:MAXBCK), X * GAMMAP(0:MAXBCK-1), P(0:MAXN), P1(0:MAXN), PWK(0:2*MAXN+1), X * ADIG(IR), XBEST(IR), ACT(IR), CB(IR) XC XC ... SPECIFICATIONS FOR EXTERNAL FUNCTIONS XC X DOUBLE PRECISION SUNORM XC XC ... EXECUTABLE STATEMENTS XC XC XC ... CHOICE OF IFLAG XC XC IFLAG .EQ. 0 Y IS SET TO r BY THE SUBROUTINE XC 0 XC IFLAG .NE. 0 THE USER MUST DEFINE Y IN XC THE MAIN PROGRAM XC X IFLAG = 1 XC X WRITE (*,*) ' OUTPUT RESULTS FROM EXA5' X WRITE (*,*) X WRITE (*,*) ' DIMENSION OF THE SYSTEM = ',N X WRITE (*,*) ' NUMBER OF CALLS = ',NBC X WRITE (*,*) ' MAXBCK (MAX JUMP MK) = ',MAXBCK X WRITE (*,*) ' MAXN (MAX NK+MK) = ',MAXN X WRITE (*,*) ' IFLAG = ',IFLAG X WRITE (*,*) ' EPS (NEAR-BREAKDOWN) = ',EPS X WRITE (*,*) ' EPS1 (GAUSS) = ',EPS1 X WRITE (*,*) ' EPS2 (SOLUTION) = ',EPS2 X WRITE (*,*) X WRITE (*,*) ' ITER. NK NORM' XC X DO 20 I = 1, N X X(I) = 0.0D0 X Y(I) = -(-1.0D0)**I X DO 10 J = 1, N X A(I,J) = 0.0D0 X 10 CONTINUE X 20 CONTINUE X Y(1) = 0.0D0 X Y(2) = 0.0D0 X Y(3) = 0.0D0 XC XC ... SET THE LAST COMPONENT OF Y XC EQUAL TO 1.0D-8 OR 0.0D0 XC X Y(N) = 1.0D-8 XC Y(N) = 0.0D0 XC X DO 30 I = 2, N X A(I-1,I) = 1.0D0 X IF (I .NE. 2) A(I,I-2) = 1.0D0 X A(I,I) = 2.0D0 X 30 CONTINUE X A(1,1) = 2.0D0 X AB(1) = 3.0D0 X AB(2) = AB(1) X AB(N) = AB(1) X CB(1) = AB(1) X CB(2) = AB(2) X CB(N) = AB(N) X DO 40 I = 3, N-1 X AB(I) = 4.0D0 X CB(I) = AB(I) X 40 CONTINUE XC X INIT = 0 XC X DO 60 K = 0, NBC X CALL BSMRZS (N, A, AB, Y, X, R, NORM, MAXBCK, MAXN, V, W, X * U, LMAT, RMAT, ETA, ETAP, GAMMA, GAMMAP, D, C, P, P1, X * PWK, IR, NK, EPS, EPS1, EPS2, INIT, IFLAG, IER) X IF (K .EQ. 0 .OR. NORM .LT. RBEST) THEN X RBEST = NORM X KSAV = K X NKSAV = NK X DO 50 I = 1, N X XBEST(I) = X(I) X 50 CONTINUE X END IF X WRITE (*,100) K, NK, NORM X IF (IER .EQ. 400) THEN X WRITE (*,*) X WRITE (*,*) ' SOLUTION NOT OBTAINED AFTER REACHING' X WRITE (*,*) ' THE DIMENSION N OF THE SYSTEM. THE' X WRITE (*,*) ' COMPUTATION CONTINUES UNTIL MAXN AT MOST' X WRITE (*,*) X IER = 0 X ELSE IF (IER .GE. 700) THEN X WRITE (*,*) X WRITE (*,*) ' JUMP FOR BREAKDOWN OR NEAR-BREAKDOWN ',IER X WRITE (*,*) X IER = 0 X END IF X IF (IER .NE. 0) THEN X WRITE (*,*) X IF (IER .EQ. 200) THEN X WRITE (*,*) ' THE SOLUTION HAS BEEN OBTAINED' X ELSE IF (IER .EQ. 300) THEN X WRITE (*,*) ' SOLUTION NOT OBTAINED AFTER' X WRITE (*,*) ' REACHING THE VALUE OF MAXN' X ELSE IF (IER .EQ. 600) THEN X WRITE (*,*) ' THE VALUE OF THE JUMP MK' X WRITE (*,*) ' IS GREATER THAN MAXBCK' X ELSE X WRITE (*,*) ' STOP DUE TO ERROR ', IER X STOP X END IF X GO TO 70 X END IF X 60 CONTINUE X WRITE (*,*) X WRITE (*,300) ' SOLUTION NOT OBTAINED AFTER', NBC X WRITE (*,*) ' ITERATIONS OF THE SUBROUTINE' XC X 70 CONTINUE X WRITE (*,*) X DO 80 I= 1, N X TMP = DABS(X(I)-1.0D0) X IF (TMP .EQ. 0.0D0) THEN X ADIG(I) = 999.00 X ELSE X ADIG(I) = -DLOG10(TMP/1.0D0) X END IF X 80 CONTINUE XC XC ... COMPUTE THE ACTUAL RESIDUAL XC X CALL MATVEC (N, A, X, 1, ACT, 1, IR, 0) X DO 90 I = 1, N X ACT(I) = ACT(I)-CB(I) X 90 CONTINUE X WRITE(*,500) ' ACTUAL RESIDUAL',SUNORM(N,ACT) X WRITE (*,*) XC X IF (IER .EQ. 200) THEN X WRITE (*,*) X * ' FINAL SOLUTION N. OF SIGN. DIG. ' X WRITE (*,200) (X(I),ADIG(I),I=1,N) X ELSE X WRITE (*,*) X * ' LAST SOLUTION N. OF SIGN. DIG. ' X WRITE (*,200) (X(I),ADIG(I),I=1,N) X WRITE (*,*) X WRITE (*,*) ' SOLUTION OF MINIMAL RESIDUAL OBTAINED AT' X WRITE (*,*) ' ITER. NK NORM' X WRITE (*,100) KSAV, NKSAV, RBEST X WRITE (*,*) X WRITE (*,400) (XBEST(I),I=1,N) X END IF X STOP X 100 FORMAT (I5,I8,D28.15) X 200 FORMAT (D25.15,F15.2) X 300 FORMAT (T2,A,I4) X 400 FORMAT (D25.15) X 500 FORMAT (T2,A,D24.15) X END SHAR_EOF if test 9016 -ne "`wc -c exa5.f`" then echo shar: error transmitting exa5.f '(should have been 9016 characters)' fi echo shar: extracting fh.f sed 's/^X//' << \SHAR_EOF > fh.f XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC FUNCTION NAME - FH + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC DESCRIPTION: + XC + XC DOUBLE PRECISION FUNCTION: + XC COMPUTES THE COEFFICIENT f (OR h ) (NEEDED IN THE SYSTEM TO + XC i,j i,j + XC BE SOLVED FOR THE COMPUTATION OF THE ALPHA , ALPHA', BETA , BETA') + XC i i i i + XC AS A LINEAR COMBINATION OF THE d (OR c ) AND OF THE COEFFICIENTS + XC (1) i i + XC OF THE POLYNOMIAL P . + XC k + XC + XC USAGE: + XC FH (I, J, CD, P1, NK, IFLAG) + XC + XC ARGUMENTS: + XC + XC I - INPUT INTEGER STATING THAT D(I+1) (OR C(I)) IS THE FIRST + XC VALUE TO BE CHOSEN IN THE LINEAR COMBINATION. + XC + XC J - INPUT INTEGER GIVING THE NUMBER OF TERMS MINUS ONE + XC IN THE LINEAR COMBINATION. + XC + XC CD - INPUT REAL VECTOR CONTAINING THE d 'S (OR c 'S). + XC i i + XC + XC P1 - INPUT REAL VECTOR CONTAINING THE COEFFICIENTS OF THE + XC (1) + XC POLYNOMIAL P OF DEGREE NK. + XC k + XC + XC NK - INPUT INTEGER EQUAL TO THE DIMENSION OF THE PREVIOUS + XC MATRIX OF THE SYSTEM. + XC + XC IFLAG - INPUT INTEGER FLAG. + XC IF IFLAG = 0, h IS COMPUTED + XC i,j + XC IF IFLAG = 1, f IS COMPUTED + XC i,j + XC + XC + XC REMARKS: + XC + XC - ALL THE REAL ARGUMENTS PASSED FROM THE CALLING PROGRAM MUST BE + XC IN DOUBLE PRECISION. + XC - IN THE CALLING PROCEDURE THE DIMENSIONS FOR THE VECTOR ARGUMENTS + XC ARE: + XC CD(0:*) + XC P1(0:*) + XC + XC AUTHORS: + XC + XC C. BREZINSKI + XC UNIVERSITY OF LILLE - FRANCE + XC M. REDIVO ZAGLIA + XC UNIVERSITY OF PADOVA - ITALY + XC + XC REFERENCE: + XC + XC - TREATMENT OF NEAR-BREAKDOWN IN THE CGS ALGORITHM + XC NUMERICAL ALGORITHMS 7 (1994) 33-73 + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC X DOUBLE PRECISION FUNCTION FH (I, J, CD, P1, NK, IFLAG) XC XC ... SPECIFICATIONS FOR ARGUMENTS XC X INTEGER I, IFLAG, J, NK X DOUBLE PRECISION CD(0:1), P1(0:1) XC XC ... SPECIFICATIONS FOR LOCAL VARIABLES XC X INTEGER L X DOUBLE PRECISION ACC XC XC ... EXECUTABLE STATEMENTS XC X ACC = 0.0D0 X IF (J .NE. 0) THEN X DO 10 L = 0, J-1 X IF (NK .GE. J-L) ACC = ACC + CD(I+L+IFLAG) * P1(NK-J+L) X 10 CONTINUE X END IF XC XC ... SET THE RESULT VALUE XC X FH = ACC + CD(I+J+IFLAG) X RETURN X END SHAR_EOF if test 6229 -ne "`wc -c fh.f`" then echo shar: error transmitting fh.f '(should have been 6229 characters)' fi echo shar: extracting gameta.f sed 's/^X//' << \SHAR_EOF > gameta.f XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC FUNCTION NAME - GAMETA + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC DESCRIPTION: + XC + XC DOUBLE PRECISION FUNCTION: + XC COMPUTES THE VALUE OF THE COEFFICIENTS GAMMA OR GAMMA' OR ETA OR + XC j j j + XC ETA' (OF THE POLYNOMIALS w OR v OR q OR t ) AS A LINEAR + XC j k k k k + XC COMBINATION OF THE BETA OR OF THE BETA' OR OF THE ALPHA OR OF + XC i i i + XC (1) + XC THE ALPHA' , WITH THE COEFFICIENTS OF THE POLYNOMIAL P . + XC i k + XC + XC USAGE: + XC GAMETA (I, J, ALBE, P1, NK) + XC + XC ARGUMENTS: + XC + XC I - INPUT INTEGER NEEDED TO DEFINE THE NUMBER OF TERMS IN THE + XC LINEAR COMBINATION. + XC + XC J - INPUT INTEGER DEFINING THE INDEX OF THE COEFFICIENT TO BE + XC COMPUTED. + XC + XC ALBE - INPUT REAL VECTOR CONTAINING THE VALUES OF ALPHA OR OF + XC i + XC ALPHA' OR OF BETA OR OF BETA'. + XC i i i + XC + XC P1 - INPUT REAL VECTOR CONTAINING THE COEFFICIENTS OF THE + XC (1) + XC POLYNOMIAL P OF DEGREE NK. + XC k + XC + XC NK - INPUT INTEGER EQUAL TO THE DIMENSION OF THE PREVIOUS + XC MATRIX OF THE SYSTEM. + XC + XC REMARKS: + XC + XC - ALL THE REAL ARGUMENTS PASSED FROM THE CALLING PROGRAM MUST BE + XC IN DOUBLE PRECISION. + XC - IN THE CALLING PROCEDURE THE DIMENSIONS FOR THE VECTOR ARGUMENTS + XC ARE: + XC ALBE(0:*) + XC P1(0:*) + XC + XC AUTHORS: + XC + XC C. BREZINSKI + XC UNIVERSITY OF LILLE - FRANCE + XC M. REDIVO ZAGLIA + XC UNIVERSITY OF PADOVA - ITALY + XC + XC REFERENCE: + XC + XC - TREATMENT OF NEAR-BREAKDOWN IN THE CGS ALGORITHM + XC NUMERICAL ALGORITHMS 7 (1994) 33-73 + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC X DOUBLE PRECISION FUNCTION GAMETA (I, J, ALBE, P1, NK) XC XC ... SPECIFICATIONS FOR ARGUMENTS XC X INTEGER I, J, NK X DOUBLE PRECISION ALBE(0:1), P1(0:1) XC XC ... SPECIFICATIONS FOR LOCAL VARIABLES XC X INTEGER L X DOUBLE PRECISION ACC XC XC ... EXECUTABLE STATEMENTS XC XC ... CASE I = J XC X IF (I .EQ. J) THEN X GAMETA = ALBE(I) X RETURN X END IF XC XC ... CASE I <> J XC X ACC = 0.0D0 X DO 10 L = 1, I-J X IF (NK .GE. L) ACC = ACC + ALBE(L+J) * P1(NK-L) X 10 CONTINUE XC XC ... SET THE RESULT VALUE XC X GAMETA = ACC + ALBE(J) X RETURN X END SHAR_EOF if test 6087 -ne "`wc -c gameta.f`" then echo shar: error transmitting gameta.f '(should have been 6087 characters)' fi echo shar: extracting gsolvd.f sed 's/^X//' << \SHAR_EOF > gsolvd.f XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC SUBPROGRAM NAME - GSOLVD + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC DESCRIPTION: + XC + XC COMPUTES THE SOLUTION OF A SYSTEM OF LINEAR EQUATIONS WITH A + XC AS COEFFICIENT MATRIX. THE MATRIX B CONTAINS THE 2 RIGHT HAND + XC SIDES. THE METHOD IS GAUSSIAN ELIMINATION. + XC IF THE MATRIX A IS NON-SINGULAR THE SOLUTIONS ARE SET IN THE + XC ROWS OF THE MATRIX B. + XC + XC USAGE: + XC CALL GSOLVD (A, B, N, EPS, INDS) + XC + XC ARGUMENTS: + XC + XC A - INPUT/WORKING REAL VECTOR OF DIMENSION (N*N) CONTAINING + XC THE COEFFICIENT MATRIX. + XC + XC B - INPUT/OUTPUT REAL ARRAY OF DIMENSION (2,N). + XC + XC N - INPUT INTEGER INDICATING THE DIMENSION OF THE SYSTEM. + XC + XC EPS - INPUT REAL VALUE USED IN TESTS FOR ZERO. + XC IF DABS(X) .LE. EPS, THEN X IS CONSIDERED TO BE ZERO. + XC + XC INDS - OUTPUT INTEGER VALUE. + XC IF INDS = 0 THEN THE MATRIX A IS NON-SINGULAR. + XC IF INDS = 1 THEN THE MATRIX A IS SINGULAR. + XC + XC REMARKS: + XC + XC - ALL THE REAL ARGUMENTS PASSED FROM THE CALLING PROGRAM MUST BE + XC IN DOUBLE PRECISION. + XC - IN THE CALLING PROCEDURE THE MINIMAL DIMENSIONS FOR THE ARGUMENTS + XC ARE: + XC A(N*N) + XC B(2,N) + XC + XC AUTHORS: + XC + XC C. BREZINSKI + XC UNIVERSITY OF LILLE - FRANCE + XC M. REDIVO ZAGLIA + XC UNIVERSITY OF PADOVA - ITALY + XC + XC REFERENCES: + XC + XC - TREATMENT OF NEAR-BREAKDOWN IN THE CGS ALGORITHMS + XC NUMERICAL ALGORITHMS 7 (1994) 33-73 + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC X SUBROUTINE GSOLVD (A, B, N, EPS, INDS) XC XC ... SPECIFICATIONS FOR ARGUMENTS XC X INTEGER INDS, N X DOUBLE PRECISION EPS X DOUBLE PRECISION A(1), B(2,1) XC XC ... SPECIFICATIONS FOR LOCAL VARIABLES XC X INTEGER I, IPIVOT, J, K, NN, NN1 X DOUBLE PRECISION TMP, PIVOT, DET XC XC ... EXECUTABLE STATEMENTS XC X INDS = 0 XC X IF (N .NE. 1) THEN XC XC ... TRIANGULARIZATION XC X DO 70 K = 1, N-1 XC XC ... SEARCH FOR THE PIVOT (PARTIAL PIVOTING) XC X PIVOT = 0.0D0 X DET = 1.0D0 X NN = N * (K - 1) X DO 10 I = K, N X TMP = DABS(A(I+NN)) X IF (TMP .GT. PIVOT) THEN X IPIVOT = I X PIVOT = TMP X END IF X 10 CONTINUE X DET = DET * PIVOT XC XC ... TEST FOR SINGULAR MATRIX XC X IF (DABS(PIVOT) .LE. EPS) THEN X INDS = 1 X RETURN X END IF XC XC ... EXCHANGE THE ROW K WITH THE ROW XC CONTAINING THE MAXIMUM PIVOT XC X DO 20 J = K, N X NN1 = N * (J - 1) X TMP = A(K+NN1) X A(K+NN1) = A(IPIVOT+NN1) X A(IPIVOT+NN1) = TMP X 20 CONTINUE X DO 30 J = 1, 2 X TMP = B(J,K) X B(J,K) = B(J,IPIVOT) X B(J,IPIVOT) = TMP X 30 CONTINUE XC XC ... APPLY THE ELEMENTS TRANSFORMATION XC X DO 60 I = K+1, N X IF (A(I+NN) .NE. 0.0D0) THEN X TMP = A(I+NN) / A(K+NN) X DO 40 J = K+1, N X NN1 = N * (J - 1) X A(I+NN1) = A(I+NN1) - TMP * A(K+NN1) X 40 CONTINUE X DO 50 J = 1, 2 X B(J,I) = B(J,I) - TMP * B(J,K) X 50 CONTINUE X END IF X 60 CONTINUE X 70 CONTINUE X END IF XC XC ... TEST FOR SINGULAR MATRIX XC (PIVOT OF THE LAST ROW) XC X DET = DET * A(N*N) X IF (DABS(DET) .LE. EPS) THEN X INDS = 1 X RETURN X END IF XC XC ... SOLVING THE TRIANGULAR SYSTEM XC FOR THE 2 RIGHT HAND SIDES XC X DO 100 K = 1, 2 X B(K,N) = B(K,N) / A(N*N) X IF (N .NE. 1) THEN X DO 90 I = N-1, 1, -1 X DO 80 J = I+1, N X B(K,I) = B(K,I) - A(I+N*(J-1)) * B(K,J) X 80 CONTINUE X B(K,I) = B(K,I) / A(I+N*(I-1)) X 90 CONTINUE X END IF X 100 CONTINUE X RETURN X END SHAR_EOF if test 7323 -ne "`wc -c gsolvd.f`" then echo shar: error transmitting gsolvd.f '(should have been 7323 characters)' fi echo shar: extracting matvec.f sed 's/^X//' << \SHAR_EOF > matvec.f XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC SUBROUTINE NAME - MATVEC + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC DESCRIPTION: + XC + XC COMPUTES THE PRODUCT BETWEEN A MATRIX (OR ITS TRANSPOSED) AND A + XC VECTOR. + XC THE VECTOR IS STORED AS A COLUMN OF THE ARRAY B. + XC THE RESULTING VECTOR IS PUT IN A PRESCRIBED COLUMN OF THE ARRAY C. + XC + XC USAGE: + XC CALL MATVEC (IP, A, B, IB, C, IC, IR, IFLAG) + XC + XC ARGUMENTS: + XC + XC + XC IP - NUMBER OF ELEMENTS TO BE CONSIDERED IN EACH VECTOR + XC STARTING FROM THE FIRST ONE. + XC + XC A - INPUT REAL MATRIX OF DIMENSION (IP,IP). + XC + XC B - INPUT REAL MATRIX CONTAINING A CERTAIN NUMBER OF COLUMN + XC VECTORS. ITS FIRST DIMENSION MUST BE IP. + XC + XC IB - INPUT COLUMN OF THE MATRIX B. SECOND VECTOR OF THE INNER + XC PRODUCT. + XC + XC C - OUTPUT REAL MATRIX CONTAINING A CERTAIN NUMBER OF COLUMN + XC VECTORS. ITS FIRST DIMENSION MUST BE IP. + XC + XC IC - COLUMN OF THE MATRIX C IN WHICH THE RESULTING VECTOR WILL + XC BE STORED. + XC + XC IR - INPUT ROWS DIMENSION EXACTLY AS SPECIFIED IN THE + XC DIMENSION STATEMENTS IN THE CALLING PROGRAM FOR THE + XC MATRIX B. + XC + XC IFLAG - IF IFLAG .EQ. 0 THE MATRIX A IS CONSIDERED + XC IF IFLAG .NE. 0 THE TRANSPOSED MATRIX A IS CONSIDERED + XC + XC REMARKS: + XC + XC - ALL THE REAL ARGUMENTS PASSED FROM THE CALLING PROGRAM MUST BE + XC IN DOUBLE PRECISION. + XC - IN THE CALLING PROCEDURE, IF THE MATRIX B (OR C) HAS ITS SECOND + XC DIMENSION WITH A LOWER BOUND NOT EQUAL TO 1, THEN THE VALUE IB = 1 + XC (OR IC = 1) CORRESPONDS TO ITS FIRST COLUMN AND NOT TO THE COLUMN + XC WHOSE INDEX IS 1. + XC - IN THE CALLING PROCEDURE THE MINIMAL DIMENSIONS FOR THE + XC MATRIX AND VECTOR ARGUMENTS (IR=IP) ARE: + XC A(IP,IP) + XC B(IP,*) + XC C(IP,*) + XC + XC AUTHORS: + XC + XC C. BREZINSKI AND H. SADOK + XC UNIVERSITY OF LILLE - FRANCE + XC M. REDIVO ZAGLIA + XC UNIVERSITY OF PADOVA - ITALY + XC + XC REFERENCES: + XC + XC - AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE ALGORITHMS, + XC NUMERICAL ALGORITHMS, 1(1991), PP. 261-284. + XC + XC - ADDENDUM TO "AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE + XC ALGORITHMS", NUMERICAL ALGORITHMS, 2(1992), PP. 133-136. + XC + XC - TREATMENT OF NEAR-BREAKDOWN IN THE CGS ALGORITHM, + XC NUMERICAL ALGORITHMS 7 (1994) 33-73 + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC X SUBROUTINE MATVEC (IP, A, B, IB, C, IC, IR, IFLAG) XC XC ... SPECIFICATIONS FOR ARGUMENTS XC X INTEGER IB, IC, IFLAG, IP, IR X DOUBLE PRECISION A(IR,1), B(1), C(1) XC XC ... SPECIFICATIONS FOR LOCAL VARIABLES XC X DOUBLE PRECISION ACC X INTEGER I, J, IBSTAR, ICSTAR XC XC ... EXECUTABLE STATEMENTS XC X IBSTAR = (IB-1)*IR X ICSTAR = (IC-1)*IR X DO 20 J = 1, IP X ACC = 0.0D0 X DO 10 I = 1, IP X IF (IFLAG .EQ. 0) THEN X ACC = ACC + A(J,I) * B(IBSTAR+I) X ELSE X ACC = ACC + A(I,J) * B(IBSTAR+I) X END IF X 10 CONTINUE XC XC ... SET THE RESULT VALUE XC X C(ICSTAR+J) = ACC X 20 CONTINUE X RETURN X END SHAR_EOF if test 6964 -ne "`wc -c matvec.f`" then echo shar: error transmitting matvec.f '(should have been 6964 characters)' fi echo shar: extracting pinner.f sed 's/^X//' << \SHAR_EOF > pinner.f XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC FUNCTION NAME - PINNER + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC DESCRIPTION: + XC + XC DOUBLE PRECISION FUNCTION: + XC COMPUTES THE INNER PRODUCT OF TWO VECTORS. + XC THE DIMENSION OF THE FIRST VECTOR IS IR. THE SECOND VECTOR + XC IS STORED AS A COLUMN OF THE ARRAY B WHOSE FIRST DIMENSION IS IR. + XC + XC USAGE: + XC PINNER (IP, A, B, IB, IR) + XC + XC ARGUMENTS: + XC + XC IP - NUMBER OF ELEMENTS TO BE CONSIDERED IN EACH VECTOR + XC STARTING FROM THE FIRST ONE. + XC + XC A - INPUT REAL VECTOR (FIRST VECTOR OF THE INNER PRODUCT). + XC + XC B - INPUT REAL MATRIX CONTAINING A CERTAIN NUMBER OF COLUMN + XC VECTORS. + XC + XC IB - INPUT INTEGER DEFINING THE NUMBER OF THE COLUMN OF THE + XC MATRIX B WHICH IS THE SECOND VECTOR OF THE INNER PRODUCT. + XC + XC IR - INPUT ROWS DIMENSION EXACTLY AS SPECIFIED IN THE + XC DIMENSION STATEMENTS IN THE CALLING PROGRAM FOR THE + XC MATRIX B. + XC + XC REMARKS: + XC + XC - ALL THE REAL ARGUMENTS PASSED FROM THE CALLING PROGRAM MUST BE + XC IN DOUBLE PRECISION. + XC - IN THE CALLING PROCEDURE THE MINIMAL DIMENSIONS FOR THE + XC MATRIX AND VECTOR ARGUMENTS (IR=IP) ARE: + XC A(IP) + XC B(IP,*) + XC - IN THE CALLING PROCEDURE, IF THE MATRIX B HAS ITS SECOND DIMENSION + XC WITH A LOWER BOUND NOT EQUAL TO 1, THEN THE VALUE IB = 1 + XC CORRESPONDS TO ITS FIRST COLUMN AND NOT TO THE COLUMN WHOSE INDEX + XC IS 1. + XC + XC AUTHORS: + XC + XC C. BREZINSKI AND H. SADOK + XC UNIVERSITY OF LILLE - FRANCE + XC M. REDIVO ZAGLIA + XC UNIVERSITY OF PADOVA - ITALY + XC + XC REFERENCES: + XC + XC - AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE ALGORITHMS, + XC NUMERICAL ALGORITHMS, 1(1991), PP. 261-284. + XC + XC - ADDENDUM TO "AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE + XC ALGORITHMS", NUMERICAL ALGORITHMS, 2(1992), PP. 133-136. + XC + XC - TREATMENT OF NEAR-BREAKDOWN IN THE CGS ALGORITHM, + XC NUMERICAL ALGORITHMS 7 (1994) 33-73 + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC X DOUBLE PRECISION FUNCTION PINNER (IP, A, B, IB, IR) XC XC ... SPECIFICATIONS FOR ARGUMENTS XC X DOUBLE PRECISION A(1), B(1) X INTEGER IB, IP, IR XC XC ... SPECIFICATIONS FOR LOCAL VARIABLES XC X DOUBLE PRECISION ACC X INTEGER I, IBSTAR XC XC ... EXECUTABLE STATEMENTS XC X IBSTAR = (IB-1)*IR X ACC = 0.0D0 X DO 10 I = 1, IP X ACC = ACC + A(I) * B(IBSTAR+I) X 10 CONTINUE XC XC ... SET THE RESULT VALUE XC X PINNER = ACC X RETURN X END SHAR_EOF if test 5838 -ne "`wc -c pinner.f`" then echo shar: error transmitting pinner.f '(should have been 5838 characters)' fi echo shar: extracting polmul.f sed 's/^X//' << \SHAR_EOF > polmul.f XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC SUBROUTINE NAME - POLMUL + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC DESCRIPTION: + XC + XC COMPUTES THE PRODUCT BETWEEN TWO POLYNOMIALS. + XC + XC USAGE: + XC CALL POLMUL (P1, IDP1, P2, IDP2, PRES, IER) + XC + XC ARGUMENTS: + XC + XC P1 - INPUT REAL VECTOR OF DIMENSION (0:IDP1) CONTAINING THE + XC COEFFICIENTS OF THE FIRST POLYNOMIAL, ORDERED FROM + XC SMALLEST TO LARGEST POWER. + XC + XC IDP1 - INPUT INTEGER VALUE. IT REPRESENTS THE DEGREE OF THE + XC FIRST POLYNOMIAL. + XC + XC P2 - INPUT REAL VECTOR OF DIMENSION (0:IDP2) CONTAINING THE + XC COEFFICIENTS OF THE SECOND POLYNOMIAL, ORDERED FROM + XC SMALLEST TO LARGEST POWER. + XC + XC IDP2 - INPUT INTEGER VALUE. IT REPRESENTS THE DEGREE OF THE + XC SECOND POLYNOMIAL. + XC + XC PRES - OUTPUT REAL VECTOR OF DIMENSION (0:IDP1+IDP2) CONTAINING + XC THE COEFFICIENTS OF THE RESULTANT POLYNOMIAL, ORDERED + XC FROM SMALLEST TO LARGEST POWER. + XC + XC IER - OUTPUT INDEX WARNING/ERROR. + XC IER = 2000 THE DEGREE IDP1 IS GREATER THAN THE DEGREE + XC IDP2. + XC + XC REMARKS: + XC + XC - ALL THE REAL ARGUMENTS PASSED FROM THE CALLING PROGRAM MUST BE + XC IN DOUBLE PRECISION. + XC - THE DEGREE OF THE FIRST POLYNOMIAL MUST BE LESS OR EQUAL THAN THE + XC DEGREE OF THE SECOND POLYNOMIAL. + XC - IN THE CALLING PROCEDURE THE DIMENSIONS FOR THE VECTOR ARGUMENTS + XC ARE: + XC P1(0:*) + XC P2(0:*) + XC PRES(0:*) + XC + XC AUTHORS: + XC + XC C. BREZINSKI + XC UNIVERSITY OF LILLE - FRANCE + XC M. REDIVO ZAGLIA + XC UNIVERSITY OF PADOVA - ITALY + XC + XC REFERENCE: + XC + XC - TREATMENT OF NEAR-BREAKDOWN IN THE CGS ALGORITHM + XC NUMERICAL ALGORITHMS 7 (1994) 33-73 + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC X SUBROUTINE POLMUL (P1, IDP1, P2, IDP2, PRES, IER) XC XC ... SPECIFICATIONS FOR ARGUMENTS XC X INTEGER IDP1, IDP2, IER X DOUBLE PRECISION P1(0:1), P2(0:1), PRES(0:1) XC XC ... SPECIFICATIONS FOR LOCAL VARIABLES XC X INTEGER I, I1, I2, IDPRES, J XC XC ... EXECUTABLE STATEMENTS XC X IF (IDP1 .GT. IDP2) THEN X IER = 2000 X RETURN X END IF XC X IDPRES = IDP1 + IDP2 XC XC ... COMPUTES THE COEFFICIENTS OF PRES XC X DO 20 I = 0, IDPRES X I1 = MAX0(0,I-IDP1) X I2 = MIN0(IDP2,I) X PRES(I) = 0.0D0 X DO 10 J = I1, I2 X PRES(I) = PRES(I) + P1(I-J) * P2(J) X 10 CONTINUE X 20 CONTINUE X RETURN X END SHAR_EOF if test 5795 -ne "`wc -c polmul.f`" then echo shar: error transmitting polmul.f '(should have been 5795 characters)' fi echo shar: extracting rxsumm.f sed 's/^X//' << \SHAR_EOF > rxsumm.f XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC FUNCTION NAME - RXSUMM + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC DESCRIPTION: + XC + XC DOUBLE PRECISION FUNCTION: + XC COMPUTES THE COMPONENT NR OF A VECTOR OBTAINED AS A LINEAR + XC COMBINATION, WITH THE COEFFICIENTS STORED IN THE VECTORS P1, P2 + XC AND P3, OF THE COLUMN VECTORS STORED IN V, W AND IN U. + XC + XC USAGE: + XC RXSUMM (P1, P2, P3, V, W, U, IR, NR, NK, MK, IFLAG) + XC + XC ARGUMENTS: + XC + XC P1 - INPUT REAL VECTOR WHICH CONTAINS THE COEFFICIENTS OF THE + XC POLYNOMIAL v(xi)*(xi*v(xi) - 2) (FOR THE COMPUTATION OF + XC THE NEW SOLUTION x) OR OF THE POLYNOMIAL (1-xi*v(xi))**2 + XC (FOR THE COMPUTATION OF THE NEW RESIDUAL r). + XC + XC P2 - INPUT REAL VECTOR WHICH CONTAINS THE COEFFICIENTS OF THE + XC POLYNOMIAL (1 - xi*v(xi)) * w(xi). + XC + XC P3 - INPUT REAL VECTOR WHICH CONTAINS THE COEFFICIENTS OF THE + XC POLYNOMIAL w(xi)**2. + XC + XC V - INPUT REAL MATRIX CONTAINING A CERTAIN NUMBER OF COLUMN + XC VECTORS. + XC + XC W - INPUT REAL MATRIX CONTAINING A CERTAIN NUMBER OF COLUMN + XC VECTORS. + XC + XC U - INPUT REAL MATRIX CONTAINING A CERTAIN NUMBER OF COLUMN + XC VECTORS. + XC + XC IR - INPUT ROWS DIMENSION EXACTLY AS SPECIFIED IN THE + XC DIMENSION STATEMENTS IN THE CALLING PROGRAM FOR THE + XC MATRICES V, W AND U. + XC + XC NR - INPUT INTEGER EQUAL TO THE COMPONENT WHICH HAS TO BE + XC USED AND COMPUTED. + XC + XC NK - INPUT INTEGER EQUAL TO THE DIMENSION OF THE PREVIOUS + XC MATRIX OF THE SYSTEM. + XC + XC MK - INPUT INTEGER EQUAL TO THE DIMENSION OF THE NEW JUMP. + XC + XC IFLAG - INPUT INTEGER FLAG. + XC IF IFLAG = 0 THE NEW SOLUTION x WILL BE COMPUTED. + XC IF IFLAG = 1 THE NEW RESIDUAL r WILL BE COMPUTED. + XC + XC REMARKS: + XC + XC - ALL THE REAL ARGUMENTS PASSED FROM THE CALLING PROGRAM MUST BE + XC IN DOUBLE PRECISION. + XC - IN THE CALLING PROCEDURE THE DIMENSIONS FOR THE ARRAY ARGUMENTS + XC ARE: + XC P1(0:*) + XC P2(0:*) + XC P3(0:*) + XC V(IR,0:*) + XC W(IR,0:*) + XC U(IR,0:*) + XC + XC AUTHORS: + XC + XC C. BREZINSKI + XC UNIVERSITY OF LILLE - FRANCE + XC M. REDIVO ZAGLIA + XC UNIVERSITY OF PADOVA - ITALY + XC + XC REFERENCE: + XC + XC - TREATMENT OF NEAR-BREAKDOWN IN THE CGS ALGORITHM + XC NUMERICAL ALGORITHMS 7 (1994) 33-73 + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC X DOUBLE PRECISION FUNCTION RXSUMM (P1, P2, P3, V, W, U, IR, NR, X * NK, MK, IFLAG) XC XC ... SPECIFICATIONS FOR ARGUMENTS XC X INTEGER IR, NK, NR, MK, IFLAG X DOUBLE PRECISION P1(0:1), P2(0:1), P3(0:1), V(IR,0:1), W(IR,0:1), X * U(IR,0:1) XC XC ... SPECIFICATIONS FOR LOCAL VARIABLES XC X DOUBLE PRECISION ACC X INTEGER I, IFL XC XC ... EXECUTABLE STATEMENTS XC XC XC ... CASE MK = 1 XC X IF (MK .EQ. 1) THEN X IF (IFLAG .EQ. 0) THEN X ACC = - 2.0D0 * P2(0) * U(NR,0) + P3(0) * V(NR,1) X ELSE X ACC = W(NR,0) - 2.0D0 * P2(0) * U(NR,1) + P3(0) * V(NR,2) X END IF X RXSUMM = ACC X RETURN X END IF XC X ACC = 0.0D0 XC XC ... CASE MK <= NK XC X IF (MK .LE. NK) THEN X DO 10 I = 0, 2*MK - 2 X IFL = I + IFLAG X ACC = ACC - 2.0D0*P2(I)*U(NR,IFL) + P3(I)*V(NR,IFL+1) X 10 CONTINUE X DO 20 I = 0, 2*MK - 3 + IFLAG X ACC = ACC + P1(I) * W(NR,I) X 20 CONTINUE XC XC ... CASE MK > NK XC X ELSE X DO 30 I = 0, 2*MK - 2 X ACC = ACC + P3(I) * V(NR,I+1+IFLAG) X 30 CONTINUE X DO 40 I = 0, NK + MK - 1 X ACC = ACC - 2.0D0 * P2(I) * U(NR,I+IFLAG) X 40 CONTINUE X IF ((IFLAG .EQ. 0 .AND. NK .NE. 0) .OR. IFLAG .EQ. 1) THEN X DO 50 I = 0, 2*NK - 1 + IFLAG X ACC = ACC + P1(I) * W(NR,I) X 50 CONTINUE X END IF X END IF XC XC ... SET THE RESULT VALUE XC X RXSUMM = ACC X RETURN X END SHAR_EOF if test 8229 -ne "`wc -c rxsumm.f`" then echo shar: error transmitting rxsumm.f '(should have been 8229 characters)' fi echo shar: extracting solshf.f sed 's/^X//' << \SHAR_EOF > solshf.f XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC SUBROUTINE NAME - SOLSHF + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC DESCRIPTION: + XC + XC COMPUTES THE SOLUTIONS OF THE AUXILIARY SYSTEMS, THAT IS THE + XC ALPHA , ALPHA' , BETA AND BETA'. + XC i i i i + XC + XC USAGE: + XC + XC CALL SOLSHF (LMAT, RMAT, C, D, P1, NK, MK, EPS, INDS) + XC + XC ARGUMENTS: + XC + XC LMAT - WORKING REAL VECTOR TO STORE THE COEFFICIENTS OF THE + XC AUXILIARY SYSTEM. + XC + XC RMAT - OUTPUT REAL MATRIX OF DIMENSION (2,0:*). + XC IF THE SYSTEM IS NON SINGULAR, IT CONTAINS IN OUTPUT: + XC ROW 1 THE BETA AND THE BETA'. + XC i i + XC ROW 2 THE ALPHA AND THE ALPHA' . + XC i i + XC + XC C - INPUT REAL VECTOR CONTAINING THE c 'S. + XC i + XC + XC D - INPUT REAL VECTOR CONTAINING THE d 'S. + XC i + XC + XC P1 - INPUT REAL VECTOR CONTAINING THE COEFFICIENTS OF THE + XC (1) + XC POLYNOMIAL P OF DEGREE NK. + XC k + XC + XC NK - INPUT INTEGER EQUAL TO THE DIMENSION OF THE PREVIOUS + XC MATRIX OF THE SYSTEM. + XC + XC MK - INPUT INTEGER EQUAL TO THE LENGTH OF THE NEW JUMP. + XC + XC EPS - INPUT REAL VALUE USED FOR TESTING THE PIVOTS IN GAUSSIAN + XC ELIMINATION. + XC IF DABS(X) .LE. EPS, THEN X IS CONSIDERED TO BE ZERO. + XC THE SUBROUTINE DOES NOT CONTROL WHETHER OR NOT EPS IS + XC A NEGATIVE OR ZERO REAL NUMBER. + XC + XC INDS - OUTPUT INTEGER VALUE. + XC IF INDS = 0 THEN THE SYSTEM IS NON-SINGULAR. + XC IF INDS = 1 THEN THE SYSTEM IS SINGULAR. + XC + XC EXTERNAL MODULES: + XC + XC - DOUBLE PRECISION FUNCTION FH + XC - SUBROUTINES GSOLVD, STOMHF, STORHF + XC + XC REMARKS: + XC + XC - ALL THE REAL ARGUMENTS PASSED FROM THE CALLING PROGRAM MUST BE + XC IN DOUBLE PRECISION. + XC - IN THE CALLING PROCEDURE THE DIMENSIONS FOR THE ARGUMENTS ARE: + XC LMAT(*) + XC RMAT(2,0:*) + XC C(0:*) + XC D(0:*) + XC P1(0:*) + XC + XC AUTHORS: + XC + XC C. BREZINSKI + XC UNIVERSITY OF LILLE - FRANCE + XC M. REDIVO ZAGLIA + XC UNIVERSITY OF PADOVA - ITALY + XC + XC REFERENCES: + XC + XC - TREATMENT OF NEAR-BREAKDOWN IN THE CGS ALGORITHMS + XC NUMERICAL ALGORITHMS 7 (1994) 33-73 + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC X SUBROUTINE SOLSHF (LMAT, RMAT, C, D, P1, NK, MK, EPS, INDS) XC XC ... SPECIFICATIONS FOR ARGUMENTS XC X INTEGER INDS, MK, NK X DOUBLE PRECISION EPS X DOUBLE PRECISION LMAT(1), RMAT(2,0:1), C(0:1), D(0:1), P1(0:1) XC XC ... SPECIFICATIONS FOR EXTERNAL MODULES XC X DOUBLE PRECISION FH XC XC ... SPECIFICATIONS FOR LOCAL VARIABLES XC X INTEGER IFLAG, NN XC XC ... EXECUTABLE STATEMENTS XC XC XC ... CASE MK = 1 XC X IF (MK .EQ. 1) THEN X INDS = 0 XC XC ... COMPUTES BETA(0) XC X RMAT(1,0) = D(0)/C(0) XC XC ... COMPUTES ALPHA(0) (CAS MK > NK = 0) XC X RMAT(2,0) = -C(1)/C(0) XC X IF (MK .LE. NK) THEN XC XC ... COMPUTES ALPHA'(0) XC X RMAT(2,1) = -C(0)/D(0) XC XC ... COMPUTES ALPHA(0) + P1(NK-1) XC (CAS MK <= NK) XC X RMAT(2,0) = (-FH(0,1,C,P1,NK,0) - D(1) * RMAT(2,1)) X * / C(0) + P1(NK-1) X END IF X RETURN X END IF XC XC ... CHOICE FOR THE SYSTEM XC X IF (MK .LE. NK) THEN X NN = MK - 1 X IFLAG = 1 X ELSE X NN = NK X IFLAG = 2 X END IF XC XC ... COMPUTATION OF BETA , BETA' AND OF XC i i XC ALPHA , ALPHA' XC i i XC XC ... STORE THE COEFFICIENTS AND THE RIGHT HAND XC SIDES XC X CALL STOMHF (LMAT, C, D, P1, NK, MK, IFLAG) X CALL STORHF (RMAT, C, D, P1, NK, MK, IFLAG) XC XC ... SOLVE THE SYSTEM XC X CALL GSOLVD (LMAT, RMAT, NN+MK, EPS, INDS) XC XC ... IF NON SINGULAR SYSTEM, STORE THE XC ADDITIONAL SOLUTION XC X IF (INDS .EQ. 0 .AND. IFLAG .EQ. 1) RMAT(2,2*MK-1) = -C(0)/D(0) XC X RETURN X END SHAR_EOF if test 8787 -ne "`wc -c solshf.f`" then echo shar: error transmitting solshf.f '(should have been 8787 characters)' fi echo shar: extracting ssumm.f sed 's/^X//' << \SHAR_EOF > ssumm.f XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC FUNCTION NAME - SSUMM + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC DESCRIPTION: + XC + XC DOUBLE PRECISION FUNCTION: + XC COMPUTES THE COMPONENT NR OF A VECTOR OBTAINED AS A LINEAR + XC COMBINATION, WITH THE COEFFICIENTS STORED IN THE VECTORS P1, P2, + XC P3 AND P4, OF THE COLUMN VECTORS STORED IN V, W AND IN U. + XC + XC USAGE: + XC SSUMM (P1, P2, P3, P4, V, W, U, IR, NR, NK, MK) + XC + XC ARGUMENTS: + XC + XC P1 - INPUT REAL VECTOR WHICH CONTAINS THE COEFFICIENTS OF THE + XC POLYNOMIAL t(xi) * (1 - xi*v(xi)). + XC + XC P2 - INPUT REAL VECTOR WHICH CONTAINS THE COEFFICIENTS OF THE + XC POLYNOMIAL (1 - xi*v(xi)) * q(xi). + XC + XC P3 - INPUT REAL VECTOR WHICH CONTAINS THE COEFFICIENTS OF THE + XC POLYNOMIAL t(xi) * w(xi). + XC + XC P4 - INPUT REAL VECTOR WHICH CONTAINS THE COEFFICIENTS OF THE + XC POLYNOMIAL w(xi) * q(xi). + XC + XC V - INPUT REAL MATRIX CONTAINING A CERTAIN NUMBER OF COLUMN + XC VECTORS. + XC + XC W - INPUT REAL MATRIX CONTAINING A CERTAIN NUMBER OF COLUMN + XC VECTORS. + XC + XC U - INPUT REAL MATRIX CONTAINING A CERTAIN NUMBER OF COLUMN + XC VECTORS. + XC + XC IR - INPUT ROWS DIMENSION EXACTLY AS SPECIFIED IN THE + XC DIMENSION STATEMENTS IN THE CALLING PROGRAM FOR THE + XC MATRICES V, W AND U. + XC + XC NR - INPUT INTEGER EQUAL TO THE COMPONENT WHICH HAS TO BE + XC USED AND COMPUTED. + XC + XC NK - INPUT INTEGER EQUAL TO THE DIMENSION OF THE PREVIOUS + XC MATRIX OF THE SYSTEM. + XC + XC MK - INPUT INTEGER EQUAL TO THE DIMENSION OF THE NEW JUMP. + XC + XC REMARKS: + XC + XC - ALL THE REAL ARGUMENTS PASSED FROM THE CALLING PROGRAM MUST BE + XC IN DOUBLE PRECISION. + XC - IN THE CALLING PROCEDURE THE DIMENSIONS FOR THE ARRAY ARGUMENTS + XC ARE: + XC P1(0:*) + XC P2(0:*) + XC P3(0:*) + XC P4(0:*) + XC V(IR,0:*) + XC W(IR,0:*) + XC U(IR,0:*) + XC + XC AUTHORS: + XC + XC C. BREZINSKI + XC UNIVERSITY OF LILLE - FRANCE + XC M. REDIVO ZAGLIA + XC UNIVERSITY OF PADOVA - ITALY + XC + XC REFERENCE: + XC + XC - TREATMENT OF NEAR-BREAKDOWN IN THE CGS ALGORITHM + XC NUMERICAL ALGORITHMS 7 (1994) 33-73 + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC X DOUBLE PRECISION FUNCTION SSUMM (P1, P2, P3, P4, V, W, U, IR, X * NR, NK, MK) XC XC ... SPECIFICATIONS FOR ARGUMENTS XC X INTEGER IR, NK, NR, MK X DOUBLE PRECISION P1(0:1), P2(0:1), P3(0:1), P4(0:1), X * V(IR,0:1), W(IR,0:1), U(IR,0:1) XC XC ... SPECIFICATIONS FOR LOCAL VARIABLES XC X DOUBLE PRECISION ACC X INTEGER I XC XC ... EXECUTABLE STATEMENTS XC XC XC ... CASE MK = 1 XC X IF (MK .EQ. 1) THEN X ACC = P2(0) * U(NR,0) + U(NR,1) - P4(0) * V(NR,1) - X * P4(1) * V(NR,2) X IF (MK .LE. NK) X * ACC = ACC + P1(0) * W(NR,0) - P3(0) * U(NR,1) X SSUMM = ACC X RETURN X END IF XC XC ... BOTH CASES MK <= NK AND MK > NK XC X ACC = 0.0D0 X DO 10 I = 0, 2*MK - 1 X ACC = ACC - P4(I) * V(NR,I+1) X 10 CONTINUE XC XC ... CASE MK <= NK XC X IF (MK .LE. NK) THEN X DO 20 I = 0, 2*MK - 2 X ACC = ACC + P1(I) * W(NR,I) + P2(I) * U(NR,I) - X * P3(I) * U(NR,I+1) X 20 CONTINUE X ACC = ACC + P2(2*MK-1) * U(NR,2*MK-1) XC XC ... CASE MK > NK XC X ELSE X DO 30 I = 0, NK + MK X ACC = ACC + P2(I) * U(NR,I) X 30 CONTINUE X IF (NK .NE. 0) THEN X DO 40 I = 0, NK + MK - 2 X ACC = ACC - P3(I) * U(NR,I+1) X 40 CONTINUE X DO 50 I = 0, 2*NK - 1 X ACC = ACC + P1(I) * W(NR,I) X 50 CONTINUE X END IF X END IF XC XC ... SET THE RESULT VALUE XC X SSUMM = ACC X RETURN X END SHAR_EOF if test 8089 -ne "`wc -c ssumm.f`" then echo shar: error transmitting ssumm.f '(should have been 8089 characters)' fi echo shar: extracting stomhf.f sed 's/^X//' << \SHAR_EOF > stomhf.f XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC SUBROUTINE NAME - STOMHF + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC DESCRIPTION: + XC + XC STORES IN THE VECTOR LMAT THE COEFFICIENTS OF THE SYSTEM TO BE + XC SOLVED FOR THE COMPUTATION OF THE ALPHA , ALPHA', BETA , BETA' . + XC i i i i + XC + XC USAGE: + XC CALL STOMAT (LMAT, C, D, P1, NK, MK, IFLAG) + XC + XC ARGUMENTS: + XC + XC LMAT - OUTPUT REAL VECTOR WHICH CONTAINS THE COEFFICIENTS OF + XC THE AUXILIARY SYSTEM. + XC + XC C - INPUT REAL VECTOR CONTAINING THE c 'S. + XC i + XC + XC D - INPUT REAL VECTOR CONTAINING THE d 'S. + XC i + XC + XC P1 - INPUT REAL VECTOR CONTAINING THE COEFFICIENTS OF THE + XC (1) + XC POLYNOMIAL P OF DEGREE NK. + XC k + XC + XC NK - INPUT INTEGER EQUAL TO THE DIMENSION OF THE PREVIOUS + XC MATRIX OF THE SYSTEM. + XC + XC MK - INPUT INTEGER EQUAL TO THE DIMENSION OF THE NEW JUMP. + XC + XC IFLAG - INPUT INTEGER FLAG. + XC IF IFLAG = 1 THE CASE MK <= NK IS CONSIDERED. + XC IF IFLAG = 2 THE CASE MK > NK IS CONSIDERED. + XC + XC EXTERNAL MODULES: + XC + XC - DOUBLE PRECISION FUNCTION FH + XC + XC REMARKS: + XC + XC - ALL THE REAL ARGUMENTS PASSED FROM THE CALLING PROGRAM MUST BE + XC IN DOUBLE PRECISION. + XC - IN THE CALLING PROCEDURE THE DIMENSIONS FOR THE VECTOR ARGUMENTS + XC ARE: + XC LMAT(*) + XC C(0:*) + XC D(0:*) + XC P1(0:*) + XC + XC AUTHORS: + XC + XC C. BREZINSKI + XC UNIVERSITY OF LILLE - FRANCE + XC M. REDIVO ZAGLIA + XC UNIVERSITY OF PADOVA - ITALY + XC + XC REFERENCE: + XC + XC - TREATMENT OF NEAR-BREAKDOWN IN THE CGS ALGORITHM + XC NUMERICAL ALGORITHMS 7 (1994) 33-73 + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC X SUBROUTINE STOMHF (LMAT, C, D, P1, NK, MK, IFLAG) XC XC ... SPECIFICATIONS FOR ARGUMENTS XC X INTEGER IFLAG, NK, MK X DOUBLE PRECISION LMAT(1), C(0:1), D(0:1), P1(0:1) XC XC ... SPECIFICATIONS FOR LOCAL VARIABLES XC X INTEGER ICOL, IROW, NI, NN, NN2, NNM XC XC ... SPECIFICATIONS FOR EXTERNAL FUNCTIONS XC X DOUBLE PRECISION FH XC XC ... EXECUTABLE STATEMENTS XC XC XC ... CASE MK <= NK XC X IF (IFLAG .EQ. 1) THEN X NN = MK - 1 X ELSE XC XC ... CASE MK > NK XC X NN = NK X END IF XC XC ... COMPUTE THE DIMENSION OF THE SYSTEM XC X NNM = NN + MK XC XC ... CLEAR THE ARRAY COMPONENTS XC X DO 10 IROW = 1, NNM * NNM X LMAT(IROW) = 0.0D0 X 10 CONTINUE XC XC ... BOTH CASES: MK <= NK AND MK > NK XC X IF (NK .NE. 0) THEN X DO 30 ICOL = NN*2, 2, -2 X NI = NN - ICOL/2 X NN2 = NNM * ICOL X DO 20 IROW = NI, NI+NN-1 X LMAT(NN2-IROW) = D(IROW-NI) X LMAT(NN2+NNM-IROW) = C(IROW-NI) X 20 CONTINUE X 30 CONTINUE X END IF X DO 50 IROW = MK-1, 0, -1 X NN2 = MK - IROW X DO 40 ICOL = 0, NN*2, 2 X NI = ICOL/2 X LMAT(NN2+NNM*ICOL) = FH(IROW,NI,C,P1,NK,0) X IF (ICOL .NE. NN*2) THEN X LMAT(NN2+NNM*(ICOL+1)) = FH(IROW,NI,D,P1,NK,1) X END IF X 40 CONTINUE X 50 CONTINUE XC XC ... CASE MK > NK XC X IF (IFLAG .EQ. 2) THEN X IF (NK .NE. 0) THEN X DO 70 ICOL = 1, MK-NN-1 X NN2 = NNM * (NN*2 + ICOL + 1) X DO 60 IROW = 0, NN-1 X LMAT(NN2-IROW) = C(IROW+ICOL) X 60 CONTINUE X 70 CONTINUE X END IF X DO 90 IROW = MK-1, 0, -1 X NN2 = NNM*NN*2 + MK - IROW X DO 80 ICOL= 1, MK-NN-1 X LMAT(NN2+NNM*ICOL) = FH(IROW,NN+ICOL,C,P1,NK,0) X 80 CONTINUE X 90 CONTINUE X END IF X RETURN X END SHAR_EOF if test 7645 -ne "`wc -c stomhf.f`" then echo shar: error transmitting stomhf.f '(should have been 7645 characters)' fi echo shar: extracting storhf.f sed 's/^X//' << \SHAR_EOF > storhf.f XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC SUBROUTINE NAME - STORHF + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC DESCRIPTION: + XC + XC STORES IN THE ARRAY RMAT THE RIGHT HAND SIDES OF THE SYSTEM TO BE + XC SOLVED FOR THE COMPUTATION OF THE ALPHA , ALPHA', BETA , BETA' . + XC i i i i + XC + XC USAGE: + XC CALL STORHF (RMAT, C, D, P1, NK, MK, IFLAG) + XC + XC ARGUMENTS: + XC + XC RMAT - OUTPUT REAL MATRIX OF DIMENSION (2,0:*) WHICH CONTAINS + XC THE 2 RIGHT HAND SIDES OF THE SYSTEM. + XC + XC C - INPUT REAL VECTOR CONTAINING THE c 'S. + XC i + XC + XC D - INPUT REAL VECTOR CONTAINING THE d 'S. + XC i + XC + XC P1 - INPUT REAL VECTOR CONTAINING THE COEFFICIENTS OF THE + XC (1) + XC POLYNOMIAL P OF DEGREE NK. + XC k + XC + XC NK - INPUT INTEGER EQUAL TO THE DIMENSION OF THE PREVIOUS + XC MATRIX OF THE SYSTEM. + XC + XC MK - INPUT INTEGER EQUAL TO THE DIMENSION OF THE NEW JUMP. + XC + XC IFLAG - INPUT INTEGER FLAG. + XC IF IFLAG = 1 THE CASE MK <= NK IS CONSIDERED. + XC IF IFLAG = 2 THE CASE MK > NK IS CONSIDERED. + XC + XC EXTERNAL MODULES: + XC + XC - DOUBLE PRECISION FUNCTION FH + XC + XC REMARKS: + XC + XC - ALL THE REAL ARGUMENTS PASSED FROM THE CALLING PROGRAM MUST BE + XC IN DOUBLE PRECISION. + XC - IN THE CALLING PROCEDURE THE DIMENSIONS FOR THE MATRIX AND VECTOR + XC ARGUMENTS ARE: + XC RMAT(2,0:*) + XC C(0:*) + XC D(0:*) + XC P1(0:*) + XC + XC AUTHORS: + XC + XC C. BREZINSKI + XC UNIVERSITY OF LILLE - FRANCE + XC M. REDIVO ZAGLIA + XC UNIVERSITY OF PADOVA - ITALY + XC + XC REFERENCE: + XC + XC - TREATMENT OF NEAR-BREAKDOWN IN THE CGS ALGORITHM + XC NUMERICAL ALGORITHMS 7 (1994) 33-73 + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC X SUBROUTINE STORHF (RMAT, C, D, P1, NK, MK, IFLAG) XC XC ... SPECIFICATIONS FOR ARGUMENTS XC X INTEGER IFLAG, NK, MK X DOUBLE PRECISION RMAT(2,0:1), C(0:1), D(0:1), P1(0:1) XC XC ... SPECIFICATIONS FOR LOCAL VARIABLES XC X INTEGER IROW, ICOL, NNM, NN2 X DOUBLE PRECISION CONST XC XC ... SPECIFICATIONS FOR EXTERNAL FUNCTIONS XC X DOUBLE PRECISION FH XC XC ... EXECUTABLE STATEMENTS XC XC XC ... COMPUTE THE MAXIMUM INDEX OF THE XC COMPONENTS IN THE ARRAY XC XC ... CASE MK <= NK XC X IF (IFLAG .EQ. 1) THEN X NNM = 2*(MK-1) X ELSE XC XC ... CASE MK > NK XC X NNM = NK + MK -1 X END IF XC XC ... CLEAR THE ARRAY COMPONENTS XC X DO 20 ICOL = 0, NNM X DO 10 IROW = 1, 2 X RMAT(IROW,ICOL) = 0.0D0 X 10 CONTINUE X 20 CONTINUE XC XC ... STORES THE RIGHT HAND SIDE FOR THE XC SYSTEM THAT COMPUTES THE BETA AND BETA' XC i i XC X IROW = 1 XC XC ... BOTH CASES: MK <= NK AND MK > NK XC X DO 30 ICOL = MK-1, 0, -1 X RMAT(IROW,MK-1-ICOL) = D(ICOL) X 30 CONTINUE XC XC ... STORES THE RIGHT HAND SIDE FOR THE XC SYSTEM THAT COMPUTES THE ALPHA AND ALPHA' XC i i XC X IROW = 2 XC XC ... CASE MK <= NK XC X IF (IFLAG .EQ. 1) CONST = -C(0)/D(0) XC XC ... BOTH CASES: MK <= NK AND MK > NK XC X IF (NK .NE. 0) THEN X DO 40 ICOL = MK, NNM X NN2 = 2*MK-1-ICOL X RMAT(IROW,ICOL) = - C(NN2) XC XC ... CASE MK <= NK XC X IF (IFLAG .EQ. 1) X * RMAT(IROW,ICOL) = RMAT(IROW,ICOL) - CONST * D(NN2) X 40 CONTINUE X END IF XC XC XC ... BOTH CASES: MK <= NK AND MK > NK XC X DO 50 ICOL = 0, MK-1 X NN2 = MK-1-ICOL X RMAT(IROW,ICOL) = - FH(NN2,MK,C,P1,NK,0) XC XC ... CASE MK <= NK XC X IF (IFLAG .EQ. 1) X * RMAT(IROW,ICOL) = RMAT(IROW,ICOL)-CONST*FH(NN2,MK-1,D,P1,NK,1) X 50 CONTINUE X RETURN X END SHAR_EOF if test 7912 -ne "`wc -c storhf.f`" then echo shar: error transmitting storhf.f '(should have been 7912 characters)' fi echo shar: extracting sunorm.f sed 's/^X//' << \SHAR_EOF > sunorm.f XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC FUNCTION NAME - SUNORM + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC DESCRIPTION: + XC + XC DOUBLE PRECISION FUNCTION: + XC COMPUTES THE MAXIMUM NORM OF A VECTOR. + XC + XC USAGE: + XC SUNORM (IP, VEC) + XC + XC ARGUMENTS: + XC + XC IP - INPUT DIMENSION OF THE VECTOR SPACE (NUMBER OF + XC ELEMENTS OF THE VECTOR TO BE CONSIDERED). + XC + XC VEC - INPUT REAL VECTOR OF DIMENSION IR >= IP. + XC + XC REMARKS: + XC + XC - THE REAL VECTOR VEC, PASSED FROM THE CALLING PROGRAM, MUST BE + XC IN DOUBLE PRECISION. + XC - IN THE CALLING PROCEDURE THE MINIMAL DIMENSION FOR THE + XC VECTOR ARGUMENT (IR=IP) IS: + XC VEC(IP) + XC + XC AUTHORS: + XC + XC C. BREZINSKI + XC UNIVERSITY OF LILLE - FRANCE + XC M. REDIVO ZAGLIA + XC UNIVERSITY OF PADOVA - ITALY + XC + XC REFERENCE: + XC + XC - TREATMENT OF NEAR-BREAKDOWN IN THE CGS ALGORITHM + XC NUMERICAL ALGORITHMS 7 (1994) 33-73 + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC X DOUBLE PRECISION FUNCTION SUNORM (IP, VEC) XC XC ... SPECIFICATIONS FOR ARGUMENTS XC X DOUBLE PRECISION VEC(1) X INTEGER IP XC XC ... SPECIFICATIONS FOR LOCAL VARIABLES XC X DOUBLE PRECISION ANORM, AMAXN X INTEGER I XC XC ... EXECUTABLE STATEMENTS XC X AMAXN = 0.0D0 X DO 10 I = 1, IP X ANORM = DABS(VEC(I)) X IF (ANORM .GT. AMAXN) AMAXN = ANORM X 10 CONTINUE XC XC ... SET THE RESULT VALUE XC X SUNORM = AMAXN X RETURN X END SHAR_EOF if test 4047 -ne "`wc -c sunorm.f`" then echo shar: error transmitting sunorm.f '(should have been 4047 characters)' fi echo shar: extracting zsumm.f sed 's/^X//' << \SHAR_EOF > zsumm.f XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC FUNCTION NAME - ZSUMM + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC + XC DESCRIPTION: + XC + XC DOUBLE PRECISION FUNCTION: + XC COMPUTES THE COMPONENT NR OF A VECTOR OBTAINED AS A LINEAR + XC COMBINATION, WITH THE COEFFICIENTS STORED IN THE VECTORS P1, P2 + XC AND P3, OF THE COLUMN VECTORS STORED IN V, W AND IN U. + XC + XC USAGE: + XC ZSUMM (P1, P2, P3, V, W, U, IR, NR, NK, MK) + XC + XC ARGUMENTS: + XC + XC P1 - INPUT REAL VECTOR WHICH CONTAINS THE COEFFICIENTS OF THE + XC POLYNOMIAL t(xi)**2. + XC + XC P2 - INPUT REAL VECTOR WHICH CONTAINS THE COEFFICIENTS OF THE + XC POLYNOMIAL t(xi) * q(xi). + XC + XC P3 - INPUT REAL VECTOR WHICH CONTAINS THE COEFFICIENTS OF THE + XC POLYNOMIAL q(xi)**2. + XC + XC V - INPUT REAL MATRIX CONTAINING A CERTAIN NUMBER OF COLUMN + XC VECTORS. + XC + XC W - INPUT REAL MATRIX CONTAINING A CERTAIN NUMBER OF COLUMN + XC VECTORS. + XC + XC U - INPUT REAL MATRIX CONTAINING A CERTAIN NUMBER OF COLUMN + XC VECTORS. + XC + XC IR - INPUT ROWS DIMENSION EXACTLY AS SPECIFIED IN THE + XC DIMENSION STATEMENTS IN THE CALLING PROGRAM FOR THE + XC MATRICES V, W AND U. + XC + XC NR - INPUT INTEGER EQUAL TO THE COMPONENT WHICH HAS TO BE + XC USED AND COMPUTED. + XC + XC NK - INPUT INTEGER EQUAL TO THE DIMENSION OF THE PREVIOUS + XC MATRIX OF THE SYSTEM. + XC + XC MK - INPUT INTEGER EQUAL TO THE DIMENSION OF THE NEW JUMP. + XC + XC REMARKS: + XC + XC - ALL THE REAL ARGUMENTS PASSED FROM THE CALLING PROGRAM MUST BE + XC IN DOUBLE PRECISION. + XC - IN THE CALLING PROCEDURE THE DIMENSIONS FOR THE ARRAY ARGUMENTS + XC ARE: + XC P1(0:*) + XC P2(0:*) + XC P3(0:*) + XC V(IR,0:*) + XC W(IR,0:*) + XC U(IR,0:*) + XC + XC AUTHORS: + XC + XC C. BREZINSKI + XC UNIVERSITY OF LILLE - FRANCE + XC M. REDIVO ZAGLIA + XC UNIVERSITY OF PADOVA - ITALY + XC + XC REFERENCE: + XC + XC - TREATMENT OF NEAR-BREAKDOWN IN THE CGS ALGORITHM + XC NUMERICAL ALGORITHMS 7 (1994) 33-73 + XC + XC+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ XC X DOUBLE PRECISION FUNCTION ZSUMM (P1, P2, P3, V, W, U, IR, NR, X * NK, MK) XC XC ... SPECIFICATIONS FOR ARGUMENTS XC X INTEGER IR, NK, NR, MK X DOUBLE PRECISION P1(0:1), P2(0:1), P3(0:2), V(IR,0:1), W(IR,0:1), X * U(IR,0:1) XC XC ... SPECIFICATIONS FOR LOCAL VARIABLES XC X DOUBLE PRECISION ACC X INTEGER I XC XC ... EXECUTABLE STATEMENTS XC XC XC ... CASE MK = 1 XC X IF (MK .EQ. 1) THEN X ACC = P3(0)*V(NR,0) + P3(1)*V(NR,1) + P3(2)*V(NR,2) X IF (MK .LE. NK) X * ACC = ACC + 2.0D0 * (P2(0)*U(NR,0) + P2(1)*U(NR,1)) X * + P1(0) * W(NR,0) X ZSUMM = ACC X RETURN X END IF XC X ACC = 0.0D0 XC XC ... BOTH CASES MK <= NK AND MK > NK XC X DO 10 I = 0, 2*MK X ACC = ACC + P3(I) * V(NR,I) X 10 CONTINUE XC XC ... CASE MK <= NK XC X IF (MK .LE. NK) THEN X DO 20 I = 0, 2*MK - 2 X ACC = ACC + P1(I) * W(NR,I) + 2.0D0 * P2(I) * U(NR,I) X 20 CONTINUE X ACC = ACC + 2.0D0 * P2(2*MK-1) * U(NR,2*MK-1) XC XC ... CASE MK > NK XC X ELSE IF (NK .NE. 0) THEN X DO 30 I = 0, 2*NK - 2 X ACC = ACC + P1(I) * W(NR,I) X 30 CONTINUE X DO 40 I = 0, NK + MK - 1 X ACC = ACC + 2.0D0 * P2(I) * U(NR,I) X 40 CONTINUE X END IF XC XC ... SET THE RESULT VALUE XC X ZSUMM = ACC X RETURN X END SHAR_EOF if test 7603 -ne "`wc -c zsumm.f`" then echo shar: error transmitting zsumm.f '(should have been 7603 characters)' fi # End of shell archive exit 0