1 /* =========================================================================
2 Copyright (c) 2010-2014, Institute for Microelectronics,
3 Institute for Analysis and Scientific Computing,
4 TU Wien.
5 Portions of this software are copyright by UChicago Argonne, LLC.
6
7 -----------------
8 ViennaCL - The Vienna Computing Library
9 -----------------
10
11 Project Head: Karl Rupp rupp@iue.tuwien.ac.at
12
13 (A list of authors and contributors can be found in the PDF manual)
14
15 License: MIT (X11), see file LICENSE in the base directory
16 ============================================================================= */
17
18 // include necessary system headers
19 #include <iostream>
20
21 #include "viennacl.hpp"
22 #include "viennacl_private.hpp"
23
24 #include "init_vector.hpp"
25 #include "init_matrix.hpp"
26
27 //include basic scalar and vector types of ViennaCL
28 #include "viennacl/scalar.hpp"
29 #include "viennacl/vector.hpp"
30 #include "viennacl/matrix.hpp"
31 #include "viennacl/linalg/direct_solve.hpp"
32 #include "viennacl/linalg/prod.hpp"
33
34 // GEMV
35
ViennaCLgemv(ViennaCLHostScalar alpha,ViennaCLMatrix A,ViennaCLVector x,ViennaCLHostScalar beta,ViennaCLVector y)36 VIENNACL_EXPORTED_FUNCTION ViennaCLStatus ViennaCLgemv(ViennaCLHostScalar alpha, ViennaCLMatrix A, ViennaCLVector x, ViennaCLHostScalar beta, ViennaCLVector y)
37 {
38 viennacl::backend::mem_handle v1_handle;
39 viennacl::backend::mem_handle v2_handle;
40 viennacl::backend::mem_handle A_handle;
41
42 if (init_vector(v1_handle, x) != ViennaCLSuccess)
43 return ViennaCLGenericFailure;
44
45 if (init_vector(v2_handle, y) != ViennaCLSuccess)
46 return ViennaCLGenericFailure;
47
48 if (init_matrix(A_handle, A) != ViennaCLSuccess)
49 return ViennaCLGenericFailure;
50
51 switch (x->precision)
52 {
53 case ViennaCLFloat:
54 {
55 typedef viennacl::vector_base<float>::size_type size_type;
56 typedef viennacl::vector_base<float>::size_type difference_type;
57
58 viennacl::vector_base<float> v1(v1_handle, size_type(x->size), size_type(x->offset), difference_type(x->inc));
59 viennacl::vector_base<float> v2(v2_handle, size_type(y->size), size_type(y->offset), difference_type(y->inc));
60
61 viennacl::matrix_base<float> mat(A_handle,
62 size_type(A->size1), size_type(A->start1), difference_type(A->stride1), size_type(A->internal_size1),
63 size_type(A->size2), size_type(A->start2), difference_type(A->stride2), size_type(A->internal_size2), A->order == ViennaCLRowMajor);
64 v2 *= beta->value_float;
65 if (A->trans == ViennaCLTrans)
66 v2 += alpha->value_float * viennacl::linalg::prod(viennacl::trans(mat), v1);
67 else
68 v2 += alpha->value_float * viennacl::linalg::prod(mat, v1);
69
70 return ViennaCLSuccess;
71 }
72
73 case ViennaCLDouble:
74 {
75 typedef viennacl::vector_base<double>::size_type size_type;
76 typedef viennacl::vector_base<double>::size_type difference_type;
77
78 viennacl::vector_base<double> v1(v1_handle, size_type(x->size), size_type(x->offset), difference_type(x->inc));
79 viennacl::vector_base<double> v2(v2_handle, size_type(y->size), size_type(y->offset), difference_type(y->inc));
80
81 viennacl::matrix_base<double> mat(A_handle,
82 size_type(A->size1), size_type(A->start1), difference_type(A->stride1), size_type(A->internal_size1),
83 size_type(A->size2), size_type(A->start2), difference_type(A->stride2), size_type(A->internal_size2), A->order == ViennaCLRowMajor);
84 v2 *= beta->value_double;
85 if (A->trans == ViennaCLTrans)
86 v2 += alpha->value_double * viennacl::linalg::prod(viennacl::trans(mat), v1);
87 else
88 v2 += alpha->value_double * viennacl::linalg::prod(mat, v1);
89
90 return ViennaCLSuccess;
91 }
92
93 default:
94 return ViennaCLGenericFailure;
95 }
96 }
97
98
99 // xTRSV
100
ViennaCLtrsv(ViennaCLMatrix A,ViennaCLVector x,ViennaCLUplo uplo)101 VIENNACL_EXPORTED_FUNCTION ViennaCLStatus ViennaCLtrsv(ViennaCLMatrix A, ViennaCLVector x, ViennaCLUplo uplo)
102 {
103 viennacl::backend::mem_handle v1_handle;
104 viennacl::backend::mem_handle A_handle;
105
106 if (init_vector(v1_handle, x) != ViennaCLSuccess)
107 return ViennaCLGenericFailure;
108
109 if (init_matrix(A_handle, A) != ViennaCLSuccess)
110 return ViennaCLGenericFailure;
111
112 switch (x->precision)
113 {
114 case ViennaCLFloat:
115 {
116 typedef viennacl::vector_base<float>::size_type size_type;
117 typedef viennacl::vector_base<float>::size_type difference_type;
118
119 viennacl::vector_base<float> v1(v1_handle, size_type(x->size), size_type(x->offset), difference_type(x->inc));
120
121 viennacl::matrix_base<float> mat(A_handle,
122 size_type(A->size1), size_type(A->start1), difference_type(A->stride1), size_type(A->internal_size1),
123 size_type(A->size2), size_type(A->start2), difference_type(A->stride2), size_type(A->internal_size2), A->order == ViennaCLRowMajor);
124 if (A->trans == ViennaCLTrans)
125 {
126 if (uplo == ViennaCLUpper)
127 viennacl::linalg::inplace_solve(viennacl::trans(mat), v1, viennacl::linalg::upper_tag());
128 else
129 viennacl::linalg::inplace_solve(viennacl::trans(mat), v1, viennacl::linalg::lower_tag());
130 }
131 else
132 {
133 if (uplo == ViennaCLUpper)
134 viennacl::linalg::inplace_solve(mat, v1, viennacl::linalg::upper_tag());
135 else
136 viennacl::linalg::inplace_solve(mat, v1, viennacl::linalg::lower_tag());
137 }
138
139 return ViennaCLSuccess;
140 }
141 case ViennaCLDouble:
142 {
143 typedef viennacl::vector_base<double>::size_type size_type;
144 typedef viennacl::vector_base<double>::size_type difference_type;
145
146 viennacl::vector_base<double> v1(v1_handle, size_type(x->size), size_type(x->offset), difference_type(x->inc));
147
148 viennacl::matrix_base<double> mat(A_handle,
149 size_type(A->size1), size_type(A->start1), difference_type(A->stride1), size_type(A->internal_size1),
150 size_type(A->size2), size_type(A->start2), difference_type(A->stride2), size_type(A->internal_size2), A->order == ViennaCLRowMajor);
151 if (A->trans == ViennaCLTrans)
152 {
153 if (uplo == ViennaCLUpper)
154 viennacl::linalg::inplace_solve(viennacl::trans(mat), v1, viennacl::linalg::upper_tag());
155 else
156 viennacl::linalg::inplace_solve(viennacl::trans(mat), v1, viennacl::linalg::lower_tag());
157 }
158 else
159 {
160 if (uplo == ViennaCLUpper)
161 viennacl::linalg::inplace_solve(mat, v1, viennacl::linalg::upper_tag());
162 else
163 viennacl::linalg::inplace_solve(mat, v1, viennacl::linalg::lower_tag());
164 }
165
166 return ViennaCLSuccess;
167 }
168
169 default:
170 return ViennaCLGenericFailure;
171 }
172 }
173
174
175 // xGER
176
ViennaCLger(ViennaCLHostScalar alpha,ViennaCLVector x,ViennaCLVector y,ViennaCLMatrix A)177 VIENNACL_EXPORTED_FUNCTION ViennaCLStatus ViennaCLger(ViennaCLHostScalar alpha, ViennaCLVector x, ViennaCLVector y, ViennaCLMatrix A)
178 {
179 viennacl::backend::mem_handle v1_handle;
180 viennacl::backend::mem_handle v2_handle;
181 viennacl::backend::mem_handle A_handle;
182
183 if (init_vector(v1_handle, x) != ViennaCLSuccess)
184 return ViennaCLGenericFailure;
185
186 if (init_vector(v2_handle, y) != ViennaCLSuccess)
187 return ViennaCLGenericFailure;
188
189 if (init_matrix(A_handle, A) != ViennaCLSuccess)
190 return ViennaCLGenericFailure;
191
192 switch (x->precision)
193 {
194 case ViennaCLFloat:
195 {
196 typedef viennacl::vector_base<float>::size_type size_type;
197 typedef viennacl::vector_base<float>::size_type difference_type;
198
199 viennacl::vector_base<float> v1(v1_handle, size_type(x->size), size_type(x->offset), difference_type(x->inc));
200 viennacl::vector_base<float> v2(v2_handle, size_type(y->size), size_type(y->offset), difference_type(y->inc));
201
202 viennacl::matrix_base<float> mat(A_handle,
203 size_type(A->size1), size_type(A->start1), difference_type(A->stride1), size_type(A->internal_size1),
204 size_type(A->size2), size_type(A->start2), difference_type(A->stride2), size_type(A->internal_size2), A->order == ViennaCLRowMajor);
205
206 mat += alpha->value_float * viennacl::linalg::outer_prod(v1, v2);
207
208 return ViennaCLSuccess;
209 }
210 case ViennaCLDouble:
211 {
212 typedef viennacl::vector_base<double>::size_type size_type;
213 typedef viennacl::vector_base<double>::size_type difference_type;
214
215 viennacl::vector_base<double> v1(v1_handle, size_type(x->size), size_type(x->offset), difference_type(x->inc));
216 viennacl::vector_base<double> v2(v2_handle, size_type(y->size), size_type(y->offset), difference_type(y->inc));
217
218 viennacl::matrix_base<double> mat(A_handle,
219 size_type(A->size1), size_type(A->start1), difference_type(A->stride1), size_type(A->internal_size1),
220 size_type(A->size2), size_type(A->start2), difference_type(A->stride2), size_type(A->internal_size2), A->order == ViennaCLRowMajor);
221
222 mat += alpha->value_double * viennacl::linalg::outer_prod(v1, v2);
223
224 return ViennaCLSuccess;
225 }
226 default:
227 return ViennaCLGenericFailure;
228 }
229 }
230
231
232