# Date: Fri, 16 Sep 1994 12:16:27 +0200 # From: "Michela Redivo Zaglia -tel.+39 49 8287625" # 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