Free cookie consent management tool by TermsFeed Policy Generator

source: branches/2789_MathNetNumerics-Exploration/HeuristicLab.Algorithms.DataAnalysis.Experimental/sbart/cher.f @ 17328

Last change on this file since 17328 was 15457, checked in by gkronber, 7 years ago

#2789 added Finbarr O'Sullivan smoothing spline code

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