C ------------------- BELOW ARE BLAS-1 ROUTINES ------------------ CUT > IDAMAX.F <<'CUT HERE............' INTEGER FUNCTION IDAMAX(N,DX,INCX) C C FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE. C JACK DONGARRA, LINPACK, 3/11/78. C MODIFIED 3/93 TO RETURN IF INCX .LE. 0. C DOUBLE PRECISION DX(1),DMAX INTEGER I,INCX,IX,N C IDAMAX = 0 IF( N.LT.1 .OR. INCX.LE.0 ) RETURN IDAMAX = 1 IF(N.EQ.1)RETURN IF(INCX.EQ.1)GO TO 20 C C CODE FOR INCREMENT NOT EQUAL TO 1 C IX = 1 DMAX = DABS(DX(1)) IX = IX + INCX DO 10 I = 2,N IF(DABS(DX(IX)).LE.DMAX) GO TO 5 IDAMAX = I DMAX = DABS(DX(IX)) 5 IX = IX + INCX 10 CONTINUE RETURN C C CODE FOR INCREMENT EQUAL TO 1 C 20 DMAX = DABS(DX(1)) DO 30 I = 2,N IF(DABS(DX(I)).LE.DMAX) GO TO 30 IDAMAX = I DMAX = DABS(DX(I)) 30 CONTINUE RETURN END CUT HERE............ CAT > DSWAP.F <<'CUT HERE............' SUBROUTINE DSWAP (N,DX,INCX,DY,INCY) C C INTERCHANGES TWO VECTORS. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL ONE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(1),DY(1),DTEMP INTEGER I,INCX,INCY,IX,IY,M,MP1,N C IF(N.LE.0)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL C TO 1 C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N DTEMP = DX(IX) DX(IX) = DY(IY) DY(IY) = DTEMP IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,3) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DTEMP = DX(I) DX(I) = DY(I) DY(I) = DTEMP 30 CONTINUE IF( N .LT. 3 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,3 DTEMP = DX(I) DX(I) = DY(I) DY(I) = DTEMP DTEMP = DX(I + 1) DX(I + 1) = DY(I + 1) DY(I + 1) = DTEMP DTEMP = DX(I + 2) DX(I + 2) = DY(I + 2) DY(I + 2) = DTEMP 50 CONTINUE RETURN END CUT HERE............ CAT > DSCAL.F <<'CUT HERE............' SUBROUTINE DSCAL(N,DA,DX,INCX) C C SCALES A VECTOR BY A CONSTANT. C USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C MODIFIED 3/93 TO RETURN IF INCX .LE. 0. C DOUBLE PRECISION DA,DX(1) INTEGER I,INCX,M,MP1,N,NINCX C IF( N.LE.0 .OR. INCX.LE.0 )RETURN IF(INCX.EQ.1)GO TO 20 C C CODE FOR INCREMENT NOT EQUAL TO 1 C NINCX = N*INCX DO 10 I = 1,NINCX,INCX DX(I) = DA*DX(I) 10 CONTINUE RETURN C C CODE FOR INCREMENT EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,5) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DX(I) = DA*DX(I) 30 CONTINUE IF( N .LT. 5 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,5 DX(I) = DA*DX(I) DX(I + 1) = DA*DX(I + 1) DX(I + 2) = DA*DX(I + 2) DX(I + 3) = DA*DX(I + 3) DX(I + 4) = DA*DX(I + 4) 50 CONTINUE RETURN END CUT HERE............ CAT > DROTG.F <<'CUT HERE............' SUBROUTINE DROTG(DA,DB,C,S) C C CONSTRUCT GIVENS PLANE ROTATION. C JACK DONGARRA, LINPACK, 3/11/78. C MODIFIED 9/27/86. C DOUBLE PRECISION DA,DB,C,S,ROE,SCALE,R,Z C ROE = DB IF( DABS(DA) .GT. DABS(DB) ) ROE = DA SCALE = DABS(DA) + DABS(DB) IF( SCALE .NE. 0.0D0 ) GO TO 10 C = 1.0D0 S = 0.0D0 R = 0.0D0 GO TO 20 10 R = SCALE*DSQRT((DA/SCALE)**2 + (DB/SCALE)**2) R = DSIGN(1.0D0,ROE)*R C = DA/R S = DB/R 20 Z = S IF( DABS(C) .GT. 0.0D0 .AND. DABS(C) .LE. S ) Z = 1.0D0/C DA = R DB = Z RETURN END CUT HERE............ CAT > DROT.F <<'CUT HERE............' SUBROUTINE DROT (N,DX,INCX,DY,INCY,C,S) C C APPLIES A PLANE ROTATION. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(1),DY(1),DTEMP,C,S INTEGER I,INCX,INCY,IX,IY,N C IF(N.LE.0)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL C TO 1 C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N DTEMP = C*DX(IX) + S*DY(IY) DY(IY) = C*DY(IY) - S*DX(IX) DX(IX) = DTEMP IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C 20 DO 30 I = 1,N DTEMP = C*DX(I) + S*DY(I) DY(I) = C*DY(I) - S*DX(I) DX(I) = DTEMP 30 CONTINUE RETURN END CUT HERE............ CAT > DNRM2.F <<'CUT HERE............' DOUBLE PRECISION FUNCTION DNRM2 ( N, DX, INCX) INTEGER I, INCX, IX, J, N, NEXT DOUBLE PRECISION DX(1), CUTLO, CUTHI, HITEST, SUM, XMAX,ZERO,ONE DATA ZERO, ONE /0.0D0, 1.0D0/ C C EUCLIDEAN NORM OF THE N-VECTOR STORED IN DX() WITH STORAGE C INCREMENT INCX . C IF N .LE. 0 RETURN WITH RESULT = 0. C IF N .GE. 1 THEN INCX MUST BE .GE. 1 C C C.L.LAWSON, 1978 JAN 08 C MODIFIED TO CORRECT FAILURE TO UPDATE IX, 1/25/92. C MODIFIED 3/93 TO RETURN IF INCX .LE. 0. C C FOUR PHASE METHOD USING TWO BUILT-IN CONSTANTS THAT ARE C HOPEFULLY APPLICABLE TO ALL MACHINES. C CUTLO = MAXIMUM OF DSQRT(U/EPS) OVER ALL KNOWN MACHINES. C CUTHI = MINIMUM OF DSQRT(V) OVER ALL KNOWN MACHINES. C WHERE C EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1. C U = SMALLEST POSITIVE NO. (UNDERFLOW LIMIT) C V = LARGEST NO. (OVERFLOW LIMIT) C C BRIEF OUTLINE OF ALGORITHM.. C C PHASE 1 SCANS ZERO COMPONENTS. C MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO C MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO C MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M C WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX. C C VALUES FOR CUTLO AND CUTHI.. C FROM THE ENVIRONMENTAL PARAMETERS LISTED IN THE IMSL CONVERTER C DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS.. C CUTLO, S.P. U/EPS = 2**(-102) FOR HONEYWELL. CLOSE SECONDS ARE C UNIVAC AND DEC AT 2**(-103) C THUS CUTLO = 2**(-51) = 4.44089E-16 C CUTHI, S.P. V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC. C THUS CUTHI = 2**(63.5) = 1.30438E19 C CUTLO, D.P. U/EPS = 2**(-67) FOR HONEYWELL AND DEC. C THUS CUTLO = 2**(-33.5) = 8.23181D-11 C CUTHI, D.P. SAME AS S.P. CUTHI = 1.30438D19 C DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 / C DATA CUTLO, CUTHI / 4.441E-16, 1.304E19 / DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 / C IF(N .GT. 0 .AND. INCX.GT.0) GO TO 10 DNRM2 = ZERO GO TO 300 C 10 ASSIGN 30 TO NEXT SUM = ZERO I = 1 IX = 1 C BEGIN MAIN LOOP 20 GO TO NEXT,(30, 50, 70, 110) 30 IF( DABS(DX(I)) .GT. CUTLO) GO TO 85 ASSIGN 50 TO NEXT XMAX = ZERO C C PHASE 1. SUM IS ZERO C 50 IF( DX(I) .EQ. ZERO) GO TO 200 IF( DABS(DX(I)) .GT. CUTLO) GO TO 85 C C PREPARE FOR PHASE 2. ASSIGN 70 TO NEXT GO TO 105 C C PREPARE FOR PHASE 4. C 100 CONTINUE IX = J ASSIGN 110 TO NEXT SUM = (SUM / DX(I)) / DX(I) 105 XMAX = DABS(DX(I)) GO TO 115 C C PHASE 2. SUM IS SMALL. C SCALE TO AVOID DESTRUCTIVE UNDERFLOW. C 70 IF( DABS(DX(I)) .GT. CUTLO ) GO TO 75 C C COMMON CODE FOR PHASES 2 AND 4. C IN PHASE 4 SUM IS LARGE. SCALE TO AVOID OVERFLOW. C 110 IF( DABS(DX(I)) .LE. XMAX ) GO TO 115 SUM = ONE + SUM * (XMAX / DX(I))**2 XMAX = DABS(DX(I)) GO TO 200 C 115 SUM = SUM + (DX(I)/XMAX)**2 GO TO 200 C C C PREPARE FOR PHASE 3. C 75 SUM = (SUM * XMAX) * XMAX C C C FOR REAL OR D.P. SET HITEST = CUTHI/N C FOR COMPLEX SET HITEST = CUTHI/(2*N) C 85 HITEST = CUTHI/FLOAT( N ) C C PHASE 3. SUM IS MID-RANGE. NO SCALING. C DO 95 J = IX,N IF(DABS(DX(I)) .GE. HITEST) GO TO 100 SUM = SUM + DX(I)**2 I = I + INCX 95 CONTINUE DNRM2 = DSQRT( SUM ) GO TO 300 C 200 CONTINUE IX = IX + 1 I = I + INCX IF( IX .LE. N ) GO TO 20 C C END OF MAIN LOOP. C C COMPUTE SQUARE ROOT AND ADJUST FOR SCALING. C DNRM2 = XMAX * DSQRT(SUM) 300 CONTINUE RETURN END CUT HERE............ CAT > DMACH.F <<'CUT HERE............' DOUBLE PRECISION FUNCTION DMACH(JOB) INTEGER JOB C C SMACH COMPUTES MACHINE PARAMETERS OF FLOATING POINT C ARITHMETIC FOR USE IN TESTING ONLY. NOT REQUIRED BY C LINPACK PROPER. C C IF TROUBLE WITH AUTOMATIC COMPUTATION OF THESE QUANTITIES, C THEY CAN BE SET BY DIRECT ASSIGNMENT STATEMENTS. C ASSUME THE COMPUTER HAS C C B = BASE OF ARITHMETIC C T = NUMBER OF BASE B DIGITS C L = SMALLEST POSSIBLE EXPONENT C U = LARGEST POSSIBLE EXPONENT C C THEN C C EPS = B**(1-T) C TINY = 100.0*B**(-L+T) C HUGE = 0.01*B**(U-T) C C DMACH SAME AS SMACH EXCEPT T, L, U APPLY TO C DOUBLE PRECISION. C C CMACH SAME AS SMACH EXCEPT IF COMPLEX DIVISION C IS DONE BY C C 1/(X+I*Y) = (X-I*Y)/(X**2+Y**2) C C THEN C C TINY = SQRT(TINY) C HUGE = SQRT(HUGE) C C C JOB IS 1, 2 OR 3 FOR EPSILON, TINY AND HUGE, RESPECTIVELY. C DOUBLE PRECISION EPS,TINY,HUGE,S C EPS = 1.0D0 10 EPS = EPS/2.0D0 S = 1.0D0 + EPS IF (S .GT. 1.0D0) GO TO 10 EPS = 2.0D0*EPS C S = 1.0D0 20 TINY = S S = S/16.0D0 IF (S*1.0 .NE. 0.0D0) GO TO 20 TINY = (TINY/EPS)*100.0 HUGE = 1.0D0/TINY C IF (JOB .EQ. 1) DMACH = EPS IF (JOB .EQ. 2) DMACH = TINY IF (JOB .EQ. 3) DMACH = HUGE RETURN END CUT HERE............ CAT > DDOT.F <<'CUT HERE............' DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY) C C FORMS THE DOT PRODUCT OF TWO VECTORS. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(1),DY(1),DTEMP INTEGER I,INCX,INCY,IX,IY,M,MP1,N C DDOT = 0.0D0 DTEMP = 0.0D0 IF(N.LE.0)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N DTEMP = DTEMP + DX(IX)*DY(IY) IX = IX + INCX IY = IY + INCY 10 CONTINUE DDOT = DTEMP RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,5) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DTEMP = DTEMP + DX(I)*DY(I) 30 CONTINUE IF( N .LT. 5 ) GO TO 60 40 MP1 = M + 1 DO 50 I = MP1,N,5 DTEMP = DTEMP + DX(I)*DY(I) + DX(I + 1)*DY(I + 1) + * DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4) 50 CONTINUE 60 DDOT = DTEMP RETURN END CUT HERE............ CAT > DCOPY.F <<'CUT HERE............' SUBROUTINE DCOPY(N,DX,INCX,DY,INCY) C C COPIES A VECTOR, X, TO A VECTOR, Y. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(1),DY(1) INTEGER I,INCX,INCY,IX,IY,M,MP1,N C IF(N.LE.0)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N DY(IY) = DX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,7) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DY(I) = DX(I) 30 CONTINUE IF( N .LT. 7 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,7 DY(I) = DX(I) DY(I + 1) = DX(I + 1) DY(I + 2) = DX(I + 2) DY(I + 3) = DX(I + 3) DY(I + 4) = DX(I + 4) DY(I + 5) = DX(I + 5) DY(I + 6) = DX(I + 6) 50 CONTINUE RETURN END CUT HERE............ CAT > DAXPY.F <<'CUT HERE............' SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY) C C CONSTANT TIMES A VECTOR PLUS A VECTOR. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(1),DY(1),DA INTEGER I,INCX,INCY,IX,IY,M,MP1,N C IF(N.LE.0)RETURN IF (DA .EQ. 0.0D0) RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N DY(IY) = DY(IY) + DA*DX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,4) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DY(I) = DY(I) + DA*DX(I) 30 CONTINUE IF( N .LT. 4 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,4 DY(I) = DY(I) + DA*DX(I) DY(I + 1) = DY(I + 1) + DA*DX(I + 1) DY(I + 2) = DY(I + 2) + DA*DX(I + 2) DY(I + 3) = DY(I + 3) + DA*DX(I + 3) 50 CONTINUE RETURN END CUT HERE............ CAT > DASUM.F <<'CUT HERE............' DOUBLE PRECISION FUNCTION DASUM(N,DX,INCX) C C TAKES THE SUM OF THE ABSOLUTE VALUES. C JACK DONGARRA, LINPACK, 3/11/78. C MODIFIED 3/93 TO RETURN IF INCX .LE. 0. C DOUBLE PRECISION DX(1),DTEMP INTEGER I,INCX,M,MP1,N,NINCX C DASUM = 0.0D0 DTEMP = 0.0D0 IF( N.LE.0 .OR. INCX.LE.0 )RETURN IF(INCX.EQ.1)GO TO 20 C C CODE FOR INCREMENT NOT EQUAL TO 1 C NINCX = N*INCX DO 10 I = 1,NINCX,INCX DTEMP = DTEMP + DABS(DX(I)) 10 CONTINUE DASUM = DTEMP RETURN C C CODE FOR INCREMENT EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,6) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DTEMP = DTEMP + DABS(DX(I)) 30 CONTINUE IF( N .LT. 6 ) GO TO 60 40 MP1 = M + 1 DO 50 I = MP1,N,6 DTEMP = DTEMP + DABS(DX(I)) + DABS(DX(I + 1)) + DABS(DX(I + 2)) * + DABS(DX(I + 3)) + DABS(DX(I + 4)) + DABS(DX(I + 5)) 50 CONTINUE 60 DASUM = DTEMP RETURN END CUT HERE............ C --------------- BELOW ARE BLAS-2 ROUTINES ----------------------- CAT > DSPR2.F <<'CUT HERE............' * ************************************************************************ * SUBROUTINE DSPR2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, AP ) * .. SCALAR ARGUMENTS .. DOUBLE PRECISION ALPHA INTEGER INCX, INCY, N CHARACTER*1 UPLO * .. ARRAY ARGUMENTS .. DOUBLE PRECISION AP( * ), X( * ), Y( * ) * .. * * PURPOSE * ======= * * DSPR2 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, SUPPLIED IN PACKED FORM. * * 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' OR 'U' THE UPPER TRIANGULAR PART OF A IS * SUPPLIED IN AP. * * UPLO = 'L' OR '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 - DOUBLE PRECISION. * ON ENTRY, ALPHA SPECIFIES THE SCALAR ALPHA. * UNCHANGED ON EXIT. * * X - DOUBLE PRECISION 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. INCX MUST NOT BE ZERO. * UNCHANGED ON EXIT. * * Y - DOUBLE PRECISION 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. INCY MUST NOT BE ZERO. * UNCHANGED ON EXIT. * * AP - DOUBLE PRECISION ARRAY OF DIMENSION AT LEAST * ( ( N*( N + 1 ) )/2 ). * BEFORE ENTRY WITH UPLO = 'U' OR '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' OR '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 22-OCTOBER-1986. * JACK DONGARRA, ARGONNE NATIONAL LAB. * JEREMY DU CROZ, NAG CENTRAL OFFICE. * SVEN HAMMARLING, NAG CENTRAL OFFICE. * RICHARD HANSON, SANDIA NATIONAL LABS. * * * .. PARAMETERS .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) * .. LOCAL SCALARS .. DOUBLE PRECISION TEMP1, TEMP2 INTEGER I, INFO, IX, IY, J, JX, JY, K, KK, KX, KY * .. EXTERNAL FUNCTIONS .. LOGICAL LSAME EXTERNAL LSAME * .. EXTERNAL SUBROUTINES .. EXTERNAL XERBLA * .. * .. EXECUTABLE STATEMENTS .. * * TEST THE INPUT PARAMETERS. * INFO = 0 IF ( .NOT.LSAME( UPLO, 'U' ).AND. $ .NOT.LSAME( UPLO, 'L' ) )THEN INFO = 1 ELSE IF( N.LT.0 )THEN INFO = 2 ELSE IF( INCX.EQ.0 )THEN INFO = 5 ELSE IF( INCY.EQ.0 )THEN INFO = 7 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DSPR2 ', INFO ) RETURN END IF * * QUICK RETURN IF POSSIBLE. * IF( ( N.EQ.0 ).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 JX = KX JY = KY END IF * * START THE OPERATIONS. IN THIS VERSION THE ELEMENTS OF THE ARRAY AP * ARE ACCESSED SEQUENTIALLY WITH ONE PASS THROUGH AP. * KK = 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 ) K = KK DO 10, I = 1, J AP( K ) = AP( K ) + X( I )*TEMP1 + Y( I )*TEMP2 K = K + 1 10 CONTINUE END IF KK = KK + J 20 CONTINUE ELSE 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, K = KK, KK + J - 1 AP( K ) = AP( K ) + X( IX )*TEMP1 + Y( IY )*TEMP2 IX = IX + INCX IY = IY + INCY 30 CONTINUE END IF JX = JX + INCX JY = JY + INCY KK = KK + J 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 ) K = KK DO 50, I = J, N AP( K ) = AP( K ) + X( I )*TEMP1 + Y( I )*TEMP2 K = K + 1 50 CONTINUE END IF KK = KK + N - J + 1 60 CONTINUE ELSE 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, K = KK, KK + N - J AP( K ) = AP( K ) + X( IX )*TEMP1 + Y( IY )*TEMP2 IX = IX + INCX IY = IY + INCY 70 CONTINUE END IF JX = JX + INCX JY = JY + INCY KK = KK + N - J + 1 80 CONTINUE END IF END IF * RETURN * * END OF DSPR2 . * END CUT HERE............ CAT > DSYR2.F <<'CUT HERE............' * ************************************************************************ * SUBROUTINE DSYR2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA ) * .. SCALAR ARGUMENTS .. DOUBLE PRECISION ALPHA INTEGER INCX, INCY, LDA, N CHARACTER*1 UPLO * .. ARRAY ARGUMENTS .. DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) * .. * * PURPOSE * ======= * * DSYR2 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' OR 'U' ONLY THE UPPER TRIANGULAR PART OF A * IS TO BE REFERENCED. * * UPLO = 'L' OR '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 - DOUBLE PRECISION. * ON ENTRY, ALPHA SPECIFIES THE SCALAR ALPHA. * UNCHANGED ON EXIT. * * X - DOUBLE PRECISION 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. INCX MUST NOT BE ZERO. * UNCHANGED ON EXIT. * * Y - DOUBLE PRECISION 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. INCY MUST NOT BE ZERO. * UNCHANGED ON EXIT. * * A - DOUBLE PRECISION ARRAY OF DIMENSION ( LDA, N ). * BEFORE ENTRY WITH UPLO = 'U' OR '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' OR '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 22-OCTOBER-1986. * JACK DONGARRA, ARGONNE NATIONAL LAB. * JEREMY DU CROZ, NAG CENTRAL OFFICE. * SVEN HAMMARLING, NAG CENTRAL OFFICE. * RICHARD HANSON, SANDIA NATIONAL LABS. * * * .. PARAMETERS .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) * .. LOCAL SCALARS .. DOUBLE PRECISION TEMP1, TEMP2 INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY * .. EXTERNAL FUNCTIONS .. LOGICAL LSAME EXTERNAL LSAME * .. EXTERNAL SUBROUTINES .. EXTERNAL XERBLA * .. INTRINSIC FUNCTIONS .. INTRINSIC MAX * .. * .. EXECUTABLE STATEMENTS .. * * TEST THE INPUT PARAMETERS. * INFO = 0 IF ( .NOT.LSAME( UPLO, 'U' ).AND. $ .NOT.LSAME( UPLO, 'L' ) )THEN INFO = 1 ELSE IF( N.LT.0 )THEN INFO = 2 ELSE IF( INCX.EQ.0 )THEN INFO = 5 ELSE IF( INCY.EQ.0 )THEN INFO = 7 ELSE IF( LDA.LT.MAX( 1, N ) )THEN INFO = 9 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DSYR2 ', INFO ) RETURN END IF * * QUICK RETURN IF POSSIBLE. * IF( ( N.EQ.0 ).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 JX = KX JY = KY 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 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 LOWER 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 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 DSYR2 . * END CUT HERE............ CAT > DSPR.F <<'CUT HERE............' * ************************************************************************ * SUBROUTINE DSPR ( UPLO, N, ALPHA, X, INCX, AP ) * .. SCALAR ARGUMENTS .. DOUBLE PRECISION ALPHA INTEGER INCX, N CHARACTER*1 UPLO * .. ARRAY ARGUMENTS .. DOUBLE PRECISION AP( * ), X( * ) * .. * * PURPOSE * ======= * * DSPR 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, SUPPLIED IN PACKED FORM. * * 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' OR 'U' THE UPPER TRIANGULAR PART OF A IS * SUPPLIED IN AP. * * UPLO = 'L' OR '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 - DOUBLE PRECISION. * ON ENTRY, ALPHA SPECIFIES THE SCALAR ALPHA. * UNCHANGED ON EXIT. * * X - DOUBLE PRECISION 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. INCX MUST NOT BE ZERO. * UNCHANGED ON EXIT. * * AP - DOUBLE PRECISION ARRAY OF DIMENSION AT LEAST * ( ( N*( N + 1 ) )/2 ). * BEFORE ENTRY WITH UPLO = 'U' OR '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' OR '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 22-OCTOBER-1986. * JACK DONGARRA, ARGONNE NATIONAL LAB. * JEREMY DU CROZ, NAG CENTRAL OFFICE. * SVEN HAMMARLING, NAG CENTRAL OFFICE. * RICHARD HANSON, SANDIA NATIONAL LABS. * * * .. PARAMETERS .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) * .. LOCAL SCALARS .. DOUBLE PRECISION TEMP INTEGER I, INFO, IX, J, JX, K, KK, KX * .. EXTERNAL FUNCTIONS .. LOGICAL LSAME EXTERNAL LSAME * .. EXTERNAL SUBROUTINES .. EXTERNAL XERBLA * .. * .. EXECUTABLE STATEMENTS .. * * TEST THE INPUT PARAMETERS. * INFO = 0 IF ( .NOT.LSAME( UPLO, 'U' ).AND. $ .NOT.LSAME( UPLO, 'L' ) )THEN INFO = 1 ELSE IF( N.LT.0 )THEN INFO = 2 ELSE IF( INCX.EQ.0 )THEN INFO = 5 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DSPR ', INFO ) RETURN END IF * * QUICK RETURN IF POSSIBLE. * IF( ( N.EQ.0 ).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. * KK = 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 ) K = KK DO 10, I = 1, J AP( K ) = AP( K ) + X( I )*TEMP K = K + 1 10 CONTINUE END IF KK = KK + J 20 CONTINUE ELSE JX = KX DO 40, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) IX = KX DO 30, K = KK, KK + J - 1 AP( K ) = AP( K ) + X( IX )*TEMP IX = IX + INCX 30 CONTINUE END IF JX = JX + INCX KK = KK + J 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 ) K = KK DO 50, I = J, N AP( K ) = AP( K ) + X( I )*TEMP K = K + 1 50 CONTINUE END IF KK = KK + N - J + 1 60 CONTINUE ELSE JX = KX DO 80, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) IX = JX DO 70, K = KK, KK + N - J AP( K ) = AP( K ) + X( IX )*TEMP IX = IX + INCX 70 CONTINUE END IF JX = JX + INCX KK = KK + N - J + 1 80 CONTINUE END IF END IF * RETURN * * END OF DSPR . * END CUT HERE............ CAT > DSYR.F <<'CUT HERE............' * ************************************************************************ * SUBROUTINE DSYR ( UPLO, N, ALPHA, X, INCX, A, LDA ) * .. SCALAR ARGUMENTS .. DOUBLE PRECISION ALPHA INTEGER INCX, LDA, N CHARACTER*1 UPLO * .. ARRAY ARGUMENTS .. DOUBLE PRECISION A( LDA, * ), X( * ) * .. * * PURPOSE * ======= * * DSYR 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' OR 'U' ONLY THE UPPER TRIANGULAR PART OF A * IS TO BE REFERENCED. * * UPLO = 'L' OR '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 - DOUBLE PRECISION. * ON ENTRY, ALPHA SPECIFIES THE SCALAR ALPHA. * UNCHANGED ON EXIT. * * X - DOUBLE PRECISION 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. INCX MUST NOT BE ZERO. * UNCHANGED ON EXIT. * * A - DOUBLE PRECISION ARRAY OF DIMENSION ( LDA, N ). * BEFORE ENTRY WITH UPLO = 'U' OR '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' OR '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 22-OCTOBER-1986. * JACK DONGARRA, ARGONNE NATIONAL LAB. * JEREMY DU CROZ, NAG CENTRAL OFFICE. * SVEN HAMMARLING, NAG CENTRAL OFFICE. * RICHARD HANSON, SANDIA NATIONAL LABS. * * * .. PARAMETERS .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) * .. LOCAL SCALARS .. DOUBLE PRECISION TEMP INTEGER I, INFO, IX, J, JX, KX * .. EXTERNAL FUNCTIONS .. LOGICAL LSAME EXTERNAL LSAME * .. EXTERNAL SUBROUTINES .. EXTERNAL XERBLA * .. INTRINSIC FUNCTIONS .. INTRINSIC MAX * .. * .. EXECUTABLE STATEMENTS .. * * TEST THE INPUT PARAMETERS. * INFO = 0 IF ( .NOT.LSAME( UPLO, 'U' ).AND. $ .NOT.LSAME( UPLO, 'L' ) )THEN INFO = 1 ELSE IF( N.LT.0 )THEN INFO = 2 ELSE IF( INCX.EQ.0 )THEN INFO = 5 ELSE IF( LDA.LT.MAX( 1, N ) )THEN INFO = 7 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DSYR ', INFO ) RETURN END IF * * QUICK RETURN IF POSSIBLE. * IF( ( N.EQ.0 ).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 DSYR . * END CUT HERE............ CAT > DGER.F <<'CUT HERE............' * ************************************************************************ * SUBROUTINE DGER ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA ) * .. SCALAR ARGUMENTS .. DOUBLE PRECISION ALPHA INTEGER INCX, INCY, LDA, M, N * .. ARRAY ARGUMENTS .. DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) * .. * * PURPOSE * ======= * * DGER 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 - DOUBLE PRECISION. * ON ENTRY, ALPHA SPECIFIES THE SCALAR ALPHA. * UNCHANGED ON EXIT. * * X - DOUBLE PRECISION 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. INCX MUST NOT BE ZERO. * UNCHANGED ON EXIT. * * Y - DOUBLE PRECISION 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. INCY MUST NOT BE ZERO. * UNCHANGED ON EXIT. * * A - DOUBLE PRECISION 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 22-OCTOBER-1986. * JACK DONGARRA, ARGONNE NATIONAL LAB. * JEREMY DU CROZ, NAG CENTRAL OFFICE. * SVEN HAMMARLING, NAG CENTRAL OFFICE. * RICHARD HANSON, SANDIA NATIONAL LABS. * * * .. PARAMETERS .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) * .. LOCAL SCALARS .. DOUBLE PRECISION TEMP INTEGER I, INFO, IX, J, JY, KX * .. EXTERNAL SUBROUTINES .. EXTERNAL XERBLA * .. INTRINSIC FUNCTIONS .. INTRINSIC MAX * .. * .. EXECUTABLE STATEMENTS .. * * TEST THE INPUT PARAMETERS. * INFO = 0 IF ( M.LT.0 )THEN INFO = 1 ELSE IF( N.LT.0 )THEN INFO = 2 ELSE IF( INCX.EQ.0 )THEN INFO = 5 ELSE IF( INCY.EQ.0 )THEN INFO = 7 ELSE IF( LDA.LT.MAX( 1, M ) )THEN INFO = 9 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DGER ', INFO ) RETURN END IF * * QUICK RETURN IF POSSIBLE. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) ) $ RETURN * * START THE OPERATIONS. IN THIS VERSION THE ELEMENTS OF A ARE * ACCESSED SEQUENTIALLY WITH ONE PASS THROUGH A. * IF( INCY.GT.0 )THEN JY = 1 ELSE JY = 1 - ( N - 1 )*INCY END IF IF( INCX.EQ.1 )THEN DO 20, J = 1, N IF( Y( JY ).NE.ZERO )THEN TEMP = ALPHA*Y( JY ) DO 10, I = 1, M A( I, J ) = A( I, J ) + X( I )*TEMP 10 CONTINUE END IF JY = JY + INCY 20 CONTINUE ELSE IF( INCX.GT.0 )THEN KX = 1 ELSE KX = 1 - ( M - 1 )*INCX 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 DGER . * END CUT HERE............ CAT > DTPSV.F <<'CUT HERE............' * ************************************************************************ * SUBROUTINE DTPSV ( UPLO, TRANS, DIAG, N, AP, X, INCX ) * .. SCALAR ARGUMENTS .. INTEGER INCX, N CHARACTER*1 DIAG, TRANS, UPLO * .. ARRAY ARGUMENTS .. DOUBLE PRECISION AP( * ), X( * ) * .. * * PURPOSE * ======= * * DTPSV 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, SUPPLIED IN PACKED FORM. * * 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' OR 'U' A IS AN UPPER TRIANGULAR MATRIX. * * UPLO = 'L' OR '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' OR 'N' A*X = B. * * TRANS = 'T' OR 'T' A'*X = B. * * TRANS = 'C' OR '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' OR 'U' A IS ASSUMED TO BE UNIT TRIANGULAR. * * DIAG = 'N' OR '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 - DOUBLE PRECISION ARRAY OF DIMENSION AT LEAST * ( ( N*( N + 1 ) )/2 ). * BEFORE ENTRY WITH UPLO = 'U' OR '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' OR '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' OR 'U', THE DIAGONAL ELEMENTS OF * A ARE NOT REFERENCED, BUT ARE ASSUMED TO BE UNITY. * UNCHANGED ON EXIT. * * X - DOUBLE PRECISION 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. INCX MUST NOT BE ZERO. * UNCHANGED ON EXIT. * * * LEVEL 2 BLAS ROUTINE. * * -- WRITTEN ON 22-OCTOBER-1986. * JACK DONGARRA, ARGONNE NATIONAL LAB. * JEREMY DU CROZ, NAG CENTRAL OFFICE. * SVEN HAMMARLING, NAG CENTRAL OFFICE. * RICHARD HANSON, SANDIA NATIONAL LABS. * * * .. PARAMETERS .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) * .. LOCAL SCALARS .. DOUBLE PRECISION TEMP INTEGER I, INFO, IX, J, JX, K, KK, KX LOGICAL NOUNIT * .. EXTERNAL FUNCTIONS .. LOGICAL LSAME EXTERNAL LSAME * .. EXTERNAL SUBROUTINES .. EXTERNAL XERBLA * .. * .. EXECUTABLE STATEMENTS .. * * TEST THE INPUT PARAMETERS. * INFO = 0 IF ( .NOT.LSAME( UPLO , 'U' ).AND. $ .NOT.LSAME( UPLO , 'L' ) )THEN INFO = 1 ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND. $ .NOT.LSAME( TRANS, 'T' ).AND. $ .NOT.LSAME( TRANS, 'C' ) )THEN INFO = 2 ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND. $ .NOT.LSAME( DIAG , 'N' ) )THEN INFO = 3 ELSE IF( N.LT.0 )THEN INFO = 4 ELSE IF( INCX.EQ.0 )THEN INFO = 7 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DTPSV ', INFO ) RETURN END IF * * QUICK RETURN IF POSSIBLE. * IF( N.EQ.0 ) $ 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 KK = ( 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( KK ) TEMP = X( J ) K = KK - 1 DO 10, I = J - 1, 1, -1 X( I ) = X( I ) - TEMP*AP( K ) K = K - 1 10 CONTINUE END IF KK = KK - J 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( KK ) TEMP = X( JX ) IX = JX DO 30, K = KK - 1, KK - J + 1, -1 IX = IX - INCX X( IX ) = X( IX ) - TEMP*AP( K ) 30 CONTINUE END IF JX = JX - INCX KK = KK - J 40 CONTINUE END IF ELSE KK = 1 IF( INCX.EQ.1 )THEN DO 60, J = 1, N IF( X( J ).NE.ZERO )THEN IF( NOUNIT ) $ X( J ) = X( J )/AP( KK ) TEMP = X( J ) K = KK + 1 DO 50, I = J + 1, N X( I ) = X( I ) - TEMP*AP( K ) K = K + 1 50 CONTINUE END IF KK = KK + ( N - J + 1 ) 60 CONTINUE ELSE JX = KX DO 80, J = 1, N IF( X( JX ).NE.ZERO )THEN IF( NOUNIT ) $ X( JX ) = X( JX )/AP( KK ) TEMP = X( JX ) IX = JX DO 70, K = KK + 1, KK + N - J IX = IX + INCX X( IX ) = X( IX ) - TEMP*AP( K ) 70 CONTINUE END IF JX = JX + INCX KK = KK + ( N - J + 1 ) 80 CONTINUE END IF END IF ELSE * * FORM X := INV( A' )*X. * IF( LSAME( UPLO, 'U' ) )THEN KK = 1 IF( INCX.EQ.1 )THEN DO 100, J = 1, N TEMP = X( J ) K = KK DO 90, I = 1, J - 1 TEMP = TEMP - AP( K )*X( I ) K = K + 1 90 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/AP( KK + J - 1 ) X( J ) = TEMP KK = KK + J 100 CONTINUE ELSE JX = KX DO 120, J = 1, N TEMP = X( JX ) IX = KX DO 110, K = KK, KK + J - 2 TEMP = TEMP - AP( K )*X( IX ) IX = IX + INCX 110 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/AP( KK + J - 1 ) X( JX ) = TEMP JX = JX + INCX KK = KK + J 120 CONTINUE END IF ELSE KK = ( N*( N + 1 ) )/2 IF( INCX.EQ.1 )THEN DO 140, J = N, 1, -1 TEMP = X( J ) K = KK DO 130, I = N, J + 1, -1 TEMP = TEMP - AP( K )*X( I ) K = K - 1 130 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/AP( KK - N + J ) X( J ) = TEMP KK = KK - ( N - J + 1 ) 140 CONTINUE ELSE KX = KX + ( N - 1 )*INCX JX = KX DO 160, J = N, 1, -1 TEMP = X( JX ) IX = KX DO 150, K = KK, KK - ( N - ( J + 1 ) ), -1 TEMP = TEMP - AP( K )*X( IX ) IX = IX - INCX 150 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/AP( KK - N + J ) X( JX ) = TEMP JX = JX - INCX KK = KK - (N - J + 1 ) 160 CONTINUE END IF END IF END IF * RETURN * * END OF DTPSV . * END CUT HERE............ CAT > DTBSV.F <<'CUT HERE............' * ************************************************************************ * SUBROUTINE DTBSV ( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX ) * .. SCALAR ARGUMENTS .. INTEGER INCX, K, LDA, N CHARACTER*1 DIAG, TRANS, UPLO * .. ARRAY ARGUMENTS .. DOUBLE PRECISION A( LDA, * ), X( * ) * .. * * PURPOSE * ======= * * DTBSV 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' OR 'U' A IS AN UPPER TRIANGULAR MATRIX. * * UPLO = 'L' OR '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' OR 'N' A*X = B. * * TRANS = 'T' OR 'T' A'*X = B. * * TRANS = 'C' OR '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' OR 'U' A IS ASSUMED TO BE UNIT TRIANGULAR. * * DIAG = 'N' OR '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' OR 'U', K SPECIFIES THE NUMBER OF * SUPER-DIAGONALS OF THE MATRIX A. * ON ENTRY WITH UPLO = 'L' OR 'L', K SPECIFIES THE NUMBER OF * SUB-DIAGONALS OF THE MATRIX A. * K MUST SATISFY 0 .LE. K. * UNCHANGED ON EXIT. * * A - DOUBLE PRECISION ARRAY OF DIMENSION ( LDA, N ). * BEFORE ENTRY WITH UPLO = 'U' OR '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. * THE FOLLOWING PROGRAM SEGMENT WILL TRANSFER AN UPPER * TRIANGULAR BAND MATRIX FROM CONVENTIONAL FULL MATRIX STORAGE * TO BAND STORAGE: * * DO 20, J = 1, N * M = K + 1 - J * DO 10, I = MAX( 1, J - K ), J * A( M + I, J ) = MATRIX( I, J ) * 10 CONTINUE * 20 CONTINUE * * BEFORE ENTRY WITH UPLO = 'L' OR '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. * THE FOLLOWING PROGRAM SEGMENT WILL TRANSFER A LOWER * TRIANGULAR BAND MATRIX FROM CONVENTIONAL FULL MATRIX STORAGE * TO BAND STORAGE: * * DO 20, J = 1, N * M = 1 - J * DO 10, I = J, MIN( N, J + K ) * A( M + I, J ) = MATRIX( I, J ) * 10 CONTINUE * 20 CONTINUE * * NOTE THAT WHEN DIAG = 'U' OR '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 FIRST DIMENSION OF A AS DECLARED * IN THE CALLING (SUB) PROGRAM. LDA MUST BE AT LEAST * ( K + 1 ). * UNCHANGED ON EXIT. * * X - DOUBLE PRECISION 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. INCX MUST NOT BE ZERO. * UNCHANGED ON EXIT. * * * LEVEL 2 BLAS ROUTINE. * * -- WRITTEN ON 22-OCTOBER-1986. * JACK DONGARRA, ARGONNE NATIONAL LAB. * JEREMY DU CROZ, NAG CENTRAL OFFICE. * SVEN HAMMARLING, NAG CENTRAL OFFICE. * RICHARD HANSON, SANDIA NATIONAL LABS. * * * .. PARAMETERS .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) * .. LOCAL SCALARS .. DOUBLE PRECISION TEMP INTEGER I, INFO, IX, J, JX, KPLUS1, KX, L LOGICAL NOUNIT * .. EXTERNAL FUNCTIONS .. LOGICAL LSAME EXTERNAL LSAME * .. EXTERNAL SUBROUTINES .. EXTERNAL XERBLA * .. INTRINSIC FUNCTIONS .. INTRINSIC MAX, MIN * .. * .. EXECUTABLE STATEMENTS .. * * TEST THE INPUT PARAMETERS. * INFO = 0 IF ( .NOT.LSAME( UPLO , 'U' ).AND. $ .NOT.LSAME( UPLO , 'L' ) )THEN INFO = 1 ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND. $ .NOT.LSAME( TRANS, 'T' ).AND. $ .NOT.LSAME( TRANS, 'C' ) )THEN INFO = 2 ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND. $ .NOT.LSAME( DIAG , 'N' ) )THEN INFO = 3 ELSE IF( N.LT.0 )THEN INFO = 4 ELSE IF( K.LT.0 )THEN INFO = 5 ELSE IF( LDA.LT.( K + 1 ) )THEN INFO = 7 ELSE IF( INCX.EQ.0 )THEN INFO = 9 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DTBSV ', INFO ) RETURN END IF * * QUICK RETURN IF POSSIBLE. * IF( N.EQ.0 ) $ 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 L = KPLUS1 - J IF( NOUNIT ) $ X( J ) = X( J )/A( KPLUS1, J ) TEMP = X( J ) DO 10, I = J - 1, MAX( 1, J - K ), -1 X( I ) = X( I ) - TEMP*A( L + I, J ) 10 CONTINUE END IF 20 CONTINUE ELSE KX = KX + ( N - 1 )*INCX JX = KX DO 40, J = N, 1, -1 KX = KX - INCX IF( X( JX ).NE.ZERO )THEN IX = KX L = KPLUS1 - J IF( NOUNIT ) $ X( JX ) = X( JX )/A( KPLUS1, J ) TEMP = X( JX ) DO 30, I = J - 1, MAX( 1, J - K ), -1 X( IX ) = X( IX ) - TEMP*A( L + I, 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 L = 1 - J IF( NOUNIT ) $ X( J ) = X( J )/A( 1, J ) TEMP = X( J ) DO 50, I = J + 1, MIN( N, J + K ) X( I ) = X( I ) - TEMP*A( L + I, J ) 50 CONTINUE END IF 60 CONTINUE ELSE JX = KX DO 80, J = 1, N KX = KX + INCX IF( X( JX ).NE.ZERO )THEN IX = KX L = 1 - J IF( NOUNIT ) $ X( JX ) = X( JX )/A( 1, J ) TEMP = X( JX ) DO 70, I = J + 1, MIN( N, J + K ) X( IX ) = X( IX ) - TEMP*A( L + I, 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 TEMP = X( J ) L = KPLUS1 - J DO 90, I = MAX( 1, J - K ), J - 1 TEMP = TEMP - A( L + I, J )*X( I ) 90 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( KPLUS1, J ) X( J ) = TEMP 100 CONTINUE ELSE JX = KX DO 120, J = 1, N TEMP = X( JX ) IX = KX L = KPLUS1 - J DO 110, I = MAX( 1, J - K ), J - 1 TEMP = TEMP - A( L + I, J )*X( IX ) IX = IX + INCX 110 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( KPLUS1, J ) X( JX ) = TEMP 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 TEMP = X( J ) L = 1 - J DO 130, I = MIN( N, J + K ), J + 1, -1 TEMP = TEMP - A( L + I, J )*X( I ) 130 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( 1, J ) X( J ) = TEMP 140 CONTINUE ELSE KX = KX + ( N - 1 )*INCX JX = KX DO 160, J = N, 1, -1 TEMP = X( JX ) IX = KX L = 1 - J DO 150, I = MIN( N, J + K ), J + 1, -1 TEMP = TEMP - A( L + I, J )*X( IX ) IX = IX - INCX 150 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( 1, J ) X( JX ) = TEMP JX = JX - INCX IF( ( N - J ).GE.K ) $ KX = KX - INCX 160 CONTINUE END IF END IF END IF * RETURN * * END OF DTBSV . * END CUT HERE............ CAT > DTRSV.F <<'CUT HERE............' * ************************************************************************ * SUBROUTINE DTRSV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX ) * .. SCALAR ARGUMENTS .. INTEGER INCX, LDA, N CHARACTER*1 DIAG, TRANS, UPLO * .. ARRAY ARGUMENTS .. DOUBLE PRECISION A( LDA, * ), X( * ) * .. * * PURPOSE * ======= * * DTRSV 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' OR 'U' A IS AN UPPER TRIANGULAR MATRIX. * * UPLO = 'L' OR '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' OR 'N' A*X = B. * * TRANS = 'T' OR 'T' A'*X = B. * * TRANS = 'C' OR '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' OR 'U' A IS ASSUMED TO BE UNIT TRIANGULAR. * * DIAG = 'N' OR '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 - DOUBLE PRECISION ARRAY OF DIMENSION ( LDA, N ). * BEFORE ENTRY WITH UPLO = 'U' OR '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' OR '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' OR '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( 1, N ). * UNCHANGED ON EXIT. * * X - DOUBLE PRECISION 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. INCX MUST NOT BE ZERO. * UNCHANGED ON EXIT. * * * LEVEL 2 BLAS ROUTINE. * * -- WRITTEN ON 22-OCTOBER-1986. * JACK DONGARRA, ARGONNE NATIONAL LAB. * JEREMY DU CROZ, NAG CENTRAL OFFICE. * SVEN HAMMARLING, NAG CENTRAL OFFICE. * RICHARD HANSON, SANDIA NATIONAL LABS. * * * .. PARAMETERS .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) * .. LOCAL SCALARS .. DOUBLE PRECISION TEMP INTEGER I, INFO, IX, J, JX, KX LOGICAL NOUNIT * .. EXTERNAL FUNCTIONS .. LOGICAL LSAME EXTERNAL LSAME * .. EXTERNAL SUBROUTINES .. EXTERNAL XERBLA * .. INTRINSIC FUNCTIONS .. INTRINSIC MAX * .. * .. EXECUTABLE STATEMENTS .. * * TEST THE INPUT PARAMETERS. * INFO = 0 IF ( .NOT.LSAME( UPLO , 'U' ).AND. $ .NOT.LSAME( UPLO , 'L' ) )THEN INFO = 1 ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND. $ .NOT.LSAME( TRANS, 'T' ).AND. $ .NOT.LSAME( TRANS, 'C' ) )THEN INFO = 2 ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND. $ .NOT.LSAME( DIAG , 'N' ) )THEN INFO = 3 ELSE IF( N.LT.0 )THEN INFO = 4 ELSE IF( LDA.LT.MAX( 1, N ) )THEN INFO = 6 ELSE IF( INCX.EQ.0 )THEN INFO = 8 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DTRSV ', INFO ) RETURN END IF * * QUICK RETURN IF POSSIBLE. * IF( N.EQ.0 ) $ 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 ) TEMP = X( J ) DO 10, I = J - 1, 1, -1 X( I ) = X( I ) - TEMP*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 ) TEMP = X( JX ) IX = JX DO 30, I = J - 1, 1, -1 IX = IX - INCX X( IX ) = X( IX ) - TEMP*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 ) TEMP = X( J ) DO 50, I = J + 1, N X( I ) = X( 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 IF( NOUNIT ) $ X( JX ) = X( JX )/A( J, J ) TEMP = X( JX ) IX = JX DO 70, I = J + 1, N IX = IX + INCX X( IX ) = X( IX ) - TEMP*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 TEMP = X( J ) DO 90, I = 1, J - 1 TEMP = TEMP - A( I, J )*X( I ) 90 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( J, J ) X( J ) = TEMP 100 CONTINUE ELSE JX = KX DO 120, J = 1, N TEMP = X( JX ) IX = KX DO 110, I = 1, J - 1 TEMP = TEMP - A( I, J )*X( IX ) IX = IX + INCX 110 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( J, J ) X( JX ) = TEMP JX = JX + INCX 120 CONTINUE END IF ELSE IF( INCX.EQ.1 )THEN DO 140, J = N, 1, -1 TEMP = X( J ) DO 130, I = N, J + 1, -1 TEMP = TEMP - A( I, J )*X( I ) 130 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( J, J ) X( J ) = TEMP 140 CONTINUE ELSE KX = KX + ( N - 1 )*INCX JX = KX DO 160, J = N, 1, -1 TEMP = X( JX ) IX = KX DO 150, I = N, J + 1, -1 TEMP = TEMP - A( I, J )*X( IX ) IX = IX - INCX 150 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( J, J ) X( JX ) = TEMP JX = JX - INCX 160 CONTINUE END IF END IF END IF * RETURN * * END OF DTRSV . * END CUT HERE............ CAT > DTPMV.F <<'CUT HERE............' * ************************************************************************ * SUBROUTINE DTPMV ( UPLO, TRANS, DIAG, N, AP, X, INCX ) * .. SCALAR ARGUMENTS .. INTEGER INCX, N CHARACTER*1 DIAG, TRANS, UPLO * .. ARRAY ARGUMENTS .. DOUBLE PRECISION AP( * ), X( * ) * .. * * PURPOSE * ======= * * DTPMV PERFORMS ONE OF THE MATRIX-VECTOR OPERATIONS * * X := A*X, OR X := A'*X, * * WHERE X IS AN N ELEMENT VECTOR AND A IS AN N BY N UNIT, OR NON-UNIT, * UPPER OR LOWER TRIANGULAR MATRIX, SUPPLIED IN PACKED FORM. * * PARAMETERS * ========== * * UPLO - CHARACTER*1. * ON ENTRY, UPLO SPECIFIES WHETHER THE MATRIX IS AN UPPER OR * LOWER TRIANGULAR MATRIX AS FOLLOWS: * * UPLO = 'U' OR 'U' A IS AN UPPER TRIANGULAR MATRIX. * * UPLO = 'L' OR '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' OR 'N' X := A*X. * * TRANS = 'T' OR 'T' X := A'*X. * * TRANS = 'C' OR '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' OR 'U' A IS ASSUMED TO BE UNIT TRIANGULAR. * * DIAG = 'N' OR '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 - DOUBLE PRECISION ARRAY OF DIMENSION AT LEAST * ( ( N*( N + 1 ) )/2 ). * BEFORE ENTRY WITH UPLO = 'U' OR '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' OR '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' OR 'U', THE DIAGONAL ELEMENTS OF * A ARE NOT REFERENCED, BUT ARE ASSUMED TO BE UNITY. * UNCHANGED ON EXIT. * * X - DOUBLE PRECISION 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. INCX MUST NOT BE ZERO. * UNCHANGED ON EXIT. * * * LEVEL 2 BLAS ROUTINE. * * -- WRITTEN ON 22-OCTOBER-1986. * JACK DONGARRA, ARGONNE NATIONAL LAB. * JEREMY DU CROZ, NAG CENTRAL OFFICE. * SVEN HAMMARLING, NAG CENTRAL OFFICE. * RICHARD HANSON, SANDIA NATIONAL LABS. * * * .. PARAMETERS .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) * .. LOCAL SCALARS .. DOUBLE PRECISION TEMP INTEGER I, INFO, IX, J, JX, K, KK, KX LOGICAL NOUNIT * .. EXTERNAL FUNCTIONS .. LOGICAL LSAME EXTERNAL LSAME * .. EXTERNAL SUBROUTINES .. EXTERNAL XERBLA * .. * .. EXECUTABLE STATEMENTS .. * * TEST THE INPUT PARAMETERS. * INFO = 0 IF ( .NOT.LSAME( UPLO , 'U' ).AND. $ .NOT.LSAME( UPLO , 'L' ) )THEN INFO = 1 ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND. $ .NOT.LSAME( TRANS, 'T' ).AND. $ .NOT.LSAME( TRANS, 'C' ) )THEN INFO = 2 ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND. $ .NOT.LSAME( DIAG , 'N' ) )THEN INFO = 3 ELSE IF( N.LT.0 )THEN INFO = 4 ELSE IF( INCX.EQ.0 )THEN INFO = 7 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DTPMV ', INFO ) RETURN END IF * * QUICK RETURN IF POSSIBLE. * IF( N.EQ.0 ) $ 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 KK =1 IF( INCX.EQ.1 )THEN DO 20, J = 1, N IF( X( J ).NE.ZERO )THEN TEMP = X( J ) K = KK DO 10, I = 1, J - 1 X( I ) = X( I ) + TEMP*AP( K ) K = K + 1 10 CONTINUE IF( NOUNIT ) $ X( J ) = X( J )*AP( KK + J - 1 ) END IF KK = KK + J 20 CONTINUE ELSE JX = KX DO 40, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = X( JX ) IX = KX DO 30, K = KK, KK + J - 2 X( IX ) = X( IX ) + TEMP*AP( K ) IX = IX + INCX 30 CONTINUE IF( NOUNIT ) $ X( JX ) = X( JX )*AP( KK + J - 1 ) END IF JX = JX + INCX KK = KK + J 40 CONTINUE END IF ELSE KK = ( N*( N + 1 ) )/2 IF( INCX.EQ.1 )THEN DO 60, J = N, 1, -1 IF( X( J ).NE.ZERO )THEN TEMP = X( J ) K = KK DO 50, I = N, J + 1, -1 X( I ) = X( I ) + TEMP*AP( K ) K = K - 1 50 CONTINUE IF( NOUNIT ) $ X( J ) = X( J )*AP( KK - N + J ) END IF KK = KK - ( N - J + 1 ) 60 CONTINUE ELSE KX = KX + ( N - 1 )*INCX JX = KX DO 80, J = N, 1, -1 IF( X( JX ).NE.ZERO )THEN TEMP = X( JX ) IX = KX DO 70, K = KK, KK - ( N - ( J + 1 ) ), -1 X( IX ) = X( IX ) + TEMP*AP( K ) IX = IX - INCX 70 CONTINUE IF( NOUNIT ) $ X( JX ) = X( JX )*AP( KK - N + J ) END IF JX = JX - INCX KK = KK - ( N - J + 1 ) 80 CONTINUE END IF END IF ELSE * * FORM X := A'*X. * IF( LSAME( UPLO, 'U' ) )THEN KK = ( N*( N + 1 ) )/2 IF( INCX.EQ.1 )THEN DO 100, J = N, 1, -1 TEMP = X( J ) IF( NOUNIT ) $ TEMP = TEMP*AP( KK ) K = KK - 1 DO 90, I = J - 1, 1, -1 TEMP = TEMP + AP( K )*X( I ) K = K - 1 90 CONTINUE X( J ) = TEMP KK = KK - J 100 CONTINUE ELSE JX = KX + ( N - 1 )*INCX DO 120, J = N, 1, -1 TEMP = X( JX ) IX = JX IF( NOUNIT ) $ TEMP = TEMP*AP( KK ) DO 110, K = KK - 1, KK - J + 1, -1 IX = IX - INCX TEMP = TEMP + AP( K )*X( IX ) 110 CONTINUE X( JX ) = TEMP JX = JX - INCX KK = KK - J 120 CONTINUE END IF ELSE KK = 1 IF( INCX.EQ.1 )THEN DO 140, J = 1, N TEMP = X( J ) IF( NOUNIT ) $ TEMP = TEMP*AP( KK ) K = KK + 1 DO 130, I = J + 1, N TEMP = TEMP + AP( K )*X( I ) K = K + 1 130 CONTINUE X( J ) = TEMP KK = KK + ( N - J + 1 ) 140 CONTINUE ELSE JX = KX DO 160, J = 1, N TEMP = X( JX ) IX = JX IF( NOUNIT ) $ TEMP = TEMP*AP( KK ) DO 150, K = KK + 1, KK + N - J IX = IX + INCX TEMP = TEMP + AP( K )*X( IX ) 150 CONTINUE X( JX ) = TEMP JX = JX + INCX KK = KK + ( N - J + 1 ) 160 CONTINUE END IF END IF END IF * RETURN * * END OF DTPMV . * END CUT HERE............ CAT > DTBMV.F <<'CUT HERE............' * ************************************************************************ * SUBROUTINE DTBMV ( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX ) * .. SCALAR ARGUMENTS .. INTEGER INCX, K, LDA, N CHARACTER*1 DIAG, TRANS, UPLO * .. ARRAY ARGUMENTS .. DOUBLE PRECISION A( LDA, * ), X( * ) * .. * * PURPOSE * ======= * * DTBMV PERFORMS ONE OF THE MATRIX-VECTOR OPERATIONS * * X := A*X, OR X := A'*X, * * WHERE X IS AN 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' OR 'U' A IS AN UPPER TRIANGULAR MATRIX. * * UPLO = 'L' OR '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' OR 'N' X := A*X. * * TRANS = 'T' OR 'T' X := A'*X. * * TRANS = 'C' OR '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' OR 'U' A IS ASSUMED TO BE UNIT TRIANGULAR. * * DIAG = 'N' OR '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' OR 'U', K SPECIFIES THE NUMBER OF * SUPER-DIAGONALS OF THE MATRIX A. * ON ENTRY WITH UPLO = 'L' OR 'L', K SPECIFIES THE NUMBER OF * SUB-DIAGONALS OF THE MATRIX A. * K MUST SATISFY 0 .LE. K. * UNCHANGED ON EXIT. * * A - DOUBLE PRECISION ARRAY OF DIMENSION ( LDA, N ). * BEFORE ENTRY WITH UPLO = 'U' OR '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. * THE FOLLOWING PROGRAM SEGMENT WILL TRANSFER AN UPPER * TRIANGULAR BAND MATRIX FROM CONVENTIONAL FULL MATRIX STORAGE * TO BAND STORAGE: * * DO 20, J = 1, N * M = K + 1 - J * DO 10, I = MAX( 1, J - K ), J * A( M + I, J ) = MATRIX( I, J ) * 10 CONTINUE * 20 CONTINUE * * BEFORE ENTRY WITH UPLO = 'L' OR '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. * THE FOLLOWING PROGRAM SEGMENT WILL TRANSFER A LOWER * TRIANGULAR BAND MATRIX FROM CONVENTIONAL FULL MATRIX STORAGE * TO BAND STORAGE: * * DO 20, J = 1, N * M = 1 - J * DO 10, I = J, MIN( N, J + K ) * A( M + I, J ) = MATRIX( I, J ) * 10 CONTINUE * 20 CONTINUE * * NOTE THAT WHEN DIAG = 'U' OR '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 FIRST DIMENSION OF A AS DECLARED * IN THE CALLING (SUB) PROGRAM. LDA MUST BE AT LEAST * ( K + 1 ). * UNCHANGED ON EXIT. * * X - DOUBLE PRECISION 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. INCX MUST NOT BE ZERO. * UNCHANGED ON EXIT. * * * LEVEL 2 BLAS ROUTINE. * * -- WRITTEN ON 22-OCTOBER-1986. * JACK DONGARRA, ARGONNE NATIONAL LAB. * JEREMY DU CROZ, NAG CENTRAL OFFICE. * SVEN HAMMARLING, NAG CENTRAL OFFICE. * RICHARD HANSON, SANDIA NATIONAL LABS. * * * .. PARAMETERS .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) * .. LOCAL SCALARS .. DOUBLE PRECISION TEMP INTEGER I, INFO, IX, J, JX, KPLUS1, KX, L LOGICAL NOUNIT * .. EXTERNAL FUNCTIONS .. LOGICAL LSAME EXTERNAL LSAME * .. EXTERNAL SUBROUTINES .. EXTERNAL XERBLA * .. INTRINSIC FUNCTIONS .. INTRINSIC MAX, MIN * .. * .. EXECUTABLE STATEMENTS .. * * TEST THE INPUT PARAMETERS. * INFO = 0 IF ( .NOT.LSAME( UPLO , 'U' ).AND. $ .NOT.LSAME( UPLO , 'L' ) )THEN INFO = 1 ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND. $ .NOT.LSAME( TRANS, 'T' ).AND. $ .NOT.LSAME( TRANS, 'C' ) )THEN INFO = 2 ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND. $ .NOT.LSAME( DIAG , 'N' ) )THEN INFO = 3 ELSE IF( N.LT.0 )THEN INFO = 4 ELSE IF( K.LT.0 )THEN INFO = 5 ELSE IF( LDA.LT.( K + 1 ) )THEN INFO = 7 ELSE IF( INCX.EQ.0 )THEN INFO = 9 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DTBMV ', INFO ) RETURN END IF * * QUICK RETURN IF POSSIBLE. * IF( N.EQ.0 ) $ 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 TEMP = X( J ) L = KPLUS1 - J DO 10, I = MAX( 1, J - K ), J - 1 X( I ) = X( I ) + TEMP*A( L + I, J ) 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 TEMP = X( JX ) IX = KX L = KPLUS1 - J DO 30, I = MAX( 1, J - K ), J - 1 X( IX ) = X( IX ) + TEMP*A( L + I, 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 TEMP = X( J ) L = 1 - J DO 50, I = MIN( N, J + K ), J + 1, -1 X( I ) = X( I ) + TEMP*A( L + I, J ) 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 TEMP = X( JX ) IX = KX L = 1 - J DO 70, I = MIN( N, J + K ), J + 1, -1 X( IX ) = X( IX ) + TEMP*A( L + I, 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 TEMP = X( J ) L = KPLUS1 - J IF( NOUNIT ) $ TEMP = TEMP*A( KPLUS1, J ) DO 90, I = J - 1, MAX( 1, J - K ), -1 TEMP = TEMP + A( L + I, J )*X( I ) 90 CONTINUE X( J ) = TEMP 100 CONTINUE ELSE KX = KX + ( N - 1 )*INCX JX = KX DO 120, J = N, 1, -1 TEMP = X( JX ) KX = KX - INCX IX = KX L = KPLUS1 - J IF( NOUNIT ) $ TEMP = TEMP*A( KPLUS1, J ) DO 110, I = J - 1, MAX( 1, J - K ), -1 TEMP = TEMP + A( L + I, J )*X( IX ) IX = IX - INCX 110 CONTINUE X( JX ) = TEMP JX = JX - INCX 120 CONTINUE END IF ELSE IF( INCX.EQ.1 )THEN DO 140, J = 1, N TEMP = X( J ) L = 1 - J IF( NOUNIT ) $ TEMP = TEMP*A( 1, J ) DO 130, I = J + 1, MIN( N, J + K ) TEMP = TEMP + A( L + I, J )*X( I ) 130 CONTINUE X( J ) = TEMP 140 CONTINUE ELSE JX = KX DO 160, J = 1, N TEMP = X( JX ) KX = KX + INCX IX = KX L = 1 - J IF( NOUNIT ) $ TEMP = TEMP*A( 1, J ) DO 150, I = J + 1, MIN( N, J + K ) TEMP = TEMP + A( L + I, J )*X( IX ) IX = IX + INCX 150 CONTINUE X( JX ) = TEMP JX = JX + INCX 160 CONTINUE END IF END IF END IF * RETURN * * END OF DTBMV . * END CUT HERE............ CAT > DTRMV.F <<'CUT HERE............' * ************************************************************************ * SUBROUTINE DTRMV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX ) * .. SCALAR ARGUMENTS .. INTEGER INCX, LDA, N CHARACTER*1 DIAG, TRANS, UPLO * .. ARRAY ARGUMENTS .. DOUBLE PRECISION A( LDA, * ), X( * ) * .. * * PURPOSE * ======= * * DTRMV PERFORMS ONE OF THE MATRIX-VECTOR OPERATIONS * * X := A*X, OR X := A'*X, * * WHERE X IS AN 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' OR 'U' A IS AN UPPER TRIANGULAR MATRIX. * * UPLO = 'L' OR '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' OR 'N' X := A*X. * * TRANS = 'T' OR 'T' X := A'*X. * * TRANS = 'C' OR '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' OR 'U' A IS ASSUMED TO BE UNIT TRIANGULAR. * * DIAG = 'N' OR '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 - DOUBLE PRECISION ARRAY OF DIMENSION ( LDA, N ). * BEFORE ENTRY WITH UPLO = 'U' OR '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' OR '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' OR '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( 1, N ). * UNCHANGED ON EXIT. * * X - DOUBLE PRECISION 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. INCX MUST NOT BE ZERO. * UNCHANGED ON EXIT. * * * LEVEL 2 BLAS ROUTINE. * * -- WRITTEN ON 22-OCTOBER-1986. * JACK DONGARRA, ARGONNE NATIONAL LAB. * JEREMY DU CROZ, NAG CENTRAL OFFICE. * SVEN HAMMARLING, NAG CENTRAL OFFICE. * RICHARD HANSON, SANDIA NATIONAL LABS. * * * .. PARAMETERS .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) * .. LOCAL SCALARS .. DOUBLE PRECISION TEMP INTEGER I, INFO, IX, J, JX, KX LOGICAL NOUNIT * .. EXTERNAL FUNCTIONS .. LOGICAL LSAME EXTERNAL LSAME * .. EXTERNAL SUBROUTINES .. EXTERNAL XERBLA * .. INTRINSIC FUNCTIONS .. INTRINSIC MAX * .. * .. EXECUTABLE STATEMENTS .. * * TEST THE INPUT PARAMETERS. * INFO = 0 IF ( .NOT.LSAME( UPLO , 'U' ).AND. $ .NOT.LSAME( UPLO , 'L' ) )THEN INFO = 1 ELSE IF( .NOT.LSAME( TRANS, 'N' ).AND. $ .NOT.LSAME( TRANS, 'T' ).AND. $ .NOT.LSAME( TRANS, 'C' ) )THEN INFO = 2 ELSE IF( .NOT.LSAME( DIAG , 'U' ).AND. $ .NOT.LSAME( DIAG , 'N' ) )THEN INFO = 3 ELSE IF( N.LT.0 )THEN INFO = 4 ELSE IF( LDA.LT.MAX( 1, N ) )THEN INFO = 6 ELSE IF( INCX.EQ.0 )THEN INFO = 8 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DTRMV ', INFO ) RETURN END IF * * QUICK RETURN IF POSSIBLE. * IF( N.EQ.0 ) $ 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 TEMP = X( J ) DO 10, I = 1, J - 1 X( I ) = X( I ) + TEMP*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 TEMP = X( JX ) IX = KX DO 30, I = 1, J - 1 X( IX ) = X( IX ) + TEMP*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 TEMP = X( J ) DO 50, I = N, J + 1, -1 X( I ) = X( I ) + TEMP*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 TEMP = X( JX ) IX = KX DO 70, I = N, J + 1, -1 X( IX ) = X( IX ) + TEMP*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 TEMP = X( J ) IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 90, I = J - 1, 1, -1 TEMP = TEMP + A( I, J )*X( I ) 90 CONTINUE X( J ) = TEMP 100 CONTINUE ELSE JX = KX + ( N - 1 )*INCX DO 120, J = N, 1, -1 TEMP = X( JX ) IX = JX IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 110, I = J - 1, 1, -1 IX = IX - INCX TEMP = TEMP + A( I, J )*X( IX ) 110 CONTINUE X( JX ) = TEMP JX = JX - INCX 120 CONTINUE END IF ELSE IF( INCX.EQ.1 )THEN DO 140, J = 1, N TEMP = X( J ) IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 130, I = J + 1, N TEMP = TEMP + A( I, J )*X( I ) 130 CONTINUE X( J ) = TEMP 140 CONTINUE ELSE JX = KX DO 160, J = 1, N TEMP = X( JX ) IX = JX IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 150, I = J + 1, N IX = IX + INCX TEMP = TEMP + A( I, J )*X( IX ) 150 CONTINUE X( JX ) = TEMP JX = JX + INCX 160 CONTINUE END IF END IF END IF * RETURN * * END OF DTRMV . * END CUT HERE............ CAT > DSPMV.F <<'CUT HERE............' * ************************************************************************ * SUBROUTINE DSPMV ( UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY ) * .. SCALAR ARGUMENTS .. DOUBLE PRECISION ALPHA, BETA INTEGER INCX, INCY, N CHARACTER*1 UPLO * .. ARRAY ARGUMENTS .. DOUBLE PRECISION AP( * ), X( * ), Y( * ) * .. * * PURPOSE * ======= * * DSPMV 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, SUPPLIED IN PACKED FORM. * * 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' OR 'U' THE UPPER TRIANGULAR PART OF A IS * SUPPLIED IN AP. * * UPLO = 'L' OR '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 - DOUBLE PRECISION. * ON ENTRY, ALPHA SPECIFIES THE SCALAR ALPHA. * UNCHANGED ON EXIT. * * AP - DOUBLE PRECISION ARRAY OF DIMENSION AT LEAST * ( ( N*( N + 1 ) )/2 ). * BEFORE ENTRY WITH UPLO = 'U' OR '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' OR '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 - DOUBLE PRECISION 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. INCX MUST NOT BE ZERO. * UNCHANGED ON EXIT. * * BETA - DOUBLE PRECISION. * 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 - DOUBLE PRECISION 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. INCY MUST NOT BE ZERO. * UNCHANGED ON EXIT. * * * LEVEL 2 BLAS ROUTINE. * * -- WRITTEN ON 22-OCTOBER-1986. * JACK DONGARRA, ARGONNE NATIONAL LAB. * JEREMY DU CROZ, NAG CENTRAL OFFICE. * SVEN HAMMARLING, NAG CENTRAL OFFICE. * RICHARD HANSON, SANDIA NATIONAL LABS. * * * .. PARAMETERS .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. LOCAL SCALARS .. DOUBLE PRECISION TEMP1, TEMP2 INTEGER I, INFO, IX, IY, J, JX, JY, K, KK, KX, KY * .. EXTERNAL FUNCTIONS .. LOGICAL LSAME EXTERNAL LSAME * .. EXTERNAL SUBROUTINES .. EXTERNAL XERBLA * .. * .. EXECUTABLE STATEMENTS .. * * TEST THE INPUT PARAMETERS. * INFO = 0 IF ( .NOT.LSAME( UPLO, 'U' ).AND. $ .NOT.LSAME( UPLO, 'L' ) )THEN INFO = 1 ELSE IF( N.LT.0 )THEN INFO = 2 ELSE IF( INCX.EQ.0 )THEN INFO = 6 ELSE IF( INCY.EQ.0 )THEN INFO = 9 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DSPMV ', INFO ) RETURN END IF * * QUICK RETURN IF POSSIBLE. * IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * SET UP THE START POINTS IN X AND Y. * 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 * * 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. * IF( BETA.NE.ONE )THEN IF( INCY.EQ.1 )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 ELSE 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 KK = 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 K = KK 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( KK + J - 1 ) + ALPHA*TEMP2 KK = KK + J 60 CONTINUE ELSE JX = KX JY = KY DO 80, J = 1, N TEMP1 = ALPHA*X( JX ) TEMP2 = ZERO IX = KX IY = KY 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( JY ) = Y( JY ) + TEMP1*AP( KK + J - 1 ) + ALPHA*TEMP2 JX = JX + INCX JY = JY + INCY KK = KK + J 80 CONTINUE END IF ELSE * * FORM Y WHEN AP CONTAINS THE 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*AP( KK ) K = KK + 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 KK = KK + ( N - J + 1 ) 100 CONTINUE ELSE JX = KX JY = KY DO 120, J = 1, N TEMP1 = ALPHA*X( JX ) TEMP2 = ZERO Y( JY ) = Y( JY ) + TEMP1*AP( KK ) IX = JX IY = JY DO 110, K = KK + 1, KK + N - J 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 KK = KK + ( N - J + 1 ) 120 CONTINUE END IF END IF * RETURN * * END OF DSPMV . * END CUT HERE............ CAT > DSBMV.F <<'CUT HERE............' * ************************************************************************ * SUBROUTINE DSBMV ( UPLO, N, K, ALPHA, A, LDA, X, INCX, $ BETA, Y, INCY ) * .. SCALAR ARGUMENTS .. DOUBLE PRECISION ALPHA, BETA INTEGER INCX, INCY, K, LDA, N CHARACTER*1 UPLO * .. ARRAY ARGUMENTS .. DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) * .. * * PURPOSE * ======= * * DSBMV 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' OR 'U' THE UPPER TRIANGULAR PART OF A IS * BEING SUPPLIED. * * UPLO = 'L' OR '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. * UNCHANGED ON EXIT. * * ALPHA - DOUBLE PRECISION. * ON ENTRY, ALPHA SPECIFIES THE SCALAR ALPHA. * UNCHANGED ON EXIT. * * A - DOUBLE PRECISION ARRAY OF DIMENSION ( LDA, N ). * BEFORE ENTRY WITH UPLO = 'U' OR '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. * THE FOLLOWING PROGRAM SEGMENT WILL TRANSFER THE UPPER * TRIANGULAR PART OF A SYMMETRIC BAND MATRIX FROM CONVENTIONAL * FULL MATRIX STORAGE TO BAND STORAGE: * * DO 20, J = 1, N * M = K + 1 - J * DO 10, I = MAX( 1, J - K ), J * A( M + I, J ) = MATRIX( I, J ) * 10 CONTINUE * 20 CONTINUE * * BEFORE ENTRY WITH UPLO = 'L' OR '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. * THE FOLLOWING PROGRAM SEGMENT WILL TRANSFER THE LOWER * TRIANGULAR PART OF A SYMMETRIC BAND MATRIX FROM CONVENTIONAL * FULL MATRIX STORAGE TO BAND STORAGE: * * DO 20, J = 1, N * M = 1 - J * DO 10, I = J, MIN( N, J + K ) * A( M + I, J ) = MATRIX( I, J ) * 10 CONTINUE * 20 CONTINUE * * 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 * ( K + 1 ). * UNCHANGED ON EXIT. * * X - DOUBLE PRECISION 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. INCX MUST NOT BE ZERO. * UNCHANGED ON EXIT. * * BETA - DOUBLE PRECISION. * ON ENTRY, BETA SPECIFIES THE SCALAR BETA. * UNCHANGED ON EXIT. * * Y - DOUBLE PRECISION 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. INCY MUST NOT BE ZERO. * UNCHANGED ON EXIT. * * * LEVEL 2 BLAS ROUTINE. * * -- WRITTEN ON 22-OCTOBER-1986. * JACK DONGARRA, ARGONNE NATIONAL LAB. * JEREMY DU CROZ, NAG CENTRAL OFFICE. * SVEN HAMMARLING, NAG CENTRAL OFFICE. * RICHARD HANSON, SANDIA NATIONAL LABS. * * * .. PARAMETERS .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. LOCAL SCALARS .. DOUBLE PRECISION TEMP1, TEMP2 INTEGER I, INFO, IX, IY, J, JX, JY, KPLUS1, KX, KY, L * .. EXTERNAL FUNCTIONS .. LOGICAL LSAME EXTERNAL LSAME * .. EXTERNAL SUBROUTINES .. EXTERNAL XERBLA * .. INTRINSIC FUNCTIONS .. INTRINSIC MAX, MIN * .. * .. EXECUTABLE STATEMENTS .. * * TEST THE INPUT PARAMETERS. * INFO = 0 IF ( .NOT.LSAME( UPLO, 'U' ).AND. $ .NOT.LSAME( UPLO, 'L' ) )THEN INFO = 1 ELSE IF( N.LT.0 )THEN INFO = 2 ELSE IF( K.LT.0 )THEN INFO = 3 ELSE IF( LDA.LT.( K + 1 ) )THEN INFO = 6 ELSE IF( INCX.EQ.0 )THEN INFO = 8 ELSE IF( INCY.EQ.0 )THEN INFO = 11 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DSBMV ', INFO ) RETURN END IF * * QUICK RETURN IF POSSIBLE. * IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * SET UP THE START POINTS IN X AND Y. * 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 * * 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. * IF( BETA.NE.ONE )THEN IF( INCY.EQ.1 )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 ELSE 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 L = KPLUS1 - J DO 50, I = MAX( 1, J - K ), J - 1 Y( I ) = Y( I ) + TEMP1*A( L + I, J ) TEMP2 = TEMP2 + A( L + I, J )*X( I ) 50 CONTINUE Y( J ) = Y( J ) + TEMP1*A( KPLUS1, J ) + ALPHA*TEMP2 60 CONTINUE ELSE JX = KX JY = KY DO 80, J = 1, N TEMP1 = ALPHA*X( JX ) TEMP2 = ZERO IX = KX IY = KY L = KPLUS1 - J DO 70, I = MAX( 1, J - K ), J - 1 Y( IY ) = Y( IY ) + TEMP1*A( L + I, J ) TEMP2 = TEMP2 + A( L + I, J )*X( IX ) IX = IX + INCX IY = IY + INCY 70 CONTINUE Y( JY ) = Y( JY ) + TEMP1*A( KPLUS1, J ) + ALPHA*TEMP2 JX = JX + INCX JY = JY + INCY 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 ) L = 1 - J DO 90, I = J + 1, MIN( N, J + K ) Y( I ) = Y( I ) + TEMP1*A( L + I, J ) TEMP2 = TEMP2 + A( L + 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( 1, J ) L = 1 - J IX = JX IY = JY DO 110, I = J + 1, MIN( N, J + K ) IX = IX + INCX IY = IY + INCY Y( IY ) = Y( IY ) + TEMP1*A( L + I, J ) TEMP2 = TEMP2 + A( L + 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 DSBMV . * END CUT HERE............ CAT > DSYMV.F <<'CUT HERE............' * ************************************************************************ * SUBROUTINE DSYMV ( UPLO, N, ALPHA, A, LDA, X, INCX, $ BETA, Y, INCY ) * .. SCALAR ARGUMENTS .. DOUBLE PRECISION ALPHA, BETA INTEGER INCX, INCY, LDA, N CHARACTER*1 UPLO * .. ARRAY ARGUMENTS .. DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) * .. * * PURPOSE * ======= * * DSYMV 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' OR 'U' ONLY THE UPPER TRIANGULAR PART OF A * IS TO BE REFERENCED. * * UPLO = 'L' OR '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 - DOUBLE PRECISION. * ON ENTRY, ALPHA SPECIFIES THE SCALAR ALPHA. * UNCHANGED ON EXIT. * * A - DOUBLE PRECISION ARRAY OF DIMENSION ( LDA, N ). * BEFORE ENTRY WITH UPLO = 'U' OR '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' OR '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( 1, N ). * UNCHANGED ON EXIT. * * X - DOUBLE PRECISION 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. INCX MUST NOT BE ZERO. * UNCHANGED ON EXIT. * * BETA - DOUBLE PRECISION. * 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 - DOUBLE PRECISION 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. INCY MUST NOT BE ZERO. * UNCHANGED ON EXIT. * * * LEVEL 2 BLAS ROUTINE. * * -- WRITTEN ON 22-OCTOBER-1986. * JACK DONGARRA, ARGONNE NATIONAL LAB. * JEREMY DU CROZ, NAG CENTRAL OFFICE. * SVEN HAMMARLING, NAG CENTRAL OFFICE. * RICHARD HANSON, SANDIA NATIONAL LABS. * * * .. PARAMETERS .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. LOCAL SCALARS .. DOUBLE PRECISION TEMP1, TEMP2 INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY * .. EXTERNAL FUNCTIONS .. LOGICAL LSAME EXTERNAL LSAME * .. EXTERNAL SUBROUTINES .. EXTERNAL XERBLA * .. INTRINSIC FUNCTIONS .. INTRINSIC MAX * .. * .. EXECUTABLE STATEMENTS .. * * TEST THE INPUT PARAMETERS. * INFO = 0 IF ( .NOT.LSAME( UPLO, 'U' ).AND. $ .NOT.LSAME( UPLO, 'L' ) )THEN INFO = 1 ELSE IF( N.LT.0 )THEN INFO = 2 ELSE IF( LDA.LT.MAX( 1, N ) )THEN INFO = 5 ELSE IF( INCX.EQ.0 )THEN INFO = 7 ELSE IF( INCY.EQ.0 )THEN INFO = 10 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DSYMV ', INFO ) RETURN END IF * * QUICK RETURN IF POSSIBLE. * IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * SET UP THE START POINTS IN X AND Y. * 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 * * 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. * IF( BETA.NE.ONE )THEN IF( INCY.EQ.1 )TH