Hi Soufiane,
from your description it seems to me that you want to solve a
least-squares problem, while the GMRES implementation in ViennaCL is for
square systems. I suggest you use a QR-factorization as outlined in
examples/tutorial/least-squares.cpp.
If you have any good pointers for GMRES applied to least-squares
problems, feel free to share so that we can consider extending our GMRES
implementation suitably.
Best regards,
Karli
On 12/03/2013 10:54 PM, Soufiane Khiat wrote:
> Hello,
>
> I am trying to solve Ax=b
> With a Dense and Rectangular (14000 Rows and 140 Cols).
> I use as presented on the sample/documentation:
>
> viennacl::vector< Scalar > _b( m_RowsCount );
> viennacl::vector< Scalar > _x( m_ColumnsCount );
>
> viennacl::fast_copy( b.Datas, b.Datas + m_RowsCount, _b.begin() );
> _x = viennacl::linalg::solve( _A, _b, viennacl::linalg::gmres_tag() );
> viennacl::fast_copy( _x.begin(), _x.end(), x.Datas );
>
> x, b is a container of contigüe datas (float*).
> But during the solve I have an assert:
> On Matrix_Operation.hpp:
> assert( (viennacl::traits::size2(mat) == viennacl::traits::size(vec)) &&
> bool("Size check failed at v1 = prod(A, v2): size2(A) != size(v2)"));
>
> On:
> template <typename NumericT, typename F>
> void prod_impl(const matrix_base<NumericT, F> & mat,
> const vector_base<NumericT> & vec,
> vector_base<NumericT> & result)
>
> But what is strange for me is in the call stack I pass through another
> assert without anyu trouble:
> template <typename NumericT, typename F>
> vector<NumericT>
> operator-=(vector_base<NumericT> & v1,
> const viennacl::vector_expression< const matrix_base<NumericT, F>, const
> vector_base<NumericT>, viennacl::op_prod> & proxy)
> {
> assert(viennacl::traits::size1(proxy.lhs()) == v1.size() &&
> bool("Size check failed for v1 -= A * v2: size1(A) != size(v1)"));
>
> vector<NumericT> result(viennacl::traits::size1(proxy.lhs()));
> viennacl::linalg::prod_impl(proxy.lhs(), proxy.rhs(), result);
> v1 -= result;
> return v1;
> }
>
> ...
>
> I initialize _A like that:
> viennacl::matrix< Scalar > _A( m_RowsCount, m_ColumnsCount );
>
> for ( u32 i = 0; i < m_ColumnsCount; ++i )
> {
> for ( u32 j = 0; j < m_RowsCount; ++j )
> {
> _A( j, i ) = A[ i ][ j ];
> }
> }
>
> I am a new user with ViennaCL. If you see some mistake can told me where.
>
> In my opinion my system is correctly balanced:
> A : 14000x140
> b : 14000x1
> x : 140x1
>
> Cheer,
>
> Soufiane
> Software Engineer
>
>
> ------------------------------------------------------------------------------
> Sponsored by Intel(R) XDK
> Develop, test and display web and hybrid apps with a single code base.
> Download it for free now!
> http://pubads.g.doubleclick.net/gampad/clk?id=111408631&iu=/4140/ostg.clktrk
>
>
>
> _______________________________________________
> ViennaCL-devel mailing list
> ViennaCL-devel@lists.sourceforge.net
> https://lists.sourceforge.net/lists/listinfo/viennacl-devel
>
------------------------------------------------------------------------------
Sponsored by Intel(R) XDK
Develop, test and display web and hybrid apps with a single code base.
Download it for free now!
http://pubads.g.doubleclick.net/gampad/clk?id=111408631&iu=/4140/ostg.clktrk
_______________________________________________
ViennaCL-devel mailing list
ViennaCL-devel@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/viennacl-devel
