***************************************************************************
* All the software contained in this library is protected by copyright. *
* Permission to use, copy, modify, and distribute this software for any *
* purpose without fee is hereby granted, provided that this entire notice *
* is included in all copies of any software which is or includes a copy *
* or modification of this software and in all copies of the supporting *
* documentation for such software. *
***************************************************************************
* THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED *
* WARRANTY. IN NO EVENT, NEITHER THE AUTHORS, NOR THE PUBLISHER, NOR ANY *
* MEMBER OF THE EDITORIAL BOARD OF THE JOURNAL "NUMERICAL ALGORITHMS", *
* NOR ITS EDITOR-IN-CHIEF, BE LIABLE FOR ANY ERROR IN THE SOFTWARE, ANY *
* MISUSE OF IT OR ANY DAMAGE ARISING OUT OF ITS USE. THE ENTIRE RISK OF *
* USING THE SOFTWARE LIES WITH THE PARTY DOING SO. *
***************************************************************************
* ANY USE OF THE SOFTWARE CONSTITUTES ACCEPTANCE OF THE TERMS OF THE *
* ABOVE STATEMENT. *
***************************************************************************
* AUTHORS:
C. BREZINSKI AND H. SADOK
UNIVERSITY OF LILLE - FRANCE
M. REDIVO ZAGLIA
UNIVERSITY OF PADOVA - ITALY
* REFERENCES:
- AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE ALGORITHMS,
NUMERICAL ALGORITHMS, 1(1991), PP. 261-284.
- ADDENDUM TO "AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS
TYPE ALGORITHMS",
NUMERICAL ALGORITHMS, 2(1992), PP. 133-136.
* SOFTWARE REVISION DATE:
MAY, 1992
* REMARKS:
- THESE PROGRAMS CAN POSSIBLY GIVE RESULTS DIFFERENT FROM THOSE PRINTED
IN THE FIRST PAPER.
THIS IS DUE TO SOME MODIFICATIONS IN THE PROGRAMMING OF THE ALGORITHM
AND TO THE NUMERICAL INSTABILITY OF THE EXAMPLE AS SHOWN BY THE
NUMERICAL RESULTS GIVEN IN THE PAPERS.
- SOME MODULES HAVE BEEN MODIFIED IN THE CODE WITH RESPECT TO THE FIRST
VERSION PUT INTO NETLIB.
THE MODULES MODIFIED ARE:
+ MAIN PROGRAMS MMRZ, MBMRZ, MSMRZ, MBSMRZ
+ SUBROUTINES MRZ, BMRZ, SMRZ, BSMRZ, GSOLVE, MATVEC,
PINNER, SOLSYS
+ DOUBLE PRECISION FUNCTIONS NVSUMM
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C PROGRAM NAME - MMRZ +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C DESCRIPTION: +
C +
C MAIN PROGRAM TO TEST THE SUBROUTINE MRZ. +
C +
C MODULES REQUIRED: +
C +
C - SUBROUTINES MRZ, MATVEC +
C - DOUBLE PRECISION FUNCTIONS EUNORM, PINNER, SSTRIA, VSUMMA +
C +
C AUTHORS: +
C +
C C. BREZINSKI AND H. SADOK +
C UNIVERSITY OF LILLE - FRANCE +
C M. REDIVO ZAGLIA +
C UNIVERSITY OF PADOVA - ITALY +
C +
C REFERENCES: +
C +
C - AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE ALGORITHMS, +
C NUMERICAL ALGORITHMS, 1(1991), PP. 261-284. +
C +
C - ADDENDUM TO "AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE +
C ALGORITHMS", NUMERICAL ALGORITHMS, 2(1992), PP. 133-136. +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
PROGRAM MMRZ
C
C ... SPECIFICATIONS FOR VARIABLES AND CONSTANTS
C
INTEGER I, IER, IFLAG, INIT, IR, J, K, MAXBCK, N, NBC, NK
DOUBLE PRECISION EPS, TMP
C
C ... SPECIFICATIONS FOR PARAMETER CONSTANTS
C
PARAMETER (N=12, NBC=12, IR=20, MAXBCK=7, EPS=1.0D-8)
C
C ... SPECIFICATIONS FOR VECTORS AND MATRICES
C
DOUBLE PRECISION A(IR,IR), AB(IR), V(IR,0:MAXBCK), D(0:MAXBCK-1),
* B(0:MAXBCK), BETA(0:MAXBCK-1), ALPHA(0:MAXBCK), P(0:MAXBCK),
* Z(IR), Y(IR), X(IR), R(IR), ADIG(IR)
C
C ... SPECIFICATIONS FOR EXTERNAL MODULES
C
DOUBLE PRECISION EUNORM
C
C ... EXECUTABLE STATEMENTS
C
C
C ... CHOICE OF IFLAG
C
C IFLAG .EQ. 0 Y IS SET TO r BY THE SUBROUTINE
C 0
C IFLAG .NE. 0 THE USER MUST DEFINE Y IN
C THE MAIN PROGRAM
C
IFLAG = 1
C
WRITE (*,*) ' OUTPUT RESULTS FROM MMRZ '
WRITE (*,*)
WRITE (*,*) ' DIMENSION OF THE SYSTEM = ',N
WRITE (*,*) ' NUMBER OF CALLS = ',NBC
WRITE (*,*) ' MAXBCK (MAX JUMP MK) = ',MAXBCK
WRITE (*,*) ' IFLAG = ',IFLAG
WRITE (*,*) ' EPS = ',EPS
WRITE (*,*)
WRITE (*,*) ' ITER. NK NORM'
C
C ... STORE THE SYSTEM AND THE VECTORS x AND y
C 0
C
DO 20 I = 1, N
X(I) = 0.0D0
Y(I) = 1.0D0
DO 10 J = 1, N
A(I,J) = 0.0D0
10 CONTINUE
20 CONTINUE
A(1,N) = -1.0D0
AB(1) = -N
DO 30 I = 2, N
A(I,I-1) = 1.0D0
AB(I) = I-1
30 CONTINUE
C
INIT = 0
C
DO 40 K = 1, NBC
C
C ... CALL MRZ
C
CALL MRZ (N, A, AB, Y, X, R, Z, MAXBCK, V, P, B, ALPHA, BETA,
* D, IR, NK, EPS, INIT, IFLAG, IER)
C
TMP = DSQRT(EUNORM(N,R))
WRITE (*,100) K, NK, TMP
IF (IER .NE. 0) THEN
WRITE (*,*)
IF (IER .NE. 200 .AND. IER .NE. 400) THEN
WRITE (*,*) ' STOP DUE TO ERROR ', IER
STOP
ELSE IF (IER .EQ. 200) THEN
WRITE (*,*) ' THE SOLUTION HAS BEEN OBTAINED'
ELSE
WRITE (*,*) ' SOLUTION NOT OBTAINED AFTER REACHING'
WRITE (*,*) ' THE DIMENSION OF THE SYSTEM'
END IF
GO TO 50
END IF
40 CONTINUE
WRITE (*,*)
WRITE (*,300) ' SOLUTION NOT OBTAINED AFTER', NBC
WRITE (*,*) ' ITERATIONS OF THE SUBROUTINE'
C
50 CONTINUE
WRITE (*,*)
IF (IER .EQ. 200) THEN
WRITE (*,*)
* ' FINAL SOLUTION N. OF SIGN. DIG. '
ELSE
WRITE (*,*)
* ' LAST SOLUTION N. OF SIGN. DIG. '
END IF
DO 60 I= 1, N
TMP = DABS(X(I)-DBLE(I))
IF (TMP .EQ. 0.0D0) THEN
ADIG(I) = 999.00
ELSE
ADIG(I) = -DLOG10(TMP/DBLE(I))
END IF
60 CONTINUE
WRITE (*,200) (X(I),ADIG(I),I=1,N)
STOP
100 FORMAT (I5,I8,D28.15)
200 FORMAT (D25.15,F15.2)
300 FORMAT (T2,A,I4)
END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C PROGRAM NAME - MBMRZ +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C DESCRIPTION: +
C +
C MAIN PROGRAM TO TEST THE SUBROUTINE BMRZ. +
C +
C MODULES REQUIRED: +
C +
C - SUBROUTINES BMRZ, MATVEC +
C - DOUBLE PRECISION FUNCTIONS EUNORM, PINNER, SSTRIA, VSUMMA +
C +
C AUTHORS: +
C +
C C. BREZINSKI AND H. SADOK +
C UNIVERSITY OF LILLE - FRANCE +
C M. REDIVO ZAGLIA +
C UNIVERSITY OF PADOVA - ITALY +
C +
C REFERENCES: +
C +
C - AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE ALGORITHMS, +
C NUMERICAL ALGORITHMS, 1(1991), PP. 261-284. +
C +
C - ADDENDUM TO "AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE +
C ALGORITHMS", NUMERICAL ALGORITHMS, 2(1992), PP. 133-136. +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
PROGRAM MBMRZ
C
C ... SPECIFICATIONS FOR VARIABLES AND CONSTANTS
C
INTEGER I, IER, IFLAG, INIT, IR, J, K, MAXBCK, N, NBC, NK
DOUBLE PRECISION EPS, TMP
C
C ... SPECIFICATIONS FOR PARAMETER CONSTANTS
C
PARAMETER (N=6, NBC=6, IR=20, MAXBCK=4, EPS=1.0D-8)
C
C ... SPECIFICATIONS FOR VECTORS AND MATRICES
C
DOUBLE PRECISION A(IR,IR), AB(IR), V(IR,0:MAXBCK), D(0:MAXBCK-1),
* B(0:MAXBCK-1), BETA(0:MAXBCK-1), YSAV(IR), Y(IR),
* X(IR), R(IR), ADIG(IR)
C
C ... SPECIFICATIONS FOR EXTERNAL MODULES
C
DOUBLE PRECISION EUNORM
C
C ... EXECUTABLE STATEMENTS
C
C
C ... CHOICE OF IFLAG
C
C IFLAG .EQ. 0 Y IS SET TO r BY THE SUBROUTINE
C 0
C IFLAG .NE. 0 THE USER MUST DEFINE Y IN
C THE MAIN PROGRAM
C
IFLAG = 0
C
WRITE (*,*) ' OUTPUT RESULTS FROM MBMRZ '
WRITE (*,*)
WRITE (*,*) ' DIMENSION OF THE SYSTEM = ',N
WRITE (*,*) ' NUMBER OF CALLS = ',NBC
WRITE (*,*) ' MAXBCK (MAX JUMP MK) = ',MAXBCK
WRITE (*,*) ' IFLAG = ',IFLAG
WRITE (*,*) ' EPS = ',EPS
WRITE (*,*)
WRITE (*,*) ' ITER. NK NORM'
C
C ... STORE THE SYSTEM AND THE VECTORS x AND y
C 0
C
DO 20 I = 1, N
X(I) = 0.0D0
Y(I) = 1.0D0
DO 10 J = 1, N
A(I,J) = 0.0D0
10 CONTINUE
20 CONTINUE
A(1,N) = -1.0D0
AB(1) = -N
DO 30 I = 2, N
A(I,I-1) = 1.0D0
AB(I) = I-1
30 CONTINUE
C
INIT = 0
C
DO 40 K = 1, NBC
C
C ... CALL BMRZ
C
CALL BMRZ (N, A, AB, Y, X, R, YSAV, MAXBCK, V, B, BETA,
* D, IR, NK, EPS, INIT, IFLAG, IER)
C
TMP = DSQRT(EUNORM(N,R))
WRITE (*,100) K, NK, TMP
IF (IER .NE. 0) THEN
WRITE (*,*)
IF (IER .NE. 200 .AND. IER .NE. 400) THEN
WRITE (*,*) ' STOP DUE TO ERROR ', IER
STOP
ELSE IF (IER .EQ. 200) THEN
WRITE (*,*) ' THE SOLUTION HAS BEEN OBTAINED'
ELSE
WRITE (*,*) ' SOLUTION NOT OBTAINED AFTER REACHING'
WRITE (*,*) ' THE DIMENSION OF THE SYSTEM'
END IF
GO TO 50
END IF
40 CONTINUE
WRITE (*,*)
WRITE (*,300) ' SOLUTION NOT OBTAINED AFTER', NBC
WRITE (*,*) ' ITERATIONS OF THE SUBROUTINE'
C
50 CONTINUE
WRITE (*,*)
WRITE (*,*)
IF (IER .EQ. 200) THEN
WRITE (*,*)
* ' FINAL SOLUTION N. OF SIGN. DIG. '
ELSE
WRITE (*,*)
* ' LAST SOLUTION N. OF SIGN. DIG. '
END IF
DO 60 I= 1, N
TMP = DABS(X(I)-DBLE(I))
IF (TMP .EQ. 0.0D0) THEN
ADIG(I) = 999.00
ELSE
ADIG(I) = -DLOG10(TMP/DBLE(I))
END IF
60 CONTINUE
WRITE (*,200) (X(I),ADIG(I),I=1,N)
STOP
100 FORMAT (I5,I8,D28.15)
200 FORMAT (D25.15,F15.2)
300 FORMAT (T2,A,I4)
END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C PROGRAM NAME - MSMRZ +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C DESCRIPTION: +
C +
C MAIN PROGRAM TO TEST THE SUBROUTINE SMRZ. +
C +
C MODULES REQUIRED: +
C +
C - SUBROUTINES SMRZ, MATVEC +
C - DOUBLE PRECISION FUNCTIONS EUNORM, PINNER, SSTRIA, VSUMMA +
C +
C AUTHORS: +
C +
C C. BREZINSKI AND H. SADOK +
C UNIVERSITY OF LILLE - FRANCE +
C M. REDIVO ZAGLIA +
C UNIVERSITY OF PADOVA - ITALY +
C +
C REFERENCES: +
C +
C - AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE ALGORITHMS, +
C NUMERICAL ALGORITHMS, 1(1991), PP. 261-284. +
C +
C - ADDENDUM TO "AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE +
C ALGORITHMS", NUMERICAL ALGORITHMS, 2(1992), PP. 133-136. +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
PROGRAM MSMRZ
C
C ... SPECIFICATIONS FOR VARIABLES AND CONSTANTS
C
INTEGER I, IER, IFLAG, INIT, IR, J, K, MAXBCK, N, NK
DOUBLE PRECISION EPS, TMP
C
C ... SPECIFICATIONS FOR PARAMETER CONSTANTS
C
PARAMETER (N=6, NBC=6, IR=20, MAXBCK=4, EPS=1.0D-8)
C
C ... SPECIFICATIONS FOR VECTORS AND MATRICES
C
DOUBLE PRECISION A(IR,IR), AB(IR), V(IR,0:MAXBCK), D(0:MAXBCK),
* B(0:MAXBCK), BETA(0:MAXBCK-1), DCOEF(0:MAXBCK),
* RSAV(IR), Y(IR), X(IR), R(IR), ADIG(IR)
C
C ... SPECIFICATIONS FOR EXTERNAL MODULES
C
DOUBLE PRECISION EUNORM
C
C ... EXECUTABLE STATEMENTS
C
C
C ... CHOICE OF IFLAG
C
C IFLAG .EQ. 0 Y IS SET TO r BY THE SUBROUTINE
C 0
C IFLAG .NE. 0 THE USER MUST DEFINE Y IN
C THE MAIN PROGRAM
C
IFLAG = 1
C
WRITE (*,*) ' OUTPUT RESULTS FROM MSMRZ '
WRITE (*,*)
WRITE (*,*) ' DIMENSION OF THE SYSTEM = ',N
WRITE (*,*) ' NUMBER OF CALLS = ',NBC
WRITE (*,*) ' MAXBCK (MAX JUMP MK) = ',MAXBCK
WRITE (*,*) ' IFLAG = ',IFLAG
WRITE (*,*) ' EPS = ',EPS
WRITE (*,*)
WRITE (*,*) ' ITER. NK NORM'
C
C ... STORE THE SYSTEM AND THE VECTORS x AND y
C 0
C
DO 20 I = 1, N
X(I) = 0.0D0
Y(I) = 1.0D0
DO 10 J = 1, N
A(I,J) = 0.0D0
10 CONTINUE
20 CONTINUE
A(1,N) = -1.0D0
AB(1) = -N
DO 30 I = 2, N
A(I,I-1) = 1.0D0
AB(I) = I-1
30 CONTINUE
C
INIT = 0
C
DO 40 K = 1, NBC
C
C ... CALL SMRZ
C
CALL SMRZ (N, A, AB, Y, X, R, RSAV, MAXBCK, V, B, DCOEF, BETA,
* D, IR, NK, EPS, INIT, IFLAG, IER)
C
TMP = DSQRT(EUNORM(N,R))
WRITE (*,100) K, NK, TMP
IF (IER .NE. 0) THEN
WRITE (*,*)
IF (IER .NE. 200 .AND. IER .NE. 400) THEN
WRITE (*,*) ' STOP DUE TO ERROR ', IER
STOP
ELSE IF (IER .EQ. 200) THEN
WRITE (*,*) ' THE SOLUTION HAS BEEN OBTAINED'
ELSE
WRITE (*,*) ' SOLUTION NOT OBTAINED AFTER REACHING'
WRITE (*,*) ' THE DIMENSION OF THE SYSTEM'
END IF
GO TO 50
END IF
40 CONTINUE
WRITE (*,*)
WRITE (*,300) ' SOLUTION NOT OBTAINED AFTER', NBC
WRITE (*,*) ' ITERATIONS OF THE SUBROUTINE'
C
50 CONTINUE
WRITE (*,*)
IF (IER .EQ. 200) THEN
WRITE (*,*)
* ' FINAL SOLUTION N. OF SIGN. DIG. '
ELSE
WRITE (*,*)
* ' LAST SOLUTION N. OF SIGN. DIG. '
END IF
DO 60 I= 1, N
TMP = DABS(X(I)-DBLE(I))
IF (TMP .EQ. 0.0D0) THEN
ADIG(I) = 999.00
ELSE
ADIG(I) = -DLOG10(TMP/DBLE(I))
END IF
60 CONTINUE
WRITE (*,200) (X(I),ADIG(I),I=1,N)
STOP
100 FORMAT (I5,I8,D28.15)
200 FORMAT (D25.15,F15.2)
300 FORMAT (T2,A,I4)
END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C PROGRAM NAME - MBSMRZ +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C DESCRIPTION: +
C +
C MAIN PROGRAM TO TEST THE SUBROUTINE BSMRZ. +
C +
C MODULES REQUIRED: +
C +
C - SUBROUTINES BSMRZ, GSOLVE, MATVEC, SOLSYS, STOMAT, +
C STORHS +
C - DOUBLE PRECISION FUNCTIONS EUNORM, NVSUMM, PINNER +
C +
C AUTHORS: +
C +
C C. BREZINSKI AND H. SADOK +
C UNIVERSITY OF LILLE - FRANCE +
C M. REDIVO ZAGLIA +
C UNIVERSITY OF PADOVA - ITALY +
C +
C REFERENCES: +
C +
C - AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE ALGORITHMS, +
C NUMERICAL ALGORITHMS, 1(1991), PP. 261-284. +
C +
C - ADDENDUM TO "AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE +
C ALGORITHMS", NUMERICAL ALGORITHMS, 2(1992), PP. 133-136. +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
PROGRAM MBSMRZ
C
C ... SPECIFICATIONS FOR VARIABLES AND CONSTANTS
C
INTEGER I, IER, IFLAG, INIT, IR, J, K, KSAV, MAXBCK, MAXNN,
* MAXN, N, NBC, NK, NKSAV
DOUBLE PRECISION EPS, EPS1, EPS2, RBEST, TMP
C
C ... SPECIFICATIONS FOR PARAMETER CONSTANTS
C
PARAMETER (N=12, IR=15, MAXBCK=6, MAXNN=25, NBC=25,
* EPS=1.0D-8, EPS1=1.0D-14, EPS2=1.0D-8)
C
C ... SPECIFICATIONS FOR VECTORS AND MATRICES
C
DOUBLE PRECISION A(IR,IR), AB(IR), Y(IR), X(IR), R(IR), YSAV(IR),
* ADIG(IR), XBEST(IR), C(0:2*MAXBCK-1), D(0:2*MAXBCK-1),
* V(IR,0:MAXBCK), W(IR,0:MAXBCK-1), RMAT(2,0:2*MAXBCK-1),
* LMAT((2*MAXBCK-1)*(2*MAXBCK-1))
C
C ... SPECIFICATIONS FOR EXTERNAL FUNCTIONS
C
DOUBLE PRECISION EUNORM
C
C ... EXECUTABLE STATEMENTS
C
C
C ... CHOICE OF IFLAG
C
C IFLAG .EQ. 0 Y IS SET TO r BY THE SUBROUTINE
C 0
C IFLAG .NE. 0 THE USER MUST DEFINE Y IN
C THE MAIN PROGRAM
C
IFLAG = 0
C
WRITE (*,*) ' OUTPUT RESULTS FROM MBSMRZ '
WRITE (*,*)
WRITE (*,*) ' DIMENSION OF THE SYSTEM = ',N
WRITE (*,*) ' NUMBER OF CALLS = ',NBC
WRITE (*,*) ' MAXBCK (MAX JUMP MK) = ',MAXBCK
WRITE (*,*) ' MAXN (MAX NK+MK) = ',MAXNN
WRITE (*,*) ' IFLAG = ',IFLAG
WRITE (*,*) ' EPS (NEAR-BREAKDOWN) = ',EPS
WRITE (*,*) ' EPS1 (GAUSS) = ',EPS1
WRITE (*,*) ' EPS2 (SOLUTION) = ',EPS2
WRITE (*,*)
WRITE (*,*) ' ITER. NK NORM'
C
DO 20 I = 1, N
X(I) = 0.0D0
Y(I) = 1.0D0
DO 10 J = 1, N
A(I,J) = 0.0D0
10 CONTINUE
20 CONTINUE
A(1,N) = -1.0D0
AB(1) = -N
DO 30 I = 2, N
A(I,I-1) = 1.0D0
AB(I) = I-1
30 CONTINUE
C
INIT = 0
MAXN = MAXNN
C
DO 50 K = 1, NBC
CALL BSMRZ (N, A, AB, Y, X, R, YSAV, MAXBCK, MAXN, V, W, LMAT,
* RMAT, D, C, IR, NK, EPS, EPS1, EPS2, INIT, IFLAG, IER)
TMP = DSQRT(EUNORM(N,R))
IF (K .EQ. 1 .OR. TMP .LT. RBEST) THEN
RBEST = TMP
KSAV = K
NKSAV = NK
DO 40 I = 1, N
XBEST(I) = X(I)
40 CONTINUE
END IF
WRITE (*,100) K, NK, TMP
IF (IER .EQ. 400) THEN
WRITE (*,*)
WRITE (*,*) ' SOLUTION NOT OBTAINED AFTER REACHING'
WRITE (*,*) ' THE DIMENSION N OF THE SYSTEM. THE'
WRITE (*,*) ' COMPUTATION CONTINUES UNTIL MAXN AT MOST'
WRITE (*,*)
IER = 0
END IF
IF (IER .NE. 0) THEN
WRITE (*,*)
IF (IER .EQ. 200) THEN
WRITE (*,*) ' THE SOLUTION HAS BEEN OBTAINED'
ELSE IF (IER .EQ. 300) THEN
WRITE (*,*) ' SOLUTION NOT OBTAINED AFTER'
WRITE (*,*) ' REACHING THE VALUE OF MAXN'
ELSE IF (IER .EQ. 600) THEN
WRITE (*,*) ' THE VALUE OF THE JUMP MK'
WRITE (*,*) ' IS GREATER THAN MAXBCK'
ELSE
WRITE (*,*) ' STOP DUE TO ERROR ', IER
STOP
END IF
GO TO 60
END IF
50 CONTINUE
WRITE (*,*)
WRITE (*,300) ' SOLUTION NOT OBTAINED AFTER', NBC
WRITE (*,*) ' ITERATIONS OF THE SUBROUTINE'
C
60 CONTINUE
WRITE (*,*)
DO 70 I= 1, N
TMP = DABS(X(I)-DBLE(I))
IF (TMP .EQ. 0.0D0) THEN
ADIG(I) = 999.00
ELSE
ADIG(I) = -DLOG10(TMP/DBLE(I))
END IF
70 CONTINUE
IF (IER .EQ. 200) THEN
WRITE (*,*)
* ' FINAL SOLUTION N. OF SIGN. DIG. '
WRITE (*,200) (X(I),ADIG(I),I=1,N)
ELSE
WRITE (*,*)
* ' LAST SOLUTION N. OF SIGN. DIG. '
WRITE (*,200) (X(I),ADIG(I),I=1,N)
WRITE (*,*)
WRITE (*,*) ' SOLUTION OF MINIMAL RESIDUAL OBTAINED AT'
WRITE (*,*) ' ITER. NK NORM'
WRITE (*,100) KSAV, NKSAV, RBEST
WRITE (*,*)
WRITE (*,400) (XBEST(I),I=1,N)
END IF
STOP
100 FORMAT (I5,I8,D28.15)
200 FORMAT (D25.15,F15.2)
300 FORMAT (T2,A,I4)
400 FORMAT (D25.15)
END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C SUBROUTINE NAME - MRZ +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C DESCRIPTION: +
C +
C COMPUTES THE SOLUTION OF A SYSTEM OF LINEAR EQUATIONS BY +
C THE MRZ ALGORITHM. +
C +
C USAGE: +
C CALL MRZ (N, A, AB, Y, X, R, Z, MAXBCK, V, P, B, ALPHA, BETA, +
C D, IR, NK, EPS, INIT, IFLAG, IER) +
C +
C ARGUMENTS: +
C +
C N - INPUT INTEGER VALUE, DIMENSION OF THE SYSTEM. +
C +
C A - INPUT REAL MATRIX OF DIMENSION (N,N) CONTAINING THE +
C MATRIX OF THE SYSTEM. +
C +
C AB - INPUT/WORKING REAL VECTOR OF DIMENSION N CONTAINING, +
C BEFORE THE FIRST CALL, THE RIGHT HAND SIDE OF THE SYSTEM. +
C THE VECTOR IS USED AS A WORKING AREA IN THE NEXT CALLS +
C AND THEN THE INPUT VALUES ARE ALWAYS DESTROYED BY THE +
C SUBROUTINE. +
C +
C Y - INPUT/WORKING REAL VECTOR OF DIMENSION N CONTAINING, +
C BEFORE THE FIRST CALL, THE AUXILIARY VECTOR y. +
C IF IFLAG .EQ. 0 THE VECTOR Y IS SET TO r = A x - b +
C 0 0 +
C BY THE SUBROUTINE DURING THE FIRST CALL. THE VECTOR IS +
C USED AS A WORKING AREA IN THE NEXT CALLS AND THEN THE +
C INPUT VALUES ARE ALWAYS DESTROYED BY THE SUBROUTINE. +
C +
C X - INPUT/OUTPUT REAL VECTOR OF DIMENSION N CONTAINING +
C AFTER THE (k+1)-TH CALL, THE SOLUTION x . +
C k+1 +
C BEFORE THE FIRST CALL, X MUST CONTAIN x . +
C 0 +
C +
C R - OUTPUT REAL VECTOR OF DIMENSION N CONTAINING, AFTER THE +
C (k+1)-TH CALL, THE RESIDUAL VECTOR r . +
C k+1 +
C IF, AT THE FIRST CALL, THE NORM OF THE VECTOR R = r IS +
C 0 +
C LESS THAN EPS THEN THE VECTOR R CONTAINS r IN OUTPUT. +
C 0 +
C Z - WORKING REAL VECTOR OF DIMENSION N. +
C +
C MAXBCK - INPUT INTEGER VALUE. IT REPRESENTS THE MAXIMUM VALUE +
C ALLOWED FOR THE JUMP m . IT MUST BE GREATER THAN ZERO. +
C k +
C MAXBCK IS USED TO CONTROL THE DIMENSION OF THE WORKING +
C ARRAYS. +
C +
C V - WORKING REAL MATRIX OF DIMENSION (N,0:MAXBCK). +
C +
C P - WORKING REAL VECTOR OF DIMENSION (0:MAXBCK). +
C +
C B - WORKING REAL VECTOR OF DIMENSION (0:MAXBCK). +
C +
C ALPHA - WORKING REAL VECTOR OF DIMENSION (0:MAXBCK). +
C +
C BETA - WORKING REAL VECTOR OF DIMENSION (0:MAXBCK-1). +
C +
C D - WORKING REAL VECTOR OF DIMENSION (0:MAXBCK-1). +
C +
C IR - INPUT ROWS DIMENSION EXACTLY AS SPECIFIED IN THE +
C DIMENSION STATEMENTS IN THE CALLING PROGRAM FOR +
C THE MATRICES A AND V. +
C +
C NK - OUTPUT INTEGER VALUE, DIMENSION OF THE INTERMEDIATE +
C KRYLOV SUBSPACE. +
C +
C EPS - INPUT REAL VALUE USED IN TESTS FOR ZERO. +
C IF DABS(X) .LT. EPS, THEN X IS CONSIDERED TO BE ZERO. +
C THE SUBROUTINE DOES NOT CONTROL WHETHER OR NOR EPS IS +
C A NEGATIVE OR ZERO REAL NUMBER. +
C +
C INIT - INPUT/OUTPUT INTEGER TO BE SET TO ZERO BEFORE THE +
C FIRST CALL OF THE SUBROUTINE. ITS VALUE IS CHANGED +
C TO 1 BY THE SUBROUTINE DURING THE FIRST CALL. +
C FOR A NEW APPLICATION OF THE METHOD INIT MUST SET +
C AGAIN TO ZERO. +
C +
C IFLAG - INPUT INTEGER. +
C IF IFLAG .EQ. 0 THEN THE AUXILIARY VECTOR y IS CHOSEN +
C TO COINCIDE WITH THE VECTOR z = r . +
C 0 0 +
C IF IFLAG .NE. 0 THEN THE USER MUST DEFINE y IN THE MAIN +
C PROGRAM, BEFORE THE FIRST CALL. +
C +
C IER - OUTPUT INDEX WARNING/ERROR. +
C IER = 100 CALL OF THE SUBROUTINE WITH A NON ZERO IER +
C VALUE. +
C IER = 200 THE NORM OF THE RESIDUAL VECTOR IS SMALLER +
C THAN EPS. THE EXACT SOLUTION HAS BEEN +
C OBTAINED. +
C IER = 300 INCURABLE BREAKDOWN. +
C IER = 400 SOLUTION NOT OBTAINED AFTER REACHING THE +
C DIMENSION N OF THE SYSTEM, DUE TO THE +
C PRECISION OF THE COMPUTER. +
C IER = 500 THE VALUE OF m EXCEEDS THE VALUE OF MAXBCK.+
C k +
C +
C EXTERNAL MODULES: +
C +
C - DOUBLE PRECISION FUNCTIONS EUNORM, PINNER, SSTRIA, VSUMMA +
C - SUBROUTINE MATVEC +
C +
C REMARKS: +
C +
C - ALL THE REAL ARGUMENTS PASSED FROM THE CALLING PROGRAM MUST BE +
C IN DOUBLE PRECISION. +
C - THE VARIABLE NK, THE VECTORS Y, X, R, AB, Z, P, B, ALPHA, BETA AND +
C D AND THE MATRIX V MUST NOT BE MODIFIED BY THE USER BETWEEN TWO +
C CONSECUTIVE CALLS OF THE SUBROUTINE. +
C - IN THE CALLING PROCEDURE, THE ROWS DIMENSION FOR THE ARRAYS A AND +
C V MUST BE THE SAME. +
C - IN THE CALLING PROCEDURE THE MINIMAL DIMENSIONS FOR THE +
C MATRIX AND VECTOR ARGUMENTS (IR=N) ARE: +
C A(N,N) +
C AB(N) +
C Y(N) +
C X(N) +
C R(N) +
C Z(N) +
C V(N,0:MAXBCK) +
C P(0:MAXBCK) +
C B(0:MAXBCK) +
C ALPHA(0:MAXBCK) +
C BETA(0:MAXBCK-1) +
C D(0:MAXBCK-1) +
C +
C AUTHORS: +
C +
C C. BREZINSKI AND H. SADOK +
C UNIVERSITY OF LILLE - FRANCE +
C M. REDIVO ZAGLIA +
C UNIVERSITY OF PADOVA - ITALY +
C +
C REFERENCES: +
C +
C - AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE ALGORITHMS, +
C NUMERICAL ALGORITHMS, 1(1991), PP. 261-284. +
C +
C - ADDENDUM TO "AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE +
C ALGORITHMS", NUMERICAL ALGORITHMS, 2(1992), PP. 133-136. +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE MRZ (N, A, AB, Y, X, R, Z, MAXBCK, V, P, B, ALPHA,
* BETA, D, IR, NK, EPS, INIT, IFLAG, IER)
C
C ... SPECIFICATIONS FOR ARGUMENTS
C
INTEGER IER, IFLAG, INIT, IR, MAXBCK, N, NK
DOUBLE PRECISION EPS
DOUBLE PRECISION A(IR,1), B(0:1), D(0:1), ALPHA(0:1), BETA(0:1),
* P(0:1), V(IR,0:1), AB(1), R(1), X(1), Y(1), Z(1)
C
C ... SPECIFICATIONS FOR EXTERNAL MODULES
C
DOUBLE PRECISION EUNORM, PINNER, SSTRIA, VSUMMA
C
C ... SPECIFICATIONS FOR LOCAL VARIABLES
C
INTEGER I, J, IND, MK, MKSAV
DOUBLE PRECISION CREC, TMP
SAVE IND, MKSAV
C
C ... EXECUTABLE STATEMENTS
C
C
C ... FIRST CALL OF THE SUBROUTINE
C
IF (INIT .EQ. 0) THEN
IER = 0
INIT = 1
NK = 0
IND = -1
MKSAV = 0
C
C ... COMPUTATION OF THE VECTORS r AND z
C 0 0
C
CALL MATVEC (N, A, X, 1, R, 1, IR, 0)
DO 10 I = 1, N
Z(I) = 0.0D0
R(I) = R(I) - AB(I)
V(I,0) = R(I)
C
C ... TEST FOR y = z
C 0
IF (IFLAG .EQ. 0) Y(I) = R(I)
10 CONTINUE
C
C ... TEST FOR THE SOLUTION
C
IF (DSQRT(EUNORM(N,R)) .LT. EPS) THEN
IER = 200
RETURN
END IF
END IF
C
C ... FIRST AND NEXT CALLS OF THE SUBROUTINE
C
IF (IER .NE. 0) THEN
IER = 100
RETURN
END IF
C
C ... IND REPRESENTS THE (NUMBER_OF_CALLS - 1)
C
IND = IND + 1
C
C ... COMPUTATION OF m
C k
C
MK = 1
CALL MATVEC (N, A, V, 1, V, 2, IR, 0)
TMP = PINNER(N, Y, V, 2, IR)
IF (DABS(TMP) .LE. EPS) THEN
IF (NK .EQ. N-1) THEN
NK = N
IER = 300
RETURN
END IF
DO 20 MK = 2, N-NK
C
C ... CONTROL THE VALUE OF m
C k
C
IF (MK .GT. MAXBCK) THEN
NK = NK + MK - 1
IER = 500
RETURN
END IF
C
C ... COMPUTATION OF TMP=(y,A**(n + m )*z )
C k k k
C
CALL MATVEC (N, A, V, MK, V, MK+1, IR, 0)
TMP = PINNER(N, Y, V, MK+1, IR)
IF (DABS(TMP) .GT. EPS) GO TO 30
20 CONTINUE
NK = N
IER = 300
RETURN
END IF
30 CONTINUE
C
C ... THE VALUE OF m HAS BEEN OBTAINED
C k
C
C ... COMPUTATION OF D(0)=(y,A**(n )*r )
C k k
C
D(0) = PINNER(N, Y, R, 1, IR)
C
C ... STORAGE OF B(0)=(y,A**(n + m )*z )
C k k k
C
B(0) = TMP
C
C ... COMPUTATION OF BETA FOR i = m - 1, ..., 0
C i k
C
BETA(MK-1) = D(0) / B(0)
DO 50 I = 1, MK
CALL MATVEC (N, A, Y, 1, AB, 1, IR, 1)
DO 40 J = 1, N
Y(J) = AB(J)
40 CONTINUE
B(I) = PINNER(N, Y, V, MK+1, IR)
IF (I. NE. MK) THEN
D(I) = PINNER(N, Y, R, 1, IR)
BETA(MK-I-1) = SSTRIA(BETA, B, D(I), I, MK, 1)
END IF
IF (IND .NE. 0 .AND. I .GT. MKSAV) THEN
P(I) = PINNER(N, Y, Z, 1, IR)
END IF
50 CONTINUE
C
C ... COMPUTATION OF THE NEXT SOLUTION VECTOR x
C k+1
C AND OF THE NEXT RESIDUAL VECTOR r
C k+1
C
DO 60 J = 1, N
X(J) = X(J) - VSUMMA(BETA, V, IR, J, MK, 0)
R(J) = R(J) - VSUMMA(BETA, V, IR, J, MK, 1)
60 CONTINUE
C
C ... COMPUTE THE VALUE OF n
C k+1
C
NK = NK + MK
C
C ... TEST FOR THE SOLUTION
C
IF (DSQRT(EUNORM(N,R)) .LT. EPS) THEN
IER = 200
RETURN
END IF
C
C ... CONTROL THE VALUE OF n
C k
C
IF (NK .EQ. N) THEN
IER = 400
RETURN
END IF
C
C ... COMPUTATION OF C
C k+1
C
IF (IND .NE. 0) CREC = B(0) / P(0)
C
C ... COMPUTATION OF ALPHA FOR i = m - 1, ..., 0
C i k
C
ALPHA(MK) = 1.0D0
DO 70 I = 1, MK
IF (IND .NE. 0) THEN
ALPHA(MK-I) = SSTRIA(ALPHA, B, CREC * P(I), I, MK, 0)
ELSE
ALPHA(MK-I) = SSTRIA(ALPHA, B, 0.0D0, I, MK, 0)
END IF
70 CONTINUE
C
C ... SAVE THE VALUES OF B IN P
C
DO 80 I = 0, MK
P(I) = B(I)
80 CONTINUE
C
C ... COMPUTATION OF THE VECTOR z
C k+1
C
DO 90 J = 1, N
TMP = V(J,0)
V(J,0) = VSUMMA(ALPHA, V, IR, J, MK, 2) - CREC * Z(J)
Z(J) = TMP
90 CONTINUE
C
C ... SAVE THE VALUE OF m
C k
C
MKSAV = MK
RETURN
END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C SUBROUTINE NAME - BMRZ +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C DESCRIPTION: +
C +
C COMPUTES THE SOLUTION OF A SYSTEM OF LINEAR EQUATIONS BY +
C THE BMRZ ALGORITHM. +
C THIS ALGORITHM IS A VARIANT OF THE MRZ AND IT CAN BE USED +
C ONLY IF ( y, A**n * r ) IS NOT EQUAL TO ZERO. +
C k k +
C +
C USAGE: +
C CALL BMRZ (N, A, AB, Y, X, R, YSAV, MAXBCK, V, B, BETA, D, +
C IR, NK, EPS, INIT, IFLAG, IER) +
C +
C ARGUMENTS: +
C +
C N - INPUT INTEGER VALUE, DIMENSION OF THE SYSTEM. +
C +
C A - INPUT REAL MATRIX OF DIMENSION (N,N) CONTAINING THE +
C MATRIX OF THE SYSTEM. +
C +
C AB - INPUT/WORKING REAL VECTOR OF DIMENSION N CONTAINING, +
C BEFORE THE FIRST CALL, THE RIGHT HAND SIDE OF THE SYSTEM. +
C THE VECTOR IS USED AS A WORKING AREA IN THE NEXT CALLS +
C AND THEN THE INPUT VALUES ARE ALWAYS DESTROYED BY THE +
C SUBROUTINE. +
C +
C Y - INPUT/WORKING REAL VECTOR OF DIMENSION N CONTAINING, +
C BEFORE THE FIRST CALL, THE AUXILIARY VECTOR y. +
C IF IFLAG .EQ. 0 THE VECTOR Y IS SET TO r = A x - b +
C 0 0 +
C BY THE SUBROUTINE DURING THE FIRST CALL. THE VECTOR IS +
C USED AS A WORKING AREA IN THE NEXT CALLS AND THEN THE +
C INPUT VALUES ARE ALWAYS DESTROYED BY THE SUBROUTINE. +
C +
C X - INPUT/OUTPUT REAL VECTOR OF DIMENSION N CONTAINING +
C AFTER THE (k+1)-TH CALL, THE SOLUTION x . +
C k+1 +
C BEFORE THE FIRST CALL, X MUST CONTAIN x . +
C 0 +
C +
C R - OUTPUT REAL VECTOR OF DIMENSION N CONTAINING, AFTER THE +
C (k+1)-TH CALL, THE RESIDUAL VECTOR r . +
C k+1 +
C IF, AT THE FIRST CALL, THE NORM OF THE VECTOR R = r IS +
C 0 +
C LESS THAN EPS THEN THE VECTOR R CONTAINS r IN OUTPUT. +
C 0 +
C YSAV - WORKING REAL VECTOR OF DIMENSION N. +
C +
C MAXBCK - INPUT INTEGER VALUE. IT REPRESENTS THE MAXIMUM VALUE +
C ALLOWED FOR THE JUMP m . IT MUST BE GREATER THAN ZERO. +
C k +
C MAXBCK IS USED TO CONTROL THE DIMENSION OF THE WORKING +
C ARRAYS. +
C +
C V - WORKING REAL MATRIX OF DIMENSION (N,0:MAXBCK). +
C +
C B - WORKING REAL VECTOR OF DIMENSION (0:MAXBCK-1). +
C +
C BETA - WORKING REAL VECTOR OF DIMENSION (0:MAXBCK-1). +
C +
C D - WORKING REAL VECTOR OF DIMENSION (0:MAXBCK-1). +
C +
C IR - INPUT ROWS DIMENSION EXACTLY AS SPECIFIED IN THE +
C DIMENSION STATEMENTS IN THE CALLING PROGRAM FOR +
C THE MATRICES A AND V. +
C +
C NK - OUTPUT INTEGER VALUE, DIMENSION OF THE INTERMEDIATE +
C KRYLOV SUBSPACE. +
C +
C EPS - INPUT REAL VALUE USED IN TESTS FOR ZERO. +
C IF DABS(X) .LT. EPS, THEN X IS CONSIDERED TO BE ZERO. +
C THE SUBROUTINE DOES NOT CONTROL WHETHER OR NOR EPS IS +
C A NEGATIVE OR ZERO REAL NUMBER. +
C +
C INIT - INPUT/OUTPUT INTEGER TO BE SET TO ZERO BEFORE THE +
C FIRST CALL OF THE SUBROUTINE. ITS VALUE IS CHANGED +
C TO 1 BY THE SUBROUTINE DURING THE FIRST CALL. +
C FOR A NEW APPLICATION OF THE METHOD INIT MUST SET +
C AGAIN TO ZERO. +
C +
C IFLAG - INPUT INTEGER. +
C IF IFLAG .EQ. 0 THEN THE AUXILIARY VECTOR y IS CHOSEN +
C TO COINCIDE WITH THE VECTOR z = r . +
C 0 0 +
C IF IFLAG .NE. 0 THEN THE USER MUST DEFINE y IN THE MAIN +
C PROGRAM, BEFORE THE FIRST CALL. +
C +
C IER - OUTPUT INDEX WARNING/ERROR. +
C IER = 100 CALL OF THE SUBROUTINE WITH A NON ZERO IER +
C VALUE. +
C IER = 200 THE NORM OF THE RESIDUAL VECTOR IS SMALLER +
C THAN EPS. THE EXACT SOLUTION HAS BEEN +
C OBTAINED. +
C IER = 300 INCURABLE BREAKDOWN. +
C IER = 400 SOLUTION NOT OBTAINED AFTER REACHING THE +
C DIMENSION N OF THE SYSTEM, DUE TO THE +
C PRECISION OF THE COMPUTER. +
C IER = 500 THE VALUE OF m EXCEEDS THE VALUE OF MAXBCK.+
C k +
C IER = 600 IT IS IMPOSSIBLE TO USE THE BMRZ. PLEASE +
C USE THE MRZ ALGORITHM. +
C +
C EXTERNAL MODULES: +
C +
C - DOUBLE PRECISION FUNCTIONS EUNORM, PINNER, SSTRIA, VSUMMA +
C - SUBROUTINE MATVEC +
C +
C REMARKS: +
C +
C - ALL THE REAL ARGUMENTS PASSED FROM THE CALLING PROGRAM MUST BE +
C IN DOUBLE PRECISION. +
C - THE VARIABLE NK, THE VECTORS Y, X, R, AB, YSAV, B, BETA AND D +
C AND THE MATRIX V MUST NOT BE MODIFIED BY THE USER BETWEEN TWO +
C CONSECUTIVE CALLS OF THE SUBROUTINE. +
C - IN THE CALLING PROCEDURE, THE ROWS DIMENSION FOR THE ARRAYS A AND +
C V MUST BE THE SAME. +
C - IN THE CALLING PROCEDURE THE MINIMAL DIMENSIONS FOR THE +
C MATRIX AND VECTOR ARGUMENTS (IR=N) ARE: +
C A(N,N) +
C AB(N) +
C Y(N) +
C X(N) +
C R(N) +
C YSAV(N) +
C V(N,0:MAXBCK) +
C B(0:MAXBCK-1) +
C BETA(0:MAXBCK-1) +
C D(0:MAXBCK-1) +
C +
C AUTHORS: +
C +
C C. BREZINSKI AND H. SADOK +
C UNIVERSITY OF LILLE - FRANCE +
C M. REDIVO ZAGLIA +
C UNIVERSITY OF PADOVA - ITALY +
C +
C REFERENCES: +
C +
C - AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE ALGORITHMS, +
C NUMERICAL ALGORITHMS, 1(1991), PP. 261-284. +
C +
C - ADDENDUM TO "AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE +
C ALGORITHMS", NUMERICAL ALGORITHMS, 2(1992), PP. 133-136. +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE BMRZ (N, A, AB, Y, X, R, YSAV, MAXBCK, V, B, BETA,
* D, IR, NK, EPS, INIT, IFLAG, IER)
C
C ... SPECIFICATIONS FOR ARGUMENTS
C
INTEGER IER, IFLAG, INIT, IR, MAXBCK, N, NK
DOUBLE PRECISION EPS
DOUBLE PRECISION A(IR,1), B(0:1), D(0:1), BETA(0:1), V(IR,0:1),
* YSAV(1), AB(1), R(1), X(1), Y(1)
C
C ... SPECIFICATIONS FOR EXTERNAL MODULES
C
DOUBLE PRECISION EUNORM, PINNER, SSTRIA, VSUMMA
C
C ... SPECIFICATIONS FOR LOCAL VARIABLES
C
INTEGER I, J, MK
DOUBLE PRECISION TMP
C
C ... EXECUTABLE STATEMENTS
C
C
C ... FIRST CALL OF THE SUBROUTINE
C
IF (INIT .EQ. 0) THEN
IER = 0
INIT = 1
NK = 0
C
C ... COMPUTATION OF THE VECTORS r AND z
C 0 0
C
CALL MATVEC (N, A, X, 1, R, 1, IR, 0)
DO 10 I = 1, N
R(I) = R(I) - AB(I)
V(I,0) = R(I)
C
C ... TEST FOR y = z
C 0
IF (IFLAG .EQ. 0) Y(I) = R(I)
10 CONTINUE
C
C ... TEST FOR THE SOLUTION
C
IF (DSQRT(EUNORM(N,R)) .LT. EPS) THEN
IER = 200
RETURN
END IF
END IF
C
C ... FIRST AND NEXT CALLS OF THE SUBROUTINE
C
IF (IER .NE. 0) THEN
IER = 100
RETURN
END IF
C
C ... COMPUTATION OF m
C k
C
MK = 1
CALL MATVEC (N, A, V, 1, V, 2, IR, 0)
TMP = PINNER(N, Y, V, 2, IR)
IF (DABS(TMP) .LE. EPS) THEN
IF (NK .EQ. N-1) THEN
NK = N
IER = 300
RETURN
END IF
DO 20 MK = 2, N-NK
C
C ... CONTROL THE VALUE OF m
C k
C
IF (MK .GT. MAXBCK) THEN
NK = NK + MK - 1
IER = 500
RETURN
END IF
C
C ... COMPUTATION OF TMP=(y,A**(n + m )*z )
C k k k
C
CALL MATVEC (N, A, V, MK, V, MK+1, IR, 0)
TMP = PINNER(N, Y, V, MK+1, IR)
IF (DABS(TMP) .GT. EPS) GO TO 30
20 CONTINUE
NK = N
IER = 300
RETURN
END IF
30 CONTINUE
C
C ... THE VALUE OF m HAS BEEN OBTAINED
C k
C
C ... COMPUTATION OF D(0)=(y,A**(n )*r )
C k k
C
D(0) = PINNER(N, Y, R, 1, IR)
C
C ... CONTROL IF THE BMRZ CAN BE USED
C
IF (DABS(D(0)) .LT. EPS) THEN
NK = NK + MK
IER = 600
RETURN
END IF
C
C ... SAVE FOR THE VECTOR y
C
DO 40 J = 1, N
YSAV(J) = Y(J)
40 CONTINUE
C
C ... STORAGE OF B(0)=(y,A**(n + m )*z )
C k k k
C
B(0) = TMP
C
C ... COMPUTATION OF BETA FOR i = m - 1, ..., 0
C i k
C
BETA(MK-1) = D(0) / B(0)
DO 60 I = 1, MK
CALL MATVEC (N, A, Y, 1, AB, 1, IR, 1)
DO 50 J = 1, N
Y(J) = AB(J)
50 CONTINUE
IF (I. NE. MK) THEN
B(I) = PINNER(N, Y, V, MK+1, IR)
D(I) = PINNER(N, Y, R, 1, IR)
BETA(MK-I-1) = SSTRIA(BETA, B, D(I), I, MK, 1)
END IF
60 CONTINUE
C
C ... COMPUTATION OF THE NEXT SOLUTION VECTOR x
C k+1
C AND OF THE NEXT RESIDUAL VECTOR r
C k+1
C
DO 70 J = 1, N
X(J) = X(J) - VSUMMA(BETA, V, IR, J, MK, 0)
R(J) = R(J) - VSUMMA(BETA, V, IR, J, MK, 1)
70 CONTINUE
C
C ... COMPUTE THE VALUE OF n
C k+1
C
NK = NK + MK
C
C ... TEST FOR THE SOLUTION
C
IF (DSQRT(EUNORM(N,R)) .LT. EPS) THEN
IER = 200
RETURN
END IF
C
C ... CONTROL THE VALUE OF n
C k
C
IF (NK .EQ. N) THEN
IER = 400
RETURN
END IF
C
C ... COMPUTATION OF TMP=(y,A**(n + m )*r )
C k k k+1
C
DO 90 I = 1, MK
CALL MATVEC (N, A, YSAV, 1, AB, 1, IR, 1)
DO 80 J = 1, N
YSAV(J) = AB(J)
80 CONTINUE
90 CONTINUE
TMP = PINNER(N, YSAV, R, 1, IR)
C
C ... COMPUTATION OF THE VECTOR z
C k+1
C
DO 100 J = 1, N
V(J,0) = - R(J)/BETA(MK-1) + TMP/D(0) * V(J,0)
100 CONTINUE
C
RETURN
END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C SUBROUTINE NAME - SMRZ +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C DESCRIPTION: +
C +
C COMPUTES THE SOLUTION OF A SYSTEM OF LINEAR EQUATIONS BY +
C THE SMRZ ALGORITHM. +
C THIS ALGORITHM IS A VARIANT OF THE MRZ AND IT CAN BE USED +
C ONLY IF ( y, A**n * r ) IS NOT EQUAL TO ZERO. +
C k k +
C +
C USAGE: +
C CALL SMRZ (N, A, AB, Y, X, R, RSAV, MAXBCK, V, B, DCOEF, BETA, +
C D, IR, NK, EPS, INIT, IFLAG, IER) +
C +
C ARGUMENTS: +
C +
C N - INPUT INTEGER VALUE, DIMENSION OF THE SYSTEM. +
C +
C A - INPUT REAL MATRIX OF DIMENSION (N,N) CONTAINING THE +
C MATRIX OF THE SYSTEM. +
C +
C AB - INPUT/WORKING REAL VECTOR OF DIMENSION N CONTAINING, +
C BEFORE THE FIRST CALL, THE RIGHT HAND SIDE OF THE SYSTEM. +
C THE VECTOR IS USED AS A WORKING AREA IN THE NEXT CALLS +
C AND THEN THE INPUT VALUES ARE ALWAYS DESTROYED BY THE +
C SUBROUTINE. +
C +
C Y - INPUT/WORKING REAL VECTOR OF DIMENSION N CONTAINING, +
C BEFORE THE FIRST CALL, THE AUXILIARY VECTOR y. +
C IF IFLAG .EQ. 0 THE VECTOR Y IS SET TO r = A x - b +
C 0 0 +
C BY THE SUBROUTINE DURING THE FIRST CALL. THE VECTOR IS +
C USED AS A WORKING AREA IN THE NEXT CALLS AND THEN THE +
C INPUT VALUES ARE ALWAYS DESTROYED BY THE SUBROUTINE. +
C +
C X - INPUT/OUTPUT REAL VECTOR OF DIMENSION N CONTAINING +
C AFTER THE (k+1)-TH CALL, THE SOLUTION x . +
C k+1 +
C BEFORE THE FIRST CALL, X MUST CONTAIN x . +
C 0 +
C +
C R - OUTPUT REAL VECTOR OF DIMENSION N CONTAINING, AFTER THE +
C (k+1)-TH CALL, THE RESIDUAL VECTOR r . +
C k+1 +
C IF, AT THE FIRST CALL, THE NORM OF THE VECTOR R = r IS +
C 0 +
C LESS THAN EPS THEN THE VECTOR R CONTAINS r IN OUTPUT. +
C 0 +
C RSAV - WORKING REAL VECTOR OF DIMENSION N. +
C +
C MAXBCK - INPUT INTEGER VALUE. IT REPRESENTS THE MAXIMUM VALUE +
C ALLOWED FOR THE JUMP m . IT MUST BE GREATER THAN ZERO. +
C k +
C MAXBCK IS USED TO CONTROL THE DIMENSION OF THE WORKING +
C ARRAYS. +
C +
C V - WORKING REAL MATRIX OF DIMENSION (N,0:MAXBCK). +
C +
C B - WORKING REAL VECTOR OF DIMENSION (0:MAXBCK). +
C +
C DCOEF - WORKING REAL VECTOR OF DIMENSION (0:MAXBCK). +
C +
C BETA - WORKING REAL VECTOR OF DIMENSION (0:MAXBCK-1). +
C +
C D - WORKING REAL VECTOR OF DIMENSION (0:MAXBCK). +
C +
C IR - INPUT ROWS DIMENSION EXACTLY AS SPECIFIED IN THE +
C DIMENSION STATEMENTS IN THE CALLING PROGRAM FOR +
C THE MATRICES A AND V. +
C +
C NK - OUTPUT INTEGER VALUE, DIMENSION OF THE INTERMEDIATE +
C KRYLOV SUBSPACE. +
C +
C EPS - INPUT REAL VALUE USED IN TESTS FOR ZERO. +
C IF DABS(X) .LT. EPS, THEN X IS CONSIDERED TO BE ZERO. +
C THE SUBROUTINE DOES NOT CONTROL WHETHER OR NOR EPS IS +
C A NEGATIVE OR ZERO REAL NUMBER. +
C +
C INIT - INPUT/OUTPUT INTEGER TO BE SET TO ZERO BEFORE THE +
C FIRST CALL OF THE SUBROUTINE. ITS VALUE IS CHANGED +
C TO 1 BY THE SUBROUTINE DURING THE FIRST CALL. +
C FOR A NEW APPLICATION OF THE METHOD INIT MUST SET +
C AGAIN TO ZERO. +
C +
C IFLAG - INPUT INTEGER. +
C IF IFLAG .EQ. 0 THEN THE AUXILIARY VECTOR y IS CHOSEN +
C TO COINCIDE WITH THE VECTOR z = r . +
C 0 0 +
C IF IFLAG .NE. 0 THEN THE USER MUST DEFINE y IN THE MAIN +
C PROGRAM, BEFORE THE FIRST CALL. +
C +
C IER - OUTPUT INDEX WARNING/ERROR. +
C IER = 100 CALL OF THE SUBROUTINE WITH A NON ZERO IER +
C VALUE. +
C IER = 200 THE NORM OF THE RESIDUAL VECTOR IS SMALLER +
C THAN EPS. THE EXACT SOLUTION HAS BEEN +
C OBTAINED. +
C IER = 300 INCURABLE BREAKDOWN. +
C IER = 400 SOLUTION NOT OBTAINED AFTER REACHING THE +
C DIMENSION N OF THE SYSTEM, DUE TO THE +
C PRECISION OF THE COMPUTER. +
C IER = 500 THE VALUE OF m EXCEEDS THE VALUE OF MAXBCK.+
C k +
C IER = 600 IT IS IMPOSSIBLE TO USE THE SMRZ. PLEASE +
C USE THE MRZ ALGORITHM. +
C +
C EXTERNAL MODULES: +
C +
C - DOUBLE PRECISION FUNCTIONS EUNORM, PINNER, SSTRIA, VSUMMA +
C - SUBROUTINE MATVEC +
C +
C REMARKS: +
C +
C - ALL THE REAL ARGUMENTS PASSED FROM THE CALLING PROGRAM MUST BE +
C IN DOUBLE PRECISION. +
C - THE VARIABLE NK, THE VECTORS Y, X, R, AB, RSAV, B, DCOEF, BETA AND +
C D AND THE MATRIX V MUST NOT BE MODIFIED BY THE USER BETWEEN TWO +
C CONSECUTIVE CALLS OF THE SUBROUTINE. +
C - IN THE CALLING PROCEDURE, THE ROWS DIMENSION FOR THE ARRAYS A AND +
C V MUST BE THE SAME. +
C - IN THE CALLING PROCEDURE THE MINIMAL DIMENSIONS FOR THE +
C MATRIX AND VECTOR ARGUMENTS (IR=N) ARE: +
C A(N,N) +
C AB(N) +
C Y(N) +
C X(N) +
C R(N) +
C RSAV(N) +
C V(N,0:MAXBCK) +
C B(0:MAXBCK) +
C DCOEF(0:MAXBCK) +
C BETA(0:MAXBCK-1) +
C D(0:MAXBCK) +
C +
C AUTHORS: +
C +
C C. BREZINSKI AND H. SADOK +
C UNIVERSITY OF LILLE - FRANCE +
C M. REDIVO ZAGLIA +
C UNIVERSITY OF PADOVA - ITALY +
C +
C REFERENCES: +
C +
C - AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE ALGORITHMS, +
C NUMERICAL ALGORITHMS, 1(1991), PP. 261-284. +
C +
C - ADDENDUM TO "AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE +
C ALGORITHMS", NUMERICAL ALGORITHMS, 2(1992), PP. 133-136. +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE SMRZ (N, A, AB, Y, X, R, RSAV, MAXBCK, V, B, DCOEF,
* BETA, D, IR, NK, EPS, INIT, IFLAG, IER)
C
C ... SPECIFICATIONS FOR ARGUMENTS
C
INTEGER IER, IFLAG, INIT, IR, MAXBCK, N, NK
DOUBLE PRECISION EPS
DOUBLE PRECISION A(IR,1), B(0:1), BETA(0:1), D(0:1), DCOEF(0:1),
* V(IR,0:1), AB(1), R(1), RSAV(1), X(1), Y(1)
C
C ... SPECIFICATIONS FOR EXTERNAL MODULES
C
DOUBLE PRECISION EUNORM, PINNER, SSTRIA, VSUMMA
C
C ... SPECIFICATIONS FOR LOCAL VARIABLES
C
INTEGER I, J, MK
DOUBLE PRECISION DREC, TMP
C
C ... EXECUTABLE STATEMENTS
C
C
C ... FIRST CALL OF THE SUBROUTINE
C
IF (INIT .EQ. 0) THEN
IER = 0
INIT = 1
NK = 0
C
C ... COMPUTATION OF THE VECTORS r AND z
C 0 0
C
CALL MATVEC (N, A, X, 1, R, 1, IR, 0)
DO 10 I = 1, N
R(I) = R(I) - AB(I)
V(I,0) = R(I)
C
C ... TEST FOR y = z
C 0
C
IF (IFLAG .EQ. 0) Y(I) = R(I)
10 CONTINUE
C
C ... TEST FOR THE SOLUTION
C
IF (DSQRT(EUNORM(N,R)) .LT. EPS) THEN
IER = 200
RETURN
END IF
END IF
C
C ... FIRST AND NEXT CALLS OF THE SUBROUTINE
C
IF (IER .NE. 0) THEN
IER = 100
RETURN
END IF
C
C ... COMPUTATION OF m
C k
C
MK = 1
CALL MATVEC (N, A, V, 1, V, 2, IR, 0)
TMP = PINNER(N, Y, V, 2, IR)
IF (DABS(TMP) .LE. EPS) THEN
IF (NK .EQ. N-1) THEN
NK = N
IER = 300
RETURN
END IF
DO 20 MK = 2, N-NK
C
C ... CONTROL THE VALUE OF m
C k
C
IF (MK .GT. MAXBCK) THEN
NK = NK + MK - 1
IER = 500
RETURN
END IF
C
C ... COMPUTATION OF TMP=(y,A**(n + m )*z )
C k k k
C
CALL MATVEC (N, A, V, MK, V, MK+1, IR, 0)
TMP = PINNER(N, Y, V, MK+1, IR)
IF (DABS(TMP) .GT. EPS) GO TO 30
20 CONTINUE
NK = N
IER = 300
RETURN
END IF
30 CONTINUE
C
C ... THE VALUE OF m HAS BEEN OBTAINED
C k
C
C ... COMPUTATION OF D(0)=(y,A**(n )*r )
C k k
C
D(0) = PINNER(N, Y, R, 1, IR)
C
C ... CONTROL IF THE SMRZ CAN BE USED
C
IF (DABS(D(0)) .LT. EPS) THEN
NK = NK + MK
IER = 600
RETURN
END IF
C
C ... STORAGE OF B(0)=(y,A**(n + m )*z )
C k k k
C
B(0) = TMP
C
C ... COMPUTATION OF BETA FOR i = m - 1, ..., 0
C i k
C
BETA(MK-1) = D(0) / B(0)
DO 50 I = 1, MK
CALL MATVEC (N, A, Y, 1, AB, 1, IR, 1)
DO 40 J = 1, N
Y(J) = AB(J)
40 CONTINUE
B(I) = PINNER(N, Y, V, MK+1, IR)
D(I) = PINNER(N, Y, R, 1, IR)
IF (I. NE. MK) BETA(MK-I-1) = SSTRIA(BETA, B, D(I), I, MK, 1)
50 CONTINUE
C
C ... COMPUTATION OF THE NEXT SOLUTION VECTOR x
C k+1
C AND OF THE NEXT RESIDUAL VECTOR r
C k+1
C
DO 60 J = 1, N
X(J) = X(J) - VSUMMA(BETA, V, IR, J, MK, 0)
C
C ... SAVE THE VALUE OF r
C k
RSAV(J) = R(J)
C
R(J) = R(J) - VSUMMA(BETA, V, IR, J, MK, 1)
60 CONTINUE
C
C ... COMPUTE THE VALUE OF n
C k+1
C
NK = NK + MK
C
C ... TEST FOR THE SOLUTION
C
IF (DSQRT(EUNORM(N,R)) .LT. EPS) THEN
IER = 200
RETURN
END IF
C
C ... CONTROL THE VALUE OF n
C k
C
IF (NK .EQ. N) THEN
IER = 400
RETURN
END IF
C
C ... COMPUTATION OF D
C k+1
C
DREC = B(0) / D(0)
C
C ... COMPUTATION OF DCOEF FOR i = m - 1, ..., 0
C i k
C
DCOEF(MK) = 1.0D0
DO 70 I = 1, MK
DCOEF(MK-I) = SSTRIA(DCOEF, B, DREC * D(I), I, MK, 0)
70 CONTINUE
C
C ... COMPUTATION OF THE VECTOR z
C k+1
C
DO 80 J = 1, N
V(J,0) = VSUMMA(DCOEF, V, IR, J, MK, 2) - DREC * RSAV(J)
80 CONTINUE
C
RETURN
END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C SUBROUTINE NAME - BSMRZ +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C DESCRIPTION: +
C +
C COMPUTES THE SOLUTION OF A SYSTEM OF LINEAR EQUATIONS BY +
C THE BSMRZ ALGORITHM. +
C THIS ALGORITHM CAN BE USED ONLY IF (y, A**n * r ) IS DIFFERENT +
C FROM ZERO. k k +
C +
C USAGE: +
C CALL BSMRZ (N, A, AB, Y, X, R, YSAV, MAXBCK, MAXN, V, W, LMAT, +
C RMAT, D, C, IR, NK, EPS, EPS1, EPS2, INIT, IFLAG, IER) +
C +
C ARGUMENTS: +
C +
C N - INPUT INTEGER VALUE, DIMENSION OF THE SYSTEM. +
C +
C A - INPUT REAL MATRIX OF DIMENSION (N,N) CONTAINING THE +
C MATRIX OF THE SYSTEM. +
C +
C AB - INPUT/WORKING REAL VECTOR OF DIMENSION N CONTAINING, +
C BEFORE THE FIRST CALL, THE RIGHT HAND SIDE OF THE SYSTEM. +
C THE VECTOR IS USED AS A WORKING AREA IN THE NEXT CALLS +
C AND THEN THE INPUT VALUES ARE ALWAYS DESTROYED BY THE +
C SUBROUTINE. +
C +
C Y - INPUT/WORKING REAL VECTOR OF DIMENSION N CONTAINING, +
C BEFORE THE FIRST CALL, THE AUXILIARY VECTOR y. +
C IF IFLAG .EQ. 0 THE VECTOR Y IS SET TO r = A x - b +
C 0 0 +
C BY THE SUBROUTINE DURING THE FIRST CALL. THE VECTOR IS +
C USED AS A WORKING AREA IN THE NEXT CALLS AND THEN THE +
C INPUT VALUES ARE ALWAYS DESTROYED BY THE SUBROUTINE. +
C +
C X - INPUT/OUTPUT REAL VECTOR OF DIMENSION N CONTAINING +
C AFTER THE (k+1)-TH CALL, THE SOLUTION x . +
C k+1 +
C BEFORE THE FIRST CALL, X MUST CONTAIN x . +
C 0 +
C +
C R - OUTPUT REAL VECTOR OF DIMENSION N CONTAINING, AFTER THE +
C (k+1)-TH CALL, THE RESIDUAL VECTOR r . +
C k+1 +
C IF, AT THE FIRST CALL, THE NORM OF THE VECTOR R = r IS +
C 0 +
C LESS THAN EPS2 THEN THE VECTOR R CONTAINS r IN OUTPUT. +
C 0 +
C YSAV - WORKING REAL VECTOR OF DIMENSION N. +
C +
C MAXBCK - INPUT INTEGER VALUE. IT REPRESENTS THE MAXIMUM VALUE +
C ALLOWED FOR THE JUMP m . IT MUST GREATER THAN ZERO. +
C k +
C MAXBCK IS USED TO CONTROL THE DIMENSION OF THE WORKING +
C ARRAYS. +
C +
C MAXN - INPUT/OUTPUT INTEGER VALUE. IT REPRESENTS THE MAXIMUM +
C VALUE ALLOWED, GREATER OR EQUAL TO N, FOR n + m . +
C k k +
C IF MAXN IS LESS THAN N, IT IS SET EQUAL TO N AT THE FIRST +
C CALL BECAUSE THIS VALUE MUST BE AT LEAST EQUAL TO N. +
C +
C V - WORKING REAL MATRIX OF DIMENSION (N,0:MAXBCK). +
C +
C W - WORKING REAL MATRIX OF DIMENSION (N,0:MAXBCK-1). +
C +
C LMAT - WORKING REAL VECTOR OF DIMENSION (2*MAXBCK-1)*(2*MAXBCK-1)+
C +
C RMAT - WORKING REAL MATRIX OF DIMENSION (2,0:2*MAXBCK-1). +
C +
C D - WORKING REAL VECTOR OF DIMENSION (0:2*MAXBCK-1). +
C +
C C - WORKING REAL VECTOR OF DIMENSION (0:2*MAXBCK-1). +
C +
C IR - INPUT ROWS DIMENSION EXACTLY AS SPECIFIED IN THE +
C DIMENSION STATEMENTS IN THE CALLING PROGRAM FOR +
C THE MATRICES A, V AND W. +
C +
C NK - OUTPUT INTEGER VALUE, DIMENSION OF THE INTERMEDIATE +
C KRYLOV SUBSPACE. +
C +
C EPS - INPUT REAL VALUE USED FOR TESTING THE NEAR BREAKDOWN. +
C THE SUBROUTINE DOES NOT CONTROL WHETHER OR NOR EPS IS +
C A NEGATIVE OR ZERO REAL NUMBER. +
C +
C EPS1 - INPUT REAL VALUE USED FOR TESTING THE PIVOTS IN GAUSSIAN +
C ELIMINATION FOR SOLVING THE AUXILIARY SYSTEMS. +
C IF DABS(X) .LT. EPS1, THEN X IS CONSIDERED TO BE ZERO. +
C THE SUBROUTINE DOES NOT CONTROL WHETHER OR NOR EPS1 IS +
C A NEGATIVE OR ZERO REAL NUMBER. +
C +
C EPS2 - INPUT REAL VALUE USED FOR TESTING THE FINAL SOLUTION. +
C THE SUBROUTINE DOES NOT CONTROL WHETHER OR NOR EPS2 IS +
C A NEGATIVE OR ZERO REAL NUMBER. +
C +
C INIT - INPUT/OUTPUT INTEGER TO BE SET TO ZERO BEFORE THE +
C FIRST CALL OF THE SUBROUTINE. ITS VALUE IS CHANGED +
C TO 1 BY THE SUBROUTINE DURING THE FIRST CALL. +
C FOR A NEW APPLICATION OF THE METHOD INIT MUST SET +
C AGAIN TO ZERO. +
C +
C IFLAG - INPUT INTEGER. +
C IF IFLAG .EQ. 0 THEN THE AUXILIARY VECTOR y IS CHOSEN +
C TO COINCIDE WITH THE VECTOR z = r . +
C 0 0 +
C IF IFLAG .NE. 0 THEN THE USER MUST DEFINE y IN THE MAIN +
C PROGRAM, BEFORE THE FIRST CALL. +
C +
C IER - OUTPUT INDEX WARNING/ERROR. +
C IER = 100 CALL OF THE SUBROUTINE WITH A NON ZERO IER +
C VALUE. +
C IER = 200 THE NORM OF THE RESIDUAL VECTOR IS SMALLER +
C THAN EPS. THE EXACT SOLUTION HAS BEEN +
C OBTAINED. +
C IER = 300 SOLUTION NOT OBTAINED AFTER REACHING THE +
C VALUE NK+MK=MAXN. +
C IER = 400 SOLUTION NOT OBTAINED AFTER REACHING THE +
C DIMENSION N OF THE SYSTEM, DUE TO THE +
C PRECISION OF THE COMPUTER. +
C THE COMPUTATION CONTINUES UNTIL NK+MK=MAXN +
C AT MOST. +
C IER = 500 IT IS IMPOSSIBLE TO USE THE BSMRZ BECAUSE +
C (y, A**n * r ) IS EQUAL TO ZERO. +
C k k +
C IER = 600 THE VALUE OF m EXCEEDS THE VALUE OF MAXBCK.+
C k +
C +
C EXTERNAL MODULES: +
C +
C - DOUBLE PRECISION FUNCTIONS EUNORM, NVSUMM, PINNER +
C - SUBROUTINES MATVEC, SOLSYS +
C +
C REMARKS: +
C +
C - ALL THE REAL ARGUMENTS PASSED FROM THE CALLING PROGRAM MUST BE +
C IN DOUBLE PRECISION. +
C - THE VARIABLE NK, THE VECTORS Y, X, R, AB, D, C, LMAT AND THE ARRAYS +
C V, W AND RMAT MUST NOT BE MODIFIED BY THE USER BETWEEN TWO +
C CONSECUTIVE CALLS OF THE SUBROUTINE. +
C - IN THE CALLING PROCEDURE, THE ROWS DIMENSION FOR THE ARRAYS A, V +
C AND W MUST BE THE SAME. +
C - IN THE CALLING PROCEDURE THE MINIMAL DIMENSIONS FOR THE MATRIX AND +
C VECTOR ARGUMENTS (IR=N) ARE: +
C A(N,N) +
C AB(N) +
C Y(N) +
C X(N) +
C R(N) +
C YSAV(N) +
C V(N,0:MAXBCK) +
C W(N,0:MAXBCK-1) +
C LMAT((2*MAXBCK-1)*(2*MAXBCK-1)) +
C RMAT(2,0:2*MAXBCK-1) +
C D(0:2*MAXBCK-1) +
C C(0:2*MAXBCK-1) +
C +
C AUTHORS: +
C +
C C. BREZINSKI AND H. SADOK +
C UNIVERSITY OF LILLE - FRANCE +
C M. REDIVO ZAGLIA +
C UNIVERSITY OF PADOVA - ITALY +
C +
C REFERENCES: +
C +
C - AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE ALGORITHMS, +
C NUMERICAL ALGORITHMS, 1(1991), PP. 261-284. +
C +
C - ADDENDUM TO "AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE +
C ALGORITHMS", NUMERICAL ALGORITHMS, 2(1992), PP. 133-136. +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE BSMRZ (N, A, AB, Y, X, R, YSAV, MAXBCK, MAXN, V, W,
* LMAT, RMAT, D, C, IR, NK, EPS, EPS1, EPS2, INIT, IFLAG, IER)
C
C ... SPECIFICATIONS FOR ARGUMENTS
C
INTEGER IER, IFLAG, INIT, IR, MAXBCK, MAXN, N, NK
DOUBLE PRECISION EPS, EPS1, EPS2
DOUBLE PRECISION A(IR,1), D(0:1), W(IR,0:1), V(IR,0:1),
* Y(1), AB(1), X(1), R(1), C(0:1), YSAV(1),
* LMAT(1), RMAT(2,0:1)
C
C ... SPECIFICATIONS FOR EXTERNAL MODULES
C
DOUBLE PRECISION EUNORM, NVSUMM, PINNER
C
C ... SPECIFICATIONS FOR LOCAL VARIABLES
C
INTEGER I, INDS, J, MK, NFIND
SAVE NFIND
C
C ... EXECUTABLE STATEMENTS
C
C
C ... FIRST CALL OF THE SUBROUTINE
C
IF (INIT .EQ. 0) THEN
IER = 0
INIT = 1
NK = 0
NFIND = 0
C
C ... CONTROL THE VALUE OF MAXN
C
IF (MAXN .LT. N) MAXN = N
C
C ... COMPUTATION OF THE VECTORS r AND z
C 0 0
C
CALL MATVEC (N, A, X, 1, R, 1, IR, 0)
DO 10 I = 1, N
R(I) = R(I) - AB(I)
C
C ... STORE THE VALUES FOR V(0) AND W(0)
C
V(I,0) = R(I)
W(I,0) = R(I)
C
C ... TEST FOR y = z
C 0
IF (IFLAG .EQ. 0) Y(I) = R(I)
C
10 CONTINUE
C
C ... TEST FOR THE SOLUTION
C
IF (DSQRT(EUNORM(N,R)) .LT. EPS2) THEN
IER = 200
RETURN
END IF
C
END IF
C
C ... FIRST AND NEXT CALLS OF THE SUBROUTINE
C
IF (IER .NE. 0) THEN
IER = 100
RETURN
END IF
C
C ... COMPUTATION OF D(0)=(y,A**(n )*r )
C k k
C
D(0) = PINNER(N, Y, R, 1, IR)
C
C ... COMPUTATION OF m
C k
MK = 1
C
C ... COMPUTATION OF V(1)
C
CALL MATVEC (N, A, V, 1, V, 2, IR, 0)
C
C ... COMPUTATION OF C(0)=(y,A**(n + 1)*z )
C k k
C(0) = PINNER(N, Y, V, 2, IR)
IF (DABS(C(0)) .LE. EPS) THEN
IF (NK .EQ. MAXN-1) THEN
NK = MAXN
IER = 300
RETURN
END IF
C
DO 30 MK = 2, MAXN-NK
C
C ... CONTROL THE VALUE OF m
C k
C
IF (MK .GT. MAXBCK) THEN
NK = NK + MK - 1
IER = 600
RETURN
END IF
C
C ... COMPUTATION OF V(i), i = 2, 3, ..., MK
C
CALL MATVEC (N, A, V, MK, V, MK+1, IR, 0)
C
C T
C ... COMPUTATION OF y = A * y
C
CALL MATVEC (N, A, Y, 1, AB, 1, IR, 1)
DO 20 J = 1, N
Y(J) = AB(J)
20 CONTINUE
C
C ... COMPUTATION OF C = (y , A**(n +1+i) * z )
C i k k
C AND OF D = (y , A**(n + i) * r )
C i k k
C FOR i = 1, ... , m - 1
C k
C(MK-1) = PINNER(N, Y, V, 2, IR)
D(MK-1) = PINNER(N, Y, R, 1, IR)
C
IF (DABS(C(MK-1)) .GT. EPS) GO TO 40
30 CONTINUE
NK = MAXN
IER = 300
RETURN
END IF
40 CONTINUE
C
C ... THE VALUE OF m HAS BEEN OBTAINED
C k
C
C ... CONTROL IF THE BSMRZ CAN BE USED
C
IF (MK .LE. NK .AND. DABS(D(0)) .EQ. 0.0D0) THEN
IER = 500
RETURN
END IF
C T
C ... COMPUTATION OF Y = A * Y TO HAVE
C T NK+MK
C Y = A y , WHERE y IS THE INITIAL
C VECTOR y. ITS VALUE IS ALSO STORED IN YSAV.
C
CALL MATVEC (N, A, Y, 1, AB, 1, IR, 1)
DO 50 J = 1, N
Y(J) = AB(J)
YSAV(J) = Y(J)
50 CONTINUE
C
C
C ... COMPUTATION OF D(m ) AND OF C(m )
C k k
C
C(MK) = PINNER(N, Y, V, 2, IR)
D(MK) = PINNER(N, Y, R, 1, IR)
C
IF (MK .NE. 1) THEN
DO 70 I = 1, MK - 1
C
C ... COMPUTATION OF W(i), i = 1, 2, ..., MK - 1
C OR i = 1, 2, ..., NK
C
IF (I .LE. NK) CALL MATVEC (N, A, W, I, W, I+1, IR, 0)
C T
C ... COMPUTATION OF YSAV = A * YSAV TO HAVE
C T (NK+2MK-1)
C YSAV = A y , WHERE y IS THE
C INITIAL VECTOR y.
C
C
CALL MATVEC (N, A, YSAV, 1, AB, 1, IR, 1)
DO 60 J = 1, N
YSAV(J) = AB(J)
60 CONTINUE
C
C
C ... COMPUTATION OF D = (y , A**(n + m + i) * r )
C i k k k
C AND OF C = (y , A**(n + m + i + 1)*z )
C i k k k
C FOR i = 2, ...
C
C(MK+I) = PINNER(N, YSAV, V, 2, IR)
IF (I .LT. NK) D(MK+I) = PINNER(N, YSAV, R, 1, IR)
70 CONTINUE
END IF
C
C ... START OF THE LOOP TO COMPUTE BETA , BETA' AND
C i i
C ALPHA , ALPHA' UNTIL THE SYSTEM IS NON SINGULAR
C i i
C
80 CONTINUE
C
C ... COMPUTATION OF BETA , BETA' AND OF
C i i
C ALPHA , ALPHA'
C i i
C
CALL SOLSYS (LMAT, RMAT, C, D, NK, MK, EPS1, INDS)
C
C ... CONTROL FOR THE SINGULARITY OF THE SYSTEM
C
IF (INDS .NE. 0) THEN
C
C ... IF SINGULAR ...
C
C
C ... CONTROL THE VALUE OF MAXN
C
IF (MK .EQ. MAXN-NK) THEN
NK = MAXN
IER = 300
RETURN
END IF
C
C ... INCREASE THE VALUE OF m
C k
MK = MK + 1
C
C ... CONTROL THE VALUE OF m
C k
C
IF (MK .GT. MAXBCK) THEN
NK = NK + MK - 1
IER = 600
RETURN
END IF
C
C ... FOR THE NEW m
C k
C
C T
C ... COMPUTATION OF Y = A * Y TO HAVE
C T NK+MK
C Y = A y , WHERE y IS THE INITIAL
C VECTOR y.
C
CALL MATVEC (N, A, Y, 1, AB, 1, IR, 1)
DO 90 J = 1, N
Y(J) = AB(J)
90 CONTINUE
C
C ... COMPUTATION OF V(MK)
C
CALL MATVEC (N, A, V, MK, V, MK+1, IR, 0)
C
C ... COMPUTATION OF W(MK-1)
C ONLY IF m <= n + 1
C k k
C
IF (MK .LE. NK+1) CALL MATVEC (N, A, W, MK-1, W, MK, IR, 0)
C
DO 110 I = 2, 1, -1
C
C ... COMPUTATION OF THE NEW YSAV
C
CALL MATVEC (N, A, YSAV, 1, AB, 1, IR, 1)
DO 100 J = 1, N
YSAV(J) = AB(J)
100 CONTINUE
C
C ... COMPUTATION OF C(2m -2) AND C(2m -1),
C k k
C ... COMPUTATION OF D(2m -2) AND D(2m -1)
C k k
C IF m <= n
C k k
C ... COMPUTATION OF D(n + m - 1)
C k k
C IF m > n
C k k
C
C(2*MK-I) = PINNER(N, YSAV, V, 2, IR)
IF (MK .LE. NK) THEN
D(2*MK-I) = PINNER(N, YSAV, R, 1, IR)
ELSE IF (I .EQ. 2) THEN
D(MK+NK-1) = PINNER(N, Y, W, NK, IR)
END IF
C
110 CONTINUE
C
C ... RETURN TO SOLVE THE SYSTEM
C
GO TO 80
END IF
C
C ... IF NON SINGULAR ...
C
C ... END OF THE LOOP TO COMPUTE BETA , BETA' AND
C i i
C ALPHA , ALPHA' IF THE SYSTEM IS NON SINGULAR
C i i
C
C ... COMPUTATION OF THE NEXT SOLUTION VECTOR x
C k+1
C AND OF THE NEXT RESIDUAL VECTOR r
C k+1
C
DO 120 J = 1, N
X(J) = X(J) - NVSUMM(RMAT, V, W, IR, J, NK, MK, 0)
R(J) = R(J) - NVSUMM(RMAT, V, W, IR, J, NK, MK, 1)
120 CONTINUE
C
C ... COMPUTE THE VALUE OF n
C k+1
NK = NK + MK
C
C ... TEST FOR THE SOLUTION
C
IF (DSQRT(EUNORM(N,R)) .LT. EPS2) THEN
IER = 200
RETURN
END IF
C
C ... CONTROL THE VALUE OF n
C k
C
IF (NK .GE. N .AND. NFIND .EQ. 0) THEN
IER = 400
NFIND = 1
END IF
C
C ... CONTROL THE VALUE OF MAXN
C
IF (NK .EQ. MAXN) THEN
IER = 300
RETURN
END IF
C
C ... COMPUTATION OF THE VECTOR z
C k+1
DO 130 J = 1, N
V(J,0) = NVSUMM(RMAT, V, W, IR, J, NK-MK, MK, 2)
130 CONTINUE
C
C ... SAVE THE VALUES OF r
C k+1
DO 140 J = 1, N
W(J,0) = R(J)
140 CONTINUE
C
RETURN
END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C FUNCTION NAME - EUNORM +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C DESCRIPTION: +
C +
C DOUBLE PRECISION FUNCTION: +
C COMPUTES THE SQUARE OF THE EUCLIDEAN NORM OF A VECTOR, THAT IS +
C THE SCALAR PRODUCT (VEC,VEC). +
C +
C USAGE: +
C EUNORM (IP, VEC) +
C +
C ARGUMENTS: +
C +
C IP - INPUT DIMENSION OF THE VECTOR SPACE (NUMBER OF +
C ELEMENTS OF THE VECTOR TO BE CONSIDERED). +
C +
C VEC - INPUT REAL VECTOR OF DIMENSION IR >= IP. +
C +
C REMARKS: +
C +
C - THE REAL VECTOR VEC, PASSED FROM THE CALLING PROGRAM, MUST BE +
C IN DOUBLE PRECISION. +
C - IN THE CALLING PROCEDURE THE MINIMAL DIMENSION FOR THE +
C VECTOR ARGUMENT (IR=IP) IS: +
C VEC(IP) +
C +
C AUTHORS: +
C +
C C. BREZINSKI AND H. SADOK +
C UNIVERSITY OF LILLE - FRANCE +
C M. REDIVO ZAGLIA +
C UNIVERSITY OF PADOVA - ITALY +
C +
C REFERENCES: +
C +
C - AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE ALGORITHMS, +
C NUMERICAL ALGORITHMS, 1(1991), PP. 261-284. +
C +
C - ADDENDUM TO "AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE +
C ALGORITHMS", NUMERICAL ALGORITHMS, 2(1992), PP. 133-136. +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
DOUBLE PRECISION FUNCTION EUNORM (IP, VEC)
C
C ... SPECIFICATIONS FOR ARGUMENTS
C
DOUBLE PRECISION VEC(1)
INTEGER IP
C
C ... SPECIFICATIONS FOR LOCAL VARIABLES
C
DOUBLE PRECISION ACC
INTEGER I
C
C ... EXECUTABLE STATEMENTS
C
ACC = 0.0D0
DO 10 I = 1, IP
ACC = ACC + VEC(I)**2
10 CONTINUE
C
C ... SET THE RESULT VALUE
C
EUNORM = ACC
RETURN
END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C SUBPROGRAM NAME - GSOLVE +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C DESCRIPTION: +
C +
C COMPUTES THE SOLUTION OF A SYSTEM OF LINEAR EQUATIONS WITH A +
C AS COEFFICIENT MATRIX. THE MATRIX B CONTAINS THE 2 RIGHT HAND +
C SIDES. THE METHOD IS GAUSSIAN ELIMINATION. +
C IF THE MATRIX A IS NON-SINGULAR THE SOLUTIONS ARE SET IN THE +
C ROWS OF THE MATRIX B. +
C +
C USAGE: +
C CALL GSOLVE (A, B, N, EPS, INDS) +
C +
C ARGUMENTS: +
C +
C A - INPUT/WORKING REAL VECTOR OF DIMENSION (N*N) CONTAINING +
C THE COEFFICIENT MATRIX. +
C +
C B - INPUT/OUTPUT REAL ARRAY OF DIMENSION (2,N). +
C +
C N - INPUT INTEGER INDICATING THE DIMENSION OF THE SYSTEM. +
C +
C EPS - INPUT REAL VALUE USED IN TESTS FOR ZERO. +
C IF DABS(X) .LT. EPS, THEN X IS CONSIDERED TO BE ZERO. +
C +
C INDS - OUTPUT INTEGER VALUE. +
C IF INDS = 0 THEN THE MATRIX A IS NON-SINGULAR. +
C IF INDS = 1 THEN THE MATRIX A IS SINGULAR. +
C +
C REMARKS: +
C +
C - ALL THE REAL ARGUMENTS PASSED FROM THE CALLING PROGRAM MUST BE +
C IN DOUBLE PRECISION. +
C - IN THE CALLING PROCEDURE THE MINIMAL DIMENSIONS FOR THE ARGUMENTS +
C ARE: +
C A(N*N) +
C B(2,N) +
C +
C AUTHORS: +
C +
C C. BREZINSKI AND H. SADOK +
C UNIVERSITY OF LILLE - FRANCE +
C M. REDIVO ZAGLIA +
C UNIVERSITY OF PADOVA - ITALY +
C +
C REFERENCES: +
C +
C - AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE ALGORITHMS, +
C NUMERICAL ALGORITHMS, 1(1991), PP. 261-284. +
C +
C - ADDENDUM TO "AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE +
C ALGORITHMS", NUMERICAL ALGORITHMS, 2(1992), PP. 133-136. +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE GSOLVE (A, B, N, EPS, INDS)
C
C ... SPECIFICATIONS FOR ARGUMENTS
C
INTEGER INDS, N
DOUBLE PRECISION EPS
DOUBLE PRECISION A(1), B(2,1)
C
C ... SPECIFICATIONS FOR LOCAL VARIABLES
C
INTEGER I, IPIVOT, J, K, NN, NN1
DOUBLE PRECISION TMP, PIVOT
C
C ... EXECUTABLE STATEMENTS
C
INDS = 0
C
IF (N .NE. 1) THEN
C
C ... TRIANGULARIZATION
C
DO 70 K = 1, N-1
C
C ... SEARCH FOR THE PIVOT (PARTIAL PIVOTING)
C
PIVOT = 0.0D0
NN = N * (K - 1)
DO 10 I = K, N
TMP = DABS(A(I+NN))
IF (TMP .GE. PIVOT) THEN
IPIVOT = I
PIVOT = TMP
END IF
10 CONTINUE
C
C ... TEST FOR SINGULAR MATRIX
C
IF (DABS(PIVOT) .LT. EPS) THEN
INDS = 1
RETURN
END IF
C
C ... EXCHANGE THE ROW K WITH THE ROW
C CONTAINING THE MAXIMUM PIVOT
C
DO 20 J = K, N
NN1 = N * (J - 1)
TMP = A(K+NN1)
A(K+NN1) = A(IPIVOT+NN1)
A(IPIVOT+NN1) = TMP
20 CONTINUE
DO 30 J = 1, 2
TMP = B(J,K)
B(J,K) = B(J,IPIVOT)
B(J,IPIVOT) = TMP
30 CONTINUE
C
C ... APPLY THE ELEMENTS TRANSFORMATION
C
DO 60 I = K+1, N
IF (A(I+NN) .NE. 0.0D0) THEN
TMP = A(I+NN) / A(K+NN)
DO 40 J = K+1, N
NN1 = N * (J - 1)
A(I+NN1) = A(I+NN1) - TMP * A(K+NN1)
40 CONTINUE
DO 50 J = 1, 2
B(J,I) = B(J,I) - TMP * B(J,K)
50 CONTINUE
END IF
60 CONTINUE
70 CONTINUE
END IF
C
C ... TEST FOR SINGULAR MATRIX
C (PIVOT OF THE LAST ROW)
C
IF (DABS(A(N*N)) .LT. EPS) THEN
INDS = 1
RETURN
END IF
C
C ... SOLVING THE TRIANGULAR SYSTEM
C FOR THE 2 RIGHT HAND SIDES
C
DO 100 K = 1, 2
B(K,N) = B(K,N) / A(N*N)
IF (N .NE. 1) THEN
DO 90 I = N-1, 1, -1
DO 80 J = I+1, N
B(K,I) = B(K,I) - A(I+N*(J-1)) * B(K,J)
80 CONTINUE
B(K,I) = B(K,I) / A(I+N*(I-1))
90 CONTINUE
END IF
100 CONTINUE
RETURN
END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C SUBROUTINE NAME - MATVEC +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C DESCRIPTION: +
C +
C COMPUTES THE PRODUCT BETWEEN A MATRIX (OR ITS TRANSPOSED) AND A +
C VECTOR. +
C THE VECTOR IS STORED AS A COLUMN OF THE ARRAY B. +
C THE RESULTING VECTOR IS PUT IN A PRESCRIBED COLUMN OF THE ARRAY C. +
C +
C USAGE: +
C CALL MATVEC (IP, A, B, IB, C, IC, IR, IFLAG) +
C +
C ARGUMENTS: +
C +
C +
C IP - NUMBER OF ELEMENTS TO BE CONSIDERED IN EACH VECTOR +
C STARTING FROM THE FIRST ONE. +
C +
C A - INPUT REAL MATRIX OF DIMENSION (IP,IP). +
C +
C B - INPUT REAL MATRIX CONTAINING A CERTAIN NUMBER OF COLUMN +
C VECTORS. ITS FIRST DIMENSION MUST BE IP. +
C +
C IB - INPUT COLUMN OF THE MATRIX B. SECOND VECTOR OF THE INNER +
C PRODUCT. +
C +
C C - OUTPUT REAL MATRIX CONTAINING A CERTAIN NUMBER OF COLUMN +
C VECTORS. ITS FIRST DIMENSION MUST BE IP. +
C +
C IC - COLUMN OF THE MATRIX C IN WHICH THE RESULTING VECTOR WILL +
C BE STORED. +
C +
C IR - INPUT ROWS DIMENSION EXACTLY AS SPECIFIED IN THE +
C DIMENSION STATEMENTS IN THE CALLING PROGRAM FOR THE +
C MATRIX B. +
C +
C IFLAG - IF IFLAG .EQ. 0 THE MATRIX A IS CONSIDERED +
C IF IFLAG .NE. 0 THE TRANSPOSED MATRIX A IS CONSIDERED +
C +
C REMARKS: +
C +
C - ALL THE REAL ARGUMENTS PASSED FROM THE CALLING PROGRAM MUST BE +
C IN DOUBLE PRECISION. +
C - IN THE CALLING PROCEDURE, IF THE MATRIX B (OR C) HAS ITS SECOND +
C DIMENSION WITH A LOWER BOUND NOT EQUAL TO 1, THEN THE VALUE IB = 1 +
C (OR IC = 1) CORRESPONDS TO ITS FIRST COLUMN AND NOT TO THE COLUMN +
C WHOSE INDEX IS 1. +
C - IN THE CALLING PROCEDURE THE MINIMAL DIMENSIONS FOR THE +
C MATRIX AND VECTOR ARGUMENTS (IR=IP) ARE: +
C A(IP,IP) +
C B(IP,*) +
C C(IP,*) +
C +
C AUTHORS: +
C +
C C. BREZINSKI AND H. SADOK +
C UNIVERSITY OF LILLE - FRANCE +
C M. REDIVO ZAGLIA +
C UNIVERSITY OF PADOVA - ITALY +
C +
C REFERENCES: +
C +
C - AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE ALGORITHMS, +
C NUMERICAL ALGORITHMS, 1(1991), PP. 261-284. +
C +
C - ADDENDUM TO "AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE +
C ALGORITHMS", NUMERICAL ALGORITHMS, 2(1992), PP. 133-136. +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE MATVEC (IP, A, B, IB, C, IC, IR, IFLAG)
C
C ... SPECIFICATIONS FOR ARGUMENTS
C
INTEGER IB, IC, IFLAG, IP, IR
DOUBLE PRECISION A(IR,1), B(1), C(1)
C
C ... SPECIFICATIONS FOR LOCAL VARIABLES
C
DOUBLE PRECISION ACC
INTEGER I, J, IBSTAR, ICSTAR
C
C ... EXECUTABLE STATEMENTS
C
IBSTAR = (IB-1)*IR
ICSTAR = (IC-1)*IR
DO 20 J = 1, IP
ACC = 0.0D0
DO 10 I = 1, IP
IF (IFLAG .EQ. 0) THEN
ACC = ACC + A(J,I) * B(IBSTAR+I)
ELSE
ACC = ACC + A(I,J) * B(IBSTAR+I)
END IF
10 CONTINUE
C
C ... SET THE RESULT VALUE
C
C(ICSTAR+J) = ACC
20 CONTINUE
RETURN
END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C FUNCTION NAME - NVSUMM +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C DESCRIPTION: +
C +
C DOUBLE PRECISION FUNCTION: +
C COMPUTES THE COMPONENT NR OF A VECTOR OBTAINED AS A LINEAR +
C COMBINATION, WITH THE COEFFICIENTS STORED IN THE FIRST OR IN THE +
C SECOND ROW OF A, OF THE COLUMN VECTORS STORED IN B AND IN C. +
C +
C USAGE: +
C NVSUMM (A, B, C, IR, NR, NK, MK, IFLAG) +
C +
C ARGUMENTS: +
C +
C A - INPUT REAL MATRIX WHICH CONTAINS THE TWO ROWS OF THE +
C COEFFICIENTS OF THE LINEAR COMBINATION. +
C +
C B - INPUT REAL MATRIX CONTAINING A CERTAIN NUMBER OF COLUMN +
C VECTORS. +
C +
C C - INPUT REAL MATRIX CONTAINING A CERTAIN NUMBER OF COLUMN +
C VECTORS. +
C +
C IR - INPUT ROWS DIMENSION EXACTLY AS SPECIFIED IN THE +
C DIMENSION STATEMENTS IN THE CALLING PROGRAM FOR THE +
C MATRICES B AND C. +
C +
C NR - INPUT INTEGER INDICATING THE COMPONENT WHICH HAS TO BE +
C USED AND COMPUTED. +
C +
C NK - INPUT INTEGER INDICATING THE DIMENSION OF THE PREVIOUS +
C MATRIX OF THE SYSTEM. +
C +
C MK - INPUT INTEGER INDICATING THE DIMENSION OF THE NEW JUMP. +
C +
C IFLAG - INPUT INTEGER FLAG. +
C IF IFLAG = 0 THE NEW SOLUTION X WILL BE COMPUTED. +
C IF IFLAG = 1 THE NEW RESIDUAL R WILL BE COMPUTED. +
C IF IFLAG = 2 THE NEW VECTOR Z WILL BE COMPUTED. +
C +
C REMARKS: +
C +
C - ALL THE REAL ARGUMENTS PASSED FROM THE CALLING PROGRAM MUST BE +
C IN DOUBLE PRECISION. +
C - IN THE CALLING PROCEDURE THE DIMENSIONS FOR THE ARRAY ARGUMENTS +
C ARE: +
C A(2,0:*) +
C B(IR,0:*) +
C C(IR,0:*) +
C +
C AUTHORS: +
C +
C C. BREZINSKI AND H. SADOK +
C UNIVERSITY OF LILLE - FRANCE +
C M. REDIVO ZAGLIA +
C UNIVERSITY OF PADOVA - ITALY +
C +
C REFERENCES: +
C +
C - AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE ALGORITHMS, +
C NUMERICAL ALGORITHMS, 1(1991), PP. 261-284. +
C +
C - ADDENDUM TO "AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE +
C ALGORITHMS", NUMERICAL ALGORITHMS, 2(1992), PP. 133-136. +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
DOUBLE PRECISION FUNCTION NVSUMM (A, B, C, IR, NR, NK, MK, IFLAG)
C
C ... SPECIFICATIONS FOR ARGUMENTS
C
INTEGER IR, NK, NR, MK, IFLAG
DOUBLE PRECISION A(2,0:1), B(IR,0:1), C(IR,0:1)
C
C ... SPECIFICATIONS FOR LOCAL VARIABLES
C
DOUBLE PRECISION ACC
INTEGER I, I2
C
C ... EXECUTABLE STATEMENTS
C
C
C ... CASE MK = 1
C
IF (MK .EQ. 1) THEN
IF (IFLAG .EQ. 0) THEN
ACC = A(1,0) * B(NR,0)
ELSE IF (IFLAG .EQ. 1) THEN
ACC = A(1,0) * B(NR,1)
ELSE IF (MK .GT. NK) THEN
ACC = A(2,0) * B(NR,0) + B(NR,1)
ELSE
ACC = A(2,0) * B(NR,0) + A(2,1) * C(NR,0) + B(NR,1)
END IF
NVSUMM = ACC
RETURN
END IF
C
ACC = 0.0D0
C
C ... CASE MK <= NK
C
IF (MK .LE. NK) THEN
IF (IFLAG .EQ. 0) THEN
DO 10 I = 0, MK - 2
I2 = I * 2
ACC = ACC + A(1,I2) * B(NR,I) + A(1,I2+1) * C(NR,I)
10 CONTINUE
ACC = ACC + A(1,2*MK-2) * B(NR,MK-1)
ELSE IF (IFLAG .EQ. 1) THEN
DO 20 I = 0, MK - 2
I2 = I * 2
ACC = ACC + A(1,I2)*B(NR,I+1) + A(1,I2+1)*C(NR,I+1)
20 CONTINUE
ACC = ACC + A(1,2*MK-2) * B(NR,MK)
ELSE
DO 30 I = 0, MK - 1
I2 = I * 2
ACC = ACC + A(2,I2) * B(NR,I) + A(2,I2+1) * C(NR,I)
30 CONTINUE
ACC = ACC + B(NR,MK)
END IF
C
C ... CASE MK > NK
C
ELSE
IF (IFLAG .EQ. 0) THEN
IF (NK .NE. 0) THEN
DO 40 I = 0, NK - 1
I2 = I * 2
ACC = ACC + A(1,I2) * B(NR,I) + A(1,I2+1) * C(NR,I)
40 CONTINUE
END IF
DO 50 I = NK, MK-1
ACC = ACC + A(1,NK+I) * B(NR,I)
50 CONTINUE
ELSE IF (IFLAG .EQ. 1) THEN
IF (NK .NE. 0) THEN
DO 60 I = 0, NK - 1
I2 = I * 2
ACC = ACC + A(1,I2)*B(NR,I+1) + A(1,I2+1)*C(NR,I+1)
60 CONTINUE
END IF
DO 70 I = NK, MK-1
ACC = ACC + A(1,NK+I) * B(NR,I+1)
70 CONTINUE
ELSE
IF (NK .NE. 0) THEN
DO 80 I = 0, NK - 1
I2 = I * 2
ACC = ACC + A(2,I2) * B(NR,I) + A(2,I2+1) * C(NR,I)
80 CONTINUE
END IF
DO 90 I = NK, MK-1
ACC = ACC + A(2,NK+I) * B(NR,I)
90 CONTINUE
ACC = ACC + B(NR,MK)
END IF
END IF
C
C ... SET THE RESULT VALUE
C
NVSUMM = ACC
RETURN
END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C FUNCTION NAME - PINNER +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C DESCRIPTION: +
C +
C DOUBLE PRECISION FUNCTION: +
C COMPUTES THE INNER PRODUCT OF TWO VECTORS. +
C THE DIMENSION OF THE FIRST VECTOR IS IR. THE SECOND VECTOR +
C IS STORED AS A COLUMN OF THE ARRAY B WHOSE FIRST DIMENSION IS IR. +
C +
C USAGE: +
C PINNER (IP, A, B, IB, IR) +
C +
C ARGUMENTS: +
C +
C IP - NUMBER OF ELEMENTS TO BE CONSIDERED IN EACH VECTOR +
C STARTING FROM THE FIRST ONE. +
C +
C A - INPUT REAL VECTOR (FIRST VECTOR OF THE INNER PRODUCT). +
C +
C B - INPUT REAL MATRIX CONTAINING A CERTAIN NUMBER OF COLUMN +
C VECTORS. +
C +
C IB - INPUT INTEGER DEFINING THE NUMBER OF THE COLUMN OF THE +
C MATRIX B WHICH IS THE SECOND VECTOR OF THE INNER PRODUCT. +
C +
C IR - INPUT ROWS DIMENSION EXACTLY AS SPECIFIED IN THE +
C DIMENSION STATEMENTS IN THE CALLING PROGRAM FOR THE +
C MATRIX B. +
C +
C REMARKS: +
C +
C - ALL THE REAL ARGUMENTS PASSED FROM THE CALLING PROGRAM MUST BE +
C IN DOUBLE PRECISION. +
C - IN THE CALLING PROCEDURE THE MINIMAL DIMENSIONS FOR THE +
C MATRIX AND VECTOR ARGUMENTS (IR=IP) ARE: +
C A(IP) +
C B(IP,*) +
C - IN THE CALLING PROCEDURE, IF THE MATRIX B HAS ITS SECOND DIMENSION +
C WITH A LOWER BOUND NOT EQUAL TO 1, THEN THE VALUE IB = 1 +
C CORRESPONDS TO ITS FIRST COLUMN AND NOT TO THE COLUMN WHOSE INDEX +
C IS 1. +
C +
C AUTHORS: +
C +
C C. BREZINSKI AND H. SADOK +
C UNIVERSITY OF LILLE - FRANCE +
C M. REDIVO ZAGLIA +
C UNIVERSITY OF PADOVA - ITALY +
C +
C REFERENCES: +
C +
C - AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE ALGORITHMS, +
C NUMERICAL ALGORITHMS, 1(1991), PP. 261-284. +
C +
C - ADDENDUM TO "AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE +
C ALGORITHMS", NUMERICAL ALGORITHMS, 2(1992), PP. 133-136. +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
DOUBLE PRECISION FUNCTION PINNER (IP, A, B, IB, IR)
C
C ... SPECIFICATIONS FOR ARGUMENTS
C
DOUBLE PRECISION A(1), B(1)
INTEGER IB, IP, IR
C
C ... SPECIFICATIONS FOR LOCAL VARIABLES
C
DOUBLE PRECISION ACC
INTEGER I, IBSTAR
C
C ... EXECUTABLE STATEMENTS
C
IBSTAR = (IB-1)*IR
ACC = 0.0D0
DO 10 I = 1, IP
ACC = ACC + A(I) * B(IBSTAR+I)
10 CONTINUE
C
C ... SET THE RESULT VALUE
C
PINNER = ACC
RETURN
END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C SUBROUTINE NAME - SOLSYS +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C DESCRIPTION: +
C +
C COMPUTES THE SOLUTIONS OF THE AUXILIARY SYSTEMS THAT IS THE +
C ALPHA , ALPHA' ,BETA AND BETA' . +
C i i i i +
C +
C USAGE: +
C +
C CALL SOLSYS (LMAT, RMAT, C, D, NK, MK, EPS, INDS) +
C +
C ARGUMENTS: +
C +
C LMAT - WORKING REAL VECTOR TO STORE THE COEFFICIENTS OF THE +
C AUXILIARY SYSTEM. +
C +
C RMAT - OUTPUT REAL MATRIX OF DIMENSION (2,0:*). +
C IF THE SYSTEM IS NON SINGULAR, IT CONTAINS IN OUTPUT: +
C ROW 1 THE BETA , BETA' SOLUTIONS. +
C i i +
C ROW 2 THE ALPHA , ALPHA' SOLUTIONS. +
C i i +
C +
C C - INPUT REAL VECTOR CONTAINING THE C 'S. +
C i +
C +
C D - INPUT REAL VECTOR CONTAINING THE D 'S. +
C i +
C +
C NK - INPUT INTEGER INDICATING THE DIMENSION OF THE PREVIOUS +
C MATRIX OF THE SYSTEM. +
C +
C MK - INPUT INTEGER INDICATING THE DIMENSION OF THE NEW JUMP. +
C +
C EPS - INPUT REAL VALUE USED FOR TESTING THE PIVOTS IN GAUSSIAN +
C ELIMINATION. +
C IF DABS(X) .LT. EPS, THEN X IS CONSIDERED TO BE ZERO. +
C THE SUBROUTINE DOES NOT CONTROL WHETHER OR NOR EPS IS +
C A NEGATIVE OR ZERO REAL NUMBER. +
C +
C INDS - OUTPUT INTEGER VALUE. +
C IF INDS = 0 THEN THE SYSTEM IS NON-SINGULAR. +
C IF INDS = 1 THEN THE SYSTEM IS SINGULAR. +
C +
C EXTERNAL MODULES: +
C +
C - SUBROUTINES GSOLVE, STOMAT, STORHS +
C +
C REMARKS: +
C +
C - ALL THE REAL ARGUMENTS PASSED FROM THE CALLING PROGRAM MUST BE +
C IN DOUBLE PRECISION. +
C - IN THE CALLING PROCEDURE THE DIMENSIONS FOR THE ARGUMENTS ARE: +
C LMAT(*) +
C RMAT(2,0:*) +
C +
C AUTHORS: +
C +
C C. BREZINSKI AND H. SADOK +
C UNIVERSITY OF LILLE - FRANCE +
C M. REDIVO ZAGLIA +
C UNIVERSITY OF PADOVA - ITALY +
C +
C REFERENCES: +
C +
C - AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE ALGORITHMS, +
C NUMERICAL ALGORITHMS, 1(1991), PP. 261-284. +
C +
C - ADDENDUM TO "AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE +
C ALGORITHMS", NUMERICAL ALGORITHMS, 2(1992), PP. 133-136. +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE SOLSYS (LMAT, RMAT, C, D, NK, MK, EPS, INDS)
C
C ... SPECIFICATIONS FOR ARGUMENTS
C
INTEGER INDS, MK, NK
DOUBLE PRECISION EPS
DOUBLE PRECISION LMAT(1), RMAT(2,0:1), C(0:1), D(0:1)
C
C ... SPECIFICATIONS FOR LOCAL VARIABLES
C
INTEGER IFLAG, NN
C
C ... EXECUTABLE STATEMENTS
C
C
C ... CASE MK = 1
C
IF (MK .EQ. 1) THEN
INDS = 0
RMAT(1,0) = D(0)/C(0)
RMAT(2,0) = -C(1)
IF (MK .LE. NK) THEN
LMAT(1) = -C(0)/D(0)
RMAT(2,1) = LMAT(1)
RMAT(2,0) = RMAT(2,0) - LMAT(1) * D(1)
END IF
RMAT(2,0) = RMAT(2,0)/C(0)
RETURN
END IF
C
C ... CHOICE FOR THE SYSTEM
C
IF (MK .LE. NK) THEN
NN = MK - 1
IFLAG = 1
ELSE
NN = NK
IFLAG = 2
END IF
C
C ... COMPUTATION OF BETA , BETA' AND OF
C i i
C ALPHA , ALPHA'
C i i
C
C ... STORE THE COEFFICIENTS AND THE RIGHT HAND
C SIDES
C
CALL STOMAT (LMAT, C, D, NK, MK, IFLAG)
CALL STORHS (RMAT, C, D, NK, MK, IFLAG)
C
C ... SOLVE THE SYSTEM
C
CALL GSOLVE (LMAT, RMAT, NN+MK, EPS, INDS)
C
C ... IF NON SINGULAR SYSTEM, STORE THE
C ADDITIONAL SOLUTION
C
IF (INDS .EQ. 0 .AND. IFLAG .EQ. 1) RMAT(2,2*MK-1) = -C(0)/D(0)
C
RETURN
END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C FUNCTION NAME - SSTRIA +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C DESCRIPTION: +
C +
C DOUBLE PRECISION FUNCTION: +
C COMPUTES ONE OF THE UNKNOWNS OF THE TRIANGULAR SYSTEMS GIVING +
C THE ALPHA'S AND THE BETA'S. +
C +
C USAGE: +
C SSTRIA(A, B, ACCINI, I, MK, INC) +
C +
C ARGUMENTS: +
C +
C A - INPUT REAL VECTOR CONTAINING THE UNKNOWNS WHICH HAVE +
C ALREADY BEEN COMPUTED. +
C +
C B - INPUT REAL VECTOR WHICH CONTAINS THE COEFFICIENTS OF THE +
C TRIANGULAR MATRIX. +
C +
C ACCINI - INPUT REAL VECTOR. INITIALIZATION FOR THE SUM. +
C +
C I - INPUT INTEGER INDICATING THE INDEX TO BE USED IN THE +
C COMPUTATION. +
C +
C MK - INPUT INTEGER INDICATING THE DIMENSION OF THE NEW JUMP. +
C +
C INC - INPUT INTEGER RELATED TO THE FIRST COMPONENT OF A TO +
C BE CONSIDERED. +
C IF INC = 0 THE ALPHA'S WILL BE COMPUTED +
C IF INC = 1 THE BETA'S WILL BE COMPUTED. +
C +
C REMARKS: +
C +
C - ALL THE REAL ARGUMENTS PASSED FROM THE CALLING PROGRAM MUST BE +
C IN DOUBLE PRECISION. +
C - IN THE CALLING PROCEDURE THE MINIMAL DIMENSIONS FOR THE +
C MATRIX AND VECTOR ARGUMENTS (IR=IP) ARE: +
C A(0:*) +
C B(0:*) +
C +
C AUTHORS: +
C +
C C. BREZINSKI AND H. SADOK +
C UNIVERSITY OF LILLE - FRANCE +
C M. REDIVO ZAGLIA +
C UNIVERSITY OF PADOVA - ITALY +
C +
C REFERENCES: +
C +
C - AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE ALGORITHMS, +
C NUMERICAL ALGORITHMS, 1(1991), PP. 261-284. +
C +
C - ADDENDUM TO "AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE +
C ALGORITHMS", NUMERICAL ALGORITHMS, 2(1992), PP. 133-136. +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
DOUBLE PRECISION FUNCTION SSTRIA(A, B, ACCINI, I, MK, INC)
C
C ... SPECIFICATIONS FOR ARGUMENTS
C
DOUBLE PRECISION A(0:1), B(0:1)
DOUBLE PRECISION ACCINI
INTEGER I, INC, MK
C
C ... SPECIFICATIONS FOR LOCAL VARIABLES
C
DOUBLE PRECISION ACC
INTEGER J
C
C ... EXECUTABLE STATEMENTS
C
ACC = ACCINI
DO 10 J = 1, I
ACC = ACC - A(MK-I+J-INC) * B(J)
10 CONTINUE
C
C ... COMPUTE THE FUNCTION VALUE
C
SSTRIA = ACC/B(0)
RETURN
END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C SUBROUTINE NAME - STOMAT +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C DESCRIPTION: +
C +
C STORES IN THE VECTOR LMAT THE COEFFICIENTS OF THE SYSTEM TO BE +
C SOLVED FOR THE COMPUTATION OF THE ALPHA , ALPHA', BETA , BETA' . +
C i i i i +
C +
C USAGE: +
C CALL STOMAT (LMAT, C, D, NK, MK, IFLAG) +
C +
C ARGUMENTS: +
C +
C LMAT - OUTPUT REAL VECTOR WHICH CONTAINS THE COEFFICIENTS OF +
C THE AUXILIARY SYSTEM. +
C +
C C - INPUT REAL VECTOR CONTAINING THE C 'S. +
C i +
C +
C D - INPUT REAL VECTOR CONTAINING THE D 'S. +
C i +
C +
C NK - INPUT INTEGER INDICATING THE DIMENSION OF THE PREVIOUS +
C MATRIX OF THE SYSTEM. +
C +
C MK - INPUT INTEGER INDICATING THE DIMENSION OF THE NEW JUMP. +
C +
C IFLAG - INPUT INTEGER FLAG. +
C IF IFLAG = 1 THE CASE MK <= NK IS CONSIDERED. +
C IF IFLAG = 2 THE CASE MK > NK IS CONSIDERED. +
C +
C REMARKS: +
C +
C - ALL THE REAL ARGUMENTS PASSED FROM THE CALLING PROGRAM MUST BE +
C IN DOUBLE PRECISION. +
C - IN THE CALLING PROCEDURE THE DIMENSIONS FOR THE VECTOR ARGUMENTS +
C ARE: +
C LMAT(*) +
C B(0:*) +
C D(0:*) +
C +
C AUTHORS: +
C +
C C. BREZINSKI AND H. SADOK +
C UNIVERSITY OF LILLE - FRANCE +
C M. REDIVO ZAGLIA +
C UNIVERSITY OF PADOVA - ITALY +
C +
C REFERENCES: +
C +
C - AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE ALGORITHMS, +
C NUMERICAL ALGORITHMS, 1(1991), PP. 261-284. +
C +
C - ADDENDUM TO "AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE +
C ALGORITHMS", NUMERICAL ALGORITHMS, 2(1992), PP. 133-136. +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE STOMAT (LMAT, C, D, NK, MK, IFLAG)
C
C ... SPECIFICATIONS FOR ARGUMENTS
C
INTEGER IFLAG, NK, MK
DOUBLE PRECISION LMAT(1), C(0:1), D(0:1)
C
C ... SPECIFICATIONS FOR LOCAL VARIABLES
C
INTEGER ICOL, IROW, NI, NN, NN2, NNM
C
C ... EXECUTABLE STATEMENTS
C
C
C ... CASE MK <= NK
C
IF (IFLAG .EQ. 1) THEN
NN = MK - 1
ELSE
C
C ... CASE MK > NK
C
NN = NK
END IF
C
C ... COMPUTE THE DIMENSION OF THE SYSTEM
C
NNM = NN + MK
C
C ... CLEAR THE ARRAY COMPONENTS
C
DO 10 IROW = 1, NNM * NNM
LMAT(IROW) = 0.0D0
10 CONTINUE
C
C ... BOTH CASES: MK <= NK AND MK > NK
C
NN2 = NN * 2
DO 30 ICOL = NN2, 0, -2
NI = NN - ICOL/2
DO 20 IROW = NI, NNM-1
LMAT(NNM*(ICOL+1)-IROW) = C(IROW-NI)
IF (ICOL .NE. 0) LMAT(NNM*ICOL-IROW) = D(IROW-NI)
20 CONTINUE
30 CONTINUE
C
C ... CASE MK > NK
C
IF (IFLAG .EQ. 2) THEN
DO 50 ICOL = 1, MK-NN-1
NN2 = NN*2 + ICOL
DO 40 IROW = 0, NNM-1
LMAT(NNM*(NN2+1)-IROW) = C(IROW+ICOL)
40 CONTINUE
50 CONTINUE
END IF
RETURN
END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C SUBROUTINE NAME - STORHS +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C DESCRIPTION: +
C +
C STORES IN THE ARRAY RMAT THE RIGHT HAND SIDES OF THE SYSTEM TO BE +
C SOLVED FOR THE COMPUTATION OF THE ALPHA , ALPHA', BETA , BETA' . +
C i i i i +
C +
C USAGE: +
C CALL STORHS (RMAT, C, D, NK, MK, IFLAG) +
C +
C ARGUMENTS: +
C +
C RMAT - OUTPUT REAL MATRIX OF DIMENSION (2,0:*) WHICH CONTAINS +
C THE 2 RIGHT HAND SIDES OF THE SYSTEM. +
C +
C C - INPUT REAL VECTOR CONTAINING THE C 'S. +
C i +
C +
C D - INPUT REAL VECTOR CONTAINING THE D 'S. +
C i +
C +
C NK - INPUT INTEGER INDICATING THE DIMENSION OF THE PREVIOUS +
C MATRIX OF THE SYSTEM. +
C +
C MK - INPUT INTEGER INDICATING THE DIMENSION OF THE NEW JUMP. +
C +
C IFLAG - INPUT INTEGER FLAG. +
C IF IFLAG = 1 THE CASE MK <= NK IS CONSIDERED. +
C IF IFLAG = 2 THE CASE MK > NK IS CONSIDERED. +
C +
C REMARKS: +
C +
C - ALL THE REAL ARGUMENTS PASSED FROM THE CALLING PROGRAM MUST BE +
C IN DOUBLE PRECISION. +
C - IN THE CALLING PROCEDURE THE DIMENSIONS FOR THE MATRIX AND VECTOR +
C ARGUMENTS ARE: +
C RMAT(2,0:*) +
C B(0:*) +
C D(0:*) +
C +
C AUTHORS: +
C +
C C. BREZINSKI AND H. SADOK +
C UNIVERSITY OF LILLE - FRANCE +
C M. REDIVO ZAGLIA +
C UNIVERSITY OF PADOVA - ITALY +
C +
C REFERENCES: +
C +
C - AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE ALGORITHMS, +
C NUMERICAL ALGORITHMS, 1(1991), PP. 261-284. +
C +
C - ADDENDUM TO "AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE +
C ALGORITHMS", NUMERICAL ALGORITHMS, 2(1992), PP. 133-136. +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE STORHS (RMAT, C, D, NK, MK, IFLAG)
C
C ... SPECIFICATIONS FOR ARGUMENTS
C
INTEGER IFLAG, NK, MK
DOUBLE PRECISION RMAT(2,0:1), C(0:1), D(0:1)
C
C ... SPECIFICATIONS FOR LOCAL VARIABLES
C
INTEGER IROW, ICOL, NNM, NN2
DOUBLE PRECISION CONST
C
C ... EXECUTABLE STATEMENTS
C
C
C ... COMPUTE THE UPPER BOUND INDEX OF THE
C COMPONENTS IN THE ARRAY
C
C ... CASE MK <= NK
C
IF (IFLAG .EQ. 1) THEN
NNM = 2*(MK-1)
ELSE
C
C ... CASE MK > NK
C
NNM = NK + MK -1
END IF
C
C ... CLEAR THE ARRAY COMPONENTS
C
DO 20 ICOL = 0, NNM
DO 10 IROW = 1, 2
RMAT(IROW,ICOL) = 0.0D0
10 CONTINUE
20 CONTINUE
C
C ... STORES THE RIGHT HAND SIDE FOR THE
C SYSTEM THAT COMPUTES THE BETA AND BETA'
C i i
C
IROW = 1
C
C ... BOTH CASES: MK <= NK AND MK > NK
C
DO 30 ICOL = MK-1, 0, -1
RMAT(IROW,MK-1-ICOL) = D(ICOL)
30 CONTINUE
C
C ... STORES THE RIGHT HAND SIDE FOR THE
C SYSTEM THAT COMPUTES THE ALPHA AND ALPHA'
C i i
C
IROW = 2
C
C ... CASE MK <= NK
C
IF (IFLAG .EQ. 1) CONST = - C(0)/D(0)
C
C ... BOTH CASES: MK <= NK AND MK > NK
C
DO 40 ICOL = 0, NNM
NN2 = 2*MK-1-ICOL
RMAT(IROW,ICOL) = - C(NN2)
C
C ... CASE MK <= NK
C
IF (IFLAG .EQ. 1)
* RMAT(IROW,ICOL) = RMAT(IROW,ICOL) - CONST * D(NN2)
40 CONTINUE
RETURN
END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C FUNCTION NAME - VSUMMA +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C +
C DESCRIPTION: +
C +
C DOUBLE PRECISION FUNCTION: +
C COMPUTES THE COMPONENT NR OF A VECTOR OBTAINED AS A LINEAR +
C COMBINATION, WITH THE COEFFICIENTS STORED IN A, OF THE MK +
C COLUMN VECTORS STORED IN B. +
C +
C USAGE: +
C VSUMMA (A, B, IR, NR, MK, IFLAG) +
C +
C ARGUMENTS: +
C +
C A - INPUT REAL VECTOR WHICH CONTAINS THE COEFFICIENTS OF +
C THE LINEAR COMBINATION. +
C +
C B - INPUT REAL MATRIX CONTAINING A CERTAIN NUMBER OF COLUMN +
C VECTORS. +
C +
C IR - INPUT ROWS DIMENSION EXACTLY AS SPECIFIED IN THE +
C DIMENSION STATEMENTS IN THE CALLING PROGRAM FOR THE +
C MATRIX B. +
C +
C NR - INPUT INTEGER INDICATING THE COMPONENT WHICH HAS TO BE +
C USED AND COMPUTED. +
C +
C MK - INPUT INTEGER INDICATING THE NUMBER OF VECTORS OF THE +
C MATRIX B TO BE CONSIDERED. +
C +
C IFLAG - INPUT INTEGER INDICATING THE LOWER AND UPPER INDEXES OF +
C THE COLUMN VECTORS OF THE MATRIX B TO BE CONSIDERED. +
C IF IFLAG = 0 THE NEW SOLUTION X WILL BE COMPUTED +
C IF IFLAG = 1 THE NEW RESIDUAL R WILL BE COMPUTED +
C IF IFLAG = 2 THE NEW VECTOR Z WILL BE COMPUTED. +
C +
C REMARKS: +
C +
C - ALL THE REAL ARGUMENTS PASSED FROM THE CALLING PROGRAM MUST BE +
C IN DOUBLE PRECISION. +
C - IN THE CALLING PROCEDURE THE MINIMAL DIMENSIONS FOR THE +
C MATRIX AND VECTOR ARGUMENTS (IR=IP) ARE: +
C A(0:*) +
C B(IP,0:*) +
C +
C AUTHORS: +
C +
C C. BREZINSKI AND H. SADOK +
C UNIVERSITY OF LILLE - FRANCE +
C M. REDIVO ZAGLIA +
C UNIVERSITY OF PADOVA - ITALY +
C +
C REFERENCES: +
C +
C - AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE ALGORITHMS, +
C NUMERICAL ALGORITHMS, 1(1991), PP. 261-284. +
C +
C - ADDENDUM TO "AVOIDING BREAKDOWN AND NEAR-BREAKDOWN IN LANCZOS TYPE +
C ALGORITHMS", NUMERICAL ALGORITHMS, 2(1992), PP. 133-136. +
C +
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
DOUBLE PRECISION FUNCTION VSUMMA (A, B, IR, NR, MK, IFLAG)
C
C ... SPECIFICATIONS FOR ARGUMENTS
C
INTEGER IFLAG, IR, MK, NR
DOUBLE PRECISION A(0:1), B(IR,0:1)
C
C ... SPECIFICATIONS FOR LOCAL VARIABLES
C
DOUBLE PRECISION ACC
INTEGER I, INK, LIND, UIND
C
C ... EXECUTABLE STATEMENTS
C
ACC = 0.0D0
IF (IFLAG .EQ. 0) THEN
INK = 0
LIND = 0
UIND = MK-1
ELSE IF (IFLAG .EQ. 1) THEN
INK = 1
LIND = 1
UIND = MK
ELSE
INK = 0
LIND = 0
UIND = MK
END IF
DO 10 I = LIND, UIND
ACC = ACC + A(I-INK) * B(NR,I)
10 CONTINUE
C
C ... SET THE RESULT VALUE
C
VSUMMA = ACC
RETURN
END