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