PROGRAM SBLAS2 C THIS IS A TEST DRIVER FOR THE REAL LEVEL TWO BLAS. C REFERENCE: A PROPOSAL FOR AN EXTENDED SET OF FORTRAN BASIC C LINEAR ALGEBRA SUBPROGRAMS. C AUTHORS: J. J. DONGARRA, J. DU CROZ, S. HAMMARLING, R. J HANSON. C REVISED 860623 C REVISED YYMMDD C AUTH=R. J. HANSON, SANDIA NATIONAL LABS. C ALBUQUERQUE, NM 87185-5800 C PHASE 1. TEST THAT FORTRAN COMPILER MAKES DOT PRODUCT LOOPS RIGHT. C THIS IS A PARTIAL TEST ONLY. C PHASE 2. TEST THAT PACKAGE SUBPROGRAMS DO NOT CHANGE DATA WHEN C THEY GET WRONG OPTIONS OR INPUT. (PARTIAL TEST ONLY). C PHASE 3. TEST THAT PACKAGE SUBPROGRAMS GET APPROXIMATELY THE C RIGHT ANSWERS TO A SET OF PROBLEMS. INTEGER LX,LY,LZ,LW,N,I,J,NMESS,ISTAT,NUNIT(3),IERROR INTEGER IFATAL,NSUBS,IG C FLAGS FOR TESTING SUBPROGRAM SNAMES(I). -I SAYS NO TEST. C +I SAYS DO TEST. PARAMETER (NSUBS=16) REAL T,ZERO,ONE,HALF,SDIFF,AVIGR PARAMETER (ZERO=0.E0,HALF=.5E0,ONE=1.E0) PARAMETER (LX=2048,LY=LX,LZ=LX,LW=LX) INTEGER ITEST(NSUBS),ISNUM REAL SW(LW),SX(LX),SY(LY),SZ(LZ) LOGICAL LSE,SAME,FATAL,SFATAL C NAMES OF THE FILES USED IN THE TESTING. CHARACTER *128 DOPE(3) C NAMES OF THE SUBROUTINES TESTED. CHARACTER *6 SNAMES(NSUBS) EXTERNAL SDIFF C C THE LONGEST ALLOWED LENGTH OF FILE NAMES IS 'NCH': C NCH=128 DATA NUNIT/3*0/ C START OF PHASE 1. N = 32 DO 30 J = 1,N DO 10 I = 1,J SX(I) = I SY(I) = J - I + 1 10 CONTINUE C COMPUTE THE DOT PRODUCT OF SX(*) AND SY(*). THE ANSWER SHOULD C BE EXACTLY J(J+1)J/2 - (J+1)J(J-1)/3, FOR EACH J. T = ZERO DO 20 I = 1,J T = T + SX(I)*SY(I) 20 CONTINUE SW(J) = T C AN ERROR WILL BE FLAGGED UNLESS THE FOLLOWING EXPRESSION C IS EXACTLY REPRESENTED AS A FLOATING POINT NUMBER. SZ(J) = J* ((J+1)*J)/2 - ((J+1)*J* (J-1))/3 30 CONTINUE SAME = LSE(SW,SZ,N,1,N) IF (SAME) GO TO 40 NMESS = 1 IFATAL = 1 C DO (PRINT A MESSAGE) C SEND OUT A NOTE ABOUT A FAILURE IN THE ARITHMETIC OR COMPILER. GO TO 180 * 40 CONTINUE C COMPUTE A GRADE SCALE TO MEASURE ACCURACY. T = ONE IG = 0 50 CONTINUE IF (SDIFF(T+ONE,ONE).EQ.ZERO) GO TO 60 IG = IG + 1 T = T*HALF GO TO 50 * 60 CONTINUE C END OF PHASE 1. C START OF PHASE 2. C READ NAMES OF SUBPROGRAMS AND REQUIRED FILE NAMES. C THE NAME OF A FILE THAT CONTAINS THESE NAMES IS C READ FROM STANDARD INPUT. THIS NAME CAN BE 'NCH' C CHARACTERS IN LENGTH. NMESS = 6 IFATAL = 0 ASSIGN 70 TO IGOMES C DO (PRINT A MESSAGE) C PRINT A MESSAGE ABOUT OPENING THE NAME FILE AND UNIT NUMBERS. GO TO 180 * 70 CONTINUE READ (*,9001,END=190) DOPE(1) READ (*,*,END=190) NUNIT(1) C ISTAT=1 FOR 'OLD', =2 FOR 'NEW', =3 FOR 'UNKNOWN'. C THE RETURN FLAG IERROR=0 FOR SUCCESS, =1 FOR FAILURE. ISTAT = 1 CALL SOPEN(NUNIT(1),DOPE(1),ISTAT,IERROR) IF (IERROR.EQ.0) GO TO 90 NMESS = 3 IFATAL = 1 ASSIGN 80 TO IGOMES C DO (PRINT A MESSAGE) C PRINT AN ERROR MESSAGE ABOUT OPENING THE NAME FILE. GO TO 180 * 80 CONTINUE 90 CONTINUE C READ THE NAMES OF THE SUBPROGRAMS, THEIR TESTING FLAGS, AND C THE NAMES OF THE FILES FOR SUMMARIES. C THE NAMES AND UNIT NUMBERS IN THIS DATA FILE ARE MACHINE C DEPENDENT. SEE LINES 17.-20. THE NEXT LINE STARTS THE FILE. C 1. SGEMV 01 PUT -01 FOR NO TEST. SAME COLUMNS. C 2. SGBMV 02 PUT -02 FOR NO TEST. SAME COLUMNS. C 3. SSYMV 03 PUT -03 FOR NO TEST. SAME COLUMNS. C 4. SSBMV 04 PUT -04 FOR NO TEST. SAME COLUMNS. C 5. SSPMV 05 PUT -05 FOR NO TEST. SAME COLUMNS. C 6. STRMV 06 PUT -06 FOR NO TEST. SAME COLUMNS. C 7. STBMV 07 PUT -07 FOR NO TEST. SAME COLUMNS. C 8. STPMV 08 PUT -08 FOR NO TEST. SAME COLUMNS. C 9. STRSV 09 PUT -09 FOR NO TEST. SAME COLUMNS. C10. STBSV 10 PUT -10 FOR NO TEST. SAME COLUMNS. C11 STPSV 11 PUT -11 FOR NO TEST. SAME COLUMNS. C12. SGER 12 PUT -12 FOR NO TEST. SAME COLUMNS. C13. SSYR 13 PUT -13 FOR NO TEST. SAME COLUMNS. C14. SSPR 14 PUT -14 FOR NO TEST. SAME COLUMNS. C15. SSYR2 15 PUT -15 FOR NO TEST. SAME COLUMNS. C16. SSPR2 16 PUT -16 FOR NO TEST. SAME COLUMNS. C17. SNOOP C18. 71 NAME THEN UNIT OF SNAPSHOT FILE. PUT -NN FOR NO FILE. C19. SUMMRY C20. 72 NAME THEN UNIT OF SUMMARY FILE. PUT -NN FOR NO FILE. C21. T LOGICAL FLAG, T TO STOP ON FAILURES, F TO CONTINUE ON. C THE ABOVE LINE 21. IS THE LAST LINE OF THE FILE. DO 100 I = 1,NSUBS READ (NUNIT(1),9011,END=200) SNAMES(I),ITEST(I) 100 CONTINUE READ (NUNIT(1),9021,END=210) (DOPE(I),NUNIT(I),I=2,3) C READ THE FLAG THAT DIRECTS STOPPING ON ANY FAILURE. READ (NUNIT(1),9031,END=210) SFATAL C OPEN THE SUMMARY FILE. ISTAT = 3 CALL SOPEN(NUNIT(3),DOPE(3),ISTAT,IERROR) IF (IERROR.EQ.0) GO TO 110 NUNIT(3) = 0 NMESS = 8 IFATAL = 1 C DO (PRINT A MESSAGE) C PRINT AN ERROR MESSAGE ABOUT OPENING THE NAME FILE. GO TO 180 * 110 CONTINUE C TEST SGEMV, 01, AND SGBMV, 02. DO 120 ISNUM = 1,2 CALL SCHCK1(ITEST(ISNUM),SNAMES(ISNUM),IG,DOPE(2),NUNIT(2), . AVIGR,FATAL) IF (FATAL .AND. SFATAL) GO TO 170 120 CONTINUE C TEST SSYMV, 03, SSBMV, 04, AND SSPMV, 05. DO 130 ISNUM = 3,5 CALL SCHCK2(ITEST(ISNUM),SNAMES(ISNUM),IG,DOPE(2),NUNIT(2), . AVIGR,FATAL) IF (FATAL .AND. SFATAL) GO TO 170 130 CONTINUE C TEST STRMV, 06, STBMV, 07, STPMV, 08, C STRSV, 09, STBSV, 10, AND STPSV, 11. DO 140 ISNUM = 6,11 CALL SCHCK3(ITEST(ISNUM),SNAMES(ISNUM),IG,DOPE(2),NUNIT(2), . AVIGR,FATAL) IF (FATAL .AND. SFATAL) GO TO 170 140 CONTINUE C TEST SGER, 12. ISNUM = 12 CALL SCHCK4(ITEST(ISNUM),SNAMES(ISNUM),IG,DOPE(2),NUNIT(2),AVIGR, . FATAL) IF (FATAL .AND. SFATAL) GO TO 170 C TEST SSYR, 13, SSPR, 14, SSYR2, 15, AND SSPR2,16. DO 150 ISNUM = 13,16 CALL SCHCK5(ITEST(ISNUM),SNAMES(ISNUM),IG,DOPE(2),NUNIT(2), . AVIGR,FATAL) IF (FATAL .AND. SFATAL) GO TO 170 150 CONTINUE C GET RID OF FILE CONTAINING SNAPSHOTS OF INDIVIDUAL TESTS. C IF TESTING PROGRAM COMPLETES NORMALLY THIS IS NOT NEEDED. ISTAT = 1 CALL SOPEN(NUNIT(2),DOPE(2),ISTAT,IERROR) IF(NUNIT(2).GT.0)CLOSE (UNIT=NUNIT(2),STATUS='DELETE',ERR=160) C CLOSE FILE CONTAINING SUMMARY INFORMATION. 160 CONTINUE IF(NUNIT(3).GT.0)CLOSE (UNIT=NUNIT(3)) 170 CONTINUE STOP C PROCEDURE (PRINT A MESSAGE) 180 CONTINUE ISTAT = 3 IF (NUNIT(3).GT.0) CALL SOPEN(NUNIT(3),DOPE(3),ISTAT,IERROR) CALL SMESSG(NUNIT,3,NMESS) IF (IFATAL.EQ.1) STOP GO TO IGOMES * 190 CONTINUE NMESS = 2 IFATAL = 1 GO TO 180 * 200 CONTINUE NMESS = 4 IFATAL = 1 GO TO 180 * 210 CONTINUE NMESS = 5 IFATAL = 1 GO TO 180 * * LAST EXECUTABLE LINE OF SBLAS2. 9001 FORMAT (A) 9011 FORMAT (06X,A,T13,I3) 9021 FORMAT (06X,A,/,12X,I3) 9031 FORMAT (14X,L1) END SUBROUTINE SMESSG(NUNIT,IP,NMESS) C DEFINE THE TEXT OF ERROR MESSAGES. C THIS IS A TEST SUBPROGRAM FOR THE LEVEL TWO BLAS. C REVISED 860623 C REVISED YYMMDD C AUTH=R. J. HANSON, SANDIA NATIONAL LABS. LOGICAL EX,OP,INQ INTEGER NUNIT(3) CHARACTER *256 MESS(09) CHARACTER *256 M DATA MESS(01)/ .' THE COMPILER IS GENERATING BAD CODE FOR IN-LINE DOT PRODUCTS OR .IS INCORRECTLY EVALUATING THE ARITHMETIC EXPRESSIONS J*((J+1)*J)/2 . - (J+1)*J*(J-1)/3, J=1 THRU 32.'/ DATA MESS(02)/ .' ABNORMAL OR EARLY END-OF-FILE WHILE READING NAME OF FILE THAT CO .NTAINS THE NAMES OF THE SUBPROGRAMS AND THE SUMMARY FILES.'/ DATA MESS(03)/ .' THE ABOVE FILE NAME MUST BE PRESENT ON THE SYSTEM. IT IS NOT. .THIS FILE CONTAINS THE NAMES OF THE SUBPROGRAMS AND THE SUMMARY FI .LES.'/ DATA MESS(04)/ .' ABNORMAL OR EARLY END-OF-FILE WHILE READING NAMES OF SUBPROGRAMS . FROM THE ABOVE FILE NAME.'/ DATA MESS(05)/ .' ABNORNAL OR EARLY END-OF-FILE WHILE READING NAMES OF FILES FOR S .UMMARY OUTPUT.'/ DATA MESS(06)/ .' ENTER NAME AND UNIT NUMBER OF FILE CONTAINING NAMES OF SUBPROGRA .MS AND SUMMARY FILES. ONE ITEM PER LINE, PLEASE.'/ DATA MESS(07)/ .' THE SNAP-SHOT FILE OF ACTIVE TESTS CANNOT BE OPENED WITH ''NEW'' . STATUS OR IT CANNOT BE DELETED. THIS FILE SHOULD NOT BE PRESENT .ON THE SYSTEM.'/ DATA MESS(08)/ .' THE SUMMARY FILE OF ACTIVE TESTS CANNOT BE OPENED WITH ''UNKOWN' .' STATUS. THIS FILE SHOULD NOT BE PRESENT ON THE SYSTEM.'/ M = MESS(NMESS) NL = 256 NS = 72 INQ = .TRUE. DO 10 I = NL,1,-1 IF (ICHAR(M(I:I)).NE.ICHAR(' ')) GO TO 20 10 CONTINUE NL = 0 GO TO 30 * 20 NL = I C FOUND NS = POINTER TO LAST NONBLANK IN MESSAGE. 30 CONTINUE C NOW OUTPUT THE MESSAGE. PARSE IT SO THAT UP TO NS CHARS. PER LINE C PRINT, BUT DO NOT BREAK WORDS ACCROSS LINES. IS = 1 40 CONTINUE IE = MIN(NL,IS+NS) IF (IS.GE.IE) GO TO 70 50 CONTINUE IF (ICHAR(M(IE:IE)).EQ.ICHAR(' ') .OR. NL-IS.LT.NS) GO TO 60 IE = IE - 1 IF (IE.GT.IS) GO TO 50 60 CONTINUE IF (INQ) THEN INQUIRE (UNIT=NUNIT(IP),EXIST=EX,OPENED=OP) END IF C IF THE INTENDED UNIT IS NOT OPENED, SEND OUTPUT TO C STANDARD OUTPUT SO IT WILL BE SEEN. IF ( .NOT. OP .OR. .NOT. EX .OR. NUNIT(IP).EQ.0) THEN IF (IE.EQ.NL) THEN WRITE (*,'(A,/)') M(IS:IE) * ELSE WRITE (*,'(A)') M(IS:IE) END IF * INQ = .FALSE. * ELSE LUNIT = NUNIT(IP) WRITE (LUNIT,'(A)') M(IS:IE) END IF * IS = IE GO TO 40 * 70 CONTINUE RETURN END * SUBROUTINE SCHCK1(ISNUM,SNAME,IG,DOPE,NUNIT,AVIGR,FATAL) C THIS IS A TEST SUBPROGRAM FOR THE LEVEL TWO BLAS. C REVISED 860623 C REVISED YYMMDD C AUTH=R. J. HANSON, SANDIA NATIONAL LABS. C THIS PROGRAM HAS TWO PARTS. THE FIRST PART CHECKS TO SEE C IF ANY DATA GETS CHANGED WHEN NONE SHOULD. FOR EXAMPLE WHEN C USING AN INVALID OPTION OR NONPOSITVE PROBLEM DIMENSIONS. C THE SECOND PART CHECKS THAT THE RESULTS ARE REASONABLY ACCURATE. C DIMENSION AND PROBLEM SIZE DATA.. INTEGER INC(04),IDIM(08),NUNIT(2) REAL ALF(04),BET(04),SDIFF LOGICAL ISAME(13),LSE,FATAL,SAME,NCHNG,RESET CHARACTER *128 DOPE(2) CHARACTER *6 SNAME CHARACTER *3 ICH CHARACTER *1 ICHS,ICI INTEGER LA,LV PARAMETER (LA=4096,LV=4096,LMN=2048) REAL A(LA),AS(LA),X(LV),XS(LV) REAL Y(LV),YS(LV),YT(LMN),XT(LMN) REAL ALPHA,ALS,BETA,BLS,T,TRANSL,XN PARAMETER (ZERO=0.E0,HALF=.5E0,ONE=1.E0) COMMON /ARRAYS/AR,AS,X,XS,Y,YS,YT,XT EXTERNAL SDIFF * DATA ALF/-1.E0,2.E0,.3E0,1.E0/ DATA BET/-1.E0,0.E0,.3E0,1.E0/ DATA INC/-2,-1,1,2/ DATA IDIM/1,2,4,8,64,128,2048,0/ DATA ICH/'NT/'/ FATAL = .FALSE. C CHECK GENERAL MATRIX-VECTOR PRODUCT, Y = ALPHA*A*X+BETA*Y, NO.1-2. IF (ISNUM.LT.0) GO TO 220 NC = 0 RESET = .TRUE. AVIGR = ZERO IX = 0 10 IX = IX + 1 IF (IX.GT.4) GO TO 200 INCX = INC(IX) ALPHA = ALF(IX) IY = 0 20 IY = IY + 1 IF (IY.GT.4) GO TO 190 INCY = INC(IY) BETA = BET(IY) MM = 0 30 MM = MM + 1 IF (MM.GT.8) GO TO 180 M = IDIM(MM) NN = 0 40 NN = NN + 1 IF (NN.GT.8) GO TO 170 N = IDIM(NN) IC = 0 50 IC = IC + 1 IF (IC.GT.3) GO TO 160 IF (FATAL) GO TO 210 C SET DEFAULT BANDWIDTH SO PRINTING WILL BE OK. KL = MAX(0,M-1) KU = MAX(0,N-1) C DEFINE THE NUMBER OF ARGUMENTS AND THE Y ARGUMENT NUMBER. IF (ISNUM.EQ.1) THEN LDA = MAX(M,1) NARGS = 11 IYARG = 10 * ELSE IF (ISNUM.EQ.2) THEN NARGS = 13 IYARG = 12 C DEFINE BANDWIDTH OF MATRIX FOR TEST OF SGBMV. KL = MAX(0,MIN(M-1,M/2)) KU = MAX(0,MIN(N-1,N/2)) LDA = MAX(KL+KU+1,M) END IF * ICI = ICH(IC:IC) IF (ICHAR(ICI).EQ.ICHAR('T')) THEN ML = N NL = M INCCA = 1 INCRA = LDA * ELSE ML = M NL = N INCCA = LDA INCRA = 1 END IF * C IF NOT ENOUGH STORAGE, SKIP THIS CASE. (AVOID EXPLICT LDA*N). IF (SQRT(REAL(N))*SQRT(REAL(LDA)).GT.SQRT(REAL(LA))) GO TO 50 C DO (PREPARE NOTES FOR THIS TEST) C C PRINT A MESSAGE THAT GIVES DEBUGGING INFORMATION. THIS C MESSAGE SAYS.. C IN SUBPROGRAM XXXXX TESTS WERE ACTIVE WITH C OPTION = 'A' C M = IIII, N = IIII, C INCX = IIII, INCY = IIII, C KL = IIII, KU = IIII. C THE MAIN IDEA HERE IS THAT A SERIOUS BUG THAT OCCURS DURING THE C TESTING OF THESE SUBPROGRAMS MAY LOSE PROGRAM CONTROL. THIS C AUXILLIARY FILE CONTAINS THE DIMENSIONS THAT RESULTED IN THE LOSS C OF CONTROL. HENCE IT PROVIDES THE IMPLEMENTOR WITH MORE COMPLETE C INFORMATION ABOUT WHERE TO START TRACKING DOWN THE BUG. IF (NUNIT(1).GT.0) THEN C IF UNIT IS NOT AVAILABLE WITH 'NEW' STATUS, OPEN WITH C 'OLD' AND THEN DELETE IT. ISTAT = 1 CALL SOPEN(NUNIT(1),DOPE(1),ISTAT,IERROR) IF (IERROR.EQ.1) GO TO 60 C GET RID OF ANY OLD FILE. CLOSE (UNIT=NUNIT(1),STATUS='DELETE',ERR=60) 60 CONTINUE ISTAT = 2 C CREATE A NEW FILE FOR THE NEXT TEST. CALL SOPEN(NUNIT(1),DOPE(1),ISTAT,IERROR) IF (IERROR.EQ.0) GO TO 80 NMESS = 7 C DO (PRINT A MESSAGE) C PRINT AN ERROR MESSAGE ABOUT OPENING THE NAME FILE. CALL SMESSG(0,1,NMESS) FATAL = .TRUE. GO TO 210 * 80 CONTINUE WRITE (NUNIT(1),9001) SNAME,ICI,M,N,INCX,INCY,KL,KU C CLOSE THE FILE SO USEFUL STATUS INFORMATION IS SEALED. CLOSE (UNIT=NUNIT(1)) END IF C DO (DEFINE A SET OF PROBLEM DATA) ASSIGN 90 TO IGO3 GO TO 340 * 90 CONTINUE C DO (CALL SUBROUTINE) ASSIGN 100 TO IGO1 GO TO 280 * 100 CONTINUE IF (M.LE.0 .OR. N.LE.0 .OR. ICHAR(ICI).EQ.ICHAR('/')) THEN C DO (SEE WHAT DATA CHANGED INSIDE SUBROUTINES) ASSIGN 110 TO IGO2 GO TO 240 * 110 CONTINUE C IF DATA WAS INCORRECTLY CHANGED, MAKE NOTES AND RETURN. SAME = .TRUE. DO 120 I = 1,NARGS SAME = SAME .AND. ISAME(I) IF ( .NOT. ISAME(I)) THEN WRITE (NUNIT(2),9011) SNAME,I,ICI,M,N,INCX,INCY,KL,KU END IF * 120 CONTINUE IF ( .NOT. SAME) THEN FATAL = .TRUE. GO TO 210 * END IF * ELSE C DO (SEE WHAT DATA CHANGED INSIDE SUBROUTINES) ASSIGN 130 TO IGO2 GO TO 240 * 130 CONTINUE C IF DATA WAS INCORRECTLY CHANGED, MAKE NOTES AND RETURN. SAME = .TRUE. DO 140 I = 1,NARGS NCHNG = (I.EQ.IYARG .OR. ISAME(I)) SAME = SAME .AND. NCHNG IF ( .NOT. NCHNG) THEN WRITE (NUNIT(2),9021) SNAME,I,ICI,M,N,INCX,INCY,KL,KU END IF * 140 CONTINUE IF ( .NOT. SAME) THEN FATAL = .TRUE. GO TO 210 * END IF * NC = NC + 1 C DO (COMPUTE A CORRECT RESULT) ASSIGN 150 TO IGO4 GO TO 370 * 150 CONTINUE C IF GOT REALLY BAD ANSWER, PRINT NOTE AND GO BACK. IF (FATAL) GO TO 200 * END IF * GO TO 50 * 160 CONTINUE GO TO 40 * 170 CONTINUE GO TO 30 * 180 CONTINUE GO TO 20 * 190 CONTINUE GO TO 10 * 200 CONTINUE C REPORT ON ACCURACY OF DATA. WRITE (NUNIT(2),9031) ISNUM,SNAME,AVIGR,IG GO TO 230 * 210 CONTINUE WRITE (NUNIT(2),9041) ISNUM,SNAME GO TO 230 * 220 CONTINUE WRITE (NUNIT(2),9051) - ISNUM,SNAME 230 CONTINUE RETURN * 240 CONTINUE C PROCEDURE (SEE WHAT DATA CHANGED INSIDE SUBROUTINES) IF (ISNUM.EQ.1) THEN ISAME(1) = ICHAR(ICI) .EQ. ICHAR(ICHS) ISAME(2) = MS .EQ. M ISAME(3) = NS .EQ. N ISAME(4) = ALS .EQ. ALPHA ISAME(5) = .TRUE. IF (M.GT.0 .AND. N.GT.0) ISAME(5) = LSE(AS,A,M,N,LDA) ISAME(6) = LDAS .EQ. LDA ISAME(7) = .TRUE. IF (NL.GT.0 .AND. INCX.NE.0) ISAME(7) = LSE(XS,X,1,NL, . ABS(INCX)) ISAME(8) = INCXS .EQ. INCX ISAME(9) = BLS .EQ. BETA ISAME(10) = .TRUE. IF (ML.GT.0 .AND. INCY.NE.0) ISAME(10) = LSE(YS,Y,1,ML, . ABS(INCY)) ISAME(11) = INCYS .EQ. INCY * ELSE IF (ISNUM.EQ.2) THEN C COMPARE THE MATRIX IN THE SGBMV DATA STRUCTURE WITH C THE SAVED COPY. ISAME(1) = ICHAR(ICI) .EQ. ICHAR(ICHS) ISAME(2) = MS .EQ. M ISAME(3) = NS .EQ. N ISAME(4) = KLS .EQ. KL ISAME(5) = KUS .EQ. KU ISAME(6) = ALS .EQ. ALPHA ISAME(7) = .TRUE. IF (N.GT.0 .AND. M.GT.0) THEN DO 260 J = 1,N DO 250 I = MAX(1,J-KU),MIN(M,J+KL) IF (AS(1+ (I-1)+ (J-1)*LDA).NE. . A(1+ (KU+I-J)+ (J-1)*LDA)) THEN ISAME(7) = .FALSE. GO TO 270 * END IF * 250 CONTINUE 260 CONTINUE 270 CONTINUE END IF * ISAME(8) = LDAS .EQ. LDA ISAME(9) = .TRUE. IF (NL.GT.0 .AND. INCX.NE.0) ISAME(9) = LSE(XS,X,1,NL, . ABS(INCX)) ISAME(10) = INCXS .EQ. INCX ISAME(11) = BLS .EQ. BETA ISAME(12) = .TRUE. IF (ML.GT.0 .AND. INCY.NE.0) ISAME(12) = LSE(YS,Y,1,ML, . ABS(INCY)) ISAME(13) = INCYS .EQ. INCY END IF * GO TO IGO2 * 280 CONTINUE C PROCEDURE (CALL SUBROUTINE) C SAVE EVERY DATUM BEFORE THE CALL. ICHS = ICI MS = M NS = N KLS = KL KUS = KU ALS = ALPHA DO 290 I = 1,LDA*N AS(I) = A(I) 290 CONTINUE LDAS = LDA C SAVE COPY OF THE X AND Y VECTORS. IBX = 1 IF (INCX.LT.0) IBX = 1 + (1-NL)*INCX DO 300 J = 1,NL XS(IBX+ (J-1)*INCX) = X(IBX+ (J-1)*INCX) 300 CONTINUE INCXS = INCX BLS = BETA IBY = 1 IF (INCY.LT.0) IBY = 1 + (1-ML)*INCY DO 310 I = 1,ML YS(IBY+ (I-1)*INCY) = Y(IBY+ (I-1)*INCY) 310 CONTINUE INCYS = INCY IF (ISNUM.EQ.1) THEN CALL SGEMV(ICI,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY) * ELSE IF (ISNUM.EQ.2) THEN C TRANSFER THE MATRIX TO THE DATA STRUCTURE USED WITH SGBMV. DO 330 J = 1,N DO 320 I = MAX(1,J-KU),MIN(M,J+KL) A(1+ (KU+I-J)+ (J-1)*LDA) = AS(1+ (I-1)+ (J-1)*LDA) 320 CONTINUE 330 CONTINUE CALL SGBMV(ICI,M,N,KL,KU,ALPHA,A,LDA,X,INCX,BETA,Y,INCY) END IF * GO TO IGO1 * 340 CONTINUE C PROCEDURE (DEFINE A SET OF PROBLEM DATA) C DO NOTHING IF BOTH DIMENSIONS ARE NOT POSITIVE. IF (M.LE.0 .OR. N.LE.0) GO TO IGO3 TRANSL = ZERO CALL SMAKE(A,M,N,LDA,RESET,TRANSL) C TRIM AWAY ELEMENTS OUTSIDE THE BANDWIDTH FOR SGBMV. IF (ISNUM.EQ.2) THEN DO 360 J = 1,N DO 350 I = 1,M T = A(1+ (I-1)+ (J-1)*LDA) IF (J.GT.I .AND. J-I.GT.KU) T = ZERO IF (I.GT.J .AND. I-J.GT.KL) T = ZERO A(1+ (I-1)+ (J-1)*LDA) = T 350 CONTINUE 360 CONTINUE END IF * TRANSL = 500.E0 RESET = .FALSE. CALL SMAKE(X,1,NL,MAX(1,ABS(INCX)),RESET,TRANSL) IF (NL.GT.1 .AND. INCX.EQ.1) X(NL/2) = ZERO TRANSL = ZERO CALL SMAKE(Y,1,ML,MAX(1,ABS(INCY)),RESET,TRANSL) GO TO IGO3 * 370 CONTINUE C PROCEDURE (COMPUTE A CORRECT RESULT) C COMPUTE THE CONDITION NUMBER TO USE AS GAUGE FOR ACCURATE RESULTS. C THIS IS RETURNED IN XT(*). C COMPUTE THE APPROXIMATE CORRECT RESULT. C THIS IS RETURNED IN YT(*). IF (INCY.LT.0) THEN IBY = (1-ML)*INCY + 1 * ELSE IBY = 1 END IF * DO 390 I = 1,ML YT(I) = BETA*YS(IBY+ (I-1)*INCY) XT(I) = YS(IBY+ (I-1)*INCY)**2 IF (INCX.LT.0) THEN IBX = (1-NL)*INCX + 1 * ELSE IBX = 1 END IF * DO 380 J = 1,NL YT(I) = YT(I) + AS(1+ (I-1)*INCRA+ (J-1)*INCCA)*ALPHA* . XS(IBX+ (J-1)*INCX) XT(I) = XT(I) + AS(1+ (I-1)*INCRA+ (J-1)*INCCA)**2 380 CONTINUE XT(I) = SQRT(XT(I)) 390 CONTINUE XN = BETA**2 DO 400 J = 1,NL XN = XN + XS(IBX+ (J-1)*INCX)**2 400 CONTINUE XN = SQRT(XN) C COMPUTE THE GAUGES FOR THE RESULTS. DO 410 I = 1,ML XT(I) = XT(I)*XN 410 CONTINUE C COMPUTE THE DIFFERENCES. THEY SHOULD BE SMALL FOR CORRECT RESULTS. DO 420 I = 1,ML YT(I) = YT(I) - Y(IBY+ (I-1)*INCY) 420 CONTINUE C COMPUTE THE GRADE OF THIS RESULT. IGR = 0 T = ONE 430 CONTINUE C THIS TEST ALLOWS UP TO A LOSS OF FULL PRECISION BEFORE QUITTING. IF (IGR.GE.IG) GO TO 460 DO 440 I = 1,ML IF (SDIFF(T*ABS(YT(I))+XT(I),XT(I)).EQ.ZERO) GO TO 440 T = T*HALF IGR = IGR + 1 GO TO 430 * 440 CONTINUE C IF THE LOOP COMPLETES, ALL VALUES ARE 'SMALL.' THE VALUE IGR/IG C IS THE GRADE ASSIGNED. THE VALUE OF IGR IS MAXIMIZED OVER ALL THE C PROBLEMS. 450 CONTINUE AVIGR = MAX(AVIGR,REAL(IGR)) GO TO IGO4 * 460 CONTINUE FATAL = .TRUE. GO TO 450 * * LAST EXECUTABLE LINE OF SCHCK1 9001 FORMAT (' IN SUBPROGRAM ',A,/,' TESTS ACTIVE WITH OPTION = ',A,/, . ' M =',I4,', N = ',I4,/,' INCX = ',I2,', INCY = ',I2,/,' KL =', . I4,', KU =',I4) 9011 FORMAT (' IN SUBPROGRAM ',A,/,' ARGUMENT ',I2, . ' WAS CHANGED WITH INVALID INPUT.',/,' OPTION = ',A,', M =',I4, . ', N = ',I4,/,' INCX = ',I2,', INCY = ',I2,/,' KL =',I4, . ', KU =',I4) 9021 FORMAT (' IN SUBPROGRAM ',A,/,' ARGUMENT ',I2, . ' WAS CHANGED WHILE COMPUTING',/,' OPTION = ',A,', M =',I4, . ', N = ',I4,/,' INCX = ',I2,', INCY = ',I2,/,' KL =',I4, . ', KU =',I4) 9031 FORMAT (1X,I2,' SUBPROGRAM ',A,T24,'RECEIVED A LOSS GRADE OF ', . F5.2,' OUT OF ',I3) 9041 FORMAT (1X,I2,' SUBPROGRAM ',A,T24,'FAILED.') 9051 FORMAT (1X,I2,' SUBPROGRAM ',A,T24,'NOT TESTED.') END SUBROUTINE SCHCK2(ISNUM,SNAME,IG,DOPE,NUNIT,AVIGR,FATAL) C THIS IS A TEST SUBPROGRAM FOR THE LEVEL TWO BLAS. C TEST SSYMV, 03, SSBMV, 04, AND SSPMV, 05. C REVISED 860623 C REVISED YYMMDD C AUTH=R. J. HANSON, SANDIA NATIONAL LABS. C THIS PROGRAM HAS TWO PARTS. THE FIRST PART CHECKS TO SEE C IF ANY DATA GETS CHANGED WHEN NONE SHOULD. FOR EXAMPLE WHEN C USING AN INVALID OPTION OR NONPOSITVE PROBLEM DIMENSIONS. C THE SECOND PART CHECKS THAT THE RESULTS ARE REASONABLY ACCURATE. C DIMENSION AND PROBLEM SIZE DATA.. INTEGER INC(04),IDIM(06),NUNIT(2) REAL ALF(04),BET(04) LOGICAL ISAME(13),LSE,FATAL,SAME,NCHNG,RESET CHARACTER *128 DOPE(2) CHARACTER *6 SNAME CHARACTER *3 ICH CHARACTER *1 ICHS,ICI INTEGER LA,LV PARAMETER (LA=4096,LV=4096,LMN=2048) REAL ALPHA,ALS,BETA,BLS,T,TRANSL,XN REAL A(LA),AS(LA),X(LV),XS(LV) REAL Y(LV),YS(LV),YT(LMN),XT(LMN) PARAMETER (ZERO=0.E0,HALF=.5E0,ONE=1.E0) COMMON /ARRAYS/AR,AS,X,XS,Y,YS,YT,XT EXTERNAL SDIFF * DATA ALF/-1.E0,2.E0,.3E0,1.E0/ DATA BET/-1.E0,0.E0,.3E0,1.E0/ DATA INC/-2,-1,1,2/ DATA IDIM/1,2,4,8,64,0/ DATA ICH/'LU/'/ FATAL = .FALSE. C CHECK SYMMETRIC MATRIX-VECTOR PRODUCT, Y = ALPHA*A*X+BETA*Y, 3-5. IF (ISNUM.LT.0) GO TO 200 NC = 0 RESET = .TRUE. AVIGR = ZERO IX = 0 10 IX = IX + 1 IF (IX.GT.4) GO TO 180 INCX = INC(IX) ALPHA = ALF(IX) IY = 0 20 IY = IY + 1 IF (IY.GT.4) GO TO 170 INCY = INC(IY) BETA = BET(IY) NN = 0 30 NN = NN + 1 IF (NN.GT.6) GO TO 160 N = IDIM(NN) IC = 0 40 IC = IC + 1 IF (IC.GT.3) GO TO 150 IF (FATAL) GO TO 190 ICI = ICH(IC:IC) C DEFINE DEFAULT VALUE OF K SO PRINTING IS OK. K = MAX(0,N-1) C DEFINE THE NUMBER OF ARGUMENTS AND THE Y ARGUMENT NUMBER. LDA = MAX(N,1) IF (ISNUM.EQ.3) THEN NARGS = 10 IYARG = 09 * ELSE IF (ISNUM.EQ.4) THEN NARGS = 11 IYARG = 10 C DEFINE BANDWIDTH OF MATRIX FOR TEST OF SSBMV. K = INT(SQRT(REAL(N))+HALF) - 1 * ELSE IF (ISNUM.EQ.5) THEN NARGS = 9 IYARG = 8 END IF C DO (PREPARE NOTES FOR THIS TEST) C C PRINT A MESSAGE THAT GIVES DEBUGGING INFORMATION. THIS C MESSAGE SAYS.. C IN SUBPROGRAM XXXXX TESTS WERE ACTIVE WITH C OPTION = 'A' C N = IIII, C INCX = IIII, INCY = IIII, C K = IIII. C THE MAIN IDEA HERE IS THAT A SERIOUS BUG THAT OCCURS DURING THE C TESTING OF THESE SUBPROGRAMS MAY LOSE PROGRAM CONTROL. THIS C AUXILLIARY FILE CONTAINS THE DIMENSIONS THAT RESULTED IN THE LOSS C OF CONTROL. HENCE IT PROVIDES THE IMPLEMENTOR WITH MORE COMPLETE C INFORMATION ABOUT WHERE TO START TRACKING DOWN THE BUG. IF (NUNIT(1).GT.0) THEN C IF UNIT IS NOT AVAILABLE WITH 'NEW' STATUS, OPEN WITH C 'OLD' AND THEN DELETE IT. ISTAT = 1 CALL SOPEN(NUNIT(1),DOPE(1),ISTAT,IERROR) IF (IERROR.EQ.1) GO TO 50 C GET RID OF ANY OLD FILE. CLOSE (UNIT=NUNIT(1),STATUS='DELETE',ERR=50) 50 CONTINUE ISTAT = 2 C CREATE A NEW FILE FOR THE NEXT TEST. CALL SOPEN(NUNIT(1),DOPE(1),ISTAT,IERROR) IF (IERROR.EQ.0) GO TO 70 60 CONTINUE NMESS = 7 C DO (PRINT A MESSAGE) C PRINT AN ERROR MESSAGE ABOUT OPENING THE NAME FILE. CALL SMESSG(0,1,NMESS) FATAL = .TRUE. GO TO 190 * 70 CONTINUE WRITE (NUNIT(1),9001) SNAME,ICI,N,INCX,INCY,K C CLOSE THE FILE SO USEFUL STATUS INFORMATION IS SEALED. CLOSE (UNIT=NUNIT(1)) END IF C DO (DEFINE A SET OF PROBLEM DATA) ASSIGN 80 TO IGO3 GO TO 370 * 80 CONTINUE C DO (CALL SUBROUTINE) ASSIGN 90 TO IGO1 GO TO 290 * 90 CONTINUE IF (N.LE.0 .OR. ICHAR(ICI).EQ.ICHAR('/')) THEN C DO (SEE WHAT DATA CHANGED INSIDE SUBROUTINES) ASSIGN 100 TO IGO2 GO TO 220 * 100 CONTINUE C IF DATA WAS INCORRECTLY CHANGED, MAKE NOTES AND RETURN. SAME = .TRUE. DO 110 I = 1,NARGS SAME = SAME .AND. ISAME(I) IF ( .NOT. ISAME(I)) THEN WRITE (NUNIT(2),9011) SNAME,I,ICI,N,INCX,INCY,K END IF * 110 CONTINUE IF ( .NOT. SAME) THEN FATAL = .TRUE. GO TO 190 * END IF * ELSE C DO (SEE WHAT DATA CHANGED INSIDE SUBROUTINES) ASSIGN 120 TO IGO2 GO TO 220 * 120 CONTINUE C IF DATA WAS INCORRECTLY CHANGED, MAKE NOTES AND RETURN. SAME = .TRUE. DO 130 I = 1,NARGS NCHNG = (I.EQ.IYARG .OR. ISAME(I)) SAME = SAME .AND. NCHNG IF ( .NOT. NCHNG) THEN WRITE (NUNIT(2),9021) SNAME,I,ICI,N,INCX,INCY,K END IF * 130 CONTINUE IF ( .NOT. SAME) THEN FATAL = .TRUE. GO TO 190 * END IF * NC = NC + 1 C DO (COMPUTE A CORRECT RESULT) ASSIGN 140 TO IGO4 GO TO 420 * 140 CONTINUE C IF GOT REALLY BAD ANSWER, PRINT NOTE AND GO BACK. IF (FATAL) GO TO 180 * END IF * GO TO 40 * 150 CONTINUE GO TO 30 * 160 CONTINUE GO TO 20 * 170 CONTINUE GO TO 10 * 180 CONTINUE C REPORT ON ACCURACY OF DATA. WRITE (NUNIT(2),9031) ISNUM,SNAME,AVIGR,IG GO TO 210 * 190 CONTINUE WRITE (NUNIT(2),9041) ISNUM,SNAME GO TO 210 * 200 CONTINUE WRITE (NUNIT(2),9051) - ISNUM,SNAME 210 CONTINUE RETURN * 220 CONTINUE C PROCEDURE (SEE WHAT DATA CHANGED INSIDE SUBROUTINES) IF (ISNUM.EQ.3) THEN ISAME(1) = ICHAR(ICI) .EQ. ICHAR(ICHS) ISAME(2) = NS .EQ. N ISAME(3) = ALS .EQ. ALPHA ISAME(4) = .TRUE. IF (N.GT.0) ISAME(4) = LSE(AS,A,N,N,LDA) ISAME(5) = LDAS .EQ. LDA ISAME(6) = .TRUE. IF (N.GT.0 .AND. INCX.NE.0) ISAME(6) = LSE(XS,X,1,N,ABS(INCX)) ISAME(7) = INCXS .EQ. INCX ISAME(8) = BLS .EQ. BETA ISAME(9) = .TRUE. IF (N.GT.0 .AND. INCY.NE.0) ISAME(9) = LSE(YS,Y,1,N,ABS(INCY)) ISAME(10) = INCYS .EQ. INCY * ELSE IF (ISNUM.EQ.4) THEN C COMPARE THE MATRIX IN THE SSBMV AND SSPMV DATA STRUCTURES WITH C THE SAVED COPY. ISAME(1) = ICHAR(ICI) .EQ. ICHAR(ICHS) ISAME(2) = NS .EQ. N ISAME(3) = KS .EQ. K ISAME(4) = ALS .EQ. ALPHA ISAME(5) = .TRUE. C TEST THE MATRIX IN THE DATA STRUCTURE USED WITH SSBMV. IF (ICHAR(ICI).EQ.ICHAR('U')) THEN KOFF = K * ELSE KOFF = 0 END IF * IF (N.GT.0) THEN DO 240 J = 1,N DO 230 I = MAX(1,J-K),MIN(N,J+K) IF (AS(1+ (I-1)+ (J-1)*LDA).NE. . A(1+ (KOFF+I-J)+ (J-1)*LDA)) THEN ISAME(5) = .FALSE. GO TO 250 * END IF * 230 CONTINUE 240 CONTINUE 250 CONTINUE END IF * ISAME(6) = LDAS .EQ. LDA ISAME(7) = .TRUE. IF (N.GT.0 .AND. INCX.NE.0) ISAME(7) = LSE(XS,X,1,N,ABS(INCX)) ISAME(8) = INCXS .EQ. INCX ISAME(9) = BLS .EQ. BETA ISAME(10) = .TRUE. IF (N.GT.0 .AND. INCY.NE.0) ISAME(10) = LSE(YS,Y,1,N, . ABS(INCY)) ISAME(11) = INCYS .EQ. INCY * ELSE IF (ISNUM.EQ.5) THEN ISAME(1) = ICHAR(ICI) .EQ. ICHAR(ICHS) ISAME(2) = NS .EQ. N ISAME(3) = ALS .EQ. ALPHA ISAME(4) = .TRUE. C TEST THE MATRIX USING THE DATA STRUCTURE USED WITH SSPMV. IOFF = 0 DO 270 J = 1,N IF (ICHAR(ICI).EQ.ICHAR('U')) THEN ISTRT = 1 IEND = J * ELSE ISTRT = J IEND = N END IF * DO 260 I = ISTRT,IEND IOFF = IOFF + 1 IF (A(IOFF).NE.AS(1+ (I-1)+ (J-1)*LDA)) THEN ISAME(4) = .FALSE. GO TO 280 * END IF * 260 CONTINUE * 270 CONTINUE 280 CONTINUE ISAME(5) = .TRUE. IF (N.GT.0 .AND. INCX.NE.0) ISAME(5) = LSE(XS,X,1,N,ABS(INCX)) ISAME(6) = INCXS .EQ. INCX ISAME(7) = BLS .EQ. BETA ISAME(8) = .TRUE. IF (N.GT.0 .AND. INCY.NE.0) ISAME(8) = LSE(YS,Y,1,N,ABS(INCY)) ISAME(9) = INCYS .EQ. INCY END IF * GO TO IGO2 * 290 CONTINUE C PROCEDURE (CALL SUBROUTINE) C SAVE EVERY DATUM BEFORE THE CALL. ICHS = ICI NS = N KS = K ALS = ALPHA DO 300 I = 1,N*N AS(I) = A(I) 300 CONTINUE LDAS = LDA C SAVE COPY OF THE X AND Y VECTORS. IBX = 1 IF (INCX.LT.0) IBX = 1 + (1-N)*INCX DO 310 J = 1,N XS(IBX+ (J-1)*INCX) = X(IBX+ (J-1)*INCX) 310 CONTINUE INCXS = INCX BLS = BETA IBY = 1 IF (INCY.LT.0) IBY = 1 + (1-N)*INCY DO 320 I = 1,N YS(IBY+ (I-1)*INCY) = Y(IBY+ (I-1)*INCY) 320 CONTINUE INCYS = INCY IF (ISNUM.EQ.3) THEN CALL SSYMV(ICI,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY) * ELSE IF (ISNUM.EQ.4) THEN C TRANSFER THE MATRIX TO THE DATA STRUCTURE USED WITH SSBMV. IF (ICHAR(ICI).EQ.ICHAR('U')) THEN KOFF = K * ELSE KOFF = 0 END IF * DO 340 J = 1,N DO 330 I = MAX(1,J-K),MIN(N,J+K) A(1+ (KOFF+I-J)+ (J-1)*LDA) = AS(1+ (I-1)+ (J-1)*LDA) 330 CONTINUE 340 CONTINUE CALL SSBMV(ICI,N,K,ALPHA,A,LDA,X,INCX,BETA,Y,INCY) * ELSE IF (ISNUM.EQ.5) THEN C TRANSFER THE MATRIX TO THE DATA STRUCTURE USED WITH SSPMV. IOFF = 0 DO 360 J = 1,N IF (ICHAR(ICI).EQ.ICHAR('U')) THEN ISTRT = 1 IEND = J * ELSE ISTRT = J IEND = N END IF * DO 350 I = ISTRT,IEND IOFF = IOFF + 1 A(IOFF) = AS(1+ (I-1)+ (J-1)*LDA) 350 CONTINUE * 360 CONTINUE CALL SSPMV(ICI,N,ALPHA,A,X,INCX,BETA,Y,INCY) END IF * GO TO IGO1 * 370 CONTINUE C PROCEDURE (DEFINE A SET OF PROBLEM DATA) C DO NOTHING IF DIMENSIONS ARE NOT POSITIVE. IF (N.LE.0) GO TO IGO3 TRANSL = ZERO CALL SMAKE(A,N,N,LDA,RESET,TRANSL) C MAKE THE DATA MATRIX SYMMETRIC. DO 390 I = 1,N DO 380 J = I,N T = (A(1+ (I-1)+ (J-1)*LDA)+A(1+ (J-1)+ (I-1)*LDA))*HALF A(1+ (I-1)+ (J-1)*LDA) = T A(1+ (J-1)+ (I-1)*LDA) = T 380 CONTINUE 390 CONTINUE C TRIM AWAY ELEMENTS OUTSIDE THE BANDWIDTH FOR SSBMV. IF (ISNUM.EQ.4) THEN DO 410 J = 1,N DO 400 I = 1,N T = A(1+ (I-1)+ (J-1)*LDA) IF (J.GT.I .AND. J-I.GT.K) T = ZERO IF (I.GT.J .AND. I-J.GT.K) T = ZERO A(1+ (I-1)+ (J-1)*LDA) = T 400 CONTINUE 410 CONTINUE END IF * TRANSL = 500.E0 RESET = .FALSE. CALL SMAKE(X,1,N,MAX(1,ABS(INCX)),RESET,TRANSL) IF (N.GT.1 .AND. INCX.EQ.1) X(N/2) = ZERO TRANSL = ZERO CALL SMAKE(Y,1,N,MAX(1,ABS(INCY)),RESET,TRANSL) GO TO IGO3 * 420 CONTINUE C PROCEDURE (COMPUTE A CORRECT RESULT) C COMPUTE THE CONDITION NUMBER TO USE AS GAUGE FOR ACCURATE RESULTS. C THIS IS RETURNED IN XT(*). C COMPUTE THE APPROXIMATE CORRECT RESULT. C THIS IS RETURNED IN YT(*). IF (INCY.LT.0) THEN IBY = (1-N)*INCY + 1 * ELSE IBY = 1 END IF * DO 440 I = 1,N YT(I) = BETA*YS(IBY+ (I-1)*INCY) XT(I) = YS(IBY+ (I-1)*INCY)**2 IF (INCX.LT.0) THEN IBX = (1-N)*INCX + 1 * ELSE IBX = 1 END IF * DO 430 J = 1,N YT(I) = YT(I) + AS(1+ (I-1)+ (J-1)*LDA)*ALPHA* . XS(IBX+ (J-1)*INCX) XT(I) = XT(I) + AS(1+ (I-1)+ (J-1)*LDA)**2 430 CONTINUE XT(I) = SQRT(XT(I)) 440 CONTINUE XN = BETA**2 DO 450 J = 1,N XN = XN + XS(IBX+ (J-1)*INCX)**2 450 CONTINUE XN = SQRT(XN) C COMPUTE THE GAUGES FOR THE RESULTS. DO 460 I = 1,N XT(I) = XT(I)*XN 460 CONTINUE C COMPUTE THE DIFFERENCES. THEY SHOULD BE SMALL FOR CORRECT RESULTS. DO 470 I = 1,N YT(I) = YT(I) - Y(IBY+ (I-1)*INCY) 470 CONTINUE C COMPUTE THE GRADE OF THIS RESULT. IGR = 0 T = ONE 480 CONTINUE C THIS TEST ALLOWS UP TO A LOSS OF FULL PRECISION BEFORE QUITTING. IF (IGR.GT.IG) GO TO 510 DO 490 I = 1,N IF (SDIFF(T*ABS(YT(I))+XT(I),XT(I)).EQ.ZERO) GO TO 490 T = T*HALF IGR = IGR + 1 GO TO 480 * 490 CONTINUE C IF THE LOOP COMPLETES, ALL VALUES ARE 'SMALL.' THE VALUE IGR/IG C IS THE GRADE ASSIGNED. THE VALUE OF IGR IS MAXIMIZED OVER ALL THE C PROBLEMS. 500 CONTINUE AVIGR = MAX(AVIGR,REAL(IGR)) GO TO IGO4 * 510 CONTINUE FATAL = .TRUE. GO TO 500 * * LAST EXECUTABLE LINE OF SCHCK2 9001 FORMAT (' IN SUBPROGRAM ',A,/,' TESTS ACTIVE WITH OPTION = ',A,/, . ' N = ',I4,/,' INCX = ',I2,', INCY = ',I2,/,' K =',I4) 9011 FORMAT (' IN SUBPROGRAM ',A,/,' ARGUMENT ',I2, . ' WAS CHANGED WITH INVALID INPUT.',/,' OPTION = ',A,/,' N = ', . I4,/,' INCX = ',I2,', INCY = ',I2,/,' K = ',I4) 9021 FORMAT (' IN SUBPROGRAM ',A,/,' ARGUMENT ',I2, . ' WAS CHANGED WHILE COMPUTING',/,' OPTION = ',A,/,' N = ',I4,/, . ' INCX = ',I2,', INCY = ',I2,/,' K = ',I4) 9031 FORMAT (1X,I2,' SUBPROGRAM ',A,T24,'RECEIVED A LOSS GRADE OF ', . F5.2,' OUT OF ',I3) 9041 FORMAT (1X,I2,' SUBPROGRAM ',A,T24,'FAILED.') 9051 FORMAT (1X,I2,' SUBPROGRAM ',A,T24,'NOT TESTED.') END SUBROUTINE SCHCK3(ISNUM,SNAME,IG,DOPE,NUNIT,AVIGR,FATAL) C THIS IS A TEST SUBPROGRAM FOR THE LEVEL TWO BLAS. C TEST STRMV, 06, STBMV, 07, STPMV, 08, C STRSV, 09, STBSV, 10, AND STPSV, 11. C REVISED 860623 C REVISED YYMMDD C AUTH=R. J. HANSON, SANDIA NATIONAL LABS. C THIS PROGRAM HAS TWO PARTS. THE FIRST PART CHECKS TO SEE C IF ANY DATA GETS CHANGED WHEN NONE SHOULD. FOR EXAMPLE WHEN C USING AN INVALID OPTION OR NONPOSITVE PROBLEM DIMENSIONS. C THE SECOND PART CHECKS THAT THE RESULTS ARE REASONABLY ACCURATE. C DIMENSION AND PROBLEM SIZE DATA.. INTEGER INC(04),IDIM(06),NUNIT(2) LOGICAL ISAME(13),LSE,FATAL,SAME,NCHNG,RESET CHARACTER *128 DOPE(2) CHARACTER *6 SNAME CHARACTER *3 ICHI,ICHJ,ICHK CHARACTER *1 ICIU,ICIT,ICID CHARACTER *1 ICIUS,ICITS,ICIDS INTEGER LA,LV PARAMETER (LA=4096,LV=4096,LMN=2048) REAL A(LA),AS(LA),X(LV),XS(LV) REAL Y(LV),YS(LV),YT(LMN),XT(LMN) PARAMETER (ZERO=0.E0,HALF=.5E0,ONE=1.E0) COMMON /ARRAYS/AR,AS,X,XS,Y,YS,XT,YT EXTERNAL SDIFF * DATA INC/-2,-1,1,2/ DATA IDIM/1,2,4,8,64,0/ DATA ICHI/'LU/'/,ICHJ/'NT/'/,ICHK/'NU/'/ FATAL = .FALSE. C CHECK TRIANGULAR MATRIX-VECTOR PRODUCT, X = A*X, 6-8, C AND TRIANGULAR SOLVERS, 9-11. IF (ISNUM.LT.0) GO TO 180 NC = 0 RESET = .TRUE. AVIGR = ZERO IX = 0 10 IX = IX + 1 IF (IX.GT.4) GO TO 160 INCX = INC(IX) NN = 0 20 NN = NN + 1 IF (NN.GT.6) GO TO 150 N = IDIM(NN) IC = 0 30 IC = IC + 1 IF (IC.GT.3) GO TO 140 IF (FATAL) GO TO 170 ICIU = ICHI(IC:IC) ICIT = ICHJ(IC:IC) ICID = ICHK(IC:IC) C DEFINE DEFAULT VALUE OF K SO PRINTING IS OK. K = MAX(0,N-1) C DEFINE THE NUMBER OF ARGUMENTS AND THE X ARGUMENT NUMBER. LDA = MAX(N,1) IF (ICHAR(ICIT).EQ.ICHAR('T')) THEN INCRA = LDA INCCA = 1 * ELSE INCRA = 1 INCCA = LDA END IF * IF (ISNUM.EQ.6 .OR. ISNUM.EQ.9) THEN NARGS = 08 IXARG = 07 * ELSE IF (ISNUM.EQ.7 .OR. ISNUM.EQ.10) THEN NARGS = 09 IXARG = 08 C DEFINE BANDWIDTH OF MATRIX FOR TEST OF STBMV. K = INT(SQRT(REAL(N))+HALF) - 1 * ELSE IF (ISNUM.EQ.8 .OR. ISNUM.EQ.11) THEN NARGS = 07 IXARG = 06 END IF C DO (PREPARE NOTES FOR THIS TEST) C C PRINT A MESSAGE THAT GIVES DEBUGGING INFORMATION. THIS C MESSAGE SAYS.. C IN SUBPROGRAM XXXXX TESTS WERE ACTIVE WITH C OPTIONS = 'A' 'A' 'A' C N = IIII, C INCX = IIII C K = IIII. C THE MAIN IDEA HERE IS THAT A SERIOUS BUG THAT OCCURS DURING THE C TESTING OF THESE SUBPROGRAMS MAY LOSE PROGRAM CONTROL. THIS C AUXILLIARY FILE CONTAINS THE DIMENSIONS THAT RESULTED IN THE LOSS C OF CONTROL. HENCE IT PROVIDES THE IMPLEMENTOR WITH MORE COMPLETE C INFORMATION ABOUT WHERE TO START TRACKING DOWN THE BUG. IF (NUNIT(1).GT.0) THEN C IF UNIT IS NOT AVAILABLE WITH 'NEW' STATUS, OPEN WITH C 'OLD' AND THEN DELETE IT. ISTAT = 1 CALL SOPEN(NUNIT(1),DOPE(1),ISTAT,IERROR) IF (IERROR.EQ.1) GO TO 40 C GET RID OF ANY OLD FILE. CLOSE (UNIT=NUNIT(1),STATUS='DELETE',ERR=40) 40 CONTINUE ISTAT = 2 C CREATE A NEW FILE FOR THE NEXT TEST. CALL SOPEN(NUNIT(1),DOPE(1),ISTAT,IERROR) IF (IERROR.EQ.0) GO TO 60 50 CONTINUE NMESS = 7 C DO (PRINT A MESSAGE) C PRINT AN ERROR MESSAGE ABOUT OPENING THE NAME FILE. CALL SMESSG(0,1,NMESS) FATAL = .TRUE. GO TO 170 * 60 CONTINUE WRITE (NUNIT(1),9001) SNAME,ICIU,ICIT,ICID,N,INCX,K C CLOSE THE FILE SO USEFUL STATUS INFORMATION IS SEALED. CLOSE (UNIT=NUNIT(1)) END IF C DO (DEFINE A SET OF PROBLEM DATA) ASSIGN 70 TO IGO3 GO TO 330 * 70 CONTINUE C DO (CALL SUBROUTINE) ASSIGN 80 TO IGO1 GO TO 260 * 80 CONTINUE IF (N.LE.0 .OR. ICHAR(ICIU).EQ.ICHAR('/') .OR. ICHAR(ICIT).EQ. . ICHAR('/') .OR. ICHAR(ICID).EQ.ICHAR('/')) THEN C DO (SEE WHAT DATA CHANGED INSIDE SUBROUTINES) ASSIGN 90 TO IGO2 GO TO 200 * 90 CONTINUE C IF DATA WAS INCORRECTLY CHANGED, MAKE NOTES AND RETURN. SAME = .TRUE. DO 100 I = 1,NARGS SAME = SAME .AND. ISAME(I) IF ( .NOT. ISAME(I)) THEN WRITE (NUNIT(2),9011) SNAME,I,ICIU,ICIT,ICID,N,INCX,K END IF * 100 CONTINUE IF ( .NOT. SAME) THEN FATAL = .TRUE. GO TO 170 * END IF * ELSE C DO (SEE WHAT DATA CHANGED INSIDE SUBROUTINES) ASSIGN 110 TO IGO2 GO TO 200 * 110 CONTINUE C IF DATA WAS INCORRECTLY CHANGED, MAKE NOTES AND RETURN. SAME = .TRUE. DO 120 I = 1,NARGS NCHNG = (I.EQ.IXARG .OR. ISAME(I)) SAME = SAME .AND. NCHNG IF ( .NOT. NCHNG) THEN WRITE (NUNIT(2),9021) SNAME,I,ICIU,ICIT,ICID,N,INCX,K END IF * 120 CONTINUE IF ( .NOT. SAME) THEN FATAL = .TRUE. GO TO 170 * END IF * NC = NC + 1 C DO (COMPUTE A CORRECT RESULT) ASSIGN 130 TO IGO4 GO TO 380 * 130 CONTINUE C IF GOT REALLY BAD ANSWER, PRINT NOTE AND GO BACK. IF (FATAL) GO TO 160 * END IF * GO TO 30 * 140 CONTINUE GO TO 20 * 150 CONTINUE GO TO 10 * 160 CONTINUE C REPORT ON ACCURACY OF DATA. WRITE (NUNIT(2),9031) ISNUM,SNAME,AVIGR,IG GO TO 190 * 170 CONTINUE WRITE (NUNIT(2),9041) ISNUM,SNAME GO TO 190 * 180 CONTINUE WRITE (NUNIT(2),9051) - ISNUM,SNAME 190 CONTINUE RETURN * 200 CONTINUE C PROCEDURE (SEE WHAT DATA CHANGED INSIDE SUBROUTINES) ISAME(1) = ICHAR(ICIU) .EQ. ICHAR(ICIUS) ISAME(2) = ICHAR(ICIT) .EQ. ICHAR(ICITS) ISAME(3) = ICHAR(ICID) .EQ. ICHAR(ICIDS) ISAME(4) = NS .EQ. N IF (ISNUM.EQ.6 .OR. ISNUM.EQ.9) THEN ISAME(5) = .TRUE. IF (N.GT.0) ISAME(5) = LSE(AS,A,N,N,LDA) ISAME(6) = LDAS .EQ. LDA ISAME(7) = .TRUE. IF (N.GT.0) ISAME(7) = LSE(XS,X,1,N,ABS(INCX)) ISAME(8) = INCXS .EQ. INCX * ELSE IF (ISNUM.EQ.7 .OR. ISNUM.EQ.10) THEN C COMPARE THE MATRIX IN THE STBMV AND STPMV DATA STRUCTURES WITH C THE SAVED COPY. ISAME(5) = KS .EQ. K ISAME(6) = .TRUE. IF (N.GT.0) THEN DO 220 J = 1,N IF (ICHAR(ICIU).EQ.ICHAR('U')) THEN ISTRT = MAX(1,J-K) IEND = J * ELSE ISTRT = J IEND = MIN(N,J+K) END IF * DO 210 I = ISTRT,IEND IF (AS(1+ (I-1)+ (J-1)*LDA).NE. . A(1+ (KOFF+I-J)+ (J-1)*LDA)) THEN ISAME(6) = .FALSE. GO TO 230 * END IF * 210 CONTINUE 220 CONTINUE 230 CONTINUE END IF * ISAME(7) = LDAS .EQ. LDA ISAME(8) = .TRUE. IF (N.GT.0) ISAME(8) = LSE(XS,X,1,N,ABS(INCX)) ISAME(9) = INCXS .EQ. INCX * ELSE IF (ISNUM.EQ.8 .OR. ISNUM.EQ.11) THEN ISAME(5) = .TRUE. IOFF = 0 DO 250 J = 1,N IF (ICHAR(ICIU).EQ.ICHAR('U')) THEN ISTRT = 1 IEND = J * ELSE ISTRT = J IEND = N END IF * DO 240 I = ISTRT,IEND IOFF = IOFF + 1 IF (A(IOFF).NE.AS(1+ (I-1)+ (J-1)* . LDA)) ISAME(5) = .FALSE. 240 CONTINUE * 250 CONTINUE ISAME(6) = .TRUE. IF (N.GT.0) ISAME(6) = LSE(XS,X,1,N,ABS(INCX)) ISAME(7) = INCXS .EQ. INCX END IF * GO TO IGO2 * 260 CONTINUE C PROCEDURE (CALL SUBROUTINE) C SAVE EVERY DATUM BEFORE THE CALL. ICIUS = ICIU ICITS = ICIT ICIDS = ICID NS = N KS = K DO 270 I = 1,N*N AS(I) = A(I) 270 CONTINUE LDAS = LDA C SAVE COPY OF THE X VECTOR. IBX = 1 IF (INCX.LT.0) IBX = 1 + (1-N)*INCX DO 280 J = 1,N XS(IBX+ (J-1)*INCX) = X(IBX+ (J-1)*INCX) 280 CONTINUE INCXS = INCX IF (ISNUM.EQ.6) THEN CALL STRMV(ICIU,ICIT,ICID,N,A,LDA,X,INCX) * ELSE IF (ISNUM.EQ.9) THEN CALL STRSV(ICIU,ICIT,ICID,N,A,LDA,X,INCX) * ELSE IF (ISNUM.EQ.7 .OR. ISNUM.EQ.10) THEN C TRANSFER THE MATRIX TO THE DATA STRUCTURE USED WITH STBMV. IF (ICHAR(ICIU).EQ.ICHAR('U')) THEN KOFF = K * ELSE KOFF = 0 END IF * DO 300 J = 1,N IF (ICHAR(ICIU).EQ.ICHAR('U')) THEN ISTRT = MAX(1,J-K) IEND = J * ELSE ISTRT = J IEND = MIN(N,J+K) END IF * DO 290 I = ISTRT,IEND A(1+ (KOFF+I-J)+ (J-1)*LDA) = AS(1+ (I-1)+ (J-1)*LDA) 290 CONTINUE 300 CONTINUE IF (ISNUM.EQ.7) CALL STBMV(ICIU,ICIT,ICID,N,K,A,LDA,X,INCX) IF (ISNUM.EQ.10) CALL STBSV(ICIU,ICIT,ICID,N,K,A,LDA,X,INCX) * ELSE IF (ISNUM.EQ.8 .OR. ISNUM.EQ.11) THEN C TRANSFER THE MATRIX TO THE DATA STRUCTURE USED WITH STPMV. IOFF = 0 DO 320 J = 1,N IF (ICHAR(ICIU).EQ.ICHAR('U')) THEN ISTRT = 1 IEND = J * ELSE ISTRT = J IEND = N END IF * DO 310 I = ISTRT,IEND IOFF = IOFF + 1 A(IOFF) = AS(1+ (I-1)+ (J-1)*LDA) 310 CONTINUE * 320 CONTINUE IF (ISNUM.EQ.8) CALL STPMV(ICIU,ICIT,ICID,N,A,X,INCX) IF (ISNUM.EQ.11) CALL STPSV(ICIU,ICIT,ICID,N,A,X,INCX) END IF * GO TO IGO1 * 330 CONTINUE C PROCEDURE (DEFINE A SET OF PROBLEM DATA) C DO NOTHING IF DIMENSIONS ARE NOT POSITIVE. IF (N.LE.0) GO TO IGO3 TRANSL = ZERO CALL SMAKE(A,N,N,LDA,RESET,TRANSL) C MAKE THE DATA MATRIX TRIANGULAR. DO 350 I = 1,N DO 340 J = 1,N T = A(1+INCRA* (I-1)+ (J-1)*INCCA) S = A(1+INCRA* (J-1)+ (I-1)*INCCA) C SCALE TERMS SO THAT UNIT MATRICES ARE WELL-CONDITIONED. S = S/1000.E0 T = T/1000.E0 IF (ICHAR(ICIU).EQ.ICHAR('L') .AND. I.LT.J) T = ZERO IF (ICHAR(ICIU).EQ.ICHAR('U') .AND. I.GT.J) S = ZERO IF (ICHAR(ICID).EQ.ICHAR('U') .AND. I.EQ.J) THEN S = ONE T = ONE END IF * A(1+INCRA* (I-1)+ (J-1)*INCCA) = T A(1+INCRA* (J-1)+ (I-1)*INCCA) = S 340 CONTINUE 350 CONTINUE C TRIM AWAY ELEMENTS OUTSIDE THE BANDWIDTH FOR STBMV. IF (ISNUM.EQ.7 .OR. ISNUM.EQ.10) THEN DO 370 I = 1,N DO 360 J = 1,N T = A(1+INCRA* (I-1)+ (J-1)*INCCA) IF (J.GT.I .AND. J-I.GT.K) T = ZERO IF (I.GT.J .AND. I-J.GT.K) T = ZERO A(1+INCRA* (I-1)+ (J-1)*INCCA) = T 360 CONTINUE 370 CONTINUE END IF * TRANSL = 500.E0 RESET = .FALSE. CALL SMAKE(X,1,N,MAX(1,ABS(INCX)),RESET,TRANSL) IF (N.GT.1 .AND. INCX.EQ.1) X(N/2) = ZERO GO TO IGO3 * 380 CONTINUE C PROCEDURE (COMPUTE A CORRECT RESULT) C COMPUTE THE CONDITION NUMBER TO USE AS GAUGE FOR ACCURATE RESULTS. C THIS IS RETURNED IN XT(*). C COMPUTE THE APPROXIMATE CORRECT RESULT. C THIS IS RETURNED IN YT(*). DO 400 I = 1,N YT(I) = ZERO XT(I) = ZERO IF (INCX.LT.0) THEN IBX = (1-N)*INCX + 1 * ELSE IBX = 1 END IF * DO 390 J = 1,N T = XS(IBX+ (J-1)*INCX) IF (ISNUM.GE.9) T = X(IBX+ (J-1)*INCX) YT(I) = YT(I) + AS(1+ (I-1)*INCRA+ (J-1)*INCCA)*T XT(I) = XT(I) + AS(1+ (I-1)*INCRA+ (J-1)*INCCA)**2 390 CONTINUE XT(I) = SQRT(XT(I)) 400 CONTINUE XN = ZERO DO 410 J = 1,N T = XS(IBX+ (J-1)*INCX) IF (ISNUM.GE.9) T = X(IBX+ (J-1)*INCX) XN = XN + T**2 410 CONTINUE XN = SQRT(XN) C COMPUTE THE GAUGES FOR THE RESULTS. DO 420 I = 1,N XT(I) = XT(I)*XN 420 CONTINUE C COMPUTE THE DIFFERENCES. THEY SHOULD BE SMALL FOR CORRECT RESULTS. DO 430 I = 1,N T = X(IBX+ (I-1)*INCX) IF (ISNUM.GE.9) T = XS(IBX+ (I-1)*INCX) YT(I) = YT(I) - T 430 CONTINUE C COMPUTE THE GRADE OF THIS RESULT. IGR = 0 T = ONE 440 CONTINUE C THIS TEST ALLOWS UP TO A LOSS OF FULL PRECISION BEFORE QUITTING. IF (IGR.GE.IG) GO TO 470 DO 450 I = 1,N IF (SDIFF(T*ABS(YT(I))+XT(I),XT(I)).EQ.ZERO) GO TO 450 T = T*HALF IGR = IGR + 1 GO TO 440 * 450 CONTINUE C IF THE LOOP COMPLETES, ALL VALUES ARE 'SMALL.' THE VALUE IGR/IG C IS THE GRADE ASSIGNED. THE VALUE OF IGR IS MAXIMIZED OVER ALL THE C PROBLEMS. 460 CONTINUE AVIGR = MAX(AVIGR,REAL(IGR)) GO TO IGO4 * 470 CONTINUE FATAL = .TRUE. GO TO 460 * * LAST EXECUTABLE LINE OF SCHCK3 9001 FORMAT (' IN SUBPROGRAM ',A,/,' TESTS ACTIVE WITH OPTIONS = ', . 3 (A,2X),/,' N = ',I4,/,' INCX = ',I2,/,' K =',I4) 9011 FORMAT (' IN SUBPROGRAM ',A,/,' ARGUMENT ',I2, . ' WAS CHANGED WITH INVALID INPUT.',/,' OPTIONS = ',3 (A,2X),/, . ' N = ',I4,/,' INCX = ',I2,/,' K = ',I4) 9021 FORMAT (' IN SUBPROGRAM ',A,/,' ARGUMENT ',I2, . ' WAS CHANGED WHILE COMPUTING',/,' OPTIONS = ',3 (A,2X),/, . ' N = ',I4,/,' INCX = ',I2,/,' K = ',I4) 9031 FORMAT (1X,I2,' SUBPROGRAM ',A,T24,'RECEIVED A LOSS GRADE OF ', . F5.2,' OUT OF ',I3) 9041 FORMAT (1X,I2,' SUBPROGRAM ',A,T24,'FAILED.') 9051 FORMAT (1X,I2,' SUBPROGRAM ',A,T24,'NOT TESTED.') END SUBROUTINE SCHCK4(ISNUM,SNAME,IG,DOPE,NUNIT,AVIGR,FATAL) C THIS IS A TEST SUBPROGRAM FOR THE LEVEL TWO BLAS. C TEST SGER, 12. C REVISED 860623 C REVISED YYMMDD C AUTH=R. J. HANSON, SANDIA NATIONAL LABS. C THIS PROGRAM HAS TWO PARTS. THE FIRST PART CHECKS TO SEE C IF ANY DATA GETS CHANGED WHEN NONE SHOULD. FOR EXAMPLE WHEN C USING AN INVALID OPTION OR NONPOSITVE PROBLEM DIMENSIONS. C THE SECOND PART CHECKS THAT THE RESULTS ARE REASONABLY ACCURATE. C DIMENSION AND PROBLEM SIZE DATA.. INTEGER INC(04),IDIM(08),NUNIT(2) REAL ALF(04),SDIFF LOGICAL ISAME(13),LSE,FATAL,SAME,NCHNG,RESET CHARACTER *128 DOPE(2) CHARACTER *6 SNAME INTEGER LA,LV PARAMETER (LA=4096,LV=4096,LMN=2048) REAL A(LA),AS(LA),X(LV),XS(LV) REAL Y(LV),YS(LV),YT(LMN),XT(LMN) PARAMETER (ZERO=0.E0,HALF=.5E0,ONE=1.E0) COMMON /ARRAYS/AR,AS,X,XS,Y,YS,YT,XT EXTERNAL SDIFF * DATA ALF/-1.E0,2.E0,.3E0,1.E0/ DATA INC/-2,-1,1,2/ DATA IDIM/1,2,4,8,64,128,2048,0/ FATAL = .FALSE. C CHECK GENERAL RANK 1 UPDATE, 12. IF (ISNUM.LT.0) GO TO 200 NC = 0 RESET = .TRUE. AVIGR = ZERO IX = 0 10 IX = IX + 1 IF (IX.GT.4) GO TO 180 INCX = INC(IX) ALPHA = ALF(IX) IY = 0 20 IY = IY + 1 IF (IY.GT.4) GO TO 170 INCY = INC(IY) MM = 0 30 MM = MM + 1 IF (MM.GT.8) GO TO 160 M = IDIM(MM) NN = 0 40 NN = NN + 1 IF (NN.GT.8) GO TO 150 N = IDIM(NN) IF (FATAL) GO TO 190 ML = N NL = M INCCA = M INCRA = 1 C DEFINE THE NUMBER OF ARGUMENTS AND THE A ARGUMENT NUMBER. LDA = MAX(M,1) NARGS = 09 IAARG = 08 C IF NOT ENOUGH STORAGE, SKIP THIS CASE. (AVOID EXPLICT M*N). IF (SQRT(REAL(N))*SQRT(REAL(M)).GT.SQRT(REAL(LA))) GO TO 40 C DO (PREPARE NOTES FOR THIS TEST) C C PRINT A MESSAGE THAT GIVES DEBUGGING INFORMATION. THIS C MESSAGE SAYS.. C IN SUBPROGRAM XXXXX TESTS WERE ACTIVE WITH C M = IIII, N = IIII, C INCX = IIII, INCY = IIII, C THE MAIN IDEA HERE IS THAT A SERIOUS BUG THAT OCCURS DURING THE C TESTING OF THESE SUBPROGRAMS MAY LOSE PROGRAM CONTROL. THIS C AUXILLIARY FILE CONTAINS THE DIMENSIONS THAT RESULTED IN THE LOSS C OF CONTROL. HENCE IT PROVIDES THE IMPLEMENTOR WITH MORE COMPLETE C INFORMATION ABOUT WHERE TO START TRACKING DOWN THE BUG. IF (NUNIT(1).GT.0) THEN C IF UNIT IS NOT AVAILABLE WITH 'NEW' STATUS, OPEN WITH C 'OLD' AND THEN DELETE IT. ISTAT = 1 CALL SOPEN(NUNIT(1),DOPE(1),ISTAT,IERROR) IF (IERROR.EQ.1) GO TO 50 C GET RID OF ANY OLD FILE. CLOSE (UNIT=NUNIT(1),STATUS='DELETE',ERR=50) 50 CONTINUE ISTAT = 2 C CREATE A NEW FILE FOR THE NEXT TEST. CALL SOPEN(NUNIT(1),DOPE(1),ISTAT,IERROR) IF (IERROR.EQ.0) GO TO 70 60 CONTINUE NMESS = 7 C DO (PRINT A MESSAGE) C PRINT AN ERROR MESSAGE ABOUT OPENING THE NAME FILE. CALL SMESSG(0,1,NMESS) FATAL = .TRUE. GO TO 190 * 70 CONTINUE WRITE (NUNIT(1),9001) SNAME,M,N,INCX,INCY C CLOSE THE FILE SO USEFUL STATUS INFORMATION IS SEALED. CLOSE (UNIT=NUNIT(1)) END IF C DO (DEFINE A SET OF PROBLEM DATA) ASSIGN 80 TO IGO3 GO TO 270 * 80 CONTINUE C DO (CALL SUBROUTINE) ASSIGN 90 TO IGO1 GO TO 230 * 90 CONTINUE IF (M.LE.0 .OR. N.LE.0) THEN C DO (SEE WHAT DATA CHANGED INSIDE SUBROUTINES) ASSIGN 100 TO IGO2 GO TO 220 * 100 CONTINUE C IF DATA WAS INCORRECTLY CHANGED, MAKE NOTES AND RETURN. SAME = .TRUE. DO 110 I = 1,NARGS SAME = SAME .AND. ISAME(I) IF ( .NOT. ISAME(I)) THEN WRITE (NUNIT(2),9011) SNAME,I,M,N,INCX,INCY END IF * 110 CONTINUE IF ( .NOT. SAME) THEN FATAL = .TRUE. GO TO 190 * END IF * ELSE C DO (SEE WHAT DATA CHANGED INSIDE SUBROUTINES) ASSIGN 120 TO IGO2 GO TO 220 * 120 CONTINUE C IF DATA WAS INCORRECTLY CHANGED, MAKE NOTES AND RETURN. SAME = .TRUE. DO 130 I = 1,NARGS NCHNG = (I.EQ.IAARG .OR. ISAME(I)) SAME = SAME .AND. NCHNG IF ( .NOT. NCHNG) THEN WRITE (NUNIT(2),9021) SNAME,I,M,N,INCX,INCY END IF * 130 CONTINUE IF ( .NOT. SAME) THEN FATAL = .TRUE. GO TO 190 * END IF * NC = NC + 1 C DO (COMPUTE A CORRECT RESULT) ASSIGN 140 TO IGO4 GO TO 280 * 140 CONTINUE C IF GOT REALLY BAD ANSWER, PRINT NOTE AND GO BACK. IF (FATAL) GO TO 180 * END IF * GO TO 40 * 150 CONTINUE GO TO 30 * 160 CONTINUE GO TO 20 * 170 CONTINUE GO TO 10 * 180 CONTINUE C REPORT ON ACCURACY OF DATA. WRITE (NUNIT(2),9031) ISNUM,SNAME,AVIGR,IG GO TO 210 * 190 CONTINUE WRITE (NUNIT(2),9041) ISNUM,SNAME GO TO 210 * 200 CONTINUE WRITE (NUNIT(2),9051) - ISNUM,SNAME 210 CONTINUE RETURN * 220 CONTINUE C PROCEDURE (SEE WHAT DATA CHANGED INSIDE SUBROUTINES) ISAME(1) = MS .EQ. M ISAME(2) = NS .EQ. N ISAME(3) = ALS .EQ. ALPHA ISAME(4) = .TRUE. IF (NL.GT.0 .AND. INCX.NE.0) ISAME(4) = LSE(XS,X,1,NL,ABS(INCX)) ISAME(5) = INCXS .EQ. INCX ISAME(6) = .TRUE. IF (ML.GT.0 .AND. INCY.NE.0) ISAME(6) = LSE(YS,Y,1,ML,ABS(INCY)) ISAME(7) = INCYS .EQ. INCY ISAME(8) = .TRUE. IF (M.GT.0 .AND. N.GT.0) ISAME(8) = LSE(AS,A,M,N,LDA) ISAME(9) = LDAS .EQ. LDA * GO TO IGO2 * 230 CONTINUE C PROCEDURE (CALL SUBROUTINE) C SAVE EVERY DATUM BEFORE THE CALL. MS = M NS = N ALS = ALPHA DO 240 I = 1,M*N AS(I) = A(I) 240 CONTINUE LDAS = LDA C SAVE COPY OF THE X AND Y VECTORS. IBX = 1 IF (INCX.LT.0) IBX = 1 + (1-NL)*INCX DO 250 J = 1,NL XS(IBX+ (J-1)*INCX) = X(IBX+ (J-1)*INCX) 250 CONTINUE INCXS = INCX IBY = 1 IF (INCY.LT.0) IBY = 1 + (1-ML)*INCY DO 260 I = 1,ML YS(IBY+ (I-1)*INCY) = Y(IBY+ (I-1)*INCY) 260 CONTINUE INCYS = INCY CALL SGER(M,N,ALPHA,X,INCX,Y,INCY,A,LDA) * GO TO IGO1 * 270 CONTINUE C PROCEDURE (DEFINE A SET OF PROBLEM DATA) C DO NOTHING IF BOTH DIMENSIONS ARE NOT POSITIVE. IF (M.LE.0 .OR. N.LE.0) GO TO IGO3 TRANSL = ZERO CALL SMAKE(A,M,N,LDA,RESET,TRANSL) * TRANSL = 500.E0 RESET = .FALSE. CALL SMAKE(X,1,NL,MAX(1,ABS(INCX)),RESET,TRANSL) IF (NL.GT.1 .AND. INCX.EQ.1) X(NL/2) = ZERO TRANSL = ZERO CALL SMAKE(Y,1,ML,MAX(1,ABS(INCY)),RESET,TRANSL) GO TO IGO3 * 280 CONTINUE C PROCEDURE (COMPUTE A CORRECT RESULT) C COMPUTE THE CONDITION NUMBER TO USE AS GAUGE FOR ACCURATE RESULTS. C THIS IS RETURNED IN XT(*). C COMPUTE THE APPROXIMATE CORRECT RESULT. C THIS IS RETURNED IN YT(*), COLUMN BY COLUMN. IF (INCY.LT.0) THEN IBY = (1-ML)*INCY + 1 * ELSE IBY = 1 END IF * DO 340 J = 1,N DO 290 I = 1,M IF (INCX.LT.0) THEN IBX = (1-NL)*INCX + 1 * ELSE IBX = 1 END IF * YT(I) = AS(1+ (I-1)*INCRA+ (J-1)*INCCA) + . ALPHA*XS(IBX+ (I-1)*INCX)*YS(IBY+ (J-1)*INCY) XT(I) = AS(1+ (I-1)*INCRA+ (J-1)*INCCA)**2 + . ALPHA**2*XS(IBX+ (I-1)*INCX)**2* . YS(IBY+ (J-1)*INCY)**2 C COMPUTE THE GAUGES FOR THE RESULTS. XT(I) = SQRT(XT(I)) 290 CONTINUE C COMPUTE THE DIFFERENCES. THEY SHOULD BE SMALL FOR CORRECT RESULTS. DO 300 I = 1,M YT(I) = YT(I) - A(1+ (I-1)*INCRA+ (J-1)*INCCA) 300 CONTINUE C COMPUTE THE GRADE OF THIS RESULT. IGR = 0 T = ONE 310 CONTINUE C THIS TEST ALLOWS UP TO A LOSS OF FULL PRECISION BEFORE QUITTING. IF (IGR.GE.IG) GO TO 360 DO 320 I = 1,M IF (SDIFF(T*ABS(YT(I))+XT(I),XT(I)).EQ.ZERO) GO TO 320 T = T*HALF IGR = IGR + 1 GO TO 310 * 320 CONTINUE C IF THE LOOP COMPLETES, ALL VALUES ARE 'SMALL.' THE VALUE IGR/IG C IS THE GRADE ASSIGNED. THE VALUE OF IGR IS MAXIMIZED OVER ALL THE C PROBLEMS. 330 CONTINUE 340 CONTINUE 350 AVIGR = MAX(AVIGR,REAL(IGR)) GO TO IGO4 * 360 CONTINUE FATAL = .TRUE. GO TO 350 * * LAST EXECUTABLE LINE OF SCHCK4 9001 FORMAT (' IN SUBPROGRAM ',A,/,' M =',I4,', N = ',I4,/,' INCX = ', . I2,', INCY = ',I2) 9011 FORMAT (' IN SUBPROGRAM ',A,/,' ARGUMENT ',I2, . ' WAS CHANGED WITH INVALID INPUT.',/,' M =',I4,', N = ',I4,/, . ' INCX = ',I2,', INCY = ',I2) 9021 FORMAT (' IN SUBPROGRAM ',A,/,' ARGUMENT ',I2, . ' WAS CHANGED WHILE COMPUTING',/,' M =',I4,', N = ',I4,/, . ' INCX = ',I2,', INCY = ',I2) 9031 FORMAT (1X,I2,' SUBPROGRAM ',A,T24,'RECEIVED A LOSS GRADE OF ', . F5.2,' OUT OF ',I3) 9041 FORMAT (1X,I2,' SUBPROGRAM ',A,T24,'FAILED.') 9051 FORMAT (1X,I2,' SUBPROGRAM ',A,T24,'NOT TESTED.') END SUBROUTINE SCHCK5(ISNUM,SNAME,IG,DOPE,NUNIT,AVIGR,FATAL) C THIS IS A TEST SUBPROGRAM FOR THE LEVEL TWO BLAS. C TEST SSYR, 13, SSPR, 14, SSYR2, 15, AND SSPR2,16. C REVISED 860623 C REVISED YYMMDD C AUTH=R. J. HANSON, SANDIA NATIONAL LABS. C THIS PROGRAM HAS TWO PARTS. THE FIRST PART CHECKS TO SEE C IF ANY DATA GETS CHANGED WHEN NONE SHOULD. FOR EXAMPLE WHEN C USING AN INVALID OPTION OR NONPOSITVE PROBLEM DIMENSIONS. C THE SECOND PART CHECKS THAT THE RESULTS ARE REASONABLY ACCURATE. C DIMENSION AND PROBLEM SIZE DATA.. INTEGER INC(04),IDIM(06),NUNIT(2) REAL ALF(04) LOGICAL ISAME(13),LSE,FATAL,SAME,NCHNG,RESET CHARACTER *128 DOPE(2) CHARACTER *6 SNAME CHARACTER *3 ICH CHARACTER *1 ICHS,ICI INTEGER LA,LV PARAMETER (LA=4096,LV=4096,LMN=2048) REAL A(LA),AS(LA),X(LV),XS(LV) REAL Y(LV),YS(LV),YT(LMN),XT(LMN) PARAMETER (ZERO=0.E0,HALF=.5E0,ONE=1.E0) COMMON /ARRAYS/AR,AS,X,XS,Y,YS,YT,XT EXTERNAL SDIFF * DATA ALF/-1.E0,2.E0,.3E0,1.E0/ DATA INC/-2,-1,1,2/ DATA IDIM/1,2,4,8,64,0/ DATA ICH/'LU/'/ FATAL = .FALSE. C CHECK SYMMETRIC MATRIX RANK 1 AND RANK 2 UPDATES. IF (ISNUM.LT.0) GO TO 200 NC = 0 RESET = .TRUE. AVIGR = ZERO IX = 0 10 IX = IX + 1 IF (IX.GT.4) GO TO 180 INCX = INC(IX) ALPHA = ALF(IX) IY = 0 20 IY = IY + 1 IF (IY.GT.4) GO TO 170 INCY = INC(IY) NN = 0 30 NN = NN + 1 IF (NN.GT.6) GO TO 160 N = IDIM(NN) IC = 0 40 IC = IC + 1 IF (IC.GT.3) GO TO 150 IF (FATAL) GO TO 190 ICI = ICH(IC:IC) C DEFINE THE NUMBER OF ARGUMENTS AND THE Y ARGUMENT NUMBER. LDA = MAX(N,1) IF (ISNUM.EQ.13) THEN NARGS = 07 IAARG = 06 * ELSE IF (ISNUM.EQ.14) THEN NARGS = 06 IAARG = 06 * ELSE IF (ISNUM.EQ.15) THEN NARGS = 9 IAARG = 8 * ELSE IF (ISNUM.EQ.16) THEN NARGS = 8 IAARG = 8 END IF C DO (PREPARE NOTES FOR THIS TEST) C C PRINT A MESSAGE THAT GIVES DEBUGGING INFORMATION. THIS C MESSAGE SAYS.. C IN SUBPROGRAM XXXXX TESTS WERE ACTIVE WITH C OPTION = 'A' C N = IIII, C INCX = IIII, INCY = IIII, C THE MAIN IDEA HERE IS THAT A SERIOUS BUG THAT OCCURS DURING THE C TESTING OF THESE SUBPROGRAMS MAY LOSE PROGRAM CONTROL. THIS C AUXILLIARY FILE CONTAINS THE DIMENSIONS THAT RESULTED IN THE LOSS C OF CONTROL. HENCE IT PROVIDES THE IMPLEMENTOR WITH MORE COMPLETE C INFORMATION ABOUT WHERE TO START TRACKING DOWN THE BUG. IF (NUNIT(1).GT.0) THEN C IF UNIT IS NOT AVAILABLE WITH 'NEW' STATUS, OPEN WITH C 'OLD' AND THEN DELETE IT. ISTAT = 1 CALL SOPEN(NUNIT(1),DOPE(1),ISTAT,IERROR) IF (IERROR.EQ.1) GO TO 50 C GET RID OF ANY OLD FILE. CLOSE (UNIT=NUNIT(1),STATUS='DELETE',ERR=50) 50 CONTINUE ISTAT = 2 C CREATE A NEW FILE FOR THE NEXT TEST. CALL SOPEN(NUNIT(1),DOPE(1),ISTAT,IERROR) IF (IERROR.EQ.0) GO TO 70 60 CONTINUE NMESS = 7 C DO (PRINT A MESSAGE) C PRINT AN ERROR MESSAGE ABOUT OPENING THE NAME FILE. CALL SMESSG(0,1,NMESS) FATAL = .TRUE. GO TO 190 * 70 CONTINUE WRITE (NUNIT(1),9001) SNAME,ICI,N,INCX,INCY C CLOSE THE FILE SO USEFUL STATUS INFORMATION IS SEALED. CLOSE (UNIT=NUNIT(1)) END IF C DO (DEFINE A SET OF PROBLEM DATA) ASSIGN 80 TO IGO3 GO TO 370 * 80 CONTINUE C DO (CALL SUBROUTINE) ASSIGN 90 TO IGO1 GO TO 290 * 90 CONTINUE IF (N.LE.0 .OR. ICHAR(ICI).EQ.ICHAR('/')) THEN C DO (SEE WHAT DATA CHANGED INSIDE SUBROUTINES) ASSIGN 100 TO IGO2 GO TO 220 * 100 CONTINUE C IF DATA WAS INCORRECTLY CHANGED, MAKE NOTES AND RETURN. SAME = .TRUE. DO 110 I = 1,NARGS SAME = SAME .AND. ISAME(I) IF ( .NOT. ISAME(I)) THEN WRITE (NUNIT(2),9011) SNAME,I,ICI,N,INCX,INCY END IF * 110 CONTINUE IF ( .NOT. SAME) THEN FATAL = .TRUE. GO TO 190 * END IF * ELSE C DO (SEE WHAT DATA CHANGED INSIDE SUBROUTINES) ASSIGN 120 TO IGO2 GO TO 220 * 120 CONTINUE C IF DATA WAS INCORRECTLY CHANGED, MAKE NOTES AND RETURN. SAME = .TRUE. DO 130 I = 1,NARGS NCHNG = (I.EQ.IAARG .OR. ISAME(I)) SAME = SAME .AND. NCHNG IF ( .NOT. NCHNG) THEN WRITE (NUNIT(2),9021) SNAME,I,ICI,N,INCX,INCY END IF * 130 CONTINUE IF ( .NOT. SAME) THEN FATAL = .TRUE. GO TO 190 * END IF * NC = NC + 1 C DO (COMPUTE A CORRECT RESULT) ASSIGN 140 TO IGO4 GO TO 400 * 140 CONTINUE C IF GOT REALLY BAD ANSWER, PRINT NOTE AND GO BACK. IF (FATAL) GO TO 180 * END IF * GO TO 40 * 150 CONTINUE GO TO 30 * 160 CONTINUE IF (ISNUM.GE.15) GO TO 20 GO TO 10 * 170 CONTINUE GO TO 10 * 180 CONTINUE C REPORT ON ACCURACY OF DATA. WRITE (NUNIT(2),9031) ISNUM,SNAME,AVIGR,IG GO TO 210 * 190 CONTINUE WRITE (NUNIT(2),9041) ISNUM,SNAME GO TO 210 * 200 CONTINUE WRITE (NUNIT(2),9051) - ISNUM,SNAME 210 CONTINUE RETURN * 220 CONTINUE C PROCEDURE (SEE WHAT DATA CHANGED INSIDE SUBROUTINES) IF (ISNUM.EQ.13) THEN ISAME(1) = ICHAR(ICI) .EQ. ICHAR(ICHS) ISAME(2) = NS .EQ. N ISAME(3) = ALS .EQ. ALPHA ISAME(4) = .TRUE. IF (N.GT.0 .AND. INCX.NE.0) ISAME(4) = LSE(XS,X,1,N,ABS(INCX)) ISAME(5) = INCXS .EQ. INCX ISAME(6) = .TRUE. IF (N.GT.0) ISAME(6) = LSE(AS,A,N,N,LDA) ISAME(7) = LDAS .EQ. LDA * ELSE IF (ISNUM.EQ.14) THEN C COMPARE THE MATRIX IN THE DATA STRUCTURES WITH THE SAVED COPY. ISAME(1) = ICHAR(ICI) .EQ. ICHAR(ICHS) ISAME(2) = NS .EQ. N ISAME(3) = ALS .EQ. ALPHA ISAME(4) = .TRUE. IF (N.GT.0 .AND. INCX.NE.0) ISAME(4) = LSE(XS,X,1,N,ABS(INCX)) ISAME(5) = INCXS .EQ. INCX ISAME(6) = .TRUE. IOFF = 0 DO 240 J = 1,N IF (ICHAR(ICI).EQ.ICHAR('U')) THEN ISTRT = 1 IEND = J * ELSE ISTRT = J IEND = N END IF * DO 230 I = ISTRT,IEND IOFF = IOFF + 1 IF (A(IOFF).NE.AS(1+ (I-1)+ (J-1)*LDA)) THEN ISAME(6) = .FALSE. GO TO 250 * END IF * 230 CONTINUE 240 CONTINUE 250 CONTINUE * ELSE IF (ISNUM.EQ.15) THEN ISAME(1) = ICHAR(ICI) .EQ. ICHAR(ICHS) ISAME(2) = NS .EQ. N ISAME(3) = ALS .EQ. ALPHA ISAME(4) = .TRUE. IF (N.GT.0 .AND. INCX.NE.0) ISAME(4) = LSE(XS,X,1,N,ABS(INCX)) ISAME(5) = INCXS .EQ. INCX ISAME(6) = .TRUE. IF (N.GT.0 .AND. INCY.NE.0) ISAME(6) = LSE(YS,Y,1,N,ABS(INCY)) ISAME(7) = INCYS .EQ. INCY ISAME(8) = .TRUE. IF (N.GT.0) ISAME(8) = LSE(AS,A,N,N,LDA) ISAME(9) = LDAS .EQ. LDA * ELSE IF (ISNUM.EQ.16) THEN ISAME(1) = ICHAR(ICI) .EQ. ICHAR(ICHS) ISAME(2) = NS .EQ. N ISAME(3) = ALS .EQ. ALPHA ISAME(4) = .TRUE. IF (N.GT.0 .AND. INCX.NE.0) ISAME(4) = LSE(XS,X,1,N,ABS(INCX)) ISAME(5) = INCXS .EQ. INCX ISAME(6) = .TRUE. IF (N.GT.0 .AND. INCY.NE.0) ISAME(6) = LSE(YS,Y,1,N,ABS(INCY)) ISAME(7) = INCYS .EQ. INCY ISAME(8) = .TRUE. IOFF = 0 DO 270 J = 1,N IF (ICHAR(ICI).EQ.ICHAR('U')) THEN ISTRT = 1 IEND = J * ELSE ISTRT = J IEND = N END IF * DO 260 I = ISTRT,IEND IOFF = IOFF + 1 IF (A(IOFF).NE.AS(1+ (I-1)+ (J-1)*LDA)) THEN ISAME(8) = .FALSE. GO TO 280 * END IF * 260 CONTINUE 270 CONTINUE 280 CONTINUE END IF * GO TO IGO2 * 290 CONTINUE C PROCEDURE (CALL SUBROUTINE) C SAVE EVERY DATUM BEFORE THE CALL. ICHS = ICI NS = N ALS = ALPHA DO 300 I = 1,N*N AS(I) = A(I) 300 CONTINUE LDAS = LDA C SAVE COPY OF THE X AND Y VECTORS. IBX = 1 IF (INCX.LT.0) IBX = 1 + (1-N)*INCX DO 310 J = 1,N XS(IBX+ (J-1)*INCX) = X(IBX+ (J-1)*INCX) 310 CONTINUE INCXS = INCX IBY = 1 IF (INCY.LT.0) IBY = 1 + (1-N)*INCY DO 320 I = 1,N YS(IBY+ (I-1)*INCY) = Y(IBY+ (I-1)*INCY) 320 CONTINUE INCYS = INCY IF (ISNUM.EQ.13) THEN CALL SSYR(ICI,N,ALPHA,X,INCX,A,LDA) * ELSE IF (ISNUM.EQ.14) THEN C TRANSFER THE MATRIX TO THE DATA STRUCTURE USED WITH SSPR. IOFF = 0 DO 340 J = 1,N IF (ICHAR(ICI).EQ.ICHAR('U')) THEN ISTRT = 1 IEND = J * ELSE ISTRT = J IEND = N END IF * DO 330 I = ISTRT,IEND IOFF = IOFF + 1 A(IOFF) = AS(1+ (I-1)+ (J-1)*LDA) 330 CONTINUE * 340 CONTINUE CALL SSPR(ICI,N,ALPHA,X,INCX,A) * ELSE IF (ISNUM.EQ.15) THEN * CALL SSYR2(ICI,N,ALPHA,X,INCX,Y,INCY,A,LDA) * ELSE IF (ISNUM.EQ.16) THEN C TRANSFER THE MATRIX TO THE DATA STRUCTURE USED WITH SSPR2. IOFF = 0 DO 360 J = 1,N IF (ICHAR(ICI).EQ.ICHAR('U')) THEN ISTRT = 1 IEND = J * ELSE ISTRT = J IEND = N END IF * DO 350 I = ISTRT,IEND IOFF = IOFF + 1 A(IOFF) = AS(1+ (I-1)+ (J-1)*LDA) 350 CONTINUE * 360 CONTINUE CALL SSPR2(ICI,N,ALPHA,X,INCX,Y,INCY,A) END IF * GO TO IGO1 * 370 CONTINUE C PROCEDURE (DEFINE A SET OF PROBLEM DATA) C DO NOTHING IF DIMENSIONS ARE NOT POSITIVE. IF (N.LE.0) GO TO IGO3 TRANSL = ZERO CALL SMAKE(A,N,N,LDA,RESET,TRANSL) C MAKE THE DATA MATRIX SYMMETRIC. DO 390 I = 1,N DO 380 J = I,N T = (A(1+ (I-1)+ (J-1)*LDA)+A(1+ (J-1)+ (I-1)*LDA))*HALF A(1+ (I-1)+ (J-1)*LDA) = T A(1+ (J-1)+ (I-1)*LDA) = T 380 CONTINUE 390 CONTINUE * TRANSL = 500.E0 RESET = .FALSE. CALL SMAKE(X,1,N,MAX(1,ABS(INCX)),RESET,TRANSL) IF (N.GT.1 .AND. INCX.EQ.1) X(N/2) = ZERO TRANSL = ZERO CALL SMAKE(Y,1,N,MAX(1,ABS(INCY)),RESET,TRANSL) GO TO IGO3 * 400 CONTINUE C PROCEDURE (COMPUTE A CORRECT RESULT) C COMPUTE THE CONDITION NUMBER TO USE AS GAUGE FOR ACCURATE RESULTS. C THIS IS RETURNED IN XT(*). C COMPUTE THE APPROXIMATE CORRECT RESULT. IF (ISNUM.EQ.13 .OR. ISNUM.EQ.14) THEN IF (INCX.LT.0) THEN IBX = (1-N)*INCX + 1 * ELSE IBX = 1 END IF * IOFF = 0 DO 450 J = 1,N IF (ICHAR(ICI).EQ.ICHAR('U')) THEN ISTRT = 1 IEND = J * ELSE ISTRT = J IEND = N END IF * DO 410 I = ISTRT,IEND YT(I) = AS(1+ (I-1)+ (J-1)*LDA) + . ALPHA*XS(IBX+ (J-1)*INCX)*XS(IBX+ (I-1)*INCX) XT(I) = AS(1+ (I-1)+ (J-1)*LDA)**2 + . ALPHA**2*XS(IBX+ (I-1)*INCX)**2* . XS(IBX+ (J-1)*INCX)**2 410 CONTINUE C COMPUTE THE DIFFERENCES. THEY SHOULD BE SMALL FOR CORRECT RESULTS. DO 420 I = ISTRT,IEND XT(I) = SQRT(XT(I)) IF (ISNUM.EQ.13) THEN YT(I) = YT(I) - A(1+ (I-1)+ (J-1)*LDA) * ELSE IF (ISNUM.EQ.14) THEN IOFF = IOFF + 1 YT(I) = YT(I) - A(IOFF) END IF * 420 CONTINUE C COMPUTE THE GRADE OF THIS RESULT. IGR = 0 T = ONE DO 440 I = ISTRT,IEND 430 CONTINUE C THIS TEST ALLOWS UP TO A LOSS OF FULL PRECISION BEFORE QUITTING. IF (IGR.GE.IG) GO TO 520 IF (SDIFF(T*ABS(YT(I))+XT(I),XT(I)).EQ.ZERO) GO TO 440 T = T*HALF IGR = IGR + 1 GO TO 430 * C IF THE LOOP COMPLETES, ALL VALUES ARE 'SMALL.' THE VALUE IGR/IG C IS THE GRADE ASSIGNED. THE VALUE OF IGR IS MAXIMIZED OVER ALL THE C PROBLEMS. 440 CONTINUE 450 CONTINUE * ELSE IF (ISNUM.EQ.15 .OR. ISNUM.EQ.16) THEN IF (INCX.LT.0) THEN IBX = (1-N)*INCX + 1 * ELSE IBX = 1 END IF * IF (INCY.LT.0) THEN IBY = (1-N)*INCY + 1 * ELSE IBY = 1 END IF * IOFF = 0 DO 500 J = 1,N IF (ICHAR(ICI).EQ.ICHAR('U')) THEN ISTRT = 1 IEND = J * ELSE ISTRT = J IEND = N END IF * DO 460 I = ISTRT,IEND YT(I) = AS(1+ (I-1)+ (J-1)*LDA) + . ALPHA*XS(IBX+ (J-1)*INCX)*YS(IBY+ (I-1)*INCY) + . ALPHA*XS(IBX+ (I-1)*INCX)*YS(IBY+ (J-1)*INCY) XT(I) = AS(1+ (I-1)+ (J-1)*LDA)**2 + . ALPHA**2*XS(IBX+ (I-1)*INCX)**2* . YS(IBY+ (J-1)*INCY)**2 + . ALPHA**2*XS(IBX+ (J-1)*INCX)**2* . YS(IBY+ (I-1)*INCY)**2 460 CONTINUE C COMPUTE THE DIFFERENCES. THEY SHOULD BE SMALL FOR CORRECT RESULTS. DO 470 I = ISTRT,IEND XT(I) = SQRT(XT(I)) IF (ISNUM.EQ.15) THEN YT(I) = YT(I) - A(1+ (I-1)+ (J-1)*LDA) * ELSE IF (ISNUM.EQ.16) THEN IOFF = IOFF + 1 YT(I) = YT(I) - A(IOFF) END IF * 470 CONTINUE C COMPUTE THE GRADE OF THIS RESULT. IGR = 0 T = ONE DO 490 I = ISTRT,IEND 480 CONTINUE C THIS TEST ALLOWS UP TO A LOSS OF FULL PRECISION BEFORE QUITTING. IF (IGR.GE.IG) GO TO 520 IF (SDIFF(T*ABS(YT(I))+XT(I),XT(I)).EQ.ZERO) GO TO 490 T = T*HALF IGR = IGR + 1 GO TO 480 * C IF THE LOOP COMPLETES, ALL VALUES ARE 'SMALL.' THE VALUE IGR/IG C IS THE GRADE ASSIGNED. THE VALUE OF IGR IS MAXIMIZED OVER ALL THE C PROBLEMS. 490 CONTINUE 500 CONTINUE END IF * 510 CONTINUE AVIGR = MAX(AVIGR,REAL(IGR)) GO TO IGO4 * 520 CONTINUE FATAL = .TRUE. GO TO 510 * * LAST EXECUTABLE LINE OF SCHCK5 9001 FORMAT (' IN SUBPROGRAM ',A,/,' TESTS ACTIVE WITH OPTION = ',A,/, . ' N = ',I4,/,' INCX = ',I2,', INCY = ',I2) 9011 FORMAT (' IN SUBPROGRAM ',A,/,' ARGUMENT ',I2, . ' WAS CHANGED WITH INVALID INPUT.',/,' OPTION = ',A,/,' N = ', . I4,/,' INCX = ',I2,', INCY = ',I2) 9021 FORMAT (' IN SUBPROGRAM ',A,/,' ARGUMENT ',I2, . ' WAS CHANGED WHILE COMPUTING',/,' OPTION = ',A,/,' N = ',I4,/, . ' INCX = ',I2,', INCY = ',I2) 9031 FORMAT (1X,I2,' SUBPROGRAM ',A,T24,'RECEIVED A LOSS GRADE OF ', . F5.2,' OUT OF ',I3) 9041 FORMAT (1X,I2,' SUBPROGRAM ',A,T24,'FAILED.') 9051 FORMAT (1X,I2,' SUBPROGRAM ',A,T24,'NOT TESTED.') END SUBROUTINE SMAKE(A,M,N,LDA,RESET,TRANS) C GENERATE VALUES FOR AN M BY N MATRIX A. C RESET THE GENERATOR IF FLAG RESET = .TRUE. C TRANSLATE THE VALUES WITH TRANS. C THIS IS A TEST SUBPROGRAM FOR THE LEVEL TWO BLAS. C REVISED 860623 C REVISED YYMMDD C AUTH=R. J. HANSON, SANDIA NATIONAL LABS. REAL A(LDA,*),TRANS,ANOISE REAL ZERO,HALF,ONE PARAMETER (ZERO=0.E0,HALF=.5E0,ONE=1.E0,THREE=3.E0) LOGICAL RESET IF (RESET) THEN ANOISE = -ONE ANOISE = SBEG(ANOISE) ANOISE = ZERO END IF * IC = 0 DO 20 I = 1,M DO 10 J = 1,N IC = IC + 1 C BREAK UP PERIODICITIES THAT ARE MULTIPLES OF 5. IF (MOD(IC,5).EQ.0) A(I,J) = SBEG(ANOISE) A(I,J) = SBEG(ANOISE) - TRANS C HERE THE PERTURBATION IN THE LAST BIT POSITION IS MADE. A(I,J) = A(I,J) + ONE/THREE ANOISE = 0.E0 10 CONTINUE 20 CONTINUE RETURN * LAST EXECUTABLE LINE OF SMAKE END SUBROUTINE SOPEN(IUNIT,NAME,ISTAT,IERROR) C OPEN UNIT IUNIT WITH FILE NAMED NAME. C ISTAT=1 FOR 'OLD', =2 FOR 'NEW', =3 FOR 'UNKNOWN'. C THE RETURN FLAG IERROR=0 FOR SUCCESS, =1 FOR FAILURE. C A BAD VALUE OF ISTAT CAN ALSO INDICATE FAILURE. C THIS IS A TEST SUBPROGRAM FOR THE LEVEL TWO BLAS. C REVISED 860623 C REVISED YYMMDD C AUTH=R. J. HANSON, SANDIA NATIONAL LABS. CHARACTER * (*) NAME IF (ISTAT.EQ.1) OPEN (UNIT=IUNIT,FILE=NAME,STATUS='OLD',ERR=10) IF (ISTAT.EQ.2) OPEN (UNIT=IUNIT,FILE=NAME,STATUS='NEW',ERR=10) IF (ISTAT.EQ.3) OPEN (UNIT=IUNIT,FILE=NAME,STATUS='UNKNOWN', . ERR=10) GO TO (20,20,20),ISTAT * 10 CONTINUE IERROR = 1 GO TO 30 * 20 CONTINUE IERROR = 0 30 CONTINUE RETURN * LAST EXECUTABLE LINE OF SOPEN END FUNCTION SDIFF(X,Y) C C.L.LAWSON AND R.J.HANSON, JET PROPULSION LABORATORY, 1973 JUNE 7 C APPEARED IN 'SOLVING LEAST SQUARES PROBLEMS', PRENTICE-HALL, 1974 C THIS IS USED AS A TEST SUBPROGRAM FOR THE LEVEL TWO BLAS. C REVISED 860623 C REVISED YYMMDD C AUTH=R. J. HANSON, SANDIA NATIONAL LABS. SDIFF = X - Y RETURN * LAST EXECUTABLE LINE OF SDIFF END * FUNCTION SBEG(ANOISE) C THIS IS A TEST SUBPROGRAM FOR THE LEVEL TWO BLAS. C REVISED 860623 C REVISED YYMMDD C AUTH=R. J. HANSON, SANDIA NATIONAL LABS. SAVE C GENERATE NUMBERS FOR CONSTRUCTION OF TEST CASES. IF (ANOISE) 10,30,20 10 MI = 891 MJ = 457 I = 7 J = 7 AJ = 0. SBEG = 0. RETURN * 20 J = J*MJ J = J - 997* (J/997) AJ = J - 498 C THE SEQUENCE OF VALUES OF I IS BOUNDED BETWEEN 1 AND 999 C IF INITIAL I = 1,2,3,6,7, OR 9, THE PERIOD WILL BE 50 C IF INITIAL I = 4 OR 8 THE PERIOD WILL BE 25 C IF INITIAL I = 5 THE PERIOD WILL BE 10 30 I = I*MI I = I - 1000* (I/1000) AI = I - 500 SBEG = AI + AJ*ANOISE RETURN * LAST EXECUTABLE LINE OF SBEG END * LOGICAL FUNCTION LSE(RI,RJ,M,N,LDI) C TEST IF TWO REAL ARRAYS ARE IDENTICAL. C THE ARRAYS ARE M BY N. C THIS IS A TEST SUBPROGRAM FOR THE LEVEL TWO BLAS. C REVISED 860623 C REVISED YYMMDD C AUTH=R. J. HANSON, SANDIA NATIONAL LABS. REAL RI(LDI,*),RJ(LDI,*) DO 20 I = 1,M DO 10 J = 1,N IF (RI(I,J).NE.RJ(I,J)) THEN LSE = .FALSE. GO TO 30 * END IF * 10 CONTINUE 20 CONTINUE LSE = .TRUE. 30 CONTINUE RETURN * LAST EXECUTABLE LINE OF LSE END * LOGICAL FUNCTION LDE(DI,DJ,M,N,LDI) C TEST IF TWO DOUBLE PRECISION ARRAYS ARE IDENTICAL. C THE ARRAYS ARE M BY N. C THIS IS A TEST SUBPROGRAM FOR THE LEVEL TWO BLAS. C REVISED 860623 C REVISED YYMMDD C AUTH=R. J. HANSON, SANDIA NATIONAL LABS. DOUBLE PRECISION DI(LDI,*),DJ(LDI,*) DO 20 I = 1,M DO 10 J = 1,N IF (DI(I,J).NE.DJ(I,J)) THEN LDE = .FALSE. GO TO 30 * END IF * 10 CONTINUE 20 CONTINUE LDE = .TRUE. 30 CONTINUE RETURN * LAST EXECUTABLE LINE OF LDE END * LOGICAL FUNCTION LCE(CI,CJ,M,N,LDI) C TEST IF TWO COMPLEX ARRAYS ARE IDENTICAL. C THE ARRAYS ARE M BY N. C THIS IS A TEST SUBPROGRAM FOR THE LEVEL TWO BLAS. C REVISED 860623 C REVISED YYMMDD C AUTH=R. J. HANSON, SANDIA NATIONAL LABS. COMPLEX CI(LDI,*),CJ(LDI,*) DO 20 I = 1,M DO 10 J = 1,N IF (REAL(CI(I,J)).NE.REAL(CJ(I,J)) .OR. AIMAG(CI(I,J)).NE. . AIMAG(CJ(I,J))) THEN LCE = .FALSE. GO TO 30 * END IF * 10 CONTINUE 20 CONTINUE LCE = .TRUE. 30 CONTINUE RETURN * LAST EXECUTABLE LINE OF LCE END C C*********************************************************************** C C File of the REAL Level 2 BLAS routines: C C SGEMV, SGBMV, SSYMV, SSBMV, SSPMV, STRMV, STBMV, STPMV, C SGER , SSYR , SSPR , C SSYR2, SSPR2, C STRSV, STBSV, STPSV. C C See: C C Dongarra J. J., Du Croz J. J., Hammarling S. and Hanson R. J.. C A proposal for an extended set of Fortran Basic Linear Algebra C Subprograms. Technical Memorandum No.41 (revision 1), C Mathematics and Computer Science Division, Argone National C Laboratory, 9700 South Cass Avenue, Argonne, Illinois 60439, C USA, or NAG Technical Report TR4/85, Numerical Algorithms Group C Inc., 1101 31st Street, Suite 100, Downers Grove, Illinois C 60606-1263, USA. C C*********************************************************************** C SUBROUTINE SGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY) CHARACTER *1 TRANS INTEGER M,N,LDA,INCX,INCY REAL ALPHA,A(LDA,*),X(*),BETA,Y(*) * * Purpose * ======= * * SGEMV performs one of the matrix-vector operations * * y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, * * where alpha and beta are scalars, x and y are vectors and A is an * m by n matrix. * * Parameters * ========== * * TRANS - CHARACTER*1. * On entry, TRANS specifies the operation to be performed as * follows: * * TRANS = 'N' y := alpha*A*x + beta*y. * * TRANS = 'T' y := alpha*A'*x + beta*y. * * TRANS = 'C' y := alpha*A'*x + beta*y *. * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of the matrix A. * M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - REAL . * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - REAL array of DIMENSION ( LDA, n ). * Before entry, the leading m by n part of the array A must * contain the matrix of coefficients. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the leading dimension of A as * declared in the calling (sub) program. LDA must be at least * max(m,1). * Unchanged on exit. * * X - REAL array of DIMENSION at least * ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' * and at least * ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. * Before entry, the incremented array X must contain the * vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. * Unchanged on exit. * * BETA - REAL * On entry, BETA specifies the scalar beta. When BETA is * supplied as zero then Y need not be set on input. * Unchanged on exit. * * Y - REAL array of DIMENSION at least * ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' * and at least * ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. * Before entry with BETA non-zero, the incremented array Y * must contain the vector y. On exit, Y is overwritten by the * updated vector y. * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. * Unchanged on exit. * * * Note that TRANS, M, N and LDA must be such that the value of the * LOGICAL variable OK in the following statement is true. * * * * * Level 2 Blas routine. * * -- Written on 30-August-1985. * Sven Hammarling, Nag Central Office. C REVISED 860623 C REVISED YYMMDD C BY R. J. HANSON, SANDIA NATIONAL LABS. * INTEGER I,IX,IY,J,JX,JY INTEGER KX,KY,LENX,LENY REAL ONE,ZERO PARAMETER (ONE=1.0E+0,ZERO=0.0E+0) REAL TEMP LOGICAL OK,LSAME OK = (LSAME(TRANS,'N') .OR. LSAME(TRANS,'T') .OR. . LSAME(TRANS,'C')) .AND. ((M.GT.0) .AND. (N.GT.0) .AND. . (LDA.GE.M)) * * Quick return if possible. * IF (((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE)) .OR. .NOT. OK) RETURN * * Set LENX and LENY, the lengths of the vectors x and y. * IF (LSAME(TRANS,'N')) THEN LENX = N LENY = M * ELSE LENX = M LENY = N END IF * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * * First form y := beta*y and set up the start points in X and Y if * the increments are not both unity. * IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN IF (BETA.NE.ONE) THEN IF (BETA.EQ.ZERO) THEN DO 10,I = 1,LENY Y(I) = ZERO 10 CONTINUE * ELSE DO 20,I = 1,LENY Y(I) = BETA*Y(I) 20 CONTINUE END IF * END IF * ELSE IF (INCX.GT.0) THEN KX = 1 * ELSE KX = 1 - (LENX-1)*INCX END IF * IF (INCY.GT.0) THEN KY = 1 * ELSE KY = 1 - (LENY-1)*INCY END IF * IF (BETA.NE.ONE) THEN IY = KY IF (BETA.EQ.ZERO) THEN DO 30,I = 1,LENY Y(IY) = ZERO IY = IY + INCY 30 CONTINUE * ELSE DO 40,I = 1,LENY Y(IY) = BETA*Y(IY) IY = IY + INCY 40 CONTINUE END IF * END IF * END IF * IF (ALPHA.EQ.ZERO) RETURN IF (LSAME(TRANS,'N')) THEN * * Form y := alpha*A*x + y. * IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN DO 60,J = 1,N IF (X(J).NE.ZERO) THEN TEMP = ALPHA*X(J) DO 50,I = 1,M Y(I) = Y(I) + TEMP*A(I,J) 50 CONTINUE END IF * 60 CONTINUE * ELSE JX = KX DO 80,J = 1,N IF (X(JX).NE.ZERO) THEN TEMP = ALPHA*X(JX) IY = KY DO 70,I = 1,M Y(IY) = Y(IY) + TEMP*A(I,J) IY = IY + INCY 70 CONTINUE END IF * JX = JX + INCX 80 CONTINUE END IF * ELSE * * Form y := alpha*A'*x + y. * IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN DO 100,J = 1,N TEMP = ZERO DO 90,I = 1,M TEMP = TEMP + A(I,J)*X(I) 90 CONTINUE Y(J) = Y(J) + ALPHA*TEMP 100 CONTINUE * ELSE JY = KY DO 120,J = 1,N TEMP = ZERO IX = KX DO 110,I = 1,M TEMP = TEMP + A(I,J)*X(IX) IX = IX + INCX 110 CONTINUE Y(JY) = Y(JY) + ALPHA*TEMP JY = JY + INCY 120 CONTINUE END IF * END IF * RETURN * * End of SGEMV . * END SUBROUTINE SGBMV(TRANS,M,N,KL,KU,ALPHA,A,LDA,X,INCX,BETA,Y,INCY) CHARACTER *1 TRANS INTEGER M,N,KL,KU,LDA,INCX,INCY REAL ALPHA,A(LDA,*),X(*),BETA,Y(*) * * Purpose * ======= * * SGBMV performs one of the matrix-vector operations * * y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, * * where alpha and beta are scalars, x and y are vectors and A is an * m by n band matrix, with kl sub-diagonals and ku super-diagonals. * * Parameters * ========== * * TRANS - CHARACTER*1. * On entry, TRANS specifies the operation to be performed as * follows: * * TRANS = 'N' y := alpha*A*x + beta*y. * * TRANS = 'T' y := alpha*A'*x + beta*y. * * TRANS = 'C' y := alpha*A'*x + beta*y. * * Unchanged on exit. * * M - INTEGER. * On entry, M specifies the number of rows of the matrix A. * M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix A. * N must be at least zero. * Unchanged on exit. * * KL - INTEGER. * On entry, KL specifies the number of sub-diagonals of the * matrix A. KL must satisfy 0 .le. KL. * Unchanged on exit. * * KU - INTEGER. * On entry, KU specifies the number of super-diagonals of the * matrix A. KU must satisfy 0 .le. KU. * Unchanged on exit. * * Users may find that efficiency of their application is enhanced by * adjusting the values of m and n so that KL .ge. max(0,m-n) and * KU .ge. max(0,n-m) or KL and KU so that KL .lt. m and KU .lt. n. * * ALPHA - REAL . * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - REAL array of DIMENSION ( LDA, n ). * Before entry, the leading ( kl + ku + 1 ) by n part of the * array A must contain the matrix of coefficients, supplied * column by column, with the leading diagonal of the matrix in * row ( ku + 1 ) of the array, the first super-diagonal * starting at position 2 in row ku, the first sub-diagonal * starting at position 1 in row ( ku + 2 ), and so on. * This placement of the data can be realized with the * following loops: * DO 20 J =1,N * K=KU+1-J * DO 10 I =MAX(1,J-KU),MIN(M,J+KL) * A(K+I,J)=matrix entry of row I, column J. * 10 CONTINUE * 20 CONTINUE * Elements in the array A that do not correspond to elements * in the band matrix (such as the top left ku by ku triangle) * are not referenced. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the leading dimension of A as * declared in the calling (sub) program. LDA must be at least * ( kl + ku + 1 ). * Unchanged on exit. * * X - REAL array of DIMENSION at least * ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' * and at least * ( 1 + ( m - 1 )*abs( INCX ) ) otherwise. * Before entry, the incremented array X must contain the * vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. * Unchanged on exit. * * BETA - REAL . * On entry, BETA specifies the scalar beta. When BETA is * supplied as zero then Y need not be set on input. * Unchanged on exit. * * Y - REAL array of DIMENSION at least * ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' * and at least * ( 1 + ( n - 1 )*abs( INCY ) ) otherwise. * Before entry, the incremented array Y must contain the * vector y. On exit, Y is overwritten by the updated vector y. * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. * Unchanged on exit. * * * * * Level 2 Blas routine. * * -- Written on 27-Sept-1985. * Sven Hammarling, Nag Central Office. C REVISED 860623 C REVISED YYMMDD C BY R. J. HANSON, SANDIA NATIONAL LABS. * INTRINSIC MAX,MIN INTEGER I,IX,IY,J,JX,JY INTEGER K,KUP1,KX,KY,LENX,LENY REAL ONE,ZERO PARAMETER (ONE=1.0E+0,ZERO=0.0E+0) REAL TEMP LOGICAL OK,LSAME OK = (LSAME(TRANS,'N') .OR. LSAME(TRANS,'T') .OR. . LSAME(TRANS,'C')) .AND. (M.GT.0) .AND. (N.GT.0) .AND. . (KL.GE.0) .AND. (KU.GE.0) .AND. . (LDA.GE. (KL+KU+1)) * * Quick return if possible. * IF ( .NOT. OK .OR. ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN * * Set LENX and LENY, the lengths of the vectors x and y. * IF (LSAME(TRANS,'N')) THEN LENX = N LENY = M * ELSE LENX = M LENY = N END IF * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through the band part of A. * * First form y := beta*y and set up the start points in X and Y * if the increments are not both unity. * IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN IF (BETA.NE.ONE) THEN IF (BETA.EQ.ZERO) THEN DO 10,I = 1,LENY Y(I) = ZERO 10 CONTINUE * ELSE DO 20,I = 1,LENY Y(I) = BETA*Y(I) 20 CONTINUE END IF * END IF * ELSE IF (INCX.GT.0) THEN KX = 1 * ELSE KX = 1 - (LENX-1)*INCX END IF * IF (INCY.GT.0) THEN KY = 1 * ELSE KY = 1 - (LENY-1)*INCY END IF * IF (BETA.NE.ONE) THEN IY = KY IF (BETA.EQ.ZERO) THEN DO 30,I = 1,LENY Y(IY) = ZERO IY = IY + INCY 30 CONTINUE * ELSE DO 40,I = 1,LENY Y(IY) = BETA*Y(IY) IY = IY + INCY 40 CONTINUE END IF * END IF * END IF * IF (ALPHA.EQ.ZERO) RETURN KUP1 = KU + 1 IF (LSAME(TRANS,'N')) THEN * * Form y := alpha*A*x + y. * IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN DO 60,J = 1,N IF (X(J).NE.ZERO) THEN TEMP = ALPHA*X(J) K = KUP1 - J DO 50,I = MAX(1,J-KU),MIN(M,J+KL) Y(I) = Y(I) + TEMP*A(K+I,J) 50 CONTINUE END IF * 60 CONTINUE * ELSE JX = KX DO 80,J = 1,N IF (X(JX).NE.ZERO) THEN TEMP = ALPHA*X(JX) IY = KY K = KUP1 - J DO 70,I = MAX(1,J-KU),MIN(M,J+KL) Y(IY) = Y(IY) + TEMP*A(K+I,J) IY = IY + INCY 70 CONTINUE END IF * JX = JX + INCX IF (J.GT.KU) KY = KY + INCY 80 CONTINUE END IF * ELSE * * Form y := alpha*A'*x + y. * IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN DO 100,J = 1,N TEMP = ZERO K = KUP1 - J DO 90,I = MAX(1,J-KU),MIN(M,J+KL) TEMP = TEMP + A(K+I,J)*X(I) 90 CONTINUE Y(J) = Y(J) + ALPHA*TEMP 100 CONTINUE * ELSE JY = KY DO 120,J = 1,N TEMP = ZERO IX = KX K = KUP1 - J DO 110,I = MAX(1,J-KU),MIN(M,J+KL) TEMP = TEMP + A(K+I,J)*X(IX) IX = IX + INCX 110 CONTINUE Y(JY) = Y(JY) + ALPHA*TEMP JY = JY + INCY IF (J.GT.KU) KX = KX + INCX 120 CONTINUE END IF * END IF * RETURN * * End of SGBMV . * END SUBROUTINE SSYMV(UPLO,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY) CHARACTER *1 UPLO INTEGER N,LDA,INCX,INCY REAL ALPHA,A(LDA,*),X(*),BETA,Y(*) * * Purpose * ======= * * SSYMV performs the matrix-vector operation * * y := alpha*A*x + beta*y, * * where alpha and beta are scalars, x and y are n element vectors and * A is an n by n symmetric matrix. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the upper or lower * triangular part of the array A is to be referenced as * follows: * * UPLO = 'U' Only the upper triangular part of A * is to be referenced. * * UPLO = 'L' Only the lower triangular part of A * is to be referenced. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - REAL . * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - REAL array of DIMENSION ( LDA, n ). * Before entry with UPLO = 'U', the leading n by n * upper triangular part of the array A must contain the upper * triangular part of the symmetric matrix and the strictly * lower triangular part of A is not referenced. * Before entry with UPLO = 'L', the leading n by n * lower triangular part of the array A must contain the lower * triangular part of the symmetric matrix and the strictly * upper triangular part of A is not referenced. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least max(n,1). * Unchanged on exit. * * X - REAL array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. * Unchanged on exit. * * BETA - REAL . * On entry, BETA specifies the scalar beta. When BETA is * supplied as zero then Y need not be set on input. * Unchanged on exit. * * Y - REAL array of dimension at least * ( 1 + ( n - 1 )*abs( INCY ) ). * Before entry, the incremented array Y must contain the n * element vector y. On exit, Y is overwritten by the updated * vector y. * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. * Unchanged on exit. * * * Level 2 Blas routine. * * -- Written on 27-Sept-1985. * Sven Hammarling, Nag Central Office. C REVISED 860623 C REVISED YYMMDD C BY R. J. HANSON, SANDIA NATIONAL LABS. * INTEGER I,IX,IY,J,JX,JY INTEGER KX,KY REAL ONE,ZERO PARAMETER (ONE=1.0E+0,ZERO=0.0E+0) REAL TEMP1,TEMP2 LOGICAL OK,LSAME OK = (LSAME(UPLO,'U') .OR. LSAME(UPLO,'L')) .AND. (N.GT.0) .AND. . (LDA.GE.N) * * Quick return if possible. * IF ( .NOT. OK .OR. ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through the triangular part * of A. * * First form y := beta*y and set up the start points in X and Y if * the increments are not both unity. * IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN IF (BETA.NE.ONE) THEN IF (BETA.EQ.ZERO) THEN DO 10,I = 1,N Y(I) = ZERO 10 CONTINUE * ELSE DO 20,I = 1,N Y(I) = BETA*Y(I) 20 CONTINUE END IF * END IF * ELSE IF (INCX.GT.0) THEN KX = 1 * ELSE KX = 1 - (N-1)*INCX END IF * IF (INCY.GT.0) THEN KY = 1 * ELSE KY = 1 - (N-1)*INCY END IF * IF (BETA.NE.ONE) THEN IY = KY IF (BETA.EQ.ZERO) THEN DO 30,I = 1,N Y(IY) = ZERO IY = IY + INCY 30 CONTINUE * ELSE DO 40,I = 1,N Y(IY) = BETA*Y(IY) IY = IY + INCY 40 CONTINUE END IF * END IF * END IF * IF (ALPHA.EQ.ZERO) RETURN IF (LSAME(UPLO,'U')) THEN * * Form y when A is stored in upper triangle. * IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN DO 60,J = 1,N TEMP1 = ALPHA*X(J) TEMP2 = ZERO DO 50,I = 1,J - 1 Y(I) = Y(I) + TEMP1*A(I,J) TEMP2 = TEMP2 + A(I,J)*X(I) 50 CONTINUE Y(J) = Y(J) + TEMP1*A(J,J) + ALPHA*TEMP2 60 CONTINUE * ELSE IX = KX - INCX DO 80,J = 1,N TEMP1 = ALPHA*X(IX+INCX) TEMP2 = ZERO IX = KX IY = KY DO 70,I = 1,J - 1 Y(IY) = Y(IY) + TEMP1*A(I,J) TEMP2 = TEMP2 + A(I,J)*X(IX) IX = IX + INCX IY = IY + INCY 70 CONTINUE Y(IY) = Y(IY) + TEMP1*A(J,J) + ALPHA*TEMP2 80 CONTINUE END IF * ELSE * * Form y when A is stored in lower triangle. * IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN DO 100,J = 1,N TEMP1 = ALPHA*X(J) TEMP2 = ZERO Y(J) = Y(J) + TEMP1*A(J,J) DO 90,I = J + 1,N Y(I) = Y(I) + TEMP1*A(I,J) TEMP2 = TEMP2 + A(I,J)*X(I) 90 CONTINUE Y(J) = Y(J) + ALPHA*TEMP2 100 CONTINUE * ELSE JX = KX JY = KY DO 120,J = 1,N TEMP1 = ALPHA*X(JX) TEMP2 = ZERO Y(JY) = Y(JY) + TEMP1*A(J,J) IX = JX IY = JY DO 110,I = J + 1,N IX = IX + INCX IY = IY + INCY Y(IY) = Y(IY) + TEMP1*A(I,J) TEMP2 = TEMP2 + A(I,J)*X(IX) 110 CONTINUE Y(JY) = Y(JY) + ALPHA*TEMP2 JX = JX + INCX JY = JY + INCY 120 CONTINUE END IF * END IF * RETURN * * End of SSYMV . * END SUBROUTINE SSBMV(UPLO,N,K,ALPHA,A,LDA,X,INCX,BETA,Y,INCY) CHARACTER *1 UPLO INTEGER N,K,LDA,INCX,INCY REAL ALPHA,A(LDA,*),X(*),BETA,Y(*) * * Purpose * ======= * * SSBMV performs the matrix-vector operation * * y := alpha*A*x + beta*y, * * where alpha and beta are scalars, x and y are n element vectors and * A is an n by n symmetric band matrix, with k super-diagonals. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the upper or lower * triangular part of the band matrix A is being supplied as * follows: * * UPLO = 'U' The upper triangular part of A is * being supplied. * * UPLO = 'L' The lower triangular part of A is * being supplied. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * K - INTEGER. * On entry, K specifies the number of super-diagonals of the * matrix A. K must satisfy 0 .le. K .lt. n. * Unchanged on exit. * * ALPHA - REAL . * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * A - REAL array of DIMENSION ( LDA, n ). * Before entry with UPLO = 'U', the leading ( k + 1 ) * by n part of the array A must contain the upper triangular * band part of the symmetric matrix, supplied column by * column, with the leading diagonal of the matrix in row * ( k + 1 ) of the array, the first super-diagonal starting at * position 2 in row k, and so on. The top left k by k triangle * of the array A is not referenced. * Before entry with UPLO = 'L', the leading ( k + 1 ) * by n part of the array A must contain the lower triangular * band part of the symmetric matrix, supplied column by * column, with the leading diagonal of the matrix in row 1 of * the array, the first sub-diagonal starting at position 1 in * row 2, and so on. The bottom right k by k triangle of the * array A is not referenced. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the leading dimension of A as * declared in the calling (sub) program. LDA must be at least * ( k + 1 ). * Unchanged on exit. * * X - REAL array of DIMENSION at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the * vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. * Unchanged on exit. * * BETA - REAL . * On entry, BETA specifies the scalar beta. * Unchanged on exit. * * Y - REAL array of DIMENSION at least * ( 1 + ( n - 1 )*abs( INCY ) ). * Before entry, the incremented array Y must contain the * vector y. On exit, Y is overwritten by the updated vector y. * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. * Unchanged on exit. * * * * Level 2 Blas routine. * * -- Written on 30-September-1985. * Sven Hammarling, Nag Central Office. C REVISED 860623 C REVISED YYMMDD C BY R. J. HANSON, SANDIA NATIONAL LABS. * INTRINSIC MAX,MIN INTEGER I,IX,IY,J,JX,JY INTEGER KPLUS1,KX,KY,L REAL ONE,ZERO PARAMETER (ONE=1.0E+0,ZERO=0.0E+0) REAL TEMP1,TEMP2 LOGICAL OK,LSAME OK = (LSAME(UPLO,'U') .OR. LSAME(UPLO,'L')) .AND. (N.GT.0) .AND. . (K.GE.0) .AND. (K.LT.N) .AND. (LDA.GE. (K+1)) * * Quick return if possible. * IF ( .NOT. OK .OR. ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN * * Start the operations. In this version the elements of the array A * are accessed sequentially with one pass through A. * * First form y := beta*y and set up the start points in X and Y if * the increments are not both unity. * IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN IF (BETA.NE.ONE) THEN IF (BETA.EQ.ZERO) THEN DO 10,I = 1,N Y(I) = ZERO 10 CONTINUE * ELSE DO 20,I = 1,N Y(I) = BETA*Y(I) 20 CONTINUE END IF * END IF * ELSE IF (INCX.GT.0) THEN KX = 1 * ELSE KX = 1 - (N-1)*INCX END IF * IF (INCY.GT.0) THEN KY = 1 * ELSE KY = 1 - (N-1)*INCY END IF * IF (BETA.NE.ONE) THEN IY = KY IF (BETA.EQ.ZERO) THEN DO 30,I = 1,N Y(IY) = ZERO IY = IY + INCY 30 CONTINUE * ELSE DO 40,I = 1,N Y(IY) = BETA*Y(IY) IY = IY + INCY 40 CONTINUE END IF * END IF * END IF * IF (ALPHA.EQ.ZERO) RETURN IF (LSAME(UPLO,'U')) THEN * * Form y when upper triangle of A is stored. * KPLUS1 = K + 1 IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN DO 60,J = 1,N TEMP1 = ALPHA*X(J) TEMP2 = ZERO I = MAX(1,J-K) DO 50,L = KPLUS1 + I - J,K Y(I) = Y(I) + TEMP1*A(L,J) TEMP2 = TEMP2 + A(L,J)*X(I) I = I + 1 50 CONTINUE Y(J) = Y(J) + TEMP1*A(KPLUS1,J) + ALPHA*TEMP2 60 CONTINUE * ELSE IX = KX - INCX DO 80,J = 1,N TEMP1 = ALPHA*X(IX+INCX) TEMP2 = ZERO IX = KX IY = KY DO 70,L = 1 + MAX(KPLUS1-J,0),K Y(IY) = Y(IY) + TEMP1*A(L,J) TEMP2 = TEMP2 + A(L,J)*X(IX) IX = IX + INCX IY = IY + INCY 70 CONTINUE Y(IY) = Y(IY) + TEMP1*A(KPLUS1,J) + ALPHA*TEMP2 IF (J.GT.K) THEN KX = KX + INCX KY = KY + INCY END IF * 80 CONTINUE END IF * ELSE * * Form y when lower triangle of A is stored. * IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN DO 100,J = 1,N TEMP1 = ALPHA*X(J) TEMP2 = ZERO Y(J) = Y(J) + TEMP1*A(1,J) I = J + 1 DO 90,L = 2,1 + MIN(K,N-J) Y(I) = Y(I) + TEMP1*A(L,J) TEMP2 = TEMP2 + A(L,J)*X(I) I = I + 1 90 CONTINUE Y(J) = Y(J) + ALPHA*TEMP2 100 CONTINUE * ELSE JX = KX JY = KY DO 120,J = 1,N TEMP1 = ALPHA*X(JX) TEMP2 = ZERO Y(JY) = Y(JY) + TEMP1*A(1,J) IX = JX IY = JY DO 110,L = 2,1 + MIN(K,N-J) IX = IX + INCX IY = IY + INCY Y(IY) = Y(IY) + TEMP1*A(L,J) TEMP2 = TEMP2 + A(L,J)*X(IX) 110 CONTINUE Y(JY) = Y(JY) + ALPHA*TEMP2 JX = JX + INCX JY = JY + INCY 120 CONTINUE END IF * END IF * RETURN * * End of SSBMV . * END SUBROUTINE SSPMV(UPLO,N,ALPHA,AP,X,INCX,BETA,Y,INCY) CHARACTER *1 UPLO INTEGER N,INCX,INCY REAL ALPHA,AP(*),X(*),BETA,Y(*) * * Purpose * ======= * * SSPMV performs the matrix-vector operation * * y := alpha*A*x + beta*y, * * where alpha and beta are scalars, x and y are n element vectors and * A is an n by n symmetric matrix. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the upper or lower * triangular part of the matrix A is supplied in the packed * array AP as follows: * * UPLO = 'U' The upper triangular part of A is * supplied in AP. * * UPLO = 'L' The lower triangular part of A is * supplied in AP. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - REAL . * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * AP - REAL array of DIMENSION at least * ( ( n*( n + 1 ) )/2 ). * Before entry with UPLO = 'U', the array AP must * contain the upper triangular part of the symmetric matrix * packed sequentially, column by column, so that AP( 1 ) * contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) * and a( 2, 2 ) respectively, and so on. * Before entry with UPLO = 'L', the array AP must * contain the lower triangular part of the symmetric matrix * packed sequentially, column by column, so that AP( 1 ) * contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) * and a( 3, 1 ) respectively, and so on. * Unchanged on exit. * * X - REAL array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. * Unchanged on exit. * * BETA - REAL . * On entry, BETA specifies the scalar beta. When BETA is * supplied as zero then Y need not be set on input. * Unchanged on exit. * * Y - REAL array of dimension at least * ( 1 + ( n - 1 )*abs( INCY ) ). * Before entry, the incremented array Y must contain the n * element vector y. On exit, Y is overwritten by the updated * vector y. * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. * Unchanged on exit. * * * * * * Level 2 Blas routine. * * -- Written on 27-Sept-1985. * Sven Hammarling, Nag Central Office. C REVISED 860623 C REVISED YYMMDD C BY R. J. HANSON, SANDIA NATIONAL LABS. * INTEGER I,IX,IY,J,JX,JY INTEGER K,KK,KX,KY REAL ONE,ZERO PARAMETER (ONE=1.0E+0,ZERO=0.0E+0) REAL TEMP1,TEMP2 LOGICAL OK,LSAME OK = (LSAME(UPLO,'U') .OR. LSAME(UPLO,'L')) .AND. (N.GT.0) * * Quick return if possible. * IF ( .NOT. OK .OR. ((ALPHA.EQ.ZERO).AND. (BETA.EQ.ONE))) RETURN * * Start the operations. In this version the elements of the array AP * are accessed sequentially with one pass through AP. * * First form y := beta*y and set up the start points in X and Y if * the increments are not both unity. * IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN IF (BETA.NE.ONE) THEN IF (BETA.EQ.ZERO) THEN DO 10,I = 1,N Y(I) = ZERO 10 CONTINUE * ELSE DO 20,I = 1,N Y(I) = BETA*Y(I) 20 CONTINUE END IF * END IF * ELSE IF (INCX.GT.0) THEN KX = 1 * ELSE KX = 1 - (N-1)*INCX END IF * IF (INCY.GT.0) THEN KY = 1 * ELSE KY = 1 - (N-1)*INCY END IF * IF (BETA.NE.ONE) THEN IY = KY IF (BETA.EQ.ZERO) THEN DO 30,I = 1,N Y(IY) = ZERO IY = IY + INCY 30 CONTINUE * ELSE DO 40,I = 1,N Y(IY) = BETA*Y(IY) IY = IY + INCY 40 CONTINUE END IF * END IF * END IF * IF (ALPHA.EQ.ZERO) RETURN K = 1 IF (LSAME(UPLO,'U')) THEN * * Form y when AP contains the upper triangle. * IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN DO 60,J = 1,N TEMP1 = ALPHA*X(J) TEMP2 = ZERO DO 50,I = 1,J - 1 Y(I) = Y(I) + TEMP1*AP(K) TEMP2 = TEMP2 + AP(K)*X(I) K = K + 1 50 CONTINUE Y(J) = Y(J) + TEMP1*AP(K) + ALPHA*TEMP2 K = K + 1 60 CONTINUE * ELSE IX = KX - INCX DO 80,J = 1,N TEMP1 = ALPHA*X(IX+INCX) TEMP2 = ZERO IX = KX IY = KY KK = K DO 70,K = KK,KK + J - 2 Y(IY) = Y(IY) + TEMP1*AP(K) TEMP2 = TEMP2 + AP(K)*X(IX) IX = IX + INCX IY = IY + INCY 70 CONTINUE Y(IY) = Y(IY) + TEMP1*AP(K) + ALPHA*TEMP2 K = K + 1 80 CONTINUE END IF * ELSE * * Form y when AP contains the upper triangle. * IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN DO 100,J = 1,N TEMP1 = ALPHA*X(J) TEMP2 = ZERO Y(J) = Y(J) + TEMP1*AP(K) K = K + 1 DO 90,I = J + 1,N Y(I) = Y(I) + TEMP1*AP(K) TEMP2 = TEMP2 + AP(K)*X(I) K = K + 1 90 CONTINUE Y(J) = Y(J) + ALPHA*TEMP2 100 CONTINUE * ELSE JX = KX JY = KY DO 120,J = 1,N TEMP1 = ALPHA*X(JX) TEMP2 = ZERO Y(JY) = Y(JY) + TEMP1*AP(K) IX = JX IY = JY KK = K + 1 DO 110,K = KK,KK + N - (J+1) IX = IX + INCX IY = IY + INCY Y(IY) = Y(IY) + TEMP1*AP(K) TEMP2 = TEMP2 + AP(K)*X(IX) 110 CONTINUE Y(JY) = Y(JY) + ALPHA*TEMP2 JX = JX + INCX JY = JY + INCY 120 CONTINUE END IF * END IF * RETURN * * End of SSPMV . * END SUBROUTINE STRMV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX) CHARACTER *1 UPLO,TRANS,DIAG INTEGER N,LDA,INCX REAL A(LDA,*),X(*) * * Purpose * ======= * * STRMV performs one of the matrix-vector operations * * x := A*x, or x := A'*x, * * where x is n element vector and A is an n by n unit, or non-unit, * upper or lower triangular matrix. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the matrix is an upper or * lower triangular matrix as follows: * * UPLO = 'U' A is an upper triangular matrix. * * UPLO = 'L' A is a lower triangular matrix. * * Unchanged on exit. * * TRANS - CHARACTER*1. * On entry, TRANS specifies the operation to be performed as * follows: * * TRANS = 'N' x := A*x. * * TRANS = 'T' x := A'*x. * * TRANS = 'C' x := A'*x. * * Unchanged on exit. * * DIAG - CHARACTER*1. * On entry, DIAG specifies whether or not A is unit * triangular as follows: * * DIAG = 'U' A is assumed to be unit triangular. * * DIAG = 'N' A is not assumed to be unit * triangular. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * A - REAL array of DIMENSION ( LDA, n ). * Before entry with UPLO = 'U', the leading n by n * upper triangular part of the array A must contain the upper * triangular matrix and the strictly lower triangular part of * A is not referenced. * Before entry with UPLO = 'L', the leading n by n * lower triangular part of the array A must contain the lower * triangular matrix and the strictly upper triangular part of * A is not referenced. * Note that when DIAG = 'U', the diagonal elements of * A are not referenced either, but are assumed to be unity. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least max(n,1). * Unchanged on exit. * * X - REAL array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element vector x. On exit, X is overwritten with the * tranformed vector x. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. * Unchanged on exit. * * * * * Level 2 Blas routine. * * -- Written on 30-September-1985. * Sven Hammarling, Nag Central Office. C REVISED 860623 C REVISED YYMMDD C BY R. J. HANSON, SANDIA NATIONAL LABS. * LOGICAL NOUNIT INTEGER I,IX,J,JX,KX REAL ZERO PARAMETER (ZERO=0.0E+0) LOGICAL OK,LSAME OK = (LSAME(UPLO,'U') .OR. LSAME(UPLO,'L')) .AND. . (LSAME(TRANS,'N') .OR. LSAME(TRANS,'T') .OR. . LSAME(TRANS,'C')) .AND. (LSAME(DIAG,'U') .OR. . LSAME(DIAG,'N')) .AND. (N.GT.0) .AND. (LDA.GE.N) * * * Quick return if possible. * IF ( .NOT. OK) RETURN NOUNIT = LSAME(DIAG,'N') * * Set up the start point in X if the increment is not unity. This * will be ( N - 1 )*INCX too small for descending loops. * IF (INCX.LE.0) THEN KX = 1 - (N-1)*INCX * ELSE IF (INCX.NE.1) THEN KX = 1 END IF * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * IF (LSAME(TRANS,'N')) THEN * * Form x := A*x. * IF (LSAME(UPLO,'U')) THEN IF (INCX.EQ.1) THEN DO 20,J = 1,N IF (X(J).NE.ZERO) THEN DO 10,I = 1,J - 1 X(I) = X(I) + X(J)*A(I,J) 10 CONTINUE IF (NOUNIT) X(J) = X(J)*A(J,J) END IF * 20 CONTINUE * ELSE JX = KX DO 40,J = 1,N IF (X(JX).NE.ZERO) THEN IX = KX DO 30,I = 1,J - 1 X(IX) = X(IX) + X(JX)*A(I,J) IX = IX + INCX 30 CONTINUE IF (NOUNIT) X(JX) = X(JX)*A(J,J) END IF * JX = JX + INCX 40 CONTINUE END IF * ELSE IF (INCX.EQ.1) THEN DO 60,J = N,1,-1 IF (X(J).NE.ZERO) THEN DO 50,I = N,J + 1,-1 X(I) = X(I) + X(J)*A(I,J) 50 CONTINUE IF (NOUNIT) X(J) = X(J)*A(J,J) END IF * 60 CONTINUE * ELSE KX = KX + (N-1)*INCX JX = KX DO 80,J = N,1,-1 IF (X(JX).NE.ZERO) THEN IX = KX DO 70,I = N,J + 1,-1 X(IX) = X(IX) + X(JX)*A(I,J) IX = IX - INCX 70 CONTINUE IF (NOUNIT) X(JX) = X(JX)*A(J,J) END IF * JX = JX - INCX 80 CONTINUE END IF * END IF * ELSE * * Form x := A'*x. * IF (LSAME(UPLO,'U')) THEN IF (INCX.EQ.1) THEN DO 100,J = N,1,-1 IF (NOUNIT) X(J) = X(J)*A(J,J) DO 90,I = J - 1,1,-1 X(J) = X(J) + A(I,J)*X(I) 90 CONTINUE 100 CONTINUE * ELSE JX = KX + (N-1)*INCX DO 120,J = N,1,-1 IX = JX IF (NOUNIT) X(JX) = X(JX)*A(J,J) DO 110,I = J - 1,1,-1 IX = IX - INCX X(JX) = X(JX) + A(I,J)*X(IX) 110 CONTINUE JX = JX - INCX 120 CONTINUE END IF * ELSE IF (INCX.EQ.1) THEN DO 140,J = 1,N IF (NOUNIT) X(J) = X(J)*A(J,J) DO 130,I = J + 1,N X(J) = X(J) + A(I,J)*X(I) 130 CONTINUE 140 CONTINUE * ELSE JX = KX DO 160,J = 1,N IX = JX IF (NOUNIT) X(JX) = X(JX)*A(J,J) DO 150,I = J + 1,N IX = IX + INCX X(JX) = X(JX) + A(I,J)*X(IX) 150 CONTINUE JX = JX + INCX 160 CONTINUE END IF * END IF * END IF * RETURN * * End of STRMV . * END SUBROUTINE STBMV(UPLO,TRANS,DIAG,N,K,A,LDA,X,INCX) CHARACTER *1 UPLO,TRANS,DIAG INTEGER N,K,LDA,INCX REAL A(LDA,*),X(*) * * Purpose * ======= * * STBMV performs one of the matrix-vector operations * * x := A*x, or x := A'*x, * * where x is n element vector and A is an n by n unit, or non-unit, * upper or lower triangular band matrix, with ( k + 1 ) diagonals. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the matrix is an upper or * lower triangular matrix as follows: * * UPLO = 'U' A is an upper triangular matrix. * * UPLO = 'L' A is a lower triangular matrix. * * Unchanged on exit. * * TRANS - CHARACTER*1. * On entry, TRANS specifies the operation to be performed as * follows: * * TRANS = 'N' x := A*x. * * TRANS = 'T' x := A'*x. * * TRANS = 'C' x := A'*x. * * Unchanged on exit. * * DIAG - CHARACTER*1. * On entry, DIAG specifies whether or not A is unit * triangular as follows: * * DIAG = 'U' A is assumed to be unit triangular. * * DIAG = 'N' A is not assumed to be unit * triangular. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * K - INTEGER. * On entry with UPLO = 'U', K specifies the number of * super-diagonals of the matrix A. * On entry with UPLO = 'L', K specifies the number of * sub-diagonals of the matrix A. * K must satisfy 0 .le. K. * Unchanged on exit. * * A - REAL array of DIMENSION ( LDA, n ). * Before entry with UPLO = 'U', the leading ( k + 1 ) * by n part of the array A must contain the upper triangular * band part of the matrix of coefficients, supplied column by * column, with the leading diagonal of the matrix in row * ( k + 1 ) of the array, the first super-diagonal starting at * position 2 in row k, and so on. The top left k by k triangle * of the array A is not referenced. * Before entry with UPLO = 'L', the leading ( k + 1 ) * by n part of the array A must contain the lower triangular * band part of the matrix of coefficients, supplied column by * column, with the leading diagonal of the matrix in row 1 of * the array, the first sub-diagonal starting at position 1 in * row 2, and so on. The bottom right k by k triangle of the * array A is not referenced. * Note that when DIAG = 'U' the elements of the array A * corresponding to the diagonal elements of the matrix are not * referenced, but are assumed to be unity. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the leading dimension of A as * declared in the calling (sub) program. LDA must be at least * ( k + 1 ). * Unchanged on exit. * * X - REAL array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element vector x. On exit, X is overwritten with the * tranformed vector x. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. * Unchanged on exit. * * * * * Level 2 Blas routine. * * -- Written on 5-November-1985. * Sven Hammarling, Nag Central Office. C REVISED 860623 C REVISED YYMMDD C BY R. J. HANSON, SANDIA NATIONAL LABS. * INTRINSIC MAX,MIN LOGICAL NOUNIT INTEGER I,IX,J,JX,KPLUS1,KX INTEGER L REAL ZERO PARAMETER (ZERO=0.0E+0) LOGICAL OK,LSAME OK = (LSAME(UPLO,'U') .OR. LSAME(UPLO,'L')) .AND. . (LSAME(TRANS,'N') .OR. LSAME(TRANS,'T') .OR. . LSAME(TRANS,'C')) .AND. (LSAME(DIAG,'U') .OR. . LSAME(DIAG,'N')) .AND. (N.GT.0) .AND. (K.GE.0) .AND. . (LDA.GE. (K+1)) * * * Quick return if possible. * IF ( .NOT. OK) RETURN NOUNIT = LSAME(DIAG,'N') * * Set up the start point in X if the increment is not unity. This * will be ( N - 1 )*INCX too small for descending loops. * IF (INCX.LE.0) THEN KX = 1 - (N-1)*INCX * ELSE IF (INCX.NE.1) THEN KX = 1 END IF * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * IF (LSAME(TRANS,'N')) THEN * * Form x := A*x. * IF (LSAME(UPLO,'U')) THEN KPLUS1 = K + 1 IF (INCX.EQ.1) THEN DO 20,J = 1,N IF (X(J).NE.ZERO) THEN I = MAX(1,J-K) DO 10,L = KPLUS1 + I - J,K X(I) = X(I) + X(J)*A(L,J) I = I + 1 10 CONTINUE IF (NOUNIT) X(J) = X(J)*A(KPLUS1,J) END IF * 20 CONTINUE * ELSE JX = KX DO 40,J = 1,N IF (X(JX).NE.ZERO) THEN IX = KX DO 30,L = 1 + MAX(KPLUS1-J,0),K X(IX) = X(IX) + X(JX)*A(L,J) IX = IX + INCX 30 CONTINUE IF (NOUNIT) X(JX) = X(JX)*A(KPLUS1,J) END IF * JX = JX + INCX IF (J.GT.K) KX = KX + INCX 40 CONTINUE END IF * ELSE IF (INCX.EQ.1) THEN DO 60,J = N,1,-1 IF (X(J).NE.ZERO) THEN I = MIN(N,J+K) DO 50,L = 1 + I - J,2,-1 X(I) = X(I) + X(J)*A(L,J) I = I - 1 50 CONTINUE IF (NOUNIT) X(J) = X(J)*A(1,J) END IF * 60 CONTINUE * ELSE KX = KX + (N-1)*INCX JX = KX DO 80,J = N,1,-1 IF (X(JX).NE.ZERO) THEN IX = KX DO 70,L = 1 + MIN(K,N-J),2,-1 X(IX) = X(IX) + X(JX)*A(L,J) IX = IX - INCX 70 CONTINUE IF (NOUNIT) X(JX) = X(JX)*A(1,J) END IF * JX = JX - INCX IF ((N-J).GE.K) KX = KX - INCX 80 CONTINUE END IF * END IF * ELSE * * Form x := A'*x. * IF (LSAME(UPLO,'U')) THEN KPLUS1 = K + 1 IF (INCX.EQ.1) THEN DO 100,J = N,1,-1 I = J IF (NOUNIT) X(J) = X(J)*A(KPLUS1,J) DO 90,L = K,1 + MAX(KPLUS1-J,0),-1 I = I - 1 X(J) = X(J) + A(L,J)*X(I) 90 CONTINUE 100 CONTINUE * ELSE KX = KX + (N-1)*INCX JX = KX DO 120,J = N,1,-1 KX = KX - INCX IX = KX IF (NOUNIT) X(JX) = X(JX)*A(KPLUS1,J) DO 110,L = K,1 + MAX(KPLUS1-J,0),-1 X(JX) = X(JX) + A(L,J)*X(IX) IX = IX - INCX 110 CONTINUE JX = JX - INCX 120 CONTINUE END IF * ELSE IF (INCX.EQ.1) THEN DO 140,J = 1,N I = J IF (NOUNIT) X(J) = X(J)*A(1,J) DO 130,L = 2,1 + MIN(K,N-J) I = I + 1 X(J) = X(J) + A(L,J)*X(I) 130 CONTINUE 140 CONTINUE * ELSE JX = KX DO 160,J = 1,N KX = KX + INCX IX = KX IF (NOUNIT) X(JX) = X(JX)*A(1,J) DO 150,L = 2,1 + MIN(K,N-J) X(JX) = X(JX) + A(L,J)*X(IX) IX = IX + INCX 150 CONTINUE JX = JX + INCX 160 CONTINUE END IF * END IF * END IF * RETURN * * End of STBMV . * END SUBROUTINE STPMV(UPLO,TRANS,DIAG,N,AP,X,INCX) CHARACTER *1 UPLO,TRANS,DIAG INTEGER N,INCX REAL AP(*),X(*) * * Purpose * ======= * * STPMV performs one of the matrix-vector operations * * x := A*x, or x := A'*x, * * where x is n element vector and A is an n by n unit, or non-unit, * upper or lower triangular matrix. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the matrix is an upper or * lower triangular matrix as follows: * * UPLO = 'U' A is an upper triangular matrix. * * UPLO = 'L' A is a lower triangular matrix. * * Unchanged on exit. * * TRANS - CHARACTER*1. * On entry, TRANS specifies the operation to be performed as * follows: * * TRANS = 'N' x := A*x. * * TRANS = 'T' x := A'*x. * * TRANS = 'C' x := A'*x. * * Unchanged on exit. * * DIAG - CHARACTER*1. * On entry, DIAG specifies whether or not A is unit * triangular as follows: * * DIAG = 'U' A is assumed to be unit triangular. * * DIAG = 'N' A is not assumed to be unit * triangular. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * AP - REAL array of DIMENSION at least * ( ( n*( n + 1 ) )/2 ). * Before entry with UPLO = 'U', the array AP must * contain the upper triangular matrix packed sequentially, * column by column, so that AP( 1 ) contains a( 1, 1 ), * AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 ) * respectively, and so on. * Before entry with UPLO = 'L', the array AP must * contain the lower triangular matrix packed sequentially, * column by column, so that AP( 1 ) contains a( 1, 1 ), * AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 ) * respectively, and so on. * Note that when DIAG = 'U', the diagonal elements of * A are not referenced, but are assumed to be unity. * Unchanged on exit. * * X - REAL array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element vector x. On exit, X is overwritten with the * tranformed vector x. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. * Unchanged on exit. * * * Note that UPLO, TRANS, DIAG and N must be such that the value of the * LOGICAL variable OK in the following statement is true. * * * * Level 2 Blas routine. * * -- Written on 2-October-1985. * Sven Hammarling, Nag Central Office. C REVISED 860623 C REVISED YYMMDD C BY R. J. HANSON, SANDIA NATIONAL LABS. * LOGICAL NOUNIT INTEGER I,IX,J,JX,K,KK INTEGER KX REAL ZERO PARAMETER (ZERO=0.0E+0) LOGICAL OK,LSAME OK = (LSAME(UPLO,'U') .OR. LSAME(UPLO,'L')) .AND. . (LSAME(TRANS,'N') .OR. LSAME(TRANS,'T') .OR. . LSAME(TRANS,'C')) .AND. (LSAME(DIAG,'U') .OR. . LSAME(DIAG,'N')) .AND. (N.GT.0) * * * Quick return if possible. * IF ( .NOT. OK) RETURN NOUNIT = LSAME(DIAG,'N') * * Set up the start point in X if the increment is not unity. This * will be ( N - 1 )*INCX too small for descending loops. * IF (INCX.LE.0) THEN KX = 1 - (N-1)*INCX * ELSE IF (INCX.NE.1) THEN KX = 1 END IF * * Start the operations. In this version the elements of AP are * accessed sequentially with one pass through AP. * IF (LSAME(TRANS,'N')) THEN * * Form x:= A*x. * IF (LSAME(UPLO,'U')) THEN K = 1 IF (INCX.EQ.1) THEN DO 20,J = 1,N IF (X(J).NE.ZERO) THEN DO 10,I = 1,J - 1 X(I) = X(I) + X(J)*AP(K) K = K + 1 10 CONTINUE IF (NOUNIT) X(J) = X(J)*AP(K) K = K + 1 * ELSE K = K + J END IF * 20 CONTINUE * ELSE JX = KX DO 40,J = 1,N IF (X(JX).NE.ZERO) THEN IX = KX KK = K DO 30,K = KK,KK + J - 2 X(IX) = X(IX) + X(JX)*AP(K) IX = IX + INCX 30 CONTINUE IF (NOUNIT) X(JX) = X(JX)*AP(K) K = K + 1 * ELSE K = K + J END IF * JX = JX + INCX 40 CONTINUE END IF * ELSE K = (N* (N+1))/2 IF (INCX.EQ.1) THEN DO 60,J = N,1,-1 IF (X(J).NE.ZERO) THEN DO 50,I = N,J + 1,-1 X(I) = X(I) + X(J)*AP(K) K = K - 1 50 CONTINUE IF (NOUNIT) X(J) = X(J)*AP(K) K = K - 1 * ELSE K = K - (N-J+1) END IF * 60 CONTINUE * ELSE KX = KX + (N-1)*INCX JX = KX DO 80,J = N,1,-1 IF (X(JX).NE.ZERO) THEN IX = KX KK = K DO 70,K = KK,KK - (N- (J+1)),-1 X(IX) = X(IX) + X(JX)*AP(K) IX = IX - INCX 70 CONTINUE IF (NOUNIT) X(JX) = X(JX)*AP(K) K = K - 1 * ELSE K = K - (N-J+1) END IF * JX = JX - INCX 80 CONTINUE END IF * END IF * ELSE * * Form x := A'*x. * IF (LSAME(UPLO,'U')) THEN K = (N* (N+1))/2 IF (INCX.EQ.1) THEN DO 100,J = N,1,-1 IF (NOUNIT) X(J) = X(J)*AP(K) K = K - 1 DO 90,I = J - 1,1,-1 X(J) = X(J) + AP(K)*X(I) K = K - 1 90 CONTINUE 100 CONTINUE * ELSE JX = KX + (N-1)*INCX DO 120,J = N,1,-1 IX = JX IF (NOUNIT) X(JX) = X(JX)*AP(K) KK = K - 1 DO 110,K = KK,KK - J + 2,-1 IX = IX - INCX X(JX) = X(JX) + AP(K)*X(IX) 110 CONTINUE JX = JX - INCX 120 CONTINUE END IF * ELSE K = 1 IF (INCX.EQ.1) THEN DO 140,J = 1,N IF (NOUNIT) X(J) = X(J)*AP(K) K = K + 1 DO 130,I = J + 1,N X(J) = X(J) + AP(K)*X(I) K = K + 1 130 CONTINUE 140 CONTINUE * ELSE JX = KX DO 160,J = 1,N IX = JX IF (NOUNIT) X(JX) = X(JX)*AP(K) KK = K + 1 DO 150,K = KK,KK + N - (J+1) IX = IX + INCX X(JX) = X(JX) + AP(K)*X(IX) 150 CONTINUE JX = JX + INCX 160 CONTINUE END IF * END IF * END IF * RETURN * * End of STPMV . * END SUBROUTINE STRSV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX) CHARACTER *1 UPLO,TRANS,DIAG INTEGER N,LDA,INCX REAL A(LDA,*),X(*) * * Purpose * ======= * * STRSV solves one of the systems of equations * * A*x = b, or A'*x = b, * * where b and x are n element vectors and A is an n by n unit, or * non-unit, upper or lower triangular matrix. * * No test for singularity or near-singularity is included in this * routine. Such tests must be performed before calling this routine. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the matrix is an upper or * lower triangular matrix as follows: * * UPLO = 'U' A is an upper triangular matrix. * * UPLO = 'L' A is a lower triangular matrix. * * Unchanged on exit. * * TRANS - CHARACTER*1. * On entry, TRANS specifies the equations to be solved as * follows: * * TRANS = 'N' A*x = b. * * TRANS = 'T' A'*x = b. * * TRANS = 'C' A'*x = b. * * Unchanged on exit. * * DIAG - CHARACTER*1. * On entry, DIAG specifies whether or not A is unit * triangular as follows: * * DIAG = 'U' A is assumed to be unit triangular. * * DIAG = 'N' A is not assumed to be unit * triangular. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * A - REAL array of DIMENSION ( LDA, n ). * Before entry with UPLO = 'U', the leading n by n * upper triangular part of the array A must contain the upper * triangular matrix and the strictly lower triangular part of * A is not referenced. * Before entry with UPLO = 'L', the leading n by n * lower triangular part of the array A must contain the lower * triangular matrix and the strictly upper triangular part of * A is not referenced. * Note that when DIAG = 'U', the diagonal elements of * A are not referenced either, but are assumed to be unity. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least max(n,1). * Unchanged on exit. * * X - REAL array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element right-hand side vector b. On exit, X is overwritten * with the solution vector x. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. * Unchanged on exit. * * * * * Level 2 Blas routine. * * -- Written on 30-September-1985. * Sven Hammarling, Nag Central Office. C REVISED 860623 C REVISED YYMMDD C BY R. J. HANSON, SANDIA NATIONAL LABS. * LOGICAL NOUNIT INTEGER I,IX,J,JX,KX REAL ZERO PARAMETER (ZERO=0.0E+0) LOGICAL OK,LSAME OK = (LSAME(UPLO,'U') .OR. LSAME(UPLO,'L')) .AND. . (LSAME(TRANS,'N') .OR. LSAME(TRANS,'T') .OR. . LSAME(TRANS,'C')) .AND. (LSAME(DIAG,'U') .OR. . LSAME(DIAG,'N')) .AND. (N.GT.0) .AND. (LDA.GE.N) * * * Quick return if possible. * IF ( .NOT. OK) RETURN NOUNIT = LSAME(DIAG,'N') * * Set up the start point in X if the increment is not unity. This * will be ( N - 1 )*INCX too small for descending loops. * IF (INCX.LE.0) THEN KX = 1 - (N-1)*INCX * ELSE IF (INCX.NE.1) THEN KX = 1 END IF * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * IF (LSAME(TRANS,'N')) THEN * * Form x := inv( A )*x. * IF (LSAME(UPLO,'U')) THEN IF (INCX.EQ.1) THEN DO 20,J = N,1,-1 IF (X(J).NE.ZERO) THEN IF (NOUNIT) X(J) = X(J)/A(J,J) DO 10,I = J - 1,1,-1 X(I) = X(I) - X(J)*A(I,J) 10 CONTINUE END IF * 20 CONTINUE * ELSE JX = KX + (N-1)*INCX DO 40,J = N,1,-1 IF (X(JX).NE.ZERO) THEN IF (NOUNIT) X(JX) = X(JX)/A(J,J) IX = JX DO 30,I = J - 1,1,-1 IX = IX - INCX X(IX) = X(IX) - X(JX)*A(I,J) 30 CONTINUE END IF * JX = JX - INCX 40 CONTINUE END IF * ELSE IF (INCX.EQ.1) THEN DO 60,J = 1,N IF (X(J).NE.ZERO) THEN IF (NOUNIT) X(J) = X(J)/A(J,J) DO 50,I = J + 1,N X(I) = X(I) - X(J)*A(I,J) 50 CONTINUE END IF * 60 CONTINUE * ELSE JX = KX DO 80,J = 1,N IF (X(JX).NE.ZERO) THEN IF (NOUNIT) X(JX) = X(JX)/A(J,J) IX = JX DO 70,I = J + 1,N IX = IX + INCX X(IX) = X(IX) - X(JX)*A(I,J) 70 CONTINUE END IF * JX = JX + INCX 80 CONTINUE END IF * END IF * ELSE * * Form x := inv( A' )*x. * IF (LSAME(UPLO,'U')) THEN IF (INCX.EQ.1) THEN DO 100,J = 1,N DO 90,I = 1,J - 1 X(J) = X(J) - A(I,J)*X(I) 90 CONTINUE IF (NOUNIT) X(J) = X(J)/A(J,J) 100 CONTINUE * ELSE JX = KX DO 120,J = 1,N IX = KX DO 110,I = 1,J - 1 X(JX) = X(JX) - A(I,J)*X(IX) IX = IX + INCX 110 CONTINUE IF (NOUNIT) X(JX) = X(JX)/A(J,J) JX = JX + INCX 120 CONTINUE END IF * ELSE IF (INCX.EQ.1) THEN DO 140,J = N,1,-1 DO 130,I = N,J + 1,-1 X(J) = X(J) - A(I,J)*X(I) 130 CONTINUE IF (NOUNIT) X(J) = X(J)/A(J,J) 140 CONTINUE * ELSE KX = KX + (N-1)*INCX JX = KX DO 160,J = N,1,-1 IX = KX DO 150,I = N,J + 1,-1 X(JX) = X(JX) - A(I,J)*X(IX) IX = IX - INCX 150 CONTINUE IF (NOUNIT) X(JX) = X(JX)/A(J,J) JX = JX - INCX 160 CONTINUE END IF * END IF * END IF * RETURN * * End of STRSV . * END SUBROUTINE STBSV(UPLO,TRANS,DIAG,N,K,A,LDA,X,INCX) CHARACTER *1 UPLO,TRANS,DIAG INTEGER N,K,LDA,INCX REAL A(LDA,*),X(*) * * Purpose * ======= * * STBSV solves one of the systems of equations * * A*x = b, or A'*x = b, * * where b and x are n element vectors and A is an n by n unit, or * non-unit, upper or lower triangular band matrix, with ( k + 1 ) * diagonals. * * No test for singularity or near-singularity is included in this * routine. Such tests must be performed before calling this routine. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the matrix is an upper or * lower triangular matrix as follows: * * UPLO = 'U' A is an upper triangular matrix. * * UPLO = 'L' A is a lower triangular matrix. * * Unchanged on exit. * * TRANS - CHARACTER*1. * On entry, TRANS specifies the equations to be solved as * follows: * * TRANS = 'N' A*x = b. * * TRANS = 'T' A'*x = b. * * TRANS = 'C' A'*x = b. * * Unchanged on exit. * * DIAG - CHARACTER*1. * On entry, DIAG specifies whether or not A is unit * triangular as follows: * * DIAG = 'U' A is assumed to be unit triangular. * * DIAG = 'N' A is not assumed to be unit * triangular. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * K - INTEGER. * On entry with UPLO = 'U', K specifies the number of * super-diagonals of the matrix A. * On entry with UPLO = 'L', K specifies the number of * sub-diagonals of the matrix A. * K must satisfy 0 .le. K. * Unchanged on exit. * * A - REAL array of DIMENSION ( LDA, n ). * Before entry with UPLO = 'U', the leading ( k + 1 ) * by n part of the array A must contain the upper triangular * band part of the matrix of coefficients, supplied column by * column, with the leading diagonal of the matrix in row * ( k + 1 ) of the array, the first super-diagonal starting at * position 2 in row k, and so on. The top left k by k triangle * of the array A is not referenced. * Before entry with UPLO = 'L', the leading ( k + 1 ) * by n part of the array A must contain the lower triangular * band part of the matrix of coefficients, supplied column by * column, with the leading diagonal of the matrix in row 1 of * the array, the first sub-diagonal starting at position 1 in * row 2, and so on. The bottom right k by k triangle of the * array A is not referenced. * Note that when DIAG = 'U' the elements of the array A * corresponding to the diagonal elements of the matrix are not * referenced, but are assumed to be unity. * Unchanged on exit. * * LDA - INTEGER. * On entry, LDA specifies the leading dimension of A as * declared in the calling (sub) program. LDA must be at least * ( k + 1 ). * Unchanged on exit. * * X - REAL array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element right-hand side vector b. On exit, X is overwritten * with the solution vector x. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. * Unchanged on exit. * * * * * * Level 2 Blas routine. * * -- Written on 7-November-1985. * Sven Hammarling, Nag Central Office. C REVISED 860623 C REVISED YYMMDD C BY R. J. HANSON, SANDIA NATIONAL LABS. * INTRINSIC MAX,MIN LOGICAL NOUNIT INTEGER I,IX,J,JX,KPLUS1,KX INTEGER L REAL ZERO PARAMETER (ZERO=0.0E+0) LOGICAL OK,LSAME OK = (LSAME(UPLO,'U') .OR. LSAME(UPLO,'L')) .AND. . (LSAME(TRANS,'N') .OR. LSAME(TRANS,'T') .OR. . LSAME(TRANS,'C')) .AND. (LSAME(DIAG,'U') .OR. . LSAME(DIAG,'N')) .AND. (N.GT.0) .AND. (K.GE.0) .AND. . (LDA.GE. (K+1)) * * * Quick return if possible. * IF ( .NOT. OK) RETURN NOUNIT = LSAME(DIAG,'N') * * Set up the start point in X if the increment is not unity. This * will be ( N - 1 )*INCX too small for descending loops. * IF (INCX.LE.0) THEN KX = 1 - (N-1)*INCX * ELSE IF (INCX.NE.1) THEN KX = 1 END IF * * Start the operations. In this version the elements of A are * accessed by sequentially with one pass through A. * IF (LSAME(TRANS,'N')) THEN * * Form x := inv( A )*x. * IF (LSAME(UPLO,'U')) THEN KPLUS1 = K + 1 IF (INCX.EQ.1) THEN DO 20,J = N,1,-1 IF (X(J).NE.ZERO) THEN IF (NOUNIT) X(J) = X(J)/A(KPLUS1,J) I = J DO 10,L = K,1 + MAX(KPLUS1-J,0),-1 I = I - 1 X(I) = X(I) - X(J)*A(L,J) 10 CONTINUE END IF * 20 CONTINUE * ELSE KX = KX + (N-1)*INCX JX = KX DO 40,J = N,1,-1 KX = KX - INCX IX = KX IF (X(JX).NE.ZERO) THEN IF (NOUNIT) X(JX) = X(JX)/A(KPLUS1,J) DO 30 L = K,1 + MAX(KPLUS1-J,0),-1 X(IX) = X(IX) - X(JX)*A(L,J) IX = IX - INCX 30 CONTINUE END IF * JX = JX - INCX 40 CONTINUE END IF * ELSE IF (INCX.EQ.1) THEN DO 60,J = 1,N IF (X(J).NE.ZERO) THEN IF (NOUNIT) X(J) = X(J)/A(1,J) I = J DO 50,L = 2,1 + MIN(K,N-J) I = I + 1 X(I) = X(I) - X(J)*A(L,J) 50 CONTINUE END IF * 60 CONTINUE * ELSE JX = KX DO 80,J = 1,N KX = KX + INCX IF (X(JX).NE.ZERO) THEN IF (NOUNIT) X(JX) = X(JX)/A(1,J) IX = KX DO 70,L = 2,1 + MIN(K,N-J) X(IX) = X(IX) - X(JX)*A(L,J) IX = IX + INCX 70 CONTINUE END IF * JX = JX + INCX 80 CONTINUE END IF * END IF * ELSE * * Form x := inv( A')*x. * IF (LSAME(UPLO,'U')) THEN KPLUS1 = K + 1 IF (INCX.EQ.1) THEN DO 100,J = 1,N I = MAX(1,J-K) DO 90,L = KPLUS1 + I - J,K X(J) = X(J) - A(L,J)*X(I) I = I + 1 90 CONTINUE IF (NOUNIT) X(J) = X(J)/A(KPLUS1,J) 100 CONTINUE * ELSE JX = KX DO 120,J = 1,N IX = KX DO 110,L = 1 + MAX(KPLUS1-J,0),K X(JX) = X(JX) - A(L,J)*X(IX) IX = IX + INCX 110 CONTINUE IF (NOUNIT) X(JX) = X(JX)/A(KPLUS1,J) JX = JX + INCX IF (J.GT.K) KX = KX + INCX 120 CONTINUE END IF * ELSE IF (INCX.EQ.1) THEN DO 140,J = N,1,-1 I = MIN(N,J+K) DO 130,L = 1 + I - J,2,-1 X(J) = X(J) - A(L,J)*X(I) I = I - 1 130 CONTINUE IF (NOUNIT) X(J) = X(J)/A(1,J) 140 CONTINUE * ELSE KX = KX + (N-1)*INCX JX = KX DO 160,J = N,1,-1 IX = KX DO 150,L = 1 + MIN(K,N-J),2,-1 X(JX) = X(JX) - A(L,J)*X(IX) IX = IX - INCX 150 CONTINUE IF (NOUNIT) X(JX) = X(JX)/A(1,J) JX = JX - INCX IF ((N-J).GE.K) KX = KX - INCX 160 CONTINUE END IF * END IF * END IF * RETURN * * End of STBSV . * END SUBROUTINE STPSV(UPLO,TRANS,DIAG,N,AP,X,INCX) CHARACTER *1 UPLO,TRANS,DIAG INTEGER N,INCX REAL AP(*),X(*) * * Purpose * ======= * * STPSV solves one of the systems of equations * * A*x = b, or A'*x = b, * * where b and x are n element vectors and A is an n by n unit, or * non-unit, upper or lower triangular matrix. * * No test for singularity or near-singularity is included in this * routine. Such tests must be performed before calling this routine. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the matrix is an upper or * lower triangular matrix as follows: * * UPLO = 'U' A is an upper triangular matrix. * * UPLO = 'L' A is a lower triangular matrix. * * Unchanged on exit. * * TRANS - CHARACTER*1. * On entry, TRANS specifies the equations to be solved as * follows: * * TRANS = 'N' A*x = b. * * TRANS = 'T' A'*x = b. * * TRANS = 'C' A'*x = b. * * Unchanged on exit. * * DIAG - CHARACTER*1. * On entry, DIAG specifies whether or not A is unit * triangular as follows: * * DIAG = 'U' A is assumed to be unit triangular. * * DIAG = 'N' A is not assumed to be unit * triangular. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * AP - REAL array of DIMENSION at least * ( ( n*( n + 1 ) )/2 ). * Before entry with UPLO = 'U', the array AP must * contain the upper triangular matrix packed sequentially, * column by column, so that AP( 1 ) contains a( 1, 1 ), * AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 ) * respectively, and so on. * Before entry with UPLO = 'L', the array AP must * contain the lower triangular matrix packed sequentially, * column by column, so that AP( 1 ) contains a( 1, 1 ), * AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 ) * respectively, and so on. * Note that when DIAG = 'U', the diagonal elements of * A are not referenced, but are assumed to be unity. * Unchanged on exit. * * X - REAL array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element right-hand side vector b. On exit, X is overwritten * with the solution vector x. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. * Unchanged on exit. * * * * * * Level 2 Blas routine. * * -- Written on 11-November-1985. * Sven Hammarling, Nag Central Office. C REVISED 860623 C REVISED YYMMDD C BY R. J. HANSON, SANDIA NATIONAL LABS. * LOGICAL NOUNIT INTEGER I,IX,J,JX,K,KK INTEGER KX REAL ZERO PARAMETER (ZERO=0.0E+0) LOGICAL OK,LSAME OK = (LSAME(UPLO,'U') .OR. LSAME(UPLO,'L')) .AND. . (LSAME(TRANS,'N') .OR. LSAME(TRANS,'T') .OR. . LSAME(TRANS,'C')) .AND. (LSAME(DIAG,'U') .OR. . LSAME(DIAG,'N')) .AND. (N.GT.0) * * Quick return if possible. * IF ( .NOT. OK) RETURN NOUNIT = LSAME(DIAG,'N') * * Set up the start point in X if the increment is not unity. This * will be ( N - 1 )*INCX too small for descending loops. * IF (INCX.LE.0) THEN KX = 1 - (N-1)*INCX * ELSE IF (INCX.NE.1) THEN KX = 1 END IF * * Start the operations. In this version the elements of AP are * accessed sequentially with one pass through AP. * IF (LSAME(TRANS,'N')) THEN * * Form x := inv( A )*x. * IF (LSAME(UPLO,'U')) THEN K = (N* (N+1))/2 IF (INCX.EQ.1) THEN DO 20,J = N,1,-1 IF (X(J).NE.ZERO) THEN IF (NOUNIT) X(J) = X(J)/AP(K) K = K - 1 DO 10,I = J - 1,1,-1 X(I) = X(I) - X(J)*AP(K) K = K - 1 10 CONTINUE * ELSE K = K - J END IF * 20 CONTINUE * ELSE JX = KX + (N-1)*INCX DO 40,J = N,1,-1 IF (X(JX).NE.ZERO) THEN IF (NOUNIT) X(JX) = X(JX)/AP(K) IX = JX KK = K - 1 DO 30,K = KK,KK - J + 2,-1 IX = IX - INCX X(IX) = X(IX) - X(JX)*AP(K) 30 CONTINUE * ELSE K = K - J END IF * JX = JX - INCX 40 CONTINUE END IF * ELSE K = 1 IF (INCX.EQ.1) THEN DO 60,J = 1,N IF (X(J).NE.ZERO) THEN IF (NOUNIT) X(J) = X(J)/AP(K) K = K + 1 DO 50,I = J + 1,N X(I) = X(I) - X(J)*AP(K) K = K + 1 50 CONTINUE * ELSE K = K + N - J + 1 END IF * 60 CONTINUE * ELSE JX = KX DO 80,J = 1,N IF (X(JX).NE.ZERO) THEN IF (NOUNIT) X(JX) = X(JX)/AP(K) IX = JX KK = K + 1 DO 70,K = KK,KK + N - (J+1) IX = IX + INCX X(IX) = X(IX) - X(JX)*AP(K) 70 CONTINUE * ELSE K = K + N - J + 1 END IF * JX = JX + INCX 80 CONTINUE END IF * END IF * ELSE * * Form x := inv( A' )*x. * IF (LSAME(UPLO,'U')) THEN K = 1 IF (INCX.EQ.1) THEN DO 100,J = 1,N DO 90,I = 1,J - 1 X(J) = X(J) - AP(K)*X(I) K = K + 1 90 CONTINUE IF (NOUNIT) X(J) = X(J)/AP(K) K = K + 1 100 CONTINUE * ELSE JX = KX DO 120,J = 1,N IX = KX KK = K DO 110,K = KK,KK + J - 2 X(JX) = X(JX) - AP(K)*X(IX) IX = IX + INCX 110 CONTINUE IF (NOUNIT) X(JX) = X(JX)/AP(K) K = K + 1 JX = JX + INCX 120 CONTINUE END IF * ELSE K = (N* (N+1))/2 IF (INCX.EQ.1) THEN DO 140,J = N,1,-1 DO 130,I = N,J + 1,-1 X(J) = X(J) - AP(K)*X(I) K = K - 1 130 CONTINUE IF (NOUNIT) X(J) = X(J)/AP(K) K = K - 1 140 CONTINUE * ELSE KX = KX + (N-1)*INCX JX = KX DO 160,J = N,1,-1 IX = KX KK = K DO 150,K = KK,KK - (N- (J+1)),-1 X(JX) = X(JX) - AP(K)*X(IX) IX = IX - INCX 150 CONTINUE IF (NOUNIT) X(JX) = X(JX)/AP(K) K = K - 1 JX = JX - INCX 160 CONTINUE END IF * END IF * END IF * RETURN * * End of STPSV . * END SUBROUTINE SGER(M,N,ALPHA,X,INCX,Y,INCY,A,LDA) INTEGER M,N,INCX,INCY,LDA REAL ALPHA,X(*),Y(*),A(LDA,*) * * Purpose * ======= * * SGER performs the rank 1 operation * * A := alpha*x*y' + A, * * where alpha is a scalar, x is an m element vector, y is an n element * vector and A is an m by n matrix. * * Parameters * ========== * * M - INTEGER. * On entry, M specifies the number of rows of the matrix A. * M must be at least zero. * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the number of columns of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - REAL . * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * X - REAL array of dimension at least * ( 1 + ( m - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the m * element vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. * Unchanged on exit. * * Y - REAL array of dimension at least * ( 1 + ( n - 1 )*abs( INCY ) ). * Before entry, the incremented array Y must contain the n * element vector y. * Unchanged on exit. * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. * Unchanged on exit. * * A - REAL array of DIMENSION ( LDA, n ). * Before entry, the leading m by n part of the array A must * contain the matrix of coefficients. On exit, A is * overwritten by the updated matrix. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least max(1,m). * Unchanged on exit. * * * * Level 2 Blas routine. * * -- Written on 30-August-1985. * Sven Hammarling, Nag Central Office. C REVISED 860623 C REVISED YYMMDD C BY R. J. HANSON, SANDIA NATIONAL LABS. * INTEGER I,IX,J,JY,KX REAL ZERO PARAMETER (ZERO=0.0E+0) REAL TEMP LOGICAL OK OK = (M.GT.0) .AND. (N.GT.0) .AND. (LDA.GE.M) * * * Quick return if possible. * IF ( .NOT. OK .OR. (ALPHA.EQ.ZERO)) RETURN * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through A. * IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN DO 20,J = 1,N IF (Y(J).NE.ZERO) THEN TEMP = ALPHA*Y(J) DO 10,I = 1,M A(I,J) = A(I,J) + X(I)*TEMP 10 CONTINUE END IF * 20 CONTINUE * ELSE IF (INCX.GT.0) THEN KX = 1 * ELSE KX = 1 - (M-1)*INCX END IF * IF (INCY.GT.0) THEN JY = 1 * ELSE JY = 1 - (N-1)*INCY END IF * DO 40,J = 1,N IF (Y(JY).NE.ZERO) THEN TEMP = ALPHA*Y(JY) IX = KX DO 30,I = 1,M A(I,J) = A(I,J) + X(IX)*TEMP IX = IX + INCX 30 CONTINUE END IF * JY = JY + INCY 40 CONTINUE END IF * RETURN * * End of SGER . * END SUBROUTINE SSYR(UPLO,N,ALPHA,X,INCX,A,LDA) CHARACTER *1 UPLO INTEGER N,INCX,LDA REAL ALPHA,X(*),A(LDA,*) * * Purpose * ======= * * SSYR performs the symmetric rank 1 operation * * A := alpha*x*x' + A, * * where alpha is a real scalar, x is an n element vector and A is an * n by n symmetric matrix. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the upper or lower * triangular part of the array A is to be referenced as * follows: * * UPLO = 'U' Only the upper triangular part of A * is to be referenced. * * UPLO = 'L' Only the lower triangular part of A * is to be referenced. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - REAL . * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * X - REAL array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. * Unchanged on exit. * * A - REAL array of DIMENSION ( LDA, n ). * Before entry with UPLO = 'U', the leading n by n * upper triangular part of the array A must contain the upper * triangular part of the symmetric matrix and the strictly * lower triangular part of A is not referenced. On exit, the * upper triangular part of the array A is overwritten by the * upper triangular part of the updated matrix. * Before entry with UPLO = 'L', the leading n by n * lower triangular part of the array A must contain the lower * triangular part of the symmetric matrix and the strictly * upper triangular part of A is not referenced. On exit, the * lower triangular part of the array A is overwritten by the * lower triangular part of the updated matrix. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least max(1,n). * Unchanged on exit. * * * * * * Level 2 Blas routine. * * -- Written on 27-September-1985. * Sven Hammarling, Nag Central Office. C REVISED 860623 C REVISED YYMMDD C BY R. J. HANSON, SANDIA NATIONAL LABS. * INTEGER I,IX,J,JX,KX REAL ZERO PARAMETER (ZERO=0.0E+0) REAL TEMP LOGICAL OK,LSAME OK = (LSAME(UPLO,'U') .OR. LSAME(UPLO,'L')) .AND. (N.GT.0) .AND. . (LDA.GE.N) * * Quick return if possible. * IF ( .NOT. OK .OR. (ALPHA.EQ.ZERO)) RETURN * * Set the start point in X if the increment is not unity. * IF (INCX.LE.0) THEN KX = 1 - (N-1)*INCX * ELSE IF (INCX.NE.1) THEN KX = 1 END IF * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through the triangular part * of A. * IF (LSAME(UPLO,'U')) THEN * * Form A when A is stored in upper triangle. * IF (INCX.EQ.1) THEN DO 20,J = 1,N IF (X(J).NE.ZERO) THEN TEMP = ALPHA*X(J) DO 10,I = 1,J A(I,J) = A(I,J) + X(I)*TEMP 10 CONTINUE END IF * 20 CONTINUE * ELSE JX = KX DO 40,J = 1,N IF (X(JX).NE.ZERO) THEN TEMP = ALPHA*X(JX) IX = KX DO 30,I = 1,J A(I,J) = A(I,J) + X(IX)*TEMP IX = IX + INCX 30 CONTINUE END IF * JX = JX + INCX 40 CONTINUE END IF * ELSE * * Form A when A is stored in lower triangle. * IF (INCX.EQ.1) THEN DO 60,J = 1,N IF (X(J).NE.ZERO) THEN TEMP = ALPHA*X(J) DO 50,I = J,N A(I,J) = A(I,J) + X(I)*TEMP 50 CONTINUE END IF * 60 CONTINUE * ELSE JX = KX DO 80,J = 1,N IF (X(JX).NE.ZERO) THEN TEMP = ALPHA*X(JX) IX = JX DO 70,I = J,N A(I,J) = A(I,J) + X(IX)*TEMP IX = IX + INCX 70 CONTINUE END IF * JX = JX + INCX 80 CONTINUE END IF * END IF * RETURN * * End of SSYR . * END SUBROUTINE SSPR(UPLO,N,ALPHA,X,INCX,AP) CHARACTER *1 UPLO INTEGER N,INCX REAL ALPHA,X(*),AP(*) * * Purpose * ======= * * SSPR performs the symmetric rank 1 operation * * A := alpha*x*x' + A, * * where alpha is a real scalar, x is an n element vector and A is an * n by n symmetric matrix. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the upper or lower * triangular part of the matrix A is supplied in the packed * array AP as follows: * * UPLO = 'U' The upper triangular part of A is * supplied in AP. * * UPLO = 'L' The lower triangular part of A is * supplied in AP. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - REAL . * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * X - REAL array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. * Unchanged on exit. * * AP - REAL array of DIMENSION at least * ( ( n*( n + 1 ) )/2 ). * Before entry with UPLO = 'U', the array AP must * contain the upper triangular part of the symmetric matrix * packed sequentially, column by column, so that AP( 1 ) * contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) * and a( 2, 2 ) respectively, and so on. On exit, the array * AP is overwritten by the upper triangular part of the * updated matrix. * Before entry with UPLO = 'L', the array AP must * contain the lower triangular part of the symmetric matrix * packed sequentially, column by column, so that AP( 1 ) * contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) * and a( 3, 1 ) respectively, and so on. On exit, the array * AP is overwritten by the lower triangular part of the * updated matrix. * * * * * Level 2 Blas routine. * * -- Written on 30-September-1985. * Sven Hammarling, Nag Central Office. C REVISED 860623 C REVISED YYMMDD C BY R. J. HANSON, SANDIA NATIONAL LABS. * INTEGER I,IX,J,JX,K,KK INTEGER KX REAL ZERO PARAMETER (ZERO=0.0E+0) REAL TEMP LOGICAL OK,LSAME OK = (LSAME(UPLO,'U') .OR. LSAME(UPLO,'L')) .AND. (N.GT.0) * * Quick return if possible. * IF ( .NOT. OK .OR. (ALPHA.EQ.ZERO)) RETURN * * Set the start point in X if the increment is not unity. * IF (INCX.LE.0) THEN KX = 1 - (N-1)*INCX * ELSE IF (INCX.NE.1) THEN KX = 1 END IF * * Start the operations. In this version the elements of the array AP * are accessed sequentially with one pass through AP. * K = 1 IF (LSAME(UPLO,'U')) THEN * * Form A when upper triangle is stored in AP. * IF (INCX.EQ.1) THEN DO 20,J = 1,N IF (X(J).NE.ZERO) THEN TEMP = ALPHA*X(J) DO 10,I = 1,J AP(K) = AP(K) + X(I)*TEMP K = K + 1 10 CONTINUE * ELSE K = K + J END IF * 20 CONTINUE * ELSE JX = KX DO 40,J = 1,N IF (X(JX).NE.ZERO) THEN TEMP = ALPHA*X(JX) IX = KX KK = K DO 30,K = KK,KK + J - 1 AP(K) = AP(K) + X(IX)*TEMP IX = IX + INCX 30 CONTINUE * ELSE K = K + J END IF * JX = JX + INCX 40 CONTINUE END IF * ELSE * * Form A when lower triangle is stored in AP. * IF (INCX.EQ.1) THEN DO 60,J = 1,N IF (X(J).NE.ZERO) THEN TEMP = ALPHA*X(J) DO 50,I = J,N AP(K) = AP(K) + X(I)*TEMP K = K + 1 50 CONTINUE * ELSE K = K + N - J + 1 END IF * 60 CONTINUE * ELSE JX = KX DO 80,J = 1,N IF (X(JX).NE.ZERO) THEN TEMP = ALPHA*X(JX) IX = JX KK = K DO 70,K = KK,KK + N - J AP(K) = AP(K) + X(IX)*TEMP IX = IX + INCX 70 CONTINUE * ELSE K = K + N - J + 1 END IF * JX = JX + INCX 80 CONTINUE END IF * END IF * RETURN * * End of SSPR . * END SUBROUTINE SSYR2(UPLO,N,ALPHA,X,INCX,Y,INCY,A,LDA) CHARACTER *1 UPLO INTEGER N,INCX,INCY,LDA REAL ALPHA,X(*),Y(*),A(LDA,*) * * Purpose * ======= * * SSYR2 performs the symmetric rank 2 operation * * A := alpha*x*y' + alpha*y*x' + A, * * where alpha is a scalar, x and y are n element vectors and A is an n * by n symmetric matrix. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the upper or lower * triangular part of the array A is to be referenced as * follows: * * UPLO = 'U' Only the upper triangular part of A * is to be referenced. * * UPLO = 'L' Only the lower triangular part of A * is to be referenced. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - REAL . * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * X - REAL array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. * Unchanged on exit. * * Y - REAL array of dimension at least * ( 1 + ( n - 1 )*abs( INCY ) ). * Before entry, the incremented array Y must contain the n * element vector y. * Unchanged on exit. * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. * Unchanged on exit. * * A - REAL array of DIMENSION ( LDA, n ). * Before entry with UPLO = 'U', the leading n by n * upper triangular part of the array A must contain the upper * triangular part of the symmetric matrix and the strictly * lower triangular part of A is not referenced. On exit, the * upper triangular part of the array A is overwritten by the * upper triangular part of the updated matrix. * Before entry with UPLO = 'L', the leading n by n * lower triangular part of the array A must contain the lower * triangular part of the symmetric matrix and the strictly * upper triangular part of A is not referenced. On exit, the * lower triangular part of the array A is overwritten by the * lower triangular part of the updated matrix. * * LDA - INTEGER. * On entry, LDA specifies the first dimension of A as declared * in the calling (sub) program. LDA must be at least max(1,n). * Unchanged on exit. * * * * * * Level 2 Blas routine. * * -- Written on 27-September-1985. * Sven Hammarling, Nag Central Office. C REVISED 860623 C REVISED YYMMDD C BY R. J. HANSON, SANDIA NATIONAL LABS. * INTEGER I,IX,IY,J,JX,JY INTEGER KX,KY REAL ZERO PARAMETER (ZERO=0.0E+0) REAL TEMP1,TEMP2 LOGICAL OK,LSAME OK = (LSAME(UPLO,'U') .OR. LSAME(UPLO,'L')) .AND. (N.GT.0) .AND. . (LDA.GE.N) * * Quick return if possible. * IF ( .NOT. OK .OR. (ALPHA.EQ.ZERO)) RETURN * * Set up the start points in X and Y if the increments are not both * unity. * IF ((INCX.NE.1) .OR. (INCY.NE.1)) THEN IF (INCX.GT.0) THEN KX = 1 * ELSE KX = 1 - (N-1)*INCX END IF * IF (INCY.GT.0) THEN KY = 1 * ELSE KY = 1 - (N-1)*INCY END IF * END IF * * Start the operations. In this version the elements of A are * accessed sequentially with one pass through the triangular part * of A. * IF (LSAME(UPLO,'U')) THEN * * Form A when A is stored in the upper triangle. * IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN DO 20,J = 1,N IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN TEMP1 = ALPHA*Y(J) TEMP2 = ALPHA*X(J) DO 10,I = 1,J A(I,J) = A(I,J) + X(I)*TEMP1 + Y(I)*TEMP2 10 CONTINUE END IF * 20 CONTINUE * ELSE JX = KX JY = KY DO 40,J = 1,N IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN TEMP1 = ALPHA*Y(JY) TEMP2 = ALPHA*X(JX) IX = KX IY = KY DO 30,I = 1,J A(I,J) = A(I,J) + X(IX)*TEMP1 + Y(IY)*TEMP2 IX = IX + INCX IY = IY + INCY 30 CONTINUE END IF * JX = JX + INCX JY = JY + INCY 40 CONTINUE END IF * ELSE * * Form A when A is stored in the upper triangle. * IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN DO 60,J = 1,N IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN TEMP1 = ALPHA*Y(J) TEMP2 = ALPHA*X(J) DO 50,I = J,N A(I,J) = A(I,J) + X(I)*TEMP1 + Y(I)*TEMP2 50 CONTINUE END IF * 60 CONTINUE * ELSE JX = KX JY = KY DO 80,J = 1,N IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN TEMP1 = ALPHA*Y(JY) TEMP2 = ALPHA*X(JX) IX = JX IY = JY DO 70,I = J,N A(I,J) = A(I,J) + X(IX)*TEMP1 + Y(IY)*TEMP2 IX = IX + INCX IY = IY + INCY 70 CONTINUE END IF * JX = JX + INCX JY = JY + INCY 80 CONTINUE END IF * END IF * RETURN * * End of SSYR2 . * END SUBROUTINE SSPR2(UPLO,N,ALPHA,X,INCX,Y,INCY,AP) CHARACTER *1 UPLO INTEGER N,INCX,INCY REAL ALPHA,X(*),Y(*),AP(*) * * Purpose * ======= * * SSPR2 performs the symmetric rank 2 operation * * A := alpha*x*y' + alpha*y*x' + A, * * where alpha is a scalar, x and y are n element vectors and A is an * n by n symmetric matrix. * * Parameters * ========== * * UPLO - CHARACTER*1. * On entry, UPLO specifies whether the upper or lower * triangular part of the matrix A is supplied in the packed * array AP as follows: * * UPLO = 'U' The upper triangular part of A is * supplied in AP. * * UPLO = 'L' The lower triangular part of A is * supplied in AP. * * Unchanged on exit. * * N - INTEGER. * On entry, N specifies the order of the matrix A. * N must be at least zero. * Unchanged on exit. * * ALPHA - REAL . * On entry, ALPHA specifies the scalar alpha. * Unchanged on exit. * * X - REAL array of dimension at least * ( 1 + ( n - 1 )*abs( INCX ) ). * Before entry, the incremented array X must contain the n * element vector x. * Unchanged on exit. * * INCX - INTEGER. * On entry, INCX specifies the increment for the elements of * X. * Unchanged on exit. * * Y - REAL array of dimension at least * ( 1 + ( n - 1 )*abs( INCY ) ). * Before entry, the incremented array Y must contain the n * element vector y. * Unchanged on exit. * * INCY - INTEGER. * On entry, INCY specifies the increment for the elements of * Y. * Unchanged on exit. * * AP - REAL array of DIMENSION at least * ( ( n*( n + 1 ) )/2 ). * Before entry with UPLO = 'U', the array AP must * contain the upper triangular part of the symmetric matrix * packed sequentially, column by column, so that AP( 1 ) * contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 ) * and a( 2, 2 ) respectively, and so on. On exit, the array * AP is overwritten by the upper triangular part of the * updated matrix. * Before entry with UPLO = 'L', the array AP must * contain the lower triangular part of the symmetric matrix * packed sequentially, column by column, so that AP( 1 ) * contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 ) * and a( 3, 1 ) respectively, and so on. On exit, the array * AP is overwritten by the lower triangular part of the * updated matrix. * * * * * * Level 2 Blas routine. * * -- Written on 30-September-1985. * Sven Hammarling, Nag Central Office. C REVISED 860623 C REVISED YYMMDD C BY R. J. HANSON, SANDIA NATIONAL LABS. * INTEGER I,IX,IY,J,JX,JY INTEGER K,KK,KX,KY REAL ZERO PARAMETER (ZERO=0.0E+0) REAL TEMP1,TEMP2 LOGICAL OK,LSAME OK = (LSAME(UPLO,'U') .OR. LSAME(UPLO,'L')) .AND. (N.GT.0) * * Quick return if possible. * IF ( .NOT. OK .OR. (ALPHA.EQ.ZERO)) RETURN * * Set up the start points in X and Y if the increments are not both * unity. * IF ((INCX.NE.1) .OR. (INCY.NE.1)) THEN IF (INCX.GT.0) THEN KX = 1 * ELSE KX = 1 - (N-1)*INCX END IF * IF (INCY.GT.0) THEN KY = 1 * ELSE KY = 1 - (N-1)*INCY END IF * END IF * * Start the operations. In this version the elements of the array AP * are accessed sequentially with one pass through AP. * K = 1 IF (LSAME(UPLO,'U')) THEN * * Form A when upper triangle is stored in AP. * IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN DO 20,J = 1,N IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN TEMP1 = ALPHA*Y(J) TEMP2 = ALPHA*X(J) DO 10,I = 1,J AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2 K = K + 1 10 CONTINUE * ELSE K = K + J END IF * 20 CONTINUE * ELSE JX = KX JY = KY DO 40,J = 1,N IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN TEMP1 = ALPHA*Y(JY) TEMP2 = ALPHA*X(JX) IX = KX IY = KY KK = K DO 30,K = KK,KK + J - 1 AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2 IX = IX + INCX IY = IY + INCY 30 CONTINUE * ELSE K = K + J END IF * JX = JX + INCX JY = JY + INCY 40 CONTINUE END IF * ELSE * * Form A when lower triangle is stored in AP. * IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN DO 60,J = 1,N IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN TEMP1 = ALPHA*Y(J) TEMP2 = ALPHA*X(J) DO 50,I = J,N AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2 K = K + 1 50 CONTINUE * ELSE K = K + N - J + 1 END IF * 60 CONTINUE * ELSE JX = KX JY = KY DO 80,J = 1,N IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN TEMP1 = ALPHA*Y(JY) TEMP2 = ALPHA*X(JX) IX = JX IY = JY KK = K DO 70,K = KK,KK + N - J AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2 IX = IX + INCX IY = IY + INCY 70 CONTINUE * ELSE K = K + N - J + 1 END IF * JX = JX + INCX JY = JY + INCY 80 CONTINUE END IF * END IF * RETURN * * End of SSPR2 . * END LOGICAL FUNCTION LSAME(CA,CB) C TEST IF TWO CHARACTERS ARE ESSENTIALLY THE SAME. C THE CHARACTER CB IS ONE OF THE FORTRAN SET. C (LOWER AND UPPER CASE LETTERS ARE EQUIVALENT.) C THIS IS A SUBPROGRAM FOR THE LEVEL TWO BLAS. C REVISED 860623 C REVISED YYMMDD C AUTH=R. J. HANSON, SANDIA NATIONAL LABS. C C THIS SUBPROGRAM IS MACHINE-DEPENDENT. C VERSION FOR CDC SYSTEMS USING 6-12 BIT REPRESENTATIONS. CHARACTER CA(*) CHARACTER *1 CB INTEGER ICIRFX DATA ICIRFX/62/ C SEE IF THE FIRST CHAR. IN STRING CA EQUALS STRING CB. LSAME = CA(1) .EQ. CB .AND. CA(1) .NE. CHAR(ICIRFX) IF (LSAME) RETURN C THE CHARS. ARE NOT IDENTICAL. NOW CHECK THEM FOR EQUIVALENCE. C LOOK FOR THE 'ESCAPE' CHARACTER, CIRCUMFLEX, FOLLOWED BY C THE LETTER. IVAL = ICHAR(CA(2)) IF (IVAL.GE.ICHAR('A') .AND. IVAL.LE.ICHAR('Z')) THEN LSAME = CA(1) .EQ. CHAR(ICIRFX) .AND. CA(2) .EQ. CB END IF * RETURN C END C LOGICAL FUNCTION LSAME(CA,CB) C TEST IF TWO CHARACTERS ARE ESSENTIALLY THE SAME. C THE CHARACTER CB IS ONE OF THE FORTRAN SET. C (LOWER AND UPPER CASE LETTERS ARE EQUIVALENT.) C THIS IS A SUBPROGRAM FOR THE LEVEL TWO BLAS. C REVISED 860623 C REVISED YYMMDD C AUTH=R. J. HANSON, SANDIA NATIONAL LABS. C C THIS SUBPROGRAM IS MACHINE-DEPENDENT. C VERSION FOR ANY ASCII MACHINE. C CHARACTER *1 CA C CHARACTER *1 CB C INTEGER IOFF C DATA IOFF/32/ C SEE IF THE CHAR. IN STRING CA EQUALS STRING CB. C LSAME = CA .EQ. CB C IF (LSAME) RETURN C THE CHARS. ARE NOT IDENTICAL. NOW CHECK THEM FOR EQUIVALENCE. C ISHIFT = ICHAR(CA) - IOFF C IF (ISHIFT.GE.ICHAR('A') .AND. ISHIFT.LE.ICHAR('Z')) THEN C LSAME = ISHIFT .EQ. ICHAR(CB) C END IF C C RETURN C END C C LOGICAL FUNCTION LSAME(CA,CB) C TEST IF TWO CHARACTERS ARE ESSENTIALLY THE SAME. C THE CHARACTER CB IS ONE OF THE FORTRAN SET. C (LOWER AND UPPER CASE LETTERS ARE EQUIVALENT.) C THIS IS A SUBPROGRAM FOR THE LEVEL TWO BLAS. C REVISED 860623 C REVISED YYMMDD C AUTH=R. J. HANSON, SANDIA NATIONAL LABS. C C THIS SUBPROGRAM IS MACHINE-DEPENDENT. C VERSION FOR ANY EBCDIC MACHINE. C CHARACTER *1 CA C CHARACTER *1 CB C INTEGER IOFF C DATA IOFF/64/ C SEE IF THE CHAR. IN STRING CA EQUALS STRING CB. C LSAME = CA .EQ. CB C IF (LSAME) RETURN C THE CHARS. ARE NOT IDENTICAL. NOW CHECK THEM FOR EQUIVALENCE. C ISHIFT = ICHAR(CA) + IOFF C IF (ISHIFT.GE.ICHAR('A') .AND. ISHIFT.LE.ICHAR('I')) THEN C LSAME = ISHIFT .EQ. ICHAR(CB) C END IF C C IF (ISHIFT.GE.ICHAR('J') .AND. ISHIFT.LE.ICHAR('R')) THEN C LSAME = ISHIFT .EQ. ICHAR(CB) C END IF C C IF (ISHIFT.GE.ICHAR('S') .AND. ISHIFT.LE.ICHAR('Z')) THEN C LSAME = ISHIFT .EQ. ICHAR(CB) C END IF C C RETURN C END END