32 #ifndef __VMML__VMMLIB_LAPACK_LINEAR_LEAST_SQUARES__HPP__
33 #define __VMML__VMMLIB_LAPACK_LINEAR_LEAST_SQUARES__HPP__
35 #include <vmmlib/matrix.hpp>
36 #include <vmmlib/vector.hpp>
37 #include <vmmlib/exception.hpp>
39 #include <vmmlib/lapack_types.hpp>
40 #include <vmmlib/lapack_includes.hpp>
73 template<
typename float_t >
88 friend std::ostream& operator << ( std::ostream& os,
97 <<
" lwork " << p.lwork
108 void dgels_(
const char *trans,
const int *M,
const int *N,
const int *nrhs,
109 double *A,
const int *lda,
double *b,
const int *ldb,
double *work,
110 const int * lwork,
int *info);
113 template<
typename float_t >
117 VMMLIB_ERROR(
"not implemented for this type.", VMMLIB_HERE );
122 llsq_call_xgels( llsq_params_xgels< float >& p )
141 llsq_call_xgels( llsq_params_xgels< double >& p )
160 template<
size_t M,
size_t N,
typename float_t >
185 template<
size_t M,
size_t N,
typename float_t >
195 llsq_call_xgels( p );
200 for(
size_t index = 0; index < N; ++index )
202 x( index ) = _b( index );
209 VMMLIB_ERROR(
"xGELS - invalid argument.", VMMLIB_HERE );
213 std::cout <<
"A\n" << A << std::endl;
214 std::cout <<
"B\n" << B << std::endl;
216 VMMLIB_ERROR(
"least squares solution could not be computed.",
224 template<
size_t M,
size_t N,
typename float_t >
225 linear_least_squares_xgels< M, N, float_t >::
226 linear_least_squares_xgels()
236 p.work =
new float_t();
240 llsq_call_xgels( p );
242 p.lwork =
static_cast< lapack_int
> ( p.work[0] );
245 p.work =
new float_t[ p.lwork ];
250 template<
size_t M,
size_t N,
typename float_t >
251 linear_least_squares_xgels< M, N, float_t >::
252 ~linear_least_squares_xgels()
265 template<
typename float_t >
277 friend std::ostream& operator << ( std::ostream& os,
282 <<
" nrhs " << p.nrhs
285 <<
" info " << p.info
294 int dgesv_(integer *n, integer *nrhs, doublereal *a, integer
295 *lda, integer *ipiv, doublereal *b, integer *ldb, integer *info);
299 template<
typename float_t >
303 VMMLIB_ERROR(
"not implemented for this type.", VMMLIB_HERE );
309 llsq_call_xgesv( llsq_params_xgesv< float >& p )
327 llsq_call_xgesv( llsq_params_xgesv< double >& p )
342 template<
size_t M,
size_t N,
typename float_t >
361 template<
size_t M,
size_t N,
typename float_t >
372 lapack::llsq_call_xgesv( p );
377 VMMLIB_ERROR(
"invalid value in input matrix", VMMLIB_HERE );
379 VMMLIB_ERROR(
"factor U is exactly singular, solution could not be computed.", VMMLIB_HERE );
385 template<
size_t M,
size_t N,
typename float_t >
386 linear_least_squares_xgesv< M, N, float_t >::
387 linear_least_squares_xgesv()
393 p.ipiv =
new lapack_int[ N ];
399 template<
size_t M,
size_t N,
typename float_t >
400 linear_least_squares_xgesv< M, N, float_t >::
401 ~linear_least_squares_xgesv()