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