1 #include <blas_extended.h>
2 #include <blas_extended_private.h>
3 #include <blas_fpu.h>
BLAS_dsymv2_s_d(enum blas_order_type order,enum blas_uplo_type uplo,int n,double alpha,const float * a,int lda,const double * x_head,const double * x_tail,int incx,double beta,double * y,int incy)4 void BLAS_dsymv2_s_d(enum blas_order_type order, enum blas_uplo_type uplo,
5 int n, double alpha, const float *a, int lda,
6 const double *x_head, const double *x_tail, int incx,
7 double beta, double *y, int incy)
8
9 /*
10 * Purpose
11 * =======
12 *
13 * This routines computes the matrix product:
14 *
15 * y <- alpha * A * (x_head + x_tail) + beta * y
16 *
17 * where A is a symmetric matrix.
18 *
19 * Arguments
20 * =========
21 *
22 * order (input) enum blas_order_type
23 * Storage format of input symmetric matrix A.
24 *
25 * uplo (input) enum blas_uplo_type
26 * Determines which half of matrix A (upper or lower triangle)
27 * is accessed.
28 *
29 * n (input) int
30 * Dimension of A and size of vectors x, y.
31 *
32 * alpha (input) double
33 *
34 * a (input) float*
35 * Matrix A.
36 *
37 * lda (input) int
38 * Leading dimension of matrix A.
39 *
40 * x_head (input) double*
41 * Vector x_head
42 *
43 * x_tail (input) double*
44 * Vector x_tail
45 *
46 * incx (input) int
47 * Stride for vector x.
48 *
49 * beta (input) double
50 *
51 * y (input) double*
52 * Vector y.
53 *
54 * incy (input) int
55 * Stride for vector y.
56 *
57 */
58 {
59 /* Routine name */
60 const char routine_name[] = "BLAS_dsymv2_s_d";
61
62 int i, j;
63 int xi, yi, xi0, yi0;
64 int aij, ai;
65 int incai;
66 int incaij, incaij2;
67
68 const float *a_i = a;
69 const double *x_head_i = x_head;
70 const double *x_tail_i = x_tail;
71 double *y_i = y;
72 double alpha_i = alpha;
73 double beta_i = beta;
74 float a_elem;
75 double x_elem;
76 double y_elem;
77 double prod1;
78 double prod2;
79 double sum1;
80 double sum2;
81 double tmp1;
82 double tmp2;
83 double tmp3;
84
85
86
87 /* Test for no-op */
88 if (n <= 0) {
89 return;
90 }
91 if (alpha_i == 0.0 && beta_i == 1.0) {
92 return;
93 }
94
95 /* Check for error conditions. */
96 if (n < 0) {
97 BLAS_error(routine_name, -3, n, NULL);
98 }
99 if (lda < n) {
100 BLAS_error(routine_name, -6, n, NULL);
101 }
102 if (incx == 0) {
103 BLAS_error(routine_name, -9, incx, NULL);
104 }
105 if (incy == 0) {
106 BLAS_error(routine_name, -12, incy, NULL);
107 }
108
109 if ((order == blas_colmajor && uplo == blas_upper) ||
110 (order == blas_rowmajor && uplo == blas_lower)) {
111 incai = lda;
112 incaij = 1;
113 incaij2 = lda;
114 } else {
115 incai = 1;
116 incaij = lda;
117 incaij2 = 1;
118 }
119
120
121
122
123
124
125 xi0 = (incx > 0) ? 0 : ((-n + 1) * incx);
126 yi0 = (incy > 0) ? 0 : ((-n + 1) * incy);
127
128
129
130 /* The most general form, y <--- alpha * A * (x_head + x_tail) + beta * y */
131 for (i = 0, yi = yi0, ai = 0; i < n; i++, yi += incy, ai += incai) {
132 sum1 = 0.0;
133 sum2 = 0.0;
134
135 for (j = 0, aij = ai, xi = xi0; j < i; j++, aij += incaij, xi += incx) {
136 a_elem = a_i[aij];
137 x_elem = x_head_i[xi];
138 prod1 = a_elem * x_elem;
139 sum1 = sum1 + prod1;
140 x_elem = x_tail_i[xi];
141 prod2 = a_elem * x_elem;
142 sum2 = sum2 + prod2;
143 }
144 for (; j < n; j++, aij += incaij2, xi += incx) {
145 a_elem = a_i[aij];
146 x_elem = x_head_i[xi];
147 prod1 = a_elem * x_elem;
148 sum1 = sum1 + prod1;
149 x_elem = x_tail_i[xi];
150 prod2 = a_elem * x_elem;
151 sum2 = sum2 + prod2;
152 }
153 sum1 = sum1 + sum2;
154 tmp1 = sum1 * alpha_i;
155 y_elem = y_i[yi];
156 tmp2 = y_elem * beta_i;
157 tmp3 = tmp1 + tmp2;
158 y_i[yi] = tmp3;
159 }
160
161
162
163 } /* end BLAS_dsymv2_s_d */
164