![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsyr2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / blas|
Return to dsyr2.f CVS log | Up to [local] / rpl / lapack / blas |
1.1 bertrand 1: SUBROUTINE DSYR2(UPLO,N,ALPHA,X,INCX,Y,INCY,A,LDA)
2: * .. Scalar Arguments ..
3: DOUBLE PRECISION ALPHA
4: INTEGER INCX,INCY,LDA,N
5: CHARACTER UPLO
6: * ..
7: * .. Array Arguments ..
8: DOUBLE PRECISION A(LDA,*),X(*),Y(*)
9: * ..
10: *
11: * Purpose
12: * =======
13: *
14: * DSYR2 performs the symmetric rank 2 operation
15: *
16: * A := alpha*x*y' + alpha*y*x' + A,
17: *
18: * where alpha is a scalar, x and y are n element vectors and A is an n
19: * by n symmetric matrix.
20: *
21: * Arguments
22: * ==========
23: *
24: * UPLO - CHARACTER*1.
25: * On entry, UPLO specifies whether the upper or lower
26: * triangular part of the array A is to be referenced as
27: * follows:
28: *
29: * UPLO = 'U' or 'u' Only the upper triangular part of A
30: * is to be referenced.
31: *
32: * UPLO = 'L' or 'l' Only the lower triangular part of A
33: * is to be referenced.
34: *
35: * Unchanged on exit.
36: *
37: * N - INTEGER.
38: * On entry, N specifies the order of the matrix A.
39: * N must be at least zero.
40: * Unchanged on exit.
41: *
42: * ALPHA - DOUBLE PRECISION.
43: * On entry, ALPHA specifies the scalar alpha.
44: * Unchanged on exit.
45: *
46: * X - DOUBLE PRECISION array of dimension at least
47: * ( 1 + ( n - 1 )*abs( INCX ) ).
48: * Before entry, the incremented array X must contain the n
49: * element vector x.
50: * Unchanged on exit.
51: *
52: * INCX - INTEGER.
53: * On entry, INCX specifies the increment for the elements of
54: * X. INCX must not be zero.
55: * Unchanged on exit.
56: *
57: * Y - DOUBLE PRECISION array of dimension at least
58: * ( 1 + ( n - 1 )*abs( INCY ) ).
59: * Before entry, the incremented array Y must contain the n
60: * element vector y.
61: * Unchanged on exit.
62: *
63: * INCY - INTEGER.
64: * On entry, INCY specifies the increment for the elements of
65: * Y. INCY must not be zero.
66: * Unchanged on exit.
67: *
68: * A - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
69: * Before entry with UPLO = 'U' or 'u', the leading n by n
70: * upper triangular part of the array A must contain the upper
71: * triangular part of the symmetric matrix and the strictly
72: * lower triangular part of A is not referenced. On exit, the
73: * upper triangular part of the array A is overwritten by the
74: * upper triangular part of the updated matrix.
75: * Before entry with UPLO = 'L' or 'l', the leading n by n
76: * lower triangular part of the array A must contain the lower
77: * triangular part of the symmetric matrix and the strictly
78: * upper triangular part of A is not referenced. On exit, the
79: * lower triangular part of the array A is overwritten by the
80: * lower triangular part of the updated matrix.
81: *
82: * LDA - INTEGER.
83: * On entry, LDA specifies the first dimension of A as declared
84: * in the calling (sub) program. LDA must be at least
85: * max( 1, n ).
86: * Unchanged on exit.
87: *
88: * Further Details
89: * ===============
90: *
91: * Level 2 Blas routine.
92: *
93: * -- Written on 22-October-1986.
94: * Jack Dongarra, Argonne National Lab.
95: * Jeremy Du Croz, Nag Central Office.
96: * Sven Hammarling, Nag Central Office.
97: * Richard Hanson, Sandia National Labs.
98: *
99: * =====================================================================
100: *
101: * .. Parameters ..
102: DOUBLE PRECISION ZERO
103: PARAMETER (ZERO=0.0D+0)
104: * ..
105: * .. Local Scalars ..
106: DOUBLE PRECISION TEMP1,TEMP2
107: INTEGER I,INFO,IX,IY,J,JX,JY,KX,KY
108: * ..
109: * .. External Functions ..
110: LOGICAL LSAME
111: EXTERNAL LSAME
112: * ..
113: * .. External Subroutines ..
114: EXTERNAL XERBLA
115: * ..
116: * .. Intrinsic Functions ..
117: INTRINSIC MAX
118: * ..
119: *
120: * Test the input parameters.
121: *
122: INFO = 0
123: IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
124: INFO = 1
125: ELSE IF (N.LT.0) THEN
126: INFO = 2
127: ELSE IF (INCX.EQ.0) THEN
128: INFO = 5
129: ELSE IF (INCY.EQ.0) THEN
130: INFO = 7
131: ELSE IF (LDA.LT.MAX(1,N)) THEN
132: INFO = 9
133: END IF
134: IF (INFO.NE.0) THEN
135: CALL XERBLA('DSYR2 ',INFO)
136: RETURN
137: END IF
138: *
139: * Quick return if possible.
140: *
141: IF ((N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN
142: *
143: * Set up the start points in X and Y if the increments are not both
144: * unity.
145: *
146: IF ((INCX.NE.1) .OR. (INCY.NE.1)) THEN
147: IF (INCX.GT.0) THEN
148: KX = 1
149: ELSE
150: KX = 1 - (N-1)*INCX
151: END IF
152: IF (INCY.GT.0) THEN
153: KY = 1
154: ELSE
155: KY = 1 - (N-1)*INCY
156: END IF
157: JX = KX
158: JY = KY
159: END IF
160: *
161: * Start the operations. In this version the elements of A are
162: * accessed sequentially with one pass through the triangular part
163: * of A.
164: *
165: IF (LSAME(UPLO,'U')) THEN
166: *
167: * Form A when A is stored in the upper triangle.
168: *
169: IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
170: DO 20 J = 1,N
171: IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN
172: TEMP1 = ALPHA*Y(J)
173: TEMP2 = ALPHA*X(J)
174: DO 10 I = 1,J
175: A(I,J) = A(I,J) + X(I)*TEMP1 + Y(I)*TEMP2
176: 10 CONTINUE
177: END IF
178: 20 CONTINUE
179: ELSE
180: DO 40 J = 1,N
181: IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN
182: TEMP1 = ALPHA*Y(JY)
183: TEMP2 = ALPHA*X(JX)
184: IX = KX
185: IY = KY
186: DO 30 I = 1,J
187: A(I,J) = A(I,J) + X(IX)*TEMP1 + Y(IY)*TEMP2
188: IX = IX + INCX
189: IY = IY + INCY
190: 30 CONTINUE
191: END IF
192: JX = JX + INCX
193: JY = JY + INCY
194: 40 CONTINUE
195: END IF
196: ELSE
197: *
198: * Form A when A is stored in the lower triangle.
199: *
200: IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
201: DO 60 J = 1,N
202: IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN
203: TEMP1 = ALPHA*Y(J)
204: TEMP2 = ALPHA*X(J)
205: DO 50 I = J,N
206: A(I,J) = A(I,J) + X(I)*TEMP1 + Y(I)*TEMP2
207: 50 CONTINUE
208: END IF
209: 60 CONTINUE
210: ELSE
211: DO 80 J = 1,N
212: IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN
213: TEMP1 = ALPHA*Y(JY)
214: TEMP2 = ALPHA*X(JX)
215: IX = JX
216: IY = JY
217: DO 70 I = J,N
218: A(I,J) = A(I,J) + X(IX)*TEMP1 + Y(IY)*TEMP2
219: IX = IX + INCX
220: IY = IY + INCY
221: 70 CONTINUE
222: END IF
223: JX = JX + INCX
224: JY = JY + INCY
225: 80 CONTINUE
226: END IF
227: END IF
228: *
229: RETURN
230: *
231: * End of DSYR2 .
232: *
233: END