*************************************************************************** * 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