Free cookie consent management tool by TermsFeed Policy Generator

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

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

#2789 added Finbarr O'Sullivan smoothing spline code

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