用 GSL 求解超定方程组及矩阵的奇异值分解(SVD) 2
2016-12-10 19:34
357 查看
接上一篇。。。
下面我们将 SVD 相关的功能封装成一个类,以方便我们提取 和 的值。
另外,当我们一个 有多组 需要求解时,也只需要计算一次 SVD 分解,用下面的类能减少很多计算量。
头文件如下:
程序文件如下:
下面是个简单的测试代码:
下面我们将 SVD 相关的功能封装成一个类,以方便我们提取 和 的值。
另外,当我们一个 有多组 需要求解时,也只需要计算一次 SVD 分解,用下面的类能减少很多计算量。
头文件如下:
#ifndef GSLSINGULARVALUEDECOMPOSITION_H #define GSLSINGULARVALUEDECOMPOSITION_H #include <gsl/gsl_matrix.h> #include <gsl/gsl_vector.h> #include <gsl/gsl_blas.h> #include <gsl/gsl_linalg.h> #include <gsl/gsl_errno.h> void linearSolve_SVD(const gsl_matrix * A, const gsl_vector * b, gsl_vector * x); class GslSVD { public: GslSVD(); ~GslSVD(); int SV_decomp(const gsl_matrix * A); int SV_decomp_mod(const gsl_matrix * A); int SV_decomp_jacobi (gsl_matrix * A); int SV_solve(const gsl_vector *b, gsl_vector *x); gsl_vector * getVectorS(); gsl_matrix * getMatrixU(); gsl_matrix * getMatrixV(); int trimVectorS(double abseps); private: gsl_vector * S; gsl_matrix * U; gsl_matrix * V; void alloc_suv(int rows, int cols); }; #endif // GSLSINGULARVALUEDECOMPOSITION_H
程序文件如下:
#include "gsl_SVD.h" void linearSolve_SVD(const gsl_matrix * A, const gsl_vector * b, gsl_vector * x) { int rows = A->size1; int cols = A->size2; gsl_vector * work = gsl_vector_alloc (cols); gsl_vector * S = gsl_vector_alloc (cols); gsl_matrix * U = gsl_matrix_alloc(rows, cols);; gsl_matrix * V = gsl_matrix_alloc(cols, cols); gsl_matrix_memcpy (U, A); // 为了不破坏 A 中原始的数据,这里全都拷贝到 U 中 gsl_linalg_SV_decomp( U, V, S, work ); gsl_linalg_SV_solve ( U, V, S, b, x ); gsl_vector_free(work); gsl_vector_free(S); gsl_matrix_free(V); gsl_matrix_free(U); } int GslSVD::trimVectorS(double abseps) { int count = 0; for(int i = 0; i < S->size; i++) { if(fabs(gsl_vector_get(S, i)) < abseps) { count ++; gsl_vector_set(S, i, 0); } } return count; } gsl_vector * GslSVD::getVectorS() { if(S == NULL) return NULL; gsl_vector * s = gsl_vector_alloc(S->size); gsl_vector_memcpy(s, S); return s; } gsl_matrix * GslSVD::getMatrixU() { if(U == NULL) return NULL; gsl_matrix * u = gsl_matrix_alloc(U->size1, U->size2); gsl_matrix_memcpy(u, U); return u; } gsl_matrix * GslSVD::getMatrixV() { if(V == NULL) return NULL; gsl_matrix * v = gsl_matrix_alloc(V->size1, V->size2); gsl_matrix_memcpy(v, V); return v; } GslSVD::GslSVD() { S = NULL; U = NULL; V = NULL; } void GslSVD::alloc_suv(int rows, int cols) { if( S != NULL ) { gsl_vector_free(S); gsl_matrix_free(U); gsl_matrix_free(V); } S = gsl_vector_alloc (cols); U = gsl_matrix_alloc(rows, cols); V = gsl_matrix_alloc(cols, cols); } int GslSVD::SV_decomp(const gsl_matrix * A) { int rows = A->size1; int cols = A->size2; gsl_vector * work = gsl_vector_alloc (cols); alloc_suv(rows, cols); gsl_matrix_memcpy (U, A); // 为了不破坏 A 中原始的数据,这里全都拷贝到 U 中 int ret = gsl_linalg_SV_decomp( U, V, S, work ); gsl_vector_free(work); return ret; } int GslSVD::SV_decomp_mod(const gsl_matrix * A) { int rows = A->size1; int cols = A->size2; gsl_vector * work = gsl_vector_alloc (cols); gsl_matrix *X = gsl_matrix_alloc(cols, cols); alloc_suv(rows, cols); gsl_matrix_memcpy (U, A); // 为了不破坏 A 中原始的数据,这里全都拷贝到 U 中 int ret = gsl_linalg_SV_decomp_mod( U, X, V, S, work ); gsl_matrix_free(X); gsl_vector_free(work); return ret; } int GslSVD::SV_decomp_jacobi (gsl_matrix * A) { int rows = A->size1; int cols = A->size2; alloc_suv(rows, cols); gsl_matrix_memcpy (U, A); // 为了不破坏 A 中原始的数据,这里全都拷贝到 U 中 int ret = gsl_linalg_SV_decomp_jacobi( U, V, S ); return ret; } int GslSVD::SV_solve(const gsl_vector *b, gsl_vector *x) { if(U != NULL) { return gsl_linalg_SV_solve (U, V, S, b, x); } return -1; } GslSVD::~GslSVD() { if(S != NULL) { gsl_vector_free(S); gsl_matrix_free(V); gsl_matrix_free(U); } }
下面是个简单的测试代码:
void test5() { double a_data[] = {1, 2, 2, 4}; gsl_matrix_view A = gsl_matrix_view_array (a_data, 2, 2); GslSVD svd; svd.SV_decomp(&A.matrix); puts("S = "); gsl_vector_fprintf (stdout, svd.getVectorS(), "%f"); puts("\nV = "); gsl_matrix_fprintf (stdout, svd.getMatrixV(), "%f"); double b_data[] = {3, 6}; gsl_vector_view b = gsl_vector_view_array (b_data, 2); gsl_vector * x = gsl_vector_alloc (2); svd.SV_solve(&b.vector, x); puts("\nx = "); gsl_vector_fprintf (stdout, x, "%f"); }
相关文章推荐
- Github仓库创建及git常用命令使用
- 第十五周oj训练—— 将字符串格式化后输出(串)(2284)
- Linux下安装JDK和Eclipse
- 第十六周 项目一 - 冒泡排序
- Linux下5种动态库运行时搜索路径的方法
- js实现全选,全不选,反选,批量删除等(Yii)
- 单例模式
- HTML5-简单触摸事件
- Exception sending context initialized event to listener instance of class org.springframework.web.co
- Entity Framework Code First ---EF Power Tool MySql
- 两个TextView控件居中显示
- 第十六周 项目一 - 基数排序
- 基数排序总结
- 1063. 计算谱半径(20)
- Android 微信第三方登录(个人笔记)
- Mac Eclipse Tomcat 权限问题
- JavaSE学习笔记之NIO.2
- C++运算符重载(二)
- 布隆过滤器
- JAVA WEB 学习笔记 Tomcat配置-附光纤光猫破解。