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 )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 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 JX = KX JY = KY DO 80, J = 1, N TEMP1 = ALPHA*X( JX ) 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( JY ) = Y( JY ) + TEMP1*A( J, J ) + ALPHA*TEMP2 JX = JX + INCX JY = JY + INCY 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 DSYMV . * END CUT HERE............ CAT > DGBMV.F <<'CUT HERE............' * ************************************************************************ * SUBROUTINE DGBMV ( TRANS, M, N, KL, KU, ALPHA, A, LDA, X, INCX, $ BETA, Y, INCY ) * .. SCALAR ARGUMENTS .. DOUBLE PRECISION ALPHA, BETA INTEGER INCX, INCY, KL, KU, LDA, M, N CHARACTER*1 TRANS * .. ARRAY ARGUMENTS .. DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) * .. * * PURPOSE * ======= * * DGBMV 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' OR 'N' Y := ALPHA*A*X + BETA*Y. * * TRANS = 'T' OR 'T' Y := ALPHA*A'*X + BETA*Y. * * TRANS = 'C' OR '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. * * ALPHA - DOUBLE PRECISION. * ON ENTRY, ALPHA SPECIFIES THE SCALAR ALPHA. * UNCHANGED ON EXIT. * * A - DOUBLE PRECISION 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. * 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. * THE FOLLOWING PROGRAM SEGMENT WILL TRANSFER A BAND MATRIX * FROM CONVENTIONAL FULL MATRIX STORAGE TO BAND STORAGE: * * 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( 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 * ( KL + KU + 1 ). * UNCHANGED ON EXIT. * * X - DOUBLE PRECISION ARRAY OF DIMENSION AT LEAST * ( 1 + ( N - 1 )*ABS( INCX ) ) WHEN TRANS = 'N' OR '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. 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 + ( M - 1 )*ABS( INCY ) ) WHEN TRANS = 'N' OR '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. 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 TEMP INTEGER I, INFO, IX, IY, J, JX, JY, K, KUP1, KX, KY, $ LENX, LENY * .. 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( TRANS, 'N' ).AND. $ .NOT.LSAME( TRANS, 'T' ).AND. $ .NOT.LSAME( TRANS, 'C' ) )THEN INFO = 1 ELSE IF( M.LT.0 )THEN INFO = 2 ELSE IF( N.LT.0 )THEN INFO = 3 ELSE IF( KL.LT.0 )THEN INFO = 4 ELSE IF( KU.LT.0 )THEN INFO = 5 ELSE IF( LDA.LT.( KL + KU + 1 ) )THEN INFO = 8 ELSE IF( INCX.EQ.0 )THEN INFO = 10 ELSE IF( INCY.EQ.0 )THEN INFO = 13 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DGBMV ', INFO ) RETURN END IF * * QUICK RETURN IF POSSIBLE. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. $ ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * SET LENX AND LENY, THE LENGTHS OF THE VECTORS X AND Y, AND SET * UP THE START POINTS IN X AND Y. * IF( LSAME( TRANS, 'N' ) )THEN LENX = N LENY = M ELSE LENX = M LENY = N END IF 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 * * 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. * IF( BETA.NE.ONE )THEN IF( INCY.EQ.1 )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 ELSE 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. * JX = KX IF( INCY.EQ.1 )THEN DO 60, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) 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 JX = JX + INCX 60 CONTINUE ELSE 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. * JY = KY IF( INCX.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( JY ) = Y( JY ) + ALPHA*TEMP JY = JY + INCY 100 CONTINUE ELSE 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 DGBMV . * END CUT HERE............ CAT > DGEMV.F <<'CUT HERE............' * ************************************************************************ * * FILE OF THE DOUBLE PRECISION LEVEL-2 BLAS. * =========================================== * * SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX, * $ BETA, Y, INCY ) * * SUBROUTINE DGBMV ( TRANS, M, N, KL, KU, ALPHA, A, LDA, X, INCX, * $ BETA, Y, INCY ) * * SUBROUTINE DSYMV ( UPLO, N, ALPHA, A, LDA, X, INCX, * $ BETA, Y, INCY ) * * SUBROUTINE DSBMV ( UPLO, N, K, ALPHA, A, LDA, X, INCX, * $ BETA, Y, INCY ) * * SUBROUTINE DSPMV ( UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY ) * * SUBROUTINE DTRMV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX ) * * SUBROUTINE DTBMV ( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX ) * * SUBROUTINE DTPMV ( UPLO, TRANS, DIAG, N, AP, X, INCX ) * * SUBROUTINE DTRSV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX ) * * SUBROUTINE DTBSV ( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX ) * * SUBROUTINE DTPSV ( UPLO, TRANS, DIAG, N, AP, X, INCX ) * * SUBROUTINE DGER ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA ) * * SUBROUTINE DSYR ( UPLO, N, ALPHA, X, INCX, A, LDA ) * * SUBROUTINE DSPR ( UPLO, N, ALPHA, X, INCX, AP ) * * SUBROUTINE DSYR2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA ) * * SUBROUTINE DSPR2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, AP ) * * SEE: * * DONGARRA J. J., DU CROZ J. J., HAMMARLING S. AND HANSON R. J.. * AN EXTENDED SET OF FORTRAN BASIC LINEAR ALGEBRA SUBPROGRAMS. * * TECHNICAL MEMORANDA NOS. 41 (REVISION 3) AND 81, MATHEMATICS * AND COMPUTER SCIENCE DIVISION, ARGONNE NATIONAL LABORATORY, * 9700 SOUTH CASS AVENUE, ARGONNE, ILLINOIS 60439, US. * * OR * * NAG TECHNICAL REPORTS TR3/87 AND TR4/87, NUMERICAL ALGORITHMS * GROUP LTD., NAG CENTRAL OFFICE, 256 BANBURY ROAD, OXFORD * OX2 7DE, UK, AND NUMERICAL ALGORITHMS GROUP INC., 1101 31ST * STREET, SUITE 100, DOWNERS GROVE, ILLINOIS 60515-1263, USA. * ************************************************************************ * SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX, $ BETA, Y, INCY ) * .. SCALAR ARGUMENTS .. DOUBLE PRECISION ALPHA, BETA INTEGER INCX, INCY, LDA, M, N CHARACTER*1 TRANS * .. ARRAY ARGUMENTS .. DOUBLE PRECISION A( LDA, * ), X( * ), Y( * ) * .. * * PURPOSE * ======= * * DGEMV 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' OR 'N' Y := ALPHA*A*X + BETA*Y. * * TRANS = 'T' OR 'T' Y := ALPHA*A'*X + BETA*Y. * * TRANS = 'C' OR '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 - DOUBLE PRECISION. * ON ENTRY, ALPHA SPECIFIES THE SCALAR ALPHA. * 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. * 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, M ). * UNCHANGED ON EXIT. * * X - DOUBLE PRECISION ARRAY OF DIMENSION AT LEAST * ( 1 + ( N - 1 )*ABS( INCX ) ) WHEN TRANS = 'N' OR '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. 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 + ( M - 1 )*ABS( INCY ) ) WHEN TRANS = 'N' OR '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. 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 TEMP INTEGER I, INFO, IX, IY, J, JX, JY, KX, KY, LENX, LENY * .. 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( TRANS, 'N' ).AND. $ .NOT.LSAME( TRANS, 'T' ).AND. $ .NOT.LSAME( TRANS, 'C' ) )THEN INFO = 1 ELSE IF( M.LT.0 )THEN INFO = 2 ELSE IF( N.LT.0 )THEN INFO = 3 ELSE IF( LDA.LT.MAX( 1, M ) )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( 'DGEMV ', INFO ) RETURN END IF * * QUICK RETURN IF POSSIBLE. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. $ ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * SET LENX AND LENY, THE LENGTHS OF THE VECTORS X AND Y, AND SET * UP THE START POINTS IN X AND Y. * IF( LSAME( TRANS, 'N' ) )THEN LENX = N LENY = M ELSE LENX = M LENY = N END IF 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 * * START THE OPERATIONS. IN THIS VERSION THE ELEMENTS OF 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, LENY Y( I ) = ZERO 10 CONTINUE ELSE DO 20, I = 1, LENY Y( I ) = BETA*Y( I ) 20 CONTINUE END IF ELSE 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. * JX = KX IF( INCY.EQ.1 )THEN DO 60, J = 1, N IF( X( JX ).NE.ZERO )THEN TEMP = ALPHA*X( JX ) DO 50, I = 1, M Y( I ) = Y( I ) + TEMP*A( I, J ) 50 CONTINUE END IF JX = JX + INCX 60 CONTINUE ELSE 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. * JY = KY IF( INCX.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( JY ) = Y( JY ) + ALPHA*TEMP JY = JY + INCY 100 CONTINUE ELSE 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 DGEMV . * END CUT HERE............ C FOLLOWING FROM V12X CODE FROM JMH SUBROUTINE XERBLA ( SRNAME, INFO ) * .. Scalar Arguments .. INTEGER INFO CHARACTER*6 SRNAME * .. * * Purpose * ======= * * XERBLA is an error handler for the Level 2 BLAS routines. * * It is called by the Level 2 BLAS routines if an input parameter is * invalid. * * Installers should consider modifying the STOP statement in order to * call system-specific exception-handling facilities. * * Parameters * ========== * * SRNAME - CHARACTER*6. * On entry, SRNAME specifies the name of the routine which * called XERBLA. * * INFO - INTEGER. * On entry, INFO specifies the position of the invalid * parameter in the parameter-list of the calling routine. * * * Auxiliary routine for Level 2 Blas. * * Written on 20-July-1986. * * .. Executable Statements .. * WRITE (*,99999) SRNAME, INFO * STOP * 99999 FORMAT ( ' ** On entry to ', A6, ' parameter number ', I2, $ ' had an illegal value' ) * * End of XERBLA. * END LOGICAL FUNCTION LSAME ( CA, CB ) * .. Scalar Arguments .. CHARACTER*1 CA, CB * .. * * Purpose * ======= * * LSAME tests if CA is the same letter as CB regardless of case. * CB is assumed to be an upper case letter. LSAME returns .TRUE. if * CA is either the same as CB or the equivalent lower case letter. * * N.B. This version of the routine is only correct for ASCII code. * Installers must modify the routine for other character-codes. * * For EBCDIC systems the constant IOFF must be changed to -64. * For CDC systems using 6-12 bit representations, the system- * specific code in comments must be activated. * * Parameters * ========== * * CA - CHARACTER*1 * CB - CHARACTER*1 * On entry, CA and CB specify characters to be compared. * Unchanged on exit. * * * Auxiliary routine for Level 2 Blas. * * -- Written on 20-July-1986 * Richard Hanson, Sandia National Labs. * Jeremy Du Croz, Nag Central Office. * * .. Parameters .. INTEGER IOFF PARAMETER ( IOFF=32 ) * .. Intrinsic Functions .. INTRINSIC ICHAR * .. Executable Statements .. * * Test if the characters are equal * LSAME = CA .EQ. CB * * Now test for equivalence * IF ( .NOT.LSAME ) THEN LSAME = ICHAR(CA) - IOFF .EQ. ICHAR(CB) END IF * RETURN * * The following comments contain code for CDC systems using 6-12 bit * representations. * * .. Parameters .. * INTEGER ICIRFX * PARAMETER ( ICIRFX=62 ) * .. Scalar Arguments .. * CHARACTER*1 CB * .. Array Arguments .. * CHARACTER*1 CA(*) * .. Local Scalars .. * INTEGER IVAL * .. Intrinsic Functions .. * INTRINSIC ICHAR, CHAR * .. Executable Statements .. * * See if the first character in string CA equals string CB. * * LSAME = CA(1) .EQ. CB .AND. CA(1) .NE. CHAR(ICIRFX) * * IF (LSAME) RETURN * * The characters are not identical. Now check them for equivalence. * Look for the 'escape' character, circumflex, followed by 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 * * End of LSAME. * END C ---------------- BELOW ARE BLAS-3 ROUTINES ------------- CUT > DGEMM.F <<'CUT HERE............' * ************************************************************************ * * FILE OF THE DOUBLE PRECISION LEVEL-3 BLAS. * ========================================== * * SUBROUTINE DGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, * $ BETA, C, LDC ) * * SUBROUTINE DSYMM ( SIDE, UPLO, M, N, ALPHA, A, LDA, B, LDB, * $ BETA, C, LDC ) * * SUBROUTINE DSYRK ( UPLO, TRANS, N, K, ALPHA, A, LDA, * $ BETA, C, LDC ) * * SUBROUTINE DSYR2K( UPLO, TRANS, N, K, ALPHA, A, LDA, B, LDB, * $ BETA, C, LDC ) * * SUBROUTINE DTRMM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, * $ B, LDB ) * * SUBROUTINE DTRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, * $ B, LDB ) * * SEE: * * DONGARRA J. J., DU CROZ J. J., DUFF I. AND HAMMARLING S. * A SET OF LEVEL 3 BASIC LINEAR ALGEBRA SUBPROGRAMS. TECHNICAL * MEMORANDUM NO.88 (REVISION 1), MATHEMATICS AND COMPUTER SCIENCE * DIVISION, ARGONNE NATIONAL LABORATORY, 9700 SOUTH CASS AVENUE, * ARGONNE, ILLINOIS 60439. * * ************************************************************************ * SUBROUTINE DGEMM ( TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, $ BETA, C, LDC ) * .. SCALAR ARGUMENTS .. CHARACTER*1 TRANSA, TRANSB INTEGER M, N, K, LDA, LDB, LDC DOUBLE PRECISION ALPHA, BETA * .. ARRAY ARGUMENTS .. DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ) * .. * * PURPOSE * ======= * * DGEMM PERFORMS ONE OF THE MATRIX-MATRIX OPERATIONS * * C := ALPHA*OP( A )*OP( B ) + BETA*C, * * WHERE OP( X ) IS ONE OF * * OP( X ) = X OR OP( X ) = X', * * ALPHA AND BETA ARE SCALARS, AND A, B AND C ARE MATRICES, WITH OP( A ) * AN M BY K MATRIX, OP( B ) A K BY N MATRIX AND C AN M BY N MATRIX. * * PARAMETERS * ========== * * TRANSA - CHARACTER*1. * ON ENTRY, TRANSA SPECIFIES THE FORM OF OP( A ) TO BE USED IN * THE MATRIX MULTIPLICATION AS FOLLOWS: * * TRANSA = 'N' OR 'N', OP( A ) = A. * * TRANSA = 'T' OR 'T', OP( A ) = A'. * * TRANSA = 'C' OR 'C', OP( A ) = A'. * * UNCHANGED ON EXIT. * * TRANSB - CHARACTER*1. * ON ENTRY, TRANSB SPECIFIES THE FORM OF OP( B ) TO BE USED IN * THE MATRIX MULTIPLICATION AS FOLLOWS: * * TRANSB = 'N' OR 'N', OP( B ) = B. * * TRANSB = 'T' OR 'T', OP( B ) = B'. * * TRANSB = 'C' OR 'C', OP( B ) = B'. * * UNCHANGED ON EXIT. * * M - INTEGER. * ON ENTRY, M SPECIFIES THE NUMBER OF ROWS OF THE MATRIX * OP( A ) AND OF THE MATRIX C. M MUST BE AT LEAST ZERO. * UNCHANGED ON EXIT. * * N - INTEGER. * ON ENTRY, N SPECIFIES THE NUMBER OF COLUMNS OF THE MATRIX * OP( B ) AND THE NUMBER OF COLUMNS OF THE MATRIX C. N MUST BE * AT LEAST ZERO. * UNCHANGED ON EXIT. * * K - INTEGER. * ON ENTRY, K SPECIFIES THE NUMBER OF COLUMNS OF THE MATRIX * OP( A ) AND THE NUMBER OF ROWS OF THE MATRIX OP( B ). K 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, KA ), WHERE KA IS * K WHEN TRANSA = 'N' OR 'N', AND IS M OTHERWISE. * BEFORE ENTRY WITH TRANSA = 'N' OR 'N', THE LEADING M BY K * PART OF THE ARRAY A MUST CONTAIN THE MATRIX A, OTHERWISE * THE LEADING K BY M PART OF THE ARRAY A MUST CONTAIN THE * MATRIX A. * UNCHANGED ON EXIT. * * LDA - INTEGER. * ON ENTRY, LDA SPECIFIES THE FIRST DIMENSION OF A AS DECLARED * IN THE CALLING (SUB) PROGRAM. WHEN TRANSA = 'N' OR 'N' THEN * LDA MUST BE AT LEAST MAX( 1, M ), OTHERWISE LDA MUST BE AT * LEAST MAX( 1, K ). * UNCHANGED ON EXIT. * * B - DOUBLE PRECISION ARRAY OF DIMENSION ( LDB, KB ), WHERE KB IS * N WHEN TRANSB = 'N' OR 'N', AND IS K OTHERWISE. * BEFORE ENTRY WITH TRANSB = 'N' OR 'N', THE LEADING K BY N * PART OF THE ARRAY B MUST CONTAIN THE MATRIX B, OTHERWISE * THE LEADING N BY K PART OF THE ARRAY B MUST CONTAIN THE * MATRIX B. * UNCHANGED ON EXIT. * * LDB - INTEGER. * ON ENTRY, LDB SPECIFIES THE FIRST DIMENSION OF B AS DECLARED * IN THE CALLING (SUB) PROGRAM. WHEN TRANSB = 'N' OR 'N' THEN * LDB MUST BE AT LEAST MAX( 1, K ), OTHERWISE LDB MUST BE AT * LEAST MAX( 1, N ). * UNCHANGED ON EXIT. * * BETA - DOUBLE PRECISION. * ON ENTRY, BETA SPECIFIES THE SCALAR BETA. WHEN BETA IS * SUPPLIED AS ZERO THEN C NEED NOT BE SET ON INPUT. * UNCHANGED ON EXIT. * * C - DOUBLE PRECISION ARRAY OF DIMENSION ( LDC, N ). * BEFORE ENTRY, THE LEADING M BY N PART OF THE ARRAY C MUST * CONTAIN THE MATRIX C, EXCEPT WHEN BETA IS ZERO, IN WHICH * CASE C NEED NOT BE SET ON ENTRY. * ON EXIT, THE ARRAY C IS OVERWRITTEN BY THE M BY N MATRIX * ( ALPHA*OP( A )*OP( B ) + BETA*C ). * * LDC - INTEGER. * ON ENTRY, LDC SPECIFIES THE FIRST DIMENSION OF C AS DECLARED * IN THE CALLING (SUB) PROGRAM. LDC MUST BE AT LEAST * MAX( 1, M ). * UNCHANGED ON EXIT. * * * LEVEL 3 BLAS ROUTINE. * * -- WRITTEN ON 8-FEBRUARY-1989. * JACK DONGARRA, ARGONNE NATIONAL LABORATORY. * IAIN DUFF, AERE HARWELL. * JEREMY DU CROZ, NUMERICAL ALGORITHMS GROUP LTD. * SVEN HAMMARLING, NUMERICAL ALGORITHMS GROUP LTD. * * * .. EXTERNAL FUNCTIONS .. LOGICAL LSAME EXTERNAL LSAME * .. EXTERNAL SUBROUTINES .. EXTERNAL XERBLA * .. INTRINSIC FUNCTIONS .. INTRINSIC MAX * .. LOCAL SCALARS .. LOGICAL NOTA, NOTB INTEGER I, INFO, J, L, NCOLA, NROWA, NROWB DOUBLE PRECISION TEMP * .. PARAMETERS .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. EXECUTABLE STATEMENTS .. * * SET NOTA AND NOTB AS TRUE IF A AND B RESPECTIVELY ARE NOT * TRANSPOSED AND SET NROWA, NCOLA AND NROWB AS THE NUMBER OF ROWS * AND COLUMNS OF A AND THE NUMBER OF ROWS OF B RESPECTIVELY. * NOTA = LSAME( TRANSA, 'N' ) NOTB = LSAME( TRANSB, 'N' ) IF( NOTA )THEN NROWA = M NCOLA = K ELSE NROWA = K NCOLA = M END IF IF( NOTB )THEN NROWB = K ELSE NROWB = N END IF * * TEST THE INPUT PARAMETERS. * INFO = 0 IF( ( .NOT.NOTA ).AND. $ ( .NOT.LSAME( TRANSA, 'C' ) ).AND. $ ( .NOT.LSAME( TRANSA, 'T' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.NOTB ).AND. $ ( .NOT.LSAME( TRANSB, 'C' ) ).AND. $ ( .NOT.LSAME( TRANSB, 'T' ) ) )THEN INFO = 2 ELSE IF( M .LT.0 )THEN INFO = 3 ELSE IF( N .LT.0 )THEN INFO = 4 ELSE IF( K .LT.0 )THEN INFO = 5 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 8 ELSE IF( LDB.LT.MAX( 1, NROWB ) )THEN INFO = 10 ELSE IF( LDC.LT.MAX( 1, M ) )THEN INFO = 13 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DGEMM ', INFO ) RETURN END IF * * QUICK RETURN IF POSSIBLE. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. $ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * AND IF ALPHA.EQ.ZERO. * IF( ALPHA.EQ.ZERO )THEN IF( BETA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, M C( I, J ) = ZERO 10 CONTINUE 20 CONTINUE ELSE DO 40, J = 1, N DO 30, I = 1, M C( I, J ) = BETA*C( I, J ) 30 CONTINUE 40 CONTINUE END IF RETURN END IF * * START THE OPERATIONS. * IF( NOTB )THEN IF( NOTA )THEN * * FORM C := ALPHA*A*B + BETA*C. * DO 90, J = 1, N IF( BETA.EQ.ZERO )THEN DO 50, I = 1, M C( I, J ) = ZERO 50 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 60, I = 1, M C( I, J ) = BETA*C( I, J ) 60 CONTINUE END IF DO 80, L = 1, K IF( B( L, J ).NE.ZERO )THEN TEMP = ALPHA*B( L, J ) DO 70, I = 1, M C( I, J ) = C( I, J ) + TEMP*A( I, L ) 70 CONTINUE END IF 80 CONTINUE 90 CONTINUE ELSE * * FORM C := ALPHA*A'*B + BETA*C * DO 120, J = 1, N DO 110, I = 1, M TEMP = ZERO DO 100, L = 1, K TEMP = TEMP + A( L, I )*B( L, J ) 100 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 110 CONTINUE 120 CONTINUE END IF ELSE IF( NOTA )THEN * * FORM C := ALPHA*A*B' + BETA*C * DO 170, J = 1, N IF( BETA.EQ.ZERO )THEN DO 130, I = 1, M C( I, J ) = ZERO 130 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 140, I = 1, M C( I, J ) = BETA*C( I, J ) 140 CONTINUE END IF DO 160, L = 1, K IF( B( J, L ).NE.ZERO )THEN TEMP = ALPHA*B( J, L ) DO 150, I = 1, M C( I, J ) = C( I, J ) + TEMP*A( I, L ) 150 CONTINUE END IF 160 CONTINUE 170 CONTINUE ELSE * * FORM C := ALPHA*A'*B' + BETA*C * DO 200, J = 1, N DO 190, I = 1, M TEMP = ZERO DO 180, L = 1, K TEMP = TEMP + A( L, I )*B( J, L ) 180 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 190 CONTINUE 200 CONTINUE END IF END IF * RETURN * * END OF DGEMM . * END CUT HERE............ CAT > DSYR2K.F <<'CUT HERE............' * ************************************************************************ * SUBROUTINE DSYR2K( UPLO, TRANS, N, K, ALPHA, A, LDA, B, LDB, $ BETA, C, LDC ) * .. SCALAR ARGUMENTS .. CHARACTER*1 UPLO, TRANS INTEGER N, K, LDA, LDB, LDC DOUBLE PRECISION ALPHA, BETA * .. ARRAY ARGUMENTS .. DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ) * .. * * PURPOSE * ======= * * DSYR2K PERFORMS ONE OF THE SYMMETRIC RANK 2K OPERATIONS * * C := ALPHA*A*B' + ALPHA*B*A' + BETA*C, * * OR * * C := ALPHA*A'*B + ALPHA*B'*A + BETA*C, * * WHERE ALPHA AND BETA ARE SCALARS, C IS AN N BY N SYMMETRIC MATRIX * AND A AND B ARE N BY K MATRICES IN THE FIRST CASE AND K BY N * MATRICES IN THE SECOND CASE. * * PARAMETERS * ========== * * UPLO - CHARACTER*1. * ON ENTRY, UPLO SPECIFIES WHETHER THE UPPER OR LOWER * TRIANGULAR PART OF THE ARRAY C IS TO BE REFERENCED AS * FOLLOWS: * * UPLO = 'U' OR 'U' ONLY THE UPPER TRIANGULAR PART OF C * IS TO BE REFERENCED. * * UPLO = 'L' OR 'L' ONLY THE LOWER TRIANGULAR PART OF C * IS TO BE REFERENCED. * * UNCHANGED ON EXIT. * * TRANS - CHARACTER*1. * ON ENTRY, TRANS SPECIFIES THE OPERATION TO BE PERFORMED AS * FOLLOWS: * * TRANS = 'N' OR 'N' C := ALPHA*A*B' + ALPHA*B*A' + * BETA*C. * * TRANS = 'T' OR 'T' C := ALPHA*A'*B + ALPHA*B'*A + * BETA*C. * * TRANS = 'C' OR 'C' C := ALPHA*A'*B + ALPHA*B'*A + * BETA*C. * * UNCHANGED ON EXIT. * * N - INTEGER. * ON ENTRY, N SPECIFIES THE ORDER OF THE MATRIX C. N MUST BE * AT LEAST ZERO. * UNCHANGED ON EXIT. * * K - INTEGER. * ON ENTRY WITH TRANS = 'N' OR 'N', K SPECIFIES THE NUMBER * OF COLUMNS OF THE MATRICES A AND B, AND ON ENTRY WITH * TRANS = 'T' OR 'T' OR 'C' OR 'C', K SPECIFIES THE NUMBER * OF ROWS OF THE MATRICES A AND B. K 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, KA ), WHERE KA IS * K WHEN TRANS = 'N' OR 'N', AND IS N OTHERWISE. * BEFORE ENTRY WITH TRANS = 'N' OR 'N', THE LEADING N BY K * PART OF THE ARRAY A MUST CONTAIN THE MATRIX A, OTHERWISE * THE LEADING K BY N PART OF THE ARRAY A MUST CONTAIN THE * MATRIX A. * UNCHANGED ON EXIT. * * LDA - INTEGER. * ON ENTRY, LDA SPECIFIES THE FIRST DIMENSION OF A AS DECLARED * IN THE CALLING (SUB) PROGRAM. WHEN TRANS = 'N' OR 'N' * THEN LDA MUST BE AT LEAST MAX( 1, N ), OTHERWISE LDA MUST * BE AT LEAST MAX( 1, K ). * UNCHANGED ON EXIT. * * B - DOUBLE PRECISION ARRAY OF DIMENSION ( LDB, KB ), WHERE KB IS * K WHEN TRANS = 'N' OR 'N', AND IS N OTHERWISE. * BEFORE ENTRY WITH TRANS = 'N' OR 'N', THE LEADING N BY K * PART OF THE ARRAY B MUST CONTAIN THE MATRIX B, OTHERWISE * THE LEADING K BY N PART OF THE ARRAY B MUST CONTAIN THE * MATRIX B. * UNCHANGED ON EXIT. * * LDB - INTEGER. * ON ENTRY, LDB SPECIFIES THE FIRST DIMENSION OF B AS DECLARED * IN THE CALLING (SUB) PROGRAM. WHEN TRANS = 'N' OR 'N' * THEN LDB MUST BE AT LEAST MAX( 1, N ), OTHERWISE LDB MUST * BE AT LEAST MAX( 1, K ). * UNCHANGED ON EXIT. * * BETA - DOUBLE PRECISION. * ON ENTRY, BETA SPECIFIES THE SCALAR BETA. * UNCHANGED ON EXIT. * * C - DOUBLE PRECISION ARRAY OF DIMENSION ( LDC, N ). * BEFORE ENTRY WITH UPLO = 'U' OR 'U', THE LEADING N BY N * UPPER TRIANGULAR PART OF THE ARRAY C MUST CONTAIN THE UPPER * TRIANGULAR PART OF THE SYMMETRIC MATRIX AND THE STRICTLY * LOWER TRIANGULAR PART OF C IS NOT REFERENCED. ON EXIT, THE * UPPER TRIANGULAR PART OF THE ARRAY C 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 C MUST CONTAIN THE LOWER * TRIANGULAR PART OF THE SYMMETRIC MATRIX AND THE STRICTLY * UPPER TRIANGULAR PART OF C IS NOT REFERENCED. ON EXIT, THE * LOWER TRIANGULAR PART OF THE ARRAY C IS OVERWRITTEN BY THE * LOWER TRIANGULAR PART OF THE UPDATED MATRIX. * * LDC - INTEGER. * ON ENTRY, LDC SPECIFIES THE FIRST DIMENSION OF C AS DECLARED * IN THE CALLING (SUB) PROGRAM. LDC MUST BE AT LEAST * MAX( 1, N ). * UNCHANGED ON EXIT. * * * LEVEL 3 BLAS ROUTINE. * * * -- WRITTEN ON 8-FEBRUARY-1989. * JACK DONGARRA, ARGONNE NATIONAL LABORATORY. * IAIN DUFF, AERE HARWELL. * JEREMY DU CROZ, NUMERICAL ALGORITHMS GROUP LTD. * SVEN HAMMARLING, NUMERICAL ALGORITHMS GROUP LTD. * * * .. EXTERNAL FUNCTIONS .. LOGICAL LSAME EXTERNAL LSAME * .. EXTERNAL SUBROUTINES .. EXTERNAL XERBLA * .. INTRINSIC FUNCTIONS .. INTRINSIC MAX * .. LOCAL SCALARS .. LOGICAL UPPER INTEGER I, INFO, J, L, NROWA DOUBLE PRECISION TEMP1, TEMP2 * .. PARAMETERS .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. EXECUTABLE STATEMENTS .. * * TEST THE INPUT PARAMETERS. * IF( LSAME( TRANS, 'N' ) )THEN NROWA = N ELSE NROWA = K END IF UPPER = LSAME( UPLO, 'U' ) * INFO = 0 IF( ( .NOT.UPPER ).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( N .LT.0 )THEN INFO = 3 ELSE IF( K .LT.0 )THEN INFO = 4 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 7 ELSE IF( LDB.LT.MAX( 1, NROWA ) )THEN INFO = 9 ELSE IF( LDC.LT.MAX( 1, N ) )THEN INFO = 12 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DSYR2K', INFO ) RETURN END IF * * QUICK RETURN IF POSSIBLE. * IF( ( N.EQ.0 ).OR. $ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * AND WHEN ALPHA.EQ.ZERO. * IF( ALPHA.EQ.ZERO )THEN IF( UPPER )THEN IF( BETA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, J C( I, J ) = ZERO 10 CONTINUE 20 CONTINUE ELSE DO 40, J = 1, N DO 30, I = 1, J C( I, J ) = BETA*C( I, J ) 30 CONTINUE 40 CONTINUE END IF ELSE IF( BETA.EQ.ZERO )THEN DO 60, J = 1, N DO 50, I = J, N C( I, J ) = ZERO 50 CONTINUE 60 CONTINUE ELSE DO 80, J = 1, N DO 70, I = J, N C( I, J ) = BETA*C( I, J ) 70 CONTINUE 80 CONTINUE END IF END IF RETURN END IF * * START THE OPERATIONS. * IF( LSAME( TRANS, 'N' ) )THEN * * FORM C := ALPHA*A*B' + ALPHA*B*A' + C. * IF( UPPER )THEN DO 130, J = 1, N IF( BETA.EQ.ZERO )THEN DO 90, I = 1, J C( I, J ) = ZERO 90 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 100, I = 1, J C( I, J ) = BETA*C( I, J ) 100 CONTINUE END IF DO 120, L = 1, K IF( ( A( J, L ).NE.ZERO ).OR. $ ( B( J, L ).NE.ZERO ) )THEN TEMP1 = ALPHA*B( J, L ) TEMP2 = ALPHA*A( J, L ) DO 110, I = 1, J C( I, J ) = C( I, J ) + $ A( I, L )*TEMP1 + B( I, L )*TEMP2 110 CONTINUE END IF 120 CONTINUE 130 CONTINUE ELSE DO 180, J = 1, N IF( BETA.EQ.ZERO )THEN DO 140, I = J, N C( I, J ) = ZERO 140 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 150, I = J, N C( I, J ) = BETA*C( I, J ) 150 CONTINUE END IF DO 170, L = 1, K IF( ( A( J, L ).NE.ZERO ).OR. $ ( B( J, L ).NE.ZERO ) )THEN TEMP1 = ALPHA*B( J, L ) TEMP2 = ALPHA*A( J, L ) DO 160, I = J, N C( I, J ) = C( I, J ) + $ A( I, L )*TEMP1 + B( I, L )*TEMP2 160 CONTINUE END IF 170 CONTINUE 180 CONTINUE END IF ELSE * * FORM C := ALPHA*A'*B + ALPHA*B'*A + C. * IF( UPPER )THEN DO 210, J = 1, N DO 200, I = 1, J TEMP1 = ZERO TEMP2 = ZERO DO 190, L = 1, K TEMP1 = TEMP1 + A( L, I )*B( L, J ) TEMP2 = TEMP2 + B( L, I )*A( L, J ) 190 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP1 + ALPHA*TEMP2 ELSE C( I, J ) = BETA *C( I, J ) + $ ALPHA*TEMP1 + ALPHA*TEMP2 END IF 200 CONTINUE 210 CONTINUE ELSE DO 240, J = 1, N DO 230, I = J, N TEMP1 = ZERO TEMP2 = ZERO DO 220, L = 1, K TEMP1 = TEMP1 + A( L, I )*B( L, J ) TEMP2 = TEMP2 + B( L, I )*A( L, J ) 220 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP1 + ALPHA*TEMP2 ELSE C( I, J ) = BETA *C( I, J ) + $ ALPHA*TEMP1 + ALPHA*TEMP2 END IF 230 CONTINUE 240 CONTINUE END IF END IF * RETURN * * END OF DSYR2K. * END CUT HERE............ CAT > DSYRK.F <<'CUT HERE............' * ************************************************************************ * SUBROUTINE DSYRK ( UPLO, TRANS, N, K, ALPHA, A, LDA, $ BETA, C, LDC ) * .. SCALAR ARGUMENTS .. CHARACTER*1 UPLO, TRANS INTEGER N, K, LDA, LDC DOUBLE PRECISION ALPHA, BETA * .. ARRAY ARGUMENTS .. DOUBLE PRECISION A( LDA, * ), C( LDC, * ) * .. * * PURPOSE * ======= * * DSYRK PERFORMS ONE OF THE SYMMETRIC RANK K OPERATIONS * * C := ALPHA*A*A' + BETA*C, * * OR * * C := ALPHA*A'*A + BETA*C, * * WHERE ALPHA AND BETA ARE SCALARS, C IS AN N BY N SYMMETRIC MATRIX * AND A IS AN N BY K MATRIX IN THE FIRST CASE AND A K BY N MATRIX * IN THE SECOND CASE. * * PARAMETERS * ========== * * UPLO - CHARACTER*1. * ON ENTRY, UPLO SPECIFIES WHETHER THE UPPER OR LOWER * TRIANGULAR PART OF THE ARRAY C IS TO BE REFERENCED AS * FOLLOWS: * * UPLO = 'U' OR 'U' ONLY THE UPPER TRIANGULAR PART OF C * IS TO BE REFERENCED. * * UPLO = 'L' OR 'L' ONLY THE LOWER TRIANGULAR PART OF C * IS TO BE REFERENCED. * * UNCHANGED ON EXIT. * * TRANS - CHARACTER*1. * ON ENTRY, TRANS SPECIFIES THE OPERATION TO BE PERFORMED AS * FOLLOWS: * * TRANS = 'N' OR 'N' C := ALPHA*A*A' + BETA*C. * * TRANS = 'T' OR 'T' C := ALPHA*A'*A + BETA*C. * * TRANS = 'C' OR 'C' C := ALPHA*A'*A + BETA*C. * * UNCHANGED ON EXIT. * * N - INTEGER. * ON ENTRY, N SPECIFIES THE ORDER OF THE MATRIX C. N MUST BE * AT LEAST ZERO. * UNCHANGED ON EXIT. * * K - INTEGER. * ON ENTRY WITH TRANS = 'N' OR 'N', K SPECIFIES THE NUMBER * OF COLUMNS OF THE MATRIX A, AND ON ENTRY WITH * TRANS = 'T' OR 'T' OR 'C' OR 'C', K SPECIFIES THE NUMBER * OF ROWS OF THE MATRIX A. K 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, KA ), WHERE KA IS * K WHEN TRANS = 'N' OR 'N', AND IS N OTHERWISE. * BEFORE ENTRY WITH TRANS = 'N' OR 'N', THE LEADING N BY K * PART OF THE ARRAY A MUST CONTAIN THE MATRIX A, OTHERWISE * THE LEADING K BY N PART OF THE ARRAY A MUST CONTAIN THE * MATRIX A. * UNCHANGED ON EXIT. * * LDA - INTEGER. * ON ENTRY, LDA SPECIFIES THE FIRST DIMENSION OF A AS DECLARED * IN THE CALLING (SUB) PROGRAM. WHEN TRANS = 'N' OR 'N' * THEN LDA MUST BE AT LEAST MAX( 1, N ), OTHERWISE LDA MUST * BE AT LEAST MAX( 1, K ). * UNCHANGED ON EXIT. * * BETA - DOUBLE PRECISION. * ON ENTRY, BETA SPECIFIES THE SCALAR BETA. * UNCHANGED ON EXIT. * * C - DOUBLE PRECISION ARRAY OF DIMENSION ( LDC, N ). * BEFORE ENTRY WITH UPLO = 'U' OR 'U', THE LEADING N BY N * UPPER TRIANGULAR PART OF THE ARRAY C MUST CONTAIN THE UPPER * TRIANGULAR PART OF THE SYMMETRIC MATRIX AND THE STRICTLY * LOWER TRIANGULAR PART OF C IS NOT REFERENCED. ON EXIT, THE * UPPER TRIANGULAR PART OF THE ARRAY C 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 C MUST CONTAIN THE LOWER * TRIANGULAR PART OF THE SYMMETRIC MATRIX AND THE STRICTLY * UPPER TRIANGULAR PART OF C IS NOT REFERENCED. ON EXIT, THE * LOWER TRIANGULAR PART OF THE ARRAY C IS OVERWRITTEN BY THE * LOWER TRIANGULAR PART OF THE UPDATED MATRIX. * * LDC - INTEGER. * ON ENTRY, LDC SPECIFIES THE FIRST DIMENSION OF C AS DECLARED * IN THE CALLING (SUB) PROGRAM. LDC MUST BE AT LEAST * MAX( 1, N ). * UNCHANGED ON EXIT. * * * LEVEL 3 BLAS ROUTINE. * * -- WRITTEN ON 8-FEBRUARY-1989. * JACK DONGARRA, ARGONNE NATIONAL LABORATORY. * IAIN DUFF, AERE HARWELL. * JEREMY DU CROZ, NUMERICAL ALGORITHMS GROUP LTD. * SVEN HAMMARLING, NUMERICAL ALGORITHMS GROUP LTD. * * * .. EXTERNAL FUNCTIONS .. LOGICAL LSAME EXTERNAL LSAME * .. EXTERNAL SUBROUTINES .. EXTERNAL XERBLA * .. INTRINSIC FUNCTIONS .. INTRINSIC MAX * .. LOCAL SCALARS .. LOGICAL UPPER INTEGER I, INFO, J, L, NROWA DOUBLE PRECISION TEMP * .. PARAMETERS .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. EXECUTABLE STATEMENTS .. * * TEST THE INPUT PARAMETERS. * IF( LSAME( TRANS, 'N' ) )THEN NROWA = N ELSE NROWA = K END IF UPPER = LSAME( UPLO, 'U' ) * INFO = 0 IF( ( .NOT.UPPER ).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( N .LT.0 )THEN INFO = 3 ELSE IF( K .LT.0 )THEN INFO = 4 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 7 ELSE IF( LDC.LT.MAX( 1, N ) )THEN INFO = 10 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DSYRK ', INFO ) RETURN END IF * * QUICK RETURN IF POSSIBLE. * IF( ( N.EQ.0 ).OR. $ ( ( ( ALPHA.EQ.ZERO ).OR.( K.EQ.0 ) ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * AND WHEN ALPHA.EQ.ZERO. * IF( ALPHA.EQ.ZERO )THEN IF( UPPER )THEN IF( BETA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, J C( I, J ) = ZERO 10 CONTINUE 20 CONTINUE ELSE DO 40, J = 1, N DO 30, I = 1, J C( I, J ) = BETA*C( I, J ) 30 CONTINUE 40 CONTINUE END IF ELSE IF( BETA.EQ.ZERO )THEN DO 60, J = 1, N DO 50, I = J, N C( I, J ) = ZERO 50 CONTINUE 60 CONTINUE ELSE DO 80, J = 1, N DO 70, I = J, N C( I, J ) = BETA*C( I, J ) 70 CONTINUE 80 CONTINUE END IF END IF RETURN END IF * * START THE OPERATIONS. * IF( LSAME( TRANS, 'N' ) )THEN * * FORM C := ALPHA*A*A' + BETA*C. * IF( UPPER )THEN DO 130, J = 1, N IF( BETA.EQ.ZERO )THEN DO 90, I = 1, J C( I, J ) = ZERO 90 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 100, I = 1, J C( I, J ) = BETA*C( I, J ) 100 CONTINUE END IF DO 120, L = 1, K IF( A( J, L ).NE.ZERO )THEN TEMP = ALPHA*A( J, L ) DO 110, I = 1, J C( I, J ) = C( I, J ) + TEMP*A( I, L ) 110 CONTINUE END IF 120 CONTINUE 130 CONTINUE ELSE DO 180, J = 1, N IF( BETA.EQ.ZERO )THEN DO 140, I = J, N C( I, J ) = ZERO 140 CONTINUE ELSE IF( BETA.NE.ONE )THEN DO 150, I = J, N C( I, J ) = BETA*C( I, J ) 150 CONTINUE END IF DO 170, L = 1, K IF( A( J, L ).NE.ZERO )THEN TEMP = ALPHA*A( J, L ) DO 160, I = J, N C( I, J ) = C( I, J ) + TEMP*A( I, L ) 160 CONTINUE END IF 170 CONTINUE 180 CONTINUE END IF ELSE * * FORM C := ALPHA*A'*A + BETA*C. * IF( UPPER )THEN DO 210, J = 1, N DO 200, I = 1, J TEMP = ZERO DO 190, L = 1, K TEMP = TEMP + A( L, I )*A( L, J ) 190 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 200 CONTINUE 210 CONTINUE ELSE DO 240, J = 1, N DO 230, I = J, N TEMP = ZERO DO 220, L = 1, K TEMP = TEMP + A( L, I )*A( L, J ) 220 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = ALPHA*TEMP ELSE C( I, J ) = ALPHA*TEMP + BETA*C( I, J ) END IF 230 CONTINUE 240 CONTINUE END IF END IF * RETURN * * END OF DSYRK . * END CUT HERE............ CAT > DSYMM.F <<'CUT HERE............' * ************************************************************************ * SUBROUTINE DSYMM ( SIDE, UPLO, M, N, ALPHA, A, LDA, B, LDB, $ BETA, C, LDC ) * .. SCALAR ARGUMENTS .. CHARACTER*1 SIDE, UPLO INTEGER M, N, LDA, LDB, LDC DOUBLE PRECISION ALPHA, BETA * .. ARRAY ARGUMENTS .. DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ) * .. * * PURPOSE * ======= * * DSYMM PERFORMS ONE OF THE MATRIX-MATRIX OPERATIONS * * C := ALPHA*A*B + BETA*C, * * OR * * C := ALPHA*B*A + BETA*C, * * WHERE ALPHA AND BETA ARE SCALARS, A IS A SYMMETRIC MATRIX AND B AND * C ARE M BY N MATRICES. * * PARAMETERS * ========== * * SIDE - CHARACTER*1. * ON ENTRY, SIDE SPECIFIES WHETHER THE SYMMETRIC MATRIX A * APPEARS ON THE LEFT OR RIGHT IN THE OPERATION AS FOLLOWS: * * SIDE = 'L' OR 'L' C := ALPHA*A*B + BETA*C, * * SIDE = 'R' OR 'R' C := ALPHA*B*A + BETA*C, * * UNCHANGED ON EXIT. * * UPLO - CHARACTER*1. * ON ENTRY, UPLO SPECIFIES WHETHER THE UPPER OR LOWER * TRIANGULAR PART OF THE SYMMETRIC MATRIX A IS TO BE * REFERENCED AS FOLLOWS: * * UPLO = 'U' OR 'U' ONLY THE UPPER TRIANGULAR PART OF THE * SYMMETRIC MATRIX IS TO BE REFERENCED. * * UPLO = 'L' OR 'L' ONLY THE LOWER TRIANGULAR PART OF THE * SYMMETRIC MATRIX IS TO BE REFERENCED. * * UNCHANGED ON EXIT. * * M - INTEGER. * ON ENTRY, M SPECIFIES THE NUMBER OF ROWS OF THE MATRIX C. * M MUST BE AT LEAST ZERO. * UNCHANGED ON EXIT. * * N - INTEGER. * ON ENTRY, N SPECIFIES THE NUMBER OF COLUMNS OF THE MATRIX C. * 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, KA ), WHERE KA IS * M WHEN SIDE = 'L' OR 'L' AND IS N OTHERWISE. * BEFORE ENTRY WITH SIDE = 'L' OR 'L', THE M BY M PART OF * THE ARRAY A MUST CONTAIN THE SYMMETRIC MATRIX, SUCH THAT * WHEN UPLO = 'U' OR 'U', THE LEADING M BY M 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, AND WHEN UPLO = 'L' OR 'L', * THE LEADING M BY M 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. * BEFORE ENTRY WITH SIDE = 'R' OR 'R', THE N BY N PART OF * THE ARRAY A MUST CONTAIN THE SYMMETRIC MATRIX, SUCH THAT * WHEN 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, AND WHEN 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. WHEN SIDE = 'L' OR 'L' THEN * LDA MUST BE AT LEAST MAX( 1, M ), OTHERWISE LDA MUST BE AT * LEAST MAX( 1, N ). * UNCHANGED ON EXIT. * * B - DOUBLE PRECISION ARRAY OF DIMENSION ( LDB, N ). * BEFORE ENTRY, THE LEADING M BY N PART OF THE ARRAY B MUST * CONTAIN THE MATRIX B. * UNCHANGED ON EXIT. * * LDB - INTEGER. * ON ENTRY, LDB SPECIFIES THE FIRST DIMENSION OF B AS DECLARED * IN THE CALLING (SUB) PROGRAM. LDB MUST BE AT LEAST * MAX( 1, M ). * UNCHANGED ON EXIT. * * BETA - DOUBLE PRECISION. * ON ENTRY, BETA SPECIFIES THE SCALAR BETA. WHEN BETA IS * SUPPLIED AS ZERO THEN C NEED NOT BE SET ON INPUT. * UNCHANGED ON EXIT. * * C - DOUBLE PRECISION ARRAY OF DIMENSION ( LDC, N ). * BEFORE ENTRY, THE LEADING M BY N PART OF THE ARRAY C MUST * CONTAIN THE MATRIX C, EXCEPT WHEN BETA IS ZERO, IN WHICH * CASE C NEED NOT BE SET ON ENTRY. * ON EXIT, THE ARRAY C IS OVERWRITTEN BY THE M BY N UPDATED * MATRIX. * * LDC - INTEGER. * ON ENTRY, LDC SPECIFIES THE FIRST DIMENSION OF C AS DECLARED * IN THE CALLING (SUB) PROGRAM. LDC MUST BE AT LEAST * MAX( 1, M ). * UNCHANGED ON EXIT. * * * LEVEL 3 BLAS ROUTINE. * * -- WRITTEN ON 8-FEBRUARY-1989. * JACK DONGARRA, ARGONNE NATIONAL LABORATORY. * IAIN DUFF, AERE HARWELL. * JEREMY DU CROZ, NUMERICAL ALGORITHMS GROUP LTD. * SVEN HAMMARLING, NUMERICAL ALGORITHMS GROUP LTD. * * * .. EXTERNAL FUNCTIONS .. LOGICAL LSAME EXTERNAL LSAME * .. EXTERNAL SUBROUTINES .. EXTERNAL XERBLA * .. INTRINSIC FUNCTIONS .. INTRINSIC MAX * .. LOCAL SCALARS .. LOGICAL UPPER INTEGER I, INFO, J, K, NROWA DOUBLE PRECISION TEMP1, TEMP2 * .. PARAMETERS .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. EXECUTABLE STATEMENTS .. * * SET NROWA AS THE NUMBER OF ROWS OF A. * IF( LSAME( SIDE, 'L' ) )THEN NROWA = M ELSE NROWA = N END IF UPPER = LSAME( UPLO, 'U' ) * * TEST THE INPUT PARAMETERS. * INFO = 0 IF( ( .NOT.LSAME( SIDE, 'L' ) ).AND. $ ( .NOT.LSAME( SIDE, 'R' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.UPPER ).AND. $ ( .NOT.LSAME( UPLO, 'L' ) ) )THEN INFO = 2 ELSE IF( M .LT.0 )THEN INFO = 3 ELSE IF( N .LT.0 )THEN INFO = 4 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 7 ELSE IF( LDB.LT.MAX( 1, M ) )THEN INFO = 9 ELSE IF( LDC.LT.MAX( 1, M ) )THEN INFO = 12 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DSYMM ', INFO ) RETURN END IF * * QUICK RETURN IF POSSIBLE. * IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR. $ ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) ) $ RETURN * * AND WHEN ALPHA.EQ.ZERO. * IF( ALPHA.EQ.ZERO )THEN IF( BETA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, M C( I, J ) = ZERO 10 CONTINUE 20 CONTINUE ELSE DO 40, J = 1, N DO 30, I = 1, M C( I, J ) = BETA*C( I, J ) 30 CONTINUE 40 CONTINUE END IF RETURN END IF * * START THE OPERATIONS. * IF( LSAME( SIDE, 'L' ) )THEN * * FORM C := ALPHA*A*B + BETA*C. * IF( UPPER )THEN DO 70, J = 1, N DO 60, I = 1, M TEMP1 = ALPHA*B( I, J ) TEMP2 = ZERO DO 50, K = 1, I - 1 C( K, J ) = C( K, J ) + TEMP1 *A( K, I ) TEMP2 = TEMP2 + B( K, J )*A( K, I ) 50 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = TEMP1*A( I, I ) + ALPHA*TEMP2 ELSE C( I, J ) = BETA *C( I, J ) + $ TEMP1*A( I, I ) + ALPHA*TEMP2 END IF 60 CONTINUE 70 CONTINUE ELSE DO 100, J = 1, N DO 90, I = M, 1, -1 TEMP1 = ALPHA*B( I, J ) TEMP2 = ZERO DO 80, K = I + 1, M C( K, J ) = C( K, J ) + TEMP1 *A( K, I ) TEMP2 = TEMP2 + B( K, J )*A( K, I ) 80 CONTINUE IF( BETA.EQ.ZERO )THEN C( I, J ) = TEMP1*A( I, I ) + ALPHA*TEMP2 ELSE C( I, J ) = BETA *C( I, J ) + $ TEMP1*A( I, I ) + ALPHA*TEMP2 END IF 90 CONTINUE 100 CONTINUE END IF ELSE * * FORM C := ALPHA*B*A + BETA*C. * DO 170, J = 1, N TEMP1 = ALPHA*A( J, J ) IF( BETA.EQ.ZERO )THEN DO 110, I = 1, M C( I, J ) = TEMP1*B( I, J ) 110 CONTINUE ELSE DO 120, I = 1, M C( I, J ) = BETA*C( I, J ) + TEMP1*B( I, J ) 120 CONTINUE END IF DO 140, K = 1, J - 1 IF( UPPER )THEN TEMP1 = ALPHA*A( K, J ) ELSE TEMP1 = ALPHA*A( J, K ) END IF DO 130, I = 1, M C( I, J ) = C( I, J ) + TEMP1*B( I, K ) 130 CONTINUE 140 CONTINUE DO 160, K = J + 1, N IF( UPPER )THEN TEMP1 = ALPHA*A( J, K ) ELSE TEMP1 = ALPHA*A( K, J ) END IF DO 150, I = 1, M C( I, J ) = C( I, J ) + TEMP1*B( I, K ) 150 CONTINUE 160 CONTINUE 170 CONTINUE END IF * RETURN * * END OF DSYMM . * END CUT HERE............ CAT > DTRSM.F <<'CUT HERE............' * ************************************************************************ * SUBROUTINE DTRSM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, $ B, LDB ) * .. SCALAR ARGUMENTS .. CHARACTER*1 SIDE, UPLO, TRANSA, DIAG INTEGER M, N, LDA, LDB DOUBLE PRECISION ALPHA * .. ARRAY ARGUMENTS .. DOUBLE PRECISION A( LDA, * ), B( LDB, * ) * .. * * PURPOSE * ======= * * DTRSM SOLVES ONE OF THE MATRIX EQUATIONS * * OP( A )*X = ALPHA*B, OR X*OP( A ) = ALPHA*B, * * WHERE ALPHA IS A SCALAR, X AND B ARE M BY N MATRICES, A IS A UNIT, OR * NON-UNIT, UPPER OR LOWER TRIANGULAR MATRIX AND OP( A ) IS ONE OF * * OP( A ) = A OR OP( A ) = A'. * * THE MATRIX X IS OVERWRITTEN ON B. * * PARAMETERS * ========== * * SIDE - CHARACTER*1. * ON ENTRY, SIDE SPECIFIES WHETHER OP( A ) APPEARS ON THE LEFT * OR RIGHT OF X AS FOLLOWS: * * SIDE = 'L' OR 'L' OP( A )*X = ALPHA*B. * * SIDE = 'R' OR 'R' X*OP( A ) = ALPHA*B. * * UNCHANGED ON EXIT. * * UPLO - CHARACTER*1. * ON ENTRY, UPLO SPECIFIES WHETHER THE MATRIX A 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. * * TRANSA - CHARACTER*1. * ON ENTRY, TRANSA SPECIFIES THE FORM OF OP( A ) TO BE USED IN * THE MATRIX MULTIPLICATION AS FOLLOWS: * * TRANSA = 'N' OR 'N' OP( A ) = A. * * TRANSA = 'T' OR 'T' OP( A ) = A'. * * TRANSA = 'C' OR 'C' OP( A ) = A'. * * 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. * * M - INTEGER. * ON ENTRY, M SPECIFIES THE NUMBER OF ROWS OF B. M MUST BE AT * LEAST ZERO. * UNCHANGED ON EXIT. * * N - INTEGER. * ON ENTRY, N SPECIFIES THE NUMBER OF COLUMNS OF B. N MUST BE * AT LEAST ZERO. * UNCHANGED ON EXIT. * * ALPHA - DOUBLE PRECISION. * ON ENTRY, ALPHA SPECIFIES THE SCALAR ALPHA. WHEN ALPHA IS * ZERO THEN A IS NOT REFERENCED AND B NEED NOT BE SET BEFORE * ENTRY. * UNCHANGED ON EXIT. * * A - DOUBLE PRECISION ARRAY OF DIMENSION ( LDA, K ), WHERE K IS M * WHEN SIDE = 'L' OR 'L' AND IS N WHEN SIDE = 'R' OR 'R'. * BEFORE ENTRY WITH UPLO = 'U' OR 'U', THE LEADING K BY K * 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 K BY K * 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. WHEN SIDE = 'L' OR 'L' THEN * LDA MUST BE AT LEAST MAX( 1, M ), WHEN SIDE = 'R' OR 'R' * THEN LDA MUST BE AT LEAST MAX( 1, N ). * UNCHANGED ON EXIT. * * B - DOUBLE PRECISION ARRAY OF DIMENSION ( LDB, N ). * BEFORE ENTRY, THE LEADING M BY N PART OF THE ARRAY B MUST * CONTAIN THE RIGHT-HAND SIDE MATRIX B, AND ON EXIT IS * OVERWRITTEN BY THE SOLUTION MATRIX X. * * LDB - INTEGER. * ON ENTRY, LDB SPECIFIES THE FIRST DIMENSION OF B AS DECLARED * IN THE CALLING (SUB) PROGRAM. LDB MUST BE AT LEAST * MAX( 1, M ). * UNCHANGED ON EXIT. * * * LEVEL 3 BLAS ROUTINE. * * * -- WRITTEN ON 8-FEBRUARY-1989. * JACK DONGARRA, ARGONNE NATIONAL LABORATORY. * IAIN DUFF, AERE HARWELL. * JEREMY DU CROZ, NUMERICAL ALGORITHMS GROUP LTD. * SVEN HAMMARLING, NUMERICAL ALGORITHMS GROUP LTD. * * * .. EXTERNAL FUNCTIONS .. LOGICAL LSAME EXTERNAL LSAME * .. EXTERNAL SUBROUTINES .. EXTERNAL XERBLA * .. INTRINSIC FUNCTIONS .. INTRINSIC MAX * .. LOCAL SCALARS .. LOGICAL LSIDE, NOUNIT, UPPER INTEGER I, INFO, J, K, NROWA DOUBLE PRECISION TEMP * .. PARAMETERS .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. EXECUTABLE STATEMENTS .. * * TEST THE INPUT PARAMETERS. * LSIDE = LSAME( SIDE , 'L' ) IF( LSIDE )THEN NROWA = M ELSE NROWA = N END IF NOUNIT = LSAME( DIAG , 'N' ) UPPER = LSAME( UPLO , 'U' ) * INFO = 0 IF( ( .NOT.LSIDE ).AND. $ ( .NOT.LSAME( SIDE , 'R' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.UPPER ).AND. $ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN INFO = 2 ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND. $ ( .NOT.LSAME( TRANSA, 'T' ) ).AND. $ ( .NOT.LSAME( TRANSA, 'C' ) ) )THEN INFO = 3 ELSE IF( ( .NOT.LSAME( DIAG , 'U' ) ).AND. $ ( .NOT.LSAME( DIAG , 'N' ) ) )THEN INFO = 4 ELSE IF( M .LT.0 )THEN INFO = 5 ELSE IF( N .LT.0 )THEN INFO = 6 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 9 ELSE IF( LDB.LT.MAX( 1, M ) )THEN INFO = 11 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DTRSM ', INFO ) RETURN END IF * * QUICK RETURN IF POSSIBLE. * IF( N.EQ.0 ) $ RETURN * * AND WHEN ALPHA.EQ.ZERO. * IF( ALPHA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, M B( I, J ) = ZERO 10 CONTINUE 20 CONTINUE RETURN END IF * * START THE OPERATIONS. * IF( LSIDE )THEN IF( LSAME( TRANSA, 'N' ) )THEN * * FORM B := ALPHA*INV( A )*B. * IF( UPPER )THEN DO 60, J = 1, N IF( ALPHA.NE.ONE )THEN DO 30, I = 1, M B( I, J ) = ALPHA*B( I, J ) 30 CONTINUE END IF DO 50, K = M, 1, -1 IF( B( K, J ).NE.ZERO )THEN IF( NOUNIT ) $ B( K, J ) = B( K, J )/A( K, K ) DO 40, I = 1, K - 1 B( I, J ) = B( I, J ) - B( K, J )*A( I, K ) 40 CONTINUE END IF 50 CONTINUE 60 CONTINUE ELSE DO 100, J = 1, N IF( ALPHA.NE.ONE )THEN DO 70, I = 1, M B( I, J ) = ALPHA*B( I, J ) 70 CONTINUE END IF DO 90 K = 1, M IF( B( K, J ).NE.ZERO )THEN IF( NOUNIT ) $ B( K, J ) = B( K, J )/A( K, K ) DO 80, I = K + 1, M B( I, J ) = B( I, J ) - B( K, J )*A( I, K ) 80 CONTINUE END IF 90 CONTINUE 100 CONTINUE END IF ELSE * * FORM B := ALPHA*INV( A' )*B. * IF( UPPER )THEN DO 130, J = 1, N DO 120, I = 1, M TEMP = ALPHA*B( I, J ) DO 110, K = 1, I - 1 TEMP = TEMP - A( K, I )*B( K, J ) 110 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( I, I ) B( I, J ) = TEMP 120 CONTINUE 130 CONTINUE ELSE DO 160, J = 1, N DO 150, I = M, 1, -1 TEMP = ALPHA*B( I, J ) DO 140, K = I + 1, M TEMP = TEMP - A( K, I )*B( K, J ) 140 CONTINUE IF( NOUNIT ) $ TEMP = TEMP/A( I, I ) B( I, J ) = TEMP 150 CONTINUE 160 CONTINUE END IF END IF ELSE IF( LSAME( TRANSA, 'N' ) )THEN * * FORM B := ALPHA*B*INV( A ). * IF( UPPER )THEN DO 210, J = 1, N IF( ALPHA.NE.ONE )THEN DO 170, I = 1, M B( I, J ) = ALPHA*B( I, J ) 170 CONTINUE END IF DO 190, K = 1, J - 1 IF( A( K, J ).NE.ZERO )THEN DO 180, I = 1, M B( I, J ) = B( I, J ) - A( K, J )*B( I, K ) 180 CONTINUE END IF 190 CONTINUE IF( NOUNIT )THEN TEMP = ONE/A( J, J ) DO 200, I = 1, M B( I, J ) = TEMP*B( I, J ) 200 CONTINUE END IF 210 CONTINUE ELSE DO 260, J = N, 1, -1 IF( ALPHA.NE.ONE )THEN DO 220, I = 1, M B( I, J ) = ALPHA*B( I, J ) 220 CONTINUE END IF DO 240, K = J + 1, N IF( A( K, J ).NE.ZERO )THEN DO 230, I = 1, M B( I, J ) = B( I, J ) - A( K, J )*B( I, K ) 230 CONTINUE END IF 240 CONTINUE IF( NOUNIT )THEN TEMP = ONE/A( J, J ) DO 250, I = 1, M B( I, J ) = TEMP*B( I, J ) 250 CONTINUE END IF 260 CONTINUE END IF ELSE * * FORM B := ALPHA*B*INV( A' ). * IF( UPPER )THEN DO 310, K = N, 1, -1 IF( NOUNIT )THEN TEMP = ONE/A( K, K ) DO 270, I = 1, M B( I, K ) = TEMP*B( I, K ) 270 CONTINUE END IF DO 290, J = 1, K - 1 IF( A( J, K ).NE.ZERO )THEN TEMP = A( J, K ) DO 280, I = 1, M B( I, J ) = B( I, J ) - TEMP*B( I, K ) 280 CONTINUE END IF 290 CONTINUE IF( ALPHA.NE.ONE )THEN DO 300, I = 1, M B( I, K ) = ALPHA*B( I, K ) 300 CONTINUE END IF 310 CONTINUE ELSE DO 360, K = 1, N IF( NOUNIT )THEN TEMP = ONE/A( K, K ) DO 320, I = 1, M B( I, K ) = TEMP*B( I, K ) 320 CONTINUE END IF DO 340, J = K + 1, N IF( A( J, K ).NE.ZERO )THEN TEMP = A( J, K ) DO 330, I = 1, M B( I, J ) = B( I, J ) - TEMP*B( I, K ) 330 CONTINUE END IF 340 CONTINUE IF( ALPHA.NE.ONE )THEN DO 350, I = 1, M B( I, K ) = ALPHA*B( I, K ) 350 CONTINUE END IF 360 CONTINUE END IF END IF END IF * RETURN * * END OF DTRSM . * END CUT HERE............ CAT > DTRMM.F <<'CUT HERE............' * ************************************************************************ * SUBROUTINE DTRMM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, $ B, LDB ) * .. SCALAR ARGUMENTS .. CHARACTER*1 SIDE, UPLO, TRANSA, DIAG INTEGER M, N, LDA, LDB DOUBLE PRECISION ALPHA * .. ARRAY ARGUMENTS .. DOUBLE PRECISION A( LDA, * ), B( LDB, * ) * .. * * PURPOSE * ======= * * DTRMM PERFORMS ONE OF THE MATRIX-MATRIX OPERATIONS * * B := ALPHA*OP( A )*B, OR B := ALPHA*B*OP( A ), * * WHERE ALPHA IS A SCALAR, B IS AN M BY N MATRIX, A IS A UNIT, OR * NON-UNIT, UPPER OR LOWER TRIANGULAR MATRIX AND OP( A ) IS ONE OF * * OP( A ) = A OR OP( A ) = A'. * * PARAMETERS * ========== * * SIDE - CHARACTER*1. * ON ENTRY, SIDE SPECIFIES WHETHER OP( A ) MULTIPLIES B FROM * THE LEFT OR RIGHT AS FOLLOWS: * * SIDE = 'L' OR 'L' B := ALPHA*OP( A )*B. * * SIDE = 'R' OR 'R' B := ALPHA*B*OP( A ). * * UNCHANGED ON EXIT. * * UPLO - CHARACTER*1. * ON ENTRY, UPLO SPECIFIES WHETHER THE MATRIX A 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. * * TRANSA - CHARACTER*1. * ON ENTRY, TRANSA SPECIFIES THE FORM OF OP( A ) TO BE USED IN * THE MATRIX MULTIPLICATION AS FOLLOWS: * * TRANSA = 'N' OR 'N' OP( A ) = A. * * TRANSA = 'T' OR 'T' OP( A ) = A'. * * TRANSA = 'C' OR 'C' OP( A ) = A'. * * 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. * * M - INTEGER. * ON ENTRY, M SPECIFIES THE NUMBER OF ROWS OF B. M MUST BE AT * LEAST ZERO. * UNCHANGED ON EXIT. * * N - INTEGER. * ON ENTRY, N SPECIFIES THE NUMBER OF COLUMNS OF B. N MUST BE * AT LEAST ZERO. * UNCHANGED ON EXIT. * * ALPHA - DOUBLE PRECISION. * ON ENTRY, ALPHA SPECIFIES THE SCALAR ALPHA. WHEN ALPHA IS * ZERO THEN A IS NOT REFERENCED AND B NEED NOT BE SET BEFORE * ENTRY. * UNCHANGED ON EXIT. * * A - DOUBLE PRECISION ARRAY OF DIMENSION ( LDA, K ), WHERE K IS M * WHEN SIDE = 'L' OR 'L' AND IS N WHEN SIDE = 'R' OR 'R'. * BEFORE ENTRY WITH UPLO = 'U' OR 'U', THE LEADING K BY K * 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 K BY K * 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. WHEN SIDE = 'L' OR 'L' THEN * LDA MUST BE AT LEAST MAX( 1, M ), WHEN SIDE = 'R' OR 'R' * THEN LDA MUST BE AT LEAST MAX( 1, N ). * UNCHANGED ON EXIT. * * B - DOUBLE PRECISION ARRAY OF DIMENSION ( LDB, N ). * BEFORE ENTRY, THE LEADING M BY N PART OF THE ARRAY B MUST * CONTAIN THE MATRIX B, AND ON EXIT IS OVERWRITTEN BY THE * TRANSFORMED MATRIX. * * LDB - INTEGER. * ON ENTRY, LDB SPECIFIES THE FIRST DIMENSION OF B AS DECLARED * IN THE CALLING (SUB) PROGRAM. LDB MUST BE AT LEAST * MAX( 1, M ). * UNCHANGED ON EXIT. * * * LEVEL 3 BLAS ROUTINE. * * -- WRITTEN ON 8-FEBRUARY-1989. * JACK DONGARRA, ARGONNE NATIONAL LABORATORY. * IAIN DUFF, AERE HARWELL. * JEREMY DU CROZ, NUMERICAL ALGORITHMS GROUP LTD. * SVEN HAMMARLING, NUMERICAL ALGORITHMS GROUP LTD. * * * .. EXTERNAL FUNCTIONS .. LOGICAL LSAME EXTERNAL LSAME * .. EXTERNAL SUBROUTINES .. EXTERNAL XERBLA * .. INTRINSIC FUNCTIONS .. INTRINSIC MAX * .. LOCAL SCALARS .. LOGICAL LSIDE, NOUNIT, UPPER INTEGER I, INFO, J, K, NROWA DOUBLE PRECISION TEMP * .. PARAMETERS .. DOUBLE PRECISION ONE , ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. EXECUTABLE STATEMENTS .. * * TEST THE INPUT PARAMETERS. * LSIDE = LSAME( SIDE , 'L' ) IF( LSIDE )THEN NROWA = M ELSE NROWA = N END IF NOUNIT = LSAME( DIAG , 'N' ) UPPER = LSAME( UPLO , 'U' ) * INFO = 0 IF( ( .NOT.LSIDE ).AND. $ ( .NOT.LSAME( SIDE , 'R' ) ) )THEN INFO = 1 ELSE IF( ( .NOT.UPPER ).AND. $ ( .NOT.LSAME( UPLO , 'L' ) ) )THEN INFO = 2 ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND. $ ( .NOT.LSAME( TRANSA, 'T' ) ).AND. $ ( .NOT.LSAME( TRANSA, 'C' ) ) )THEN INFO = 3 ELSE IF( ( .NOT.LSAME( DIAG , 'U' ) ).AND. $ ( .NOT.LSAME( DIAG , 'N' ) ) )THEN INFO = 4 ELSE IF( M .LT.0 )THEN INFO = 5 ELSE IF( N .LT.0 )THEN INFO = 6 ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN INFO = 9 ELSE IF( LDB.LT.MAX( 1, M ) )THEN INFO = 11 END IF IF( INFO.NE.0 )THEN CALL XERBLA( 'DTRMM ', INFO ) RETURN END IF * * QUICK RETURN IF POSSIBLE. * IF( N.EQ.0 ) $ RETURN * * AND WHEN ALPHA.EQ.ZERO. * IF( ALPHA.EQ.ZERO )THEN DO 20, J = 1, N DO 10, I = 1, M B( I, J ) = ZERO 10 CONTINUE 20 CONTINUE RETURN END IF * * START THE OPERATIONS. * IF( LSIDE )THEN IF( LSAME( TRANSA, 'N' ) )THEN * * FORM B := ALPHA*A*B. * IF( UPPER )THEN DO 50, J = 1, N DO 40, K = 1, M IF( B( K, J ).NE.ZERO )THEN TEMP = ALPHA*B( K, J ) DO 30, I = 1, K - 1 B( I, J ) = B( I, J ) + TEMP*A( I, K ) 30 CONTINUE IF( NOUNIT ) $ TEMP = TEMP*A( K, K ) B( K, J ) = TEMP END IF 40 CONTINUE 50 CONTINUE ELSE DO 80, J = 1, N DO 70 K = M, 1, -1 IF( B( K, J ).NE.ZERO )THEN TEMP = ALPHA*B( K, J ) B( K, J ) = TEMP IF( NOUNIT ) $ B( K, J ) = B( K, J )*A( K, K ) DO 60, I = K + 1, M B( I, J ) = B( I, J ) + TEMP*A( I, K ) 60 CONTINUE END IF 70 CONTINUE 80 CONTINUE END IF ELSE * * FORM B := ALPHA*B*A'. * IF( UPPER )THEN DO 110, J = 1, N DO 100, I = M, 1, -1 TEMP = B( I, J ) IF( NOUNIT ) $ TEMP = TEMP*A( I, I ) DO 90, K = 1, I - 1 TEMP = TEMP + A( K, I )*B( K, J ) 90 CONTINUE B( I, J ) = ALPHA*TEMP 100 CONTINUE 110 CONTINUE ELSE DO 140, J = 1, N DO 130, I = 1, M TEMP = B( I, J ) IF( NOUNIT ) $ TEMP = TEMP*A( I, I ) DO 120, K = I + 1, M TEMP = TEMP + A( K, I )*B( K, J ) 120 CONTINUE B( I, J ) = ALPHA*TEMP 130 CONTINUE 140 CONTINUE END IF END IF ELSE IF( LSAME( TRANSA, 'N' ) )THEN * * FORM B := ALPHA*B*A. * IF( UPPER )THEN DO 180, J = N, 1, -1 TEMP = ALPHA IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 150, I = 1, M B( I, J ) = TEMP*B( I, J ) 150 CONTINUE DO 170, K = 1, J - 1 IF( A( K, J ).NE.ZERO )THEN TEMP = ALPHA*A( K, J ) DO 160, I = 1, M B( I, J ) = B( I, J ) + TEMP*B( I, K ) 160 CONTINUE END IF 170 CONTINUE 180 CONTINUE ELSE DO 220, J = 1, N TEMP = ALPHA IF( NOUNIT ) $ TEMP = TEMP*A( J, J ) DO 190, I = 1, M B( I, J ) = TEMP*B( I, J ) 190 CONTINUE DO 210, K = J + 1, N IF( A( K, J ).NE.ZERO )THEN TEMP = ALPHA*A( K, J ) DO 200, I = 1, M B( I, J ) = B( I, J ) + TEMP*B( I, K ) 200 CONTINUE END IF 210 CONTINUE 220 CONTINUE END IF ELSE * * FORM B := ALPHA*B*A'. * IF( UPPER )THEN DO 260, K = 1, N DO 240, J = 1, K - 1 IF( A( J, K ).NE.ZERO )THEN TEMP = ALPHA*A( J, K ) DO 230, I = 1, M B( I, J ) = B( I, J ) + TEMP*B( I, K ) 230 CONTINUE END IF 240 CONTINUE TEMP = ALPHA IF( NOUNIT ) $ TEMP = TEMP*A( K, K ) IF( TEMP.NE.ONE )THEN DO 250, I = 1, M B( I, K ) = TEMP*B( I, K ) 250 CONTINUE END IF 260 CONTINUE ELSE DO 300, K = N, 1, -1 DO 280, J = K + 1, N IF( A( J, K ).NE.ZERO )THEN TEMP = ALPHA*A( J, K ) DO 270, I = 1, M B( I, J ) = B( I, J ) + TEMP*B( I, K ) 270 CONTINUE END IF 280 CONTINUE TEMP = ALPHA IF( NOUNIT ) $ TEMP = TEMP*A( K, K ) IF( TEMP.NE.ONE )THEN DO 290, I = 1, M B( I, K ) = TEMP*B( I, K ) 290 CONTINUE END IF 300 CONTINUE END IF END IF END IF * RETURN * * END OF DTRMM . * END CUT HERE............