← Back to branch summary

~ubuntu-branches/ubuntu/vivid/atlas/vivid

« back to all changes in this revision

Viewing changes to src/blas/f77reference/dtpmv.f

  • Committer: Package Import Robot
  • Author(s): Sébastien Villemot
  • Date: 2013-06-11 15:58:16 UTC
  • mfrom: (1.1.3 upstream)
  • mto: (2.2.21 experimental)
  • mto: This revision was merged to the branch mainline in revision 26.
  • Revision ID: package-import@ubuntu.com-20130611155816-b72z8f621tuhbzn0
Tags: upstream-3.10.1
Import upstream version 3.10.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/blas/f77reference/dtpmv.f
 
1
      SUBROUTINE DTPMV(UPLO,TRANS,DIAG,N,AP,X,INCX)
 
2
*     .. Scalar Arguments ..
 
3
      INTEGER INCX,N
 
4
      CHARACTER DIAG,TRANS,UPLO
 
5
*     ..
 
6
*     .. Array Arguments ..
 
7
      DOUBLE PRECISION AP(*),X(*)
 
8
*     ..
 
9
*
 
10
*  Purpose
 
11
*  =======
 
12
*
 
13
*  DTPMV  performs one of the matrix-vector operations
 
14
*
 
15
*     x := A*x,   or   x := A'*x,
 
16
*
 
17
*  where x is an n element vector and  A is an n by n unit, or non-unit,
 
18
*  upper or lower triangular matrix, supplied in packed form.
 
19
*
 
20
*  Arguments
 
21
*  ==========
 
22
*
 
23
*  UPLO   - CHARACTER*1.
 
24
*           On entry, UPLO specifies whether the matrix is an upper or
 
25
*           lower triangular matrix as follows:
 
26
*
 
27
*              UPLO = 'U' or 'u'   A is an upper triangular matrix.
 
28
*
 
29
*              UPLO = 'L' or 'l'   A is a lower triangular matrix.
 
30
*
 
31
*           Unchanged on exit.
 
32
*
 
33
*  TRANS  - CHARACTER*1.
 
34
*           On entry, TRANS specifies the operation to be performed as
 
35
*           follows:
 
36
*
 
37
*              TRANS = 'N' or 'n'   x := A*x.
 
38
*
 
39
*              TRANS = 'T' or 't'   x := A'*x.
 
40
*
 
41
*              TRANS = 'C' or 'c'   x := A'*x.
 
42
*
 
43
*           Unchanged on exit.
 
44
*
 
45
*  DIAG   - CHARACTER*1.
 
46
*           On entry, DIAG specifies whether or not A is unit
 
47
*           triangular as follows:
 
48
*
 
49
*              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
 
50
*
 
51
*              DIAG = 'N' or 'n'   A is not assumed to be unit
 
52
*                                  triangular.
 
53
*
 
54
*           Unchanged on exit.
 
55
*
 
56
*  N      - INTEGER.
 
57
*           On entry, N specifies the order of the matrix A.
 
58
*           N must be at least zero.
 
59
*           Unchanged on exit.
 
60
*
 
61
*  AP     - DOUBLE PRECISION array of DIMENSION at least
 
62
*           ( ( n*( n + 1 ) )/2 ).
 
63
*           Before entry with  UPLO = 'U' or 'u', the array AP must
 
64
*           contain the upper triangular matrix packed sequentially,
 
65
*           column by column, so that AP( 1 ) contains a( 1, 1 ),
 
66
*           AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 )
 
67
*           respectively, and so on.
 
68
*           Before entry with UPLO = 'L' or 'l', the array AP must
 
69
*           contain the lower triangular matrix packed sequentially,
 
70
*           column by column, so that AP( 1 ) contains a( 1, 1 ),
 
71
*           AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 )
 
72
*           respectively, and so on.
 
73
*           Note that when  DIAG = 'U' or 'u', the diagonal elements of
 
74
*           A are not referenced, but are assumed to be unity.
 
75
*           Unchanged on exit.
 
76
*
 
77
*  X      - DOUBLE PRECISION array of dimension at least
 
78
*           ( 1 + ( n - 1 )*abs( INCX ) ).
 
79
*           Before entry, the incremented array X must contain the n
 
80
*           element vector x. On exit, X is overwritten with the
 
81
*           tranformed vector x.
 
82
*
 
83
*  INCX   - INTEGER.
 
84
*           On entry, INCX specifies the increment for the elements of
 
85
*           X. INCX must not be zero.
 
86
*           Unchanged on exit.
 
87
*
 
88
*
 
89
*  Level 2 Blas routine.
 
90
*
 
91
*  -- Written on 22-October-1986.
 
92
*     Jack Dongarra, Argonne National Lab.
 
93
*     Jeremy Du Croz, Nag Central Office.
 
94
*     Sven Hammarling, Nag Central Office.
 
95
*     Richard Hanson, Sandia National Labs.
 
96
*
 
97
*
 
98
*     .. Parameters ..
 
99
      DOUBLE PRECISION ZERO
 
100
      PARAMETER (ZERO=0.0D+0)
 
101
*     ..
 
102
*     .. Local Scalars ..
 
103
      DOUBLE PRECISION TEMP
 
104
      INTEGER I,INFO,IX,J,JX,K,KK,KX
 
105
      LOGICAL NOUNIT
 
106
*     ..
 
107
*     .. External Functions ..
 
108
      LOGICAL LSAME
 
109
      EXTERNAL LSAME
 
110
*     ..
 
111
*     .. External Subroutines ..
 
112
      EXTERNAL XERBLA
 
113
*     ..
 
114
*
 
115
*     Test the input parameters.
 
116
*
 
117
      INFO = 0
 
118
      IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
 
119
          INFO = 1
 
120
      ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
 
121
     +         .NOT.LSAME(TRANS,'C')) THEN
 
122
          INFO = 2
 
123
      ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN
 
124
          INFO = 3
 
125
      ELSE IF (N.LT.0) THEN
 
126
          INFO = 4
 
127
      ELSE IF (INCX.EQ.0) THEN
 
128
          INFO = 7
 
129
      END IF
 
130
      IF (INFO.NE.0) THEN
 
131
          CALL XERBLA('DTPMV ',INFO)
 
132
          RETURN
 
133
      END IF
 
134
*
 
135
*     Quick return if possible.
 
136
*
 
137
      IF (N.EQ.0) RETURN
 
138
*
 
139
      NOUNIT = LSAME(DIAG,'N')
 
140
*
 
141
*     Set up the start point in X if the increment is not unity. This
 
142
*     will be  ( N - 1 )*INCX  too small for descending loops.
 
143
*
 
144
      IF (INCX.LE.0) THEN
 
145
          KX = 1 - (N-1)*INCX
 
146
      ELSE IF (INCX.NE.1) THEN
 
147
          KX = 1
 
148
      END IF
 
149
*
 
150
*     Start the operations. In this version the elements of AP are
 
151
*     accessed sequentially with one pass through AP.
 
152
*
 
153
      IF (LSAME(TRANS,'N')) THEN
 
154
*
 
155
*        Form  x:= A*x.
 
156
*
 
157
          IF (LSAME(UPLO,'U')) THEN
 
158
              KK = 1
 
159
              IF (INCX.EQ.1) THEN
 
160
                  DO 20 J = 1,N
 
161
                      IF (X(J).NE.ZERO) THEN
 
162
                          TEMP = X(J)
 
163
                          K = KK
 
164
                          DO 10 I = 1,J - 1
 
165
                              X(I) = X(I) + TEMP*AP(K)
 
166
                              K = K + 1
 
167
   10                     CONTINUE
 
168
                          IF (NOUNIT) X(J) = X(J)*AP(KK+J-1)
 
169
                      END IF
 
170
                      KK = KK + J
 
171
   20             CONTINUE
 
172
              ELSE
 
173
                  JX = KX
 
174
                  DO 40 J = 1,N
 
175
                      IF (X(JX).NE.ZERO) THEN
 
176
                          TEMP = X(JX)
 
177
                          IX = KX
 
178
                          DO 30 K = KK,KK + J - 2
 
179
                              X(IX) = X(IX) + TEMP*AP(K)
 
180
                              IX = IX + INCX
 
181
   30                     CONTINUE
 
182
                          IF (NOUNIT) X(JX) = X(JX)*AP(KK+J-1)
 
183
                      END IF
 
184
                      JX = JX + INCX
 
185
                      KK = KK + J
 
186
   40             CONTINUE
 
187
              END IF
 
188
          ELSE
 
189
              KK = (N* (N+1))/2
 
190
              IF (INCX.EQ.1) THEN
 
191
                  DO 60 J = N,1,-1
 
192
                      IF (X(J).NE.ZERO) THEN
 
193
                          TEMP = X(J)
 
194
                          K = KK
 
195
                          DO 50 I = N,J + 1,-1
 
196
                              X(I) = X(I) + TEMP*AP(K)
 
197
                              K = K - 1
 
198
   50                     CONTINUE
 
199
                          IF (NOUNIT) X(J) = X(J)*AP(KK-N+J)
 
200
                      END IF
 
201
                      KK = KK - (N-J+1)
 
202
   60             CONTINUE
 
203
              ELSE
 
204
                  KX = KX + (N-1)*INCX
 
205
                  JX = KX
 
206
                  DO 80 J = N,1,-1
 
207
                      IF (X(JX).NE.ZERO) THEN
 
208
                          TEMP = X(JX)
 
209
                          IX = KX
 
210
                          DO 70 K = KK,KK - (N- (J+1)),-1
 
211
                              X(IX) = X(IX) + TEMP*AP(K)
 
212
                              IX = IX - INCX
 
213
   70                     CONTINUE
 
214
                          IF (NOUNIT) X(JX) = X(JX)*AP(KK-N+J)
 
215
                      END IF
 
216
                      JX = JX - INCX
 
217
                      KK = KK - (N-J+1)
 
218
   80             CONTINUE
 
219
              END IF
 
220
          END IF
 
221
      ELSE
 
222
*
 
223
*        Form  x := A'*x.
 
224
*
 
225
          IF (LSAME(UPLO,'U')) THEN
 
226
              KK = (N* (N+1))/2
 
227
              IF (INCX.EQ.1) THEN
 
228
                  DO 100 J = N,1,-1
 
229
                      TEMP = X(J)
 
230
                      IF (NOUNIT) TEMP = TEMP*AP(KK)
 
231
                      K = KK - 1
 
232
                      DO 90 I = J - 1,1,-1
 
233
                          TEMP = TEMP + AP(K)*X(I)
 
234
                          K = K - 1
 
235
   90                 CONTINUE
 
236
                      X(J) = TEMP
 
237
                      KK = KK - J
 
238
  100             CONTINUE
 
239
              ELSE
 
240
                  JX = KX + (N-1)*INCX
 
241
                  DO 120 J = N,1,-1
 
242
                      TEMP = X(JX)
 
243
                      IX = JX
 
244
                      IF (NOUNIT) TEMP = TEMP*AP(KK)
 
245
                      DO 110 K = KK - 1,KK - J + 1,-1
 
246
                          IX = IX - INCX
 
247
                          TEMP = TEMP + AP(K)*X(IX)
 
248
  110                 CONTINUE
 
249
                      X(JX) = TEMP
 
250
                      JX = JX - INCX
 
251
                      KK = KK - J
 
252
  120             CONTINUE
 
253
              END IF
 
254
          ELSE
 
255
              KK = 1
 
256
              IF (INCX.EQ.1) THEN
 
257
                  DO 140 J = 1,N
 
258
                      TEMP = X(J)
 
259
                      IF (NOUNIT) TEMP = TEMP*AP(KK)
 
260
                      K = KK + 1
 
261
                      DO 130 I = J + 1,N
 
262
                          TEMP = TEMP + AP(K)*X(I)
 
263
                          K = K + 1
 
264
  130                 CONTINUE
 
265
                      X(J) = TEMP
 
266
                      KK = KK + (N-J+1)
 
267
  140             CONTINUE
 
268
              ELSE
 
269
                  JX = KX
 
270
                  DO 160 J = 1,N
 
271
                      TEMP = X(JX)
 
272
                      IX = JX
 
273
                      IF (NOUNIT) TEMP = TEMP*AP(KK)
 
274
                      DO 150 K = KK + 1,KK + N - J
 
275
                          IX = IX + INCX
 
276
                          TEMP = TEMP + AP(K)*X(IX)
 
277
  150                 CONTINUE
 
278
                      X(JX) = TEMP
 
279
                      JX = JX + INCX
 
280
                      KK = KK + (N-J+1)
 
281
  160             CONTINUE
 
282
              END IF
 
283
          END IF
 
284
      END IF
 
285
*
 
286
      RETURN
 
287
*
 
288
*     End of DTPMV .
 
289
*
 
290
      END