C C THIS DRIVER TESTS EISPACK FOR THE CLASS OF REAL C BAND MATRICES EXHIBITING THE USE OF EISPACK TO FIND THE C SOLUTION TO THE EQUATION A*X = B . THIS DRIVER SUMMARIZES C THE FIGURES OF MERIT FOR THE PATH. C C THIS DRIVER IS CATALOGUED AS EISPDRV4(RBLSUMAR). C C THE DIMENSION OF A AND AHOLD SHOULD BE NM BY MBB. C THE DIMENSION OF X,B, AND BHOLD SHOULD BE NM BY NM. C THE DIMENSION OF RESDUL SHOULD BE NM. C THE DIMENSION OF RV SHOULD BE NM*MBB. C THE DIMENSION OF RV1 AND RV6 SHOULD BE NM. C HERE NM = 20 AND MBB = 39. C C REAL A( 20, 39),X( 20, 20),B( 20, 20),RV( 780), X RESDUL( 20),RV1( 20),RV6( 20),TCRIT( 1),EPSLON, X RESMAX,E21,DFL REAL AHOLD( 20, 39),BHOLD( 20, 20) INTEGER IERR( 1),ERROR DATA IREAD1/1/,IREADC/5/,IWRITE/6/ C OPEN(UNIT=IREAD1,FILE='FILE49') OPEN(UNIT=IREADC,FILE='FILE50') REWIND IREAD1 REWIND IREADC C NM = 20 MBB = 39 LCOUNT = 0 WRITE(IWRITE,1) 1 FORMAT(1H1,19X,57H EXPLANATION OF COLUMN ENTRIES FOR THE SUMMARY S XTATISTICS//1H ,95(1H-)/22H ORDER B.W. IP BANDV / X1H ,95(1H-)// X49H UNDER 'ORDER' IS THE ORDER OF EACH TEST MATRIX. // X52H UNDER 'B.W.' IS THE FULL BAND WIDTH OF THE MATRIX. // X59H UNDER 'IP' IS THE COLUMN DIMENSION OF THE CONSTANT MATRIX. // X62H UNDER 'BANDV' ARE TWO NUMBERS. THE FIRST NUMBER, AN INTEGER,, X29H IS THE ABSOLUTE VALUE OF THE / X92H ERROR FLAG RETURNED FROM BANDV. THE SECOND NUMBER IS THE MEA XSURE OF PERFORMANCE BASED / X46H UPON THE RESIDUAL COMPUTED FROM THE SOLUTION. /) WRITE(IWRITE,15) 15 FORMAT(1X,21HS.P. VERSION 04/15/83 ) 5 FORMAT( 80H1 TABULATION OF THE ERROR FLAG ERROR AND THE ME XASURE OF PERFORMANCE FOR /5X,13HTHE EISPACK , X73H CODES. THIS RUN DISPLAYS THESE STATISTICS FOR REAL BAND LINEA XR SYSTEMS. // X22H ORDER B.W. IP BANDV ) 10 CALL RMATIN(NM,MBB,N,MBW,A,AHOLD,BHOLD,0,0) CALL RMATIN(NM,NM,N,IP,X,AHOLD,BHOLD,1,0) READ(IREADC,20) IE21 20 FORMAT(I6) E21 = IE21 C C RBL USING BANDV C ICT = 1 DO 905 I = 1,N RV1(I) = 0.0E0 905 CONTINUE CALL BANDV(NM,N,MBW,A,E21,IP,RV1,X,ERROR,780,RV,RV6) IERR(ICT) = IABS(ERROR) CALL RMATIN(NM,MBB,N,MBW,A,AHOLD,BHOLD,0,1) CALL RMATIN(NM,NM,N,IP,B,AHOLD,BHOLD,1,1) CALL RBLRES(NM,N,MBW,IP,A,B,X,RESDUL,E21) RESMAX = 0.0E0 DO 935 I = 1,IP IF( RESDUL(I) .GT. RESMAX ) RESMAX = RESDUL(I) 935 CONTINUE DFL = N TCRIT(ICT) = RESMAX/EPSLON(DFL) C IF( MOD(LCOUNT,35) .EQ. 0 ) WRITE(IWRITE,5) LCOUNT = LCOUNT + 1 MBF = MBW IF( E21 .EQ. 1.0E0 ) MBF = MBW*2 - 1 WRITE(IWRITE,520) N,MBF,IP,IERR(1),TCRIT(1) 520 FORMAT(1X,I3,2I5,I4,F6.2) GO TO 10 END SUBROUTINE RBLRES(NM,N,MBW,IP,A,B,X,NORM,E21) C REAL A(NM,MBW),B(NM,IP),X(NM,IP),NORM(N) REAL NORMA,NORMX,NORMR,SUM,SUMA,E21 C C THIS SUBROUTINE FORMS THE 1-NORM OF THE RESIDUAL MATRIX C A*X-B WHERE A IS A REAL BAND MATRIX, WHICH MAY BE SYMMETRIC, C X IS A REAL N BY IP MATRIX, AND B IS A REAL C N BY IP MATRIX. C C THIS SUBROUTINE IS CATALOGUED AS EISPDRV4(RBLRES). C C INPUT. C NM IS THE ROW DIMENSION OF TWO-DIMENSIONAL ARRAY PARAMETERS C AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT; C C N IS THE ROW DIMENSION OF THE MATRIX A; C C MBW IS THE NUMBER OF COLUMNS OF THE ARRAY A USED TO STORE THE C BAND MATRIX. IF THE MATRIX IS SYMMETRIC, MBW IS ITS (HALF) C BAND WIDTH, DENOTED MB AND DEFINED AS THE NUMBER OF ADJACENT C DIAGONALS, INCLUDING THE PRINCIPAL DIAGONAL, REQUIRED TO C SPECIFY THE NON-ZERO PORTION OF THE LOWER TRIANGLE OF THE C MATRIX. IF THE MATRIX IS NON-SYMMETRIC IT MUST HAVE THE SAM C NUMBER OF ADJACENT DIAGONALS ABOVE THE MAIN DIAGONAL AS C BELOW. IN THIS CASE MBW = 2*MB-1; C C IP IS THE COLUMN DIMENSION OF B AND X; C C A(N,MBW) IS AN ARRAY WHICH CONTAINS IN ITS COLUMNS THE C SUBDIAGONALS AND DIAGONAL OF THE SYMMETRIC BAND C MATRIX. IF A IS NON-SYMMETRIC IT ALSO CONTAINS C THE SUPER-DIAGONALS; C C B IS A REAL FULL N BY IP MATRIX; C C X IS A REAL FULL N BY IP MATRIX; C C E21 TELLS IF THE MATRIX IS SYMMETRIC. IF E21 EQUALS ONE THE C MATRIX IS SYMMETRIC AND IF E21 EQUALS MINUS ONE THE C MATRIX IS NON-SYMMETRIC. C C OUTPUT. C C NORM(N) IS AN ARRAY SUCH THAT FOR EACH K, C NORM(K) = !!A*X(K)-B(K)!!/(!!A!!*!!X(K)!!) . C C ------------------------------------------------------------------ C IF( E21 .EQ. 1.0E0 ) GO TO 75 MBWT = MBW MB = MBW/2 + 1 NORMA = 0.0E0 C DO 35 I = 1,N L = MAX0(MB - I + 1,1) M = MIN0(MB - I + N,MBWT) SUMA = 0.0E0 C DO 34 J = L,M SUMA = SUMA + ABS(A(I,J)) 34 CONTINUE C NORMA = AMAX1(NORMA,SUMA) 35 CONTINUE C DO 70 IPP = 1,IP NORMR = 0.0E0 NORMX = 0.0E0 DO 50 I = 1,N L = MAX0(MB - I + 1,1) M = MIN0(MB - I + N,MBWT) K = MAX0(-MB + 1 + I,1) SUM = -B(I,IPP) SUMA = 0.0E0 C DO 40 J = L,M SUM = SUM + A(I,J)*X(K,IPP) K = K + 1 SUMA = SUMA + ABS(A(I,J)) 40 CONTINUE C NORMR = NORMR + ABS(SUM) NORMX = NORMX + ABS(X(I,IPP)) 50 CONTINUE C IF( NORMX .EQ. 0.0E0 ) NORMX = 1.0E0 IF( NORMA .EQ. 0.0E0 ) NORMA = 1.0E0 NORM(IPP) = NORMR/(NORMA*NORMX) 70 CONTINUE C RETURN C 75 MB = MBW NORMA = 0.0E0 MB1 = MB - 1 C DO 110 I = 1,N SUMA = 0.0E0 IF( I .EQ. 1 ) GO TO 90 LSTART = MAX0(1,MB + 1 - I) C DO 80 L = LSTART,MB1 SUMA = SUMA + ABS(A(I,L)) 80 CONTINUE C 90 LSTOP = MIN0(MB,N + 1 - I) J = I C DO 100 L = 1,LSTOP L1 = MB + 1 - L SUMA = SUMA + ABS(A(J,L1)) J = J + 1 100 CONTINUE C NORMA = AMAX1(NORMA,SUMA) 110 CONTINUE C DO 170 I = 1,IP NORMR = 0.0E0 NORMX = 0.0E0 C DO 150 L = 1,N SUM = -B(L,I) J = MAX0(0,L - MB) IF( L .EQ. 1 ) GO TO 130 KSTART = MAX0(1,MB + 1 - L) C DO 120 K = KSTART,MB1 J = J + 1 SUM = SUM + A(L,K)*X(J,I) 120 CONTINUE C 130 KSTOP = MIN0(MB,N + 1 - L) C DO 140 K = 1,KSTOP J = J + 1 K1 = MB + 1 - K SUM = SUM + A(J,K1)*X(J,I) 140 CONTINUE C NORMR = NORMR + ABS(SUM) 150 CONTINUE C DO 160 K = 1,N NORMX = NORMX + ABS(X(K,I)) 160 CONTINUE C IF( NORMX .EQ. 0.0E0 ) NORMX = 1.0E0 IF( NORMA .EQ. 0.0E0 ) NORMA = 1.0E0 NORM(I) = NORMR/(NORMA*NORMX) 170 CONTINUE RETURN C END SUBROUTINE RMATIN(NM,NN,M,N,A,AHOLD,BHOLD,AORB,INITIL) C C THIS INPUT SUBROUTINE READS A REAL MATRIX, WITH DIMENSION C M BY N, FROM SYSIN. C TO GENERATE THE MATRIX A (OR B) INITIALLY, INITIL IS C SET TO 0. TO REGENERATE THE MATRIX A (OR B) FOR THE C PURPOSE OF THE RESIDUAL CALCULATION, INITIL IS SET TO 1. C C THIS ROUTINE IS CATALOGUED AS EISPDRV4(RLREADI). C INTEGER IA( 20),NM,NN,M,N,AORB,INITIL REAL A(NM,NN),AHOLD(NM,NN),BHOLD(NM,NM) DATA IREADA/1/,IWRITE/6/ C IF( INITIL .EQ. 1 ) GO TO 90 READ(IREADA,5) M,N 5 FORMAT(I6,6X,I6) IF( M .EQ. 0 ) GO TO 170 DO 20 I = 1,M READ(IREADA,10) (IA(J),J=1,N) 10 FORMAT(6I12) DO 15 J = 1,N A(I,J) = IA(J) 15 CONTINUE 20 CONTINUE IF( AORB .EQ. 1 ) GO TO 50 DO 40 J = 1,N DO 30 I = 1,M AHOLD(I,J) = A(I,J) 30 CONTINUE 40 CONTINUE GO TO 80 50 DO 70 J = 1,N DO 60 I = 1,M BHOLD(I,J) = A(I,J) 60 CONTINUE 70 CONTINUE 80 RETURN C 90 IF( AORB .EQ. 1 ) GO TO 120 DO 110 J = 1,N DO 100 I = 1,M A(I,J) = AHOLD(I,J) 100 CONTINUE 110 CONTINUE GO TO 150 120 DO 140 J = 1,N DO 130 I = 1,M A(I,J) = BHOLD(I,J) 130 CONTINUE 140 CONTINUE 150 RETURN C 170 WRITE(IWRITE,175) 175 FORMAT(45H0END OF DATA FOR SUBROUTINE RMATIN(RLREADI). /1H1) STOP END