3309 lines
80 KiB
C++
3309 lines
80 KiB
C++
/********************************************************
|
||
* ██████╗ ██████╗████████╗██╗
|
||
* ██╔════╝ ██╔════╝╚══██╔══╝██║
|
||
* ██║ ███╗██║ ██║ ██║
|
||
* ██║ ██║██║ ██║ ██║
|
||
* ╚██████╔╝╚██████╗ ██║ ███████╗
|
||
* ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝
|
||
* Geophysical Computational Tools & Library (GCTL)
|
||
*
|
||
* Copyright (c) 2023 Yi Zhang (yizhang-geo@zju.edu.cn)
|
||
*
|
||
* GCTL is distributed under a dual licensing scheme. You can redistribute
|
||
* it and/or modify it under the terms of the GNU Lesser General Public
|
||
* License as published by the Free Software Foundation, either version 2
|
||
* of the License, or (at your option) any later version. You should have
|
||
* received a copy of the GNU Lesser General Public License along with this
|
||
* program. If not, see <http://www.gnu.org/licenses/>.
|
||
*
|
||
* If the terms and conditions of the LGPL v.2. would prevent you from using
|
||
* the GCTL, please consider the option to obtain a commercial license for a
|
||
* fee. These licenses are offered by the GCTL's original author. As a rule,
|
||
* licenses are provided "as-is", unlimited in time for a one time fee. Please
|
||
* send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget
|
||
* to include some description of your company and the realm of its activities.
|
||
* Also add information on how to contact you by electronic and paper mail.
|
||
******************************************************/
|
||
|
||
#ifndef _GCTL_SPMAT_H
|
||
#define _GCTL_SPMAT_H
|
||
|
||
#include "enum.h"
|
||
#include "array.h"
|
||
#include "matrix.h"
|
||
|
||
namespace gctl
|
||
{
|
||
template <typename MatrixValueType> class spmat;
|
||
|
||
/**
|
||
* 我们在下面定义一些常用的稀疏矩阵类型，方便我们使用。
|
||
*/
|
||
typedef spmat<int> _1i_spmat; ///< 整形浮点稀疏矩阵
|
||
typedef spmat<float> _1f_spmat; ///< 单精度浮点稀疏矩阵
|
||
typedef spmat<double> _1d_spmat; ///< 双精度浮点稀疏矩阵
|
||
typedef spmat<std::complex<float>> _1cf_spmat; ///< 单精度复浮点稀疏矩阵
|
||
typedef spmat<std::complex<double>> _1cd_spmat; ///< 双精度复浮点稀疏矩阵
|
||
|
||
/**
|
||
* @brief 稀疏矩阵的节点元素类型
|
||
*
|
||
* @tparam NodeValueType 模版变量类型
|
||
*/
|
||
template <typename NodeValueType>
|
||
struct mat_node
|
||
{
|
||
int r_id, c_id; ///< 节点的行列号
|
||
NodeValueType val; ///< 节点值
|
||
mat_node<NodeValueType> *next_right; ///< 指向同行右方节点的指针
|
||
mat_node<NodeValueType> *next_down; ///< 指向同列下方节点的指针
|
||
|
||
/**
|
||
* 构造函数
|
||
*/
|
||
mat_node()
|
||
{
|
||
r_id = c_id = -1;
|
||
next_right = next_down = nullptr;
|
||
}
|
||
|
||
mat_node(int r, int c, NodeValueType v)
|
||
{
|
||
r_id = r; c_id = c;
|
||
val = v;
|
||
next_right = next_down = nullptr;
|
||
}
|
||
|
||
NodeValueType value() const
|
||
{
|
||
return val;
|
||
}
|
||
|
||
void value(NodeValueType in_val)
|
||
{
|
||
val = in_val;
|
||
return;
|
||
}
|
||
};
|
||
|
||
/**
|
||
* @brief 基于十字指针的稀疏矩阵的模版类 定义常规的加减乘除与矩阵运算
|
||
*
|
||
* @tparam MatrixValueType 模版变量类型
|
||
*/
|
||
template <typename MatrixValueType>
|
||
class spmat
|
||
{
|
||
protected:
|
||
mat_node<MatrixValueType> **r_head, **c_head; ///< 行和列的头指针
|
||
mat_node<MatrixValueType> **r_end, **c_end; ///< 行和列的尾指针 这两个变量的使用主要是为了提高插入操作的效率
|
||
mat_node<MatrixValueType> **d_link; ///< 所有非零元素按行优先排序的指针数组，可在稀疏矩阵元素位置不变但值经常改变时快速索引并改变值
|
||
int *r_count, *c_count; ///< 行与列的非零元素的个数
|
||
int r_num, c_num; ///< 满矩阵的行列数
|
||
int d_num; ///< d_link长度
|
||
MatrixValueType zero_val; ///< 类型的0值
|
||
|
||
public:
|
||
/**
|
||
* @brief 默认构造函数
|
||
*/
|
||
spmat();
|
||
|
||
/**
|
||
* @brief 构造函数：空的行列头指针
|
||
*
|
||
* @param[in] in_rnum 矩阵行数
|
||
* @param[in] in_cnum 矩阵列数
|
||
* @param[in] null_val 矩阵零值
|
||
*/
|
||
spmat(int in_rnum, int in_cnum, MatrixValueType null_val);
|
||
|
||
/**
|
||
* @brief 拷贝构造函数
|
||
*
|
||
* @param[in] b 拷贝稀疏矩阵
|
||
*/
|
||
spmat(const spmat<MatrixValueType> &b);
|
||
|
||
/**
|
||
* @brief 析构函数
|
||
*/
|
||
virtual ~spmat();
|
||
|
||
/**
|
||
* @brief 初始化为空的行列头指针
|
||
*
|
||
* @param[in] in_rnum 矩阵行数
|
||
* @param[in] in_cnum 矩阵列数
|
||
* @param[in] null_val 矩阵零值
|
||
*/
|
||
void malloc(int in_rnum, int in_cnum, MatrixValueType null_val);
|
||
|
||
/**
|
||
* @brief 清空稀疏矩阵
|
||
*
|
||
* @note 清空整个矩阵时会删除所有元素和辅助数组，必须重新调用malloc函数再能插入元素。删除某行（列）时仅删除已有元素，可以直接再次插入元素。
|
||
*
|
||
* @param[in] index 行或者列的索引号，如果为-1则清空整个稀疏矩阵
|
||
* @param[in] mode 清空制定的行或者列，默认为行模式
|
||
*/
|
||
void clear(int index = -1, matrix_order_e mode = RowMajor);
|
||
|
||
/**
|
||
* @brief 重置稀疏矩阵，仅删除所有已有元素，可以直接插入元素。相当于对所有行调用clear函数。
|
||
*
|
||
*/
|
||
void reset();
|
||
|
||
/**
|
||
* @brief 删除所有等于零值的元素
|
||
*/
|
||
void purge();
|
||
|
||
/**
|
||
* @brief 设置所有非零值为指定值，如果输入为零值则直接调用reset函数
|
||
*
|
||
* @note 注意此函数可能会重置稀疏矩阵
|
||
*
|
||
* @param in_val 输入值
|
||
*/
|
||
void set_all(MatrixValueType in_val);
|
||
|
||
/**
|
||
* @brief 从三元组数组赋值，三元组的排列必须是按行优先方式排序，顺序增加，可有重复索引的元素（函数会检查并报错）
|
||
*
|
||
* @note 注意此函数会重置稀疏矩阵
|
||
*
|
||
* @param triplts 三元组数组
|
||
* @param val_mode 相同位置重复输入值的处理方式
|
||
*/
|
||
void set_triplts(const std::vector<mat_node<MatrixValueType> > &triplts, value_operator_e val_mode = AppendVal);
|
||
|
||
/**
|
||
* @brief 从COO格式储存的稀疏矩阵对矩阵进行赋值，元素排列必须是按行优先方式排序，顺序增加，不可有重复索引的元素（函数不检查和报错）
|
||
*
|
||
* @note 注意此函数会重置稀疏矩阵
|
||
*
|
||
* @param row_idx 非零元素的行索引数组
|
||
* @param col_idx 非零元素的列索引数组
|
||
* @param data 非零元素数组
|
||
*/
|
||
void set_coo(const array<size_t> &row_idx, const array<size_t> &col_idx, const array<MatrixValueType> &data);
|
||
|
||
/**
|
||
* @brief 从CSR格式储存的稀疏矩阵对矩阵进行赋值，元素排列必须是按行优先方式排序，顺序增加，不可有重复索引的元素（函数不检查和报错）
|
||
*
|
||
* @note 注意此函数会重置稀疏矩阵
|
||
*
|
||
* @param row_off 非零元素的行偏移数组
|
||
* @param col_idx 非零元素的列索引数组
|
||
* @param data 非零元素数组
|
||
*/
|
||
void set_csr(const array<size_t> &row_off, const array<size_t> &col_idx, const array<MatrixValueType> &data);
|
||
|
||
/**
|
||
* @brief 从二维向量对矩阵进行赋值
|
||
*
|
||
* @param[in] den_mat_ptr 二维向量
|
||
* @param[in] off_set 初始化时零值的偏差范围
|
||
*/
|
||
void set_2d_vector(const std::vector<std::vector<MatrixValueType> > &den_mat_ptr, MatrixValueType off_set);
|
||
|
||
/**
|
||
* @brief 从二维数组对矩阵进行赋值
|
||
*
|
||
* @param den_mat_ptr 二维数组的指针
|
||
* @param[in] off_set 初始化时零值的偏差范围
|
||
*/
|
||
void set_2d_array(MatrixValueType **den_mat_ptr, MatrixValueType off_set);
|
||
|
||
/**
|
||
* @brief 从matrix数组对矩阵进行赋值
|
||
*
|
||
* @param[in] den_mat_ptr matrix数组
|
||
* @param[in] off_set 初始化时零值的偏差范围
|
||
*/
|
||
void set_matrix(const matrix<MatrixValueType> &den_mat_ptr, MatrixValueType off_set);
|
||
|
||
/**
|
||
* @brief 储存非零元素的储存位置，建立快速查找索引（目前只在set_triplts函数中使用储存的位置）
|
||
*
|
||
* @note 元素数量或位置改变后失效
|
||
*
|
||
*/
|
||
void init_pattern();
|
||
|
||
/**
|
||
* @brief 如果稀疏矩阵的行列数为零，则为空的稀疏矩阵
|
||
*
|
||
* @return 矩阵是否为空
|
||
*/
|
||
bool empty() const;
|
||
|
||
/**
|
||
* @brief 返回行数
|
||
*
|
||
* @return 行数
|
||
*/
|
||
int row_size() const;
|
||
|
||
/**
|
||
* @brief 返回列数
|
||
*
|
||
* @return 列数
|
||
*/
|
||
int col_size() const;
|
||
|
||
/**
|
||
* @brief 返回非0元素个数
|
||
*
|
||
* 此函数可返回矩阵中所有或者指定行或者列的零值元素个数
|
||
*
|
||
* @param[in] index 行或列的索引
|
||
* @param[in] mode 计算行还是列的非零元素个数
|
||
*
|
||
* @return 非0元素个数
|
||
*/
|
||
int ele_size(int index = -1, matrix_order_e mode = RowMajor) const;
|
||
|
||
/**
|
||
* @brief 返回0值
|
||
*
|
||
* @return 0值
|
||
*/
|
||
MatrixValueType zero_value() const;
|
||
|
||
/**
|
||
* @brief 按矩阵形式显示所有元素，包括0元素
|
||
*/
|
||
void show_matrix(char delimiter = '\t') const;
|
||
|
||
/**
|
||
* @brief 按列表形式显示三元组
|
||
*
|
||
* @param[in] mode 按行优先或者列优先形式排序
|
||
*/
|
||
void show_list(matrix_order_e mode = RowMajor, char delimiter = ' ', bool full = false) const;
|
||
|
||
/**
|
||
* @brief 返回位于(i, j)位置的元素
|
||
*
|
||
* @param[in] i_index 行索引
|
||
* @param[in] j_index 列索引
|
||
*
|
||
* @return 矩阵元素
|
||
*/
|
||
MatrixValueType at(size_t i_index, size_t j_index) const;
|
||
|
||
/**
|
||
* @brief 输出稠密矩阵
|
||
*
|
||
* @param ret_arr2d 以引用形式返回的matrix数组
|
||
*/
|
||
void export_dense(matrix<MatrixValueType> &ret_arr2d) const;
|
||
|
||
/**
|
||
* @brief 输出矩阵中某一行或列
|
||
*
|
||
* @param[in] out_id 行或列的索引
|
||
* @param ret_arr 以引用形式返回的array数组
|
||
* @param[in] multipler 倍数，默认为1.0
|
||
* @param[in] line_mode 是否为行优先
|
||
* @param[in] val_mode 返回数组的数值操作方式
|
||
*/
|
||
void export_dense(int out_id, array<MatrixValueType> &ret_arr,
|
||
double multipler = 1.0, matrix_order_e line_mode = RowMajor,
|
||
value_operator_e val_mode = ReplaceVal) const;
|
||
|
||
/**
|
||
* @brief 使用COO格式输出稀疏矩阵
|
||
*
|
||
* @param row_idx 非零元素的行索引数组
|
||
* @param col_idx 非零元素的列索引数组
|
||
* @param data 非零元素数组
|
||
*/
|
||
void export_coo(array<size_t> &row_idx, array<size_t> &col_idx, array<MatrixValueType> &data) const;
|
||
|
||
/**
|
||
* @brief 使用CSR格式输出稀疏矩阵
|
||
*
|
||
* @param row_off 非零元素的行偏移数组
|
||
* @param col_idx 非零元素的列索引数组
|
||
* @param data 非零元素数组
|
||
*/
|
||
void export_csr(array<size_t> &row_off, array<size_t> &col_idx, array<MatrixValueType> &data) const;
|
||
|
||
/**
|
||
* @brief 添加一个元素 若插入位置已有值则替换
|
||
*
|
||
* @param[in] r 行索引
|
||
* @param[in] c 列索引
|
||
* @param[in] v 插入值
|
||
*/
|
||
void insert(int r, int c, MatrixValueType v);
|
||
|
||
/**
|
||
* @brief 删除一个元素
|
||
*
|
||
* @param[in] r 行索引
|
||
* @param[in] c 列索引
|
||
*/
|
||
void remove(int r, int c);
|
||
|
||
/**
|
||
* @brief 计算非零元素值的总和
|
||
*
|
||
* @param[in] index 行或者列的索引号，如果为-1则计算整个稀疏矩阵
|
||
* @param[in] mode 清空制定的行或者列，默认为行模式
|
||
*/
|
||
MatrixValueType sum(int index = -1, matrix_order_e mode = RowMajor);
|
||
|
||
/**
|
||
* @brief 计算行或列的模长
|
||
*
|
||
* @param index 行或者列的索引号
|
||
* @param mode 切换行与列模式
|
||
* @return 测度值
|
||
*/
|
||
MatrixValueType module(int index, matrix_order_e mode = RowMajor);
|
||
|
||
/**
|
||
* @brief 单个元素位置的加法
|
||
*
|
||
* @param[in] r 行索引
|
||
* @param[in] c 列索引
|
||
* @param[in] v 被加数
|
||
*/
|
||
void add_scalar(int r, int c, MatrixValueType v);
|
||
|
||
/**
|
||
* @brief 单个元素位置的减法
|
||
*
|
||
* @param[in] r 行索引
|
||
* @param[in] c 列索引
|
||
* @param[in] v 被减数
|
||
*/
|
||
void minus_scalar(int r, int c, MatrixValueType v);
|
||
|
||
// 矩阵运算
|
||
/**
|
||
* @brief 稀疏矩阵相加
|
||
*
|
||
* @param[in] b 被加稀疏矩阵
|
||
*/
|
||
void add_sparse_matrix(const spmat<MatrixValueType> &b);
|
||
|
||
/**
|
||
* @brief 稀疏矩阵相减
|
||
*
|
||
* @param[in] b 被减稀疏矩阵
|
||
*/
|
||
void minus_sparse_matrix(const spmat<MatrixValueType> &b);
|
||
|
||
/**
|
||
* @brief 稠密矩阵相加
|
||
*
|
||
* @param[in] arr2d 引用的形式返回的相加后矩阵
|
||
*/
|
||
void add_matrix(matrix<MatrixValueType> &arr2d) const;
|
||
|
||
/**
|
||
* @brief 稠密矩阵相减
|
||
*
|
||
* @param arr2d 引用的形式返回的相减后矩阵
|
||
* @param[in] preposition 稠密矩阵为减数
|
||
*/
|
||
void minus_matrix(matrix<MatrixValueType> &arr2d, bool preposition = false) const;
|
||
|
||
/**
|
||
* @brief 映射某行或列的数据到另一行或列
|
||
*
|
||
* @param[in] from_index 来源行或列的索引
|
||
* @param[in] to_index 目标行或列的索引
|
||
* @param[in] multipler 倍数
|
||
* @param[in] line_mode 是否为行优先
|
||
* @param[in] val_mode 返回数组的数值操作方式
|
||
*/
|
||
void map(int from_index, int to_index, MatrixValueType multipler,
|
||
matrix_order_e line_mode = RowMajor, value_operator_e val_mode = AppendVal);
|
||
|
||
/**
|
||
* @brief 计算转置矩阵
|
||
*
|
||
* @param ret_mat 转置后矩阵
|
||
*/
|
||
void transpose(spmat<MatrixValueType> &ret_mat) const;
|
||
|
||
/**
|
||
* @brief 单个元素位置的乘法
|
||
*
|
||
* @param[in] r 行索引
|
||
* @param[in] c 列索引
|
||
* @param[in] v 乘数
|
||
*/
|
||
void multiply_scalar(int r, int c, MatrixValueType v);
|
||
|
||
/**
|
||
* @brief 稀疏矩阵的乘法
|
||
*
|
||
* @param[in] b 待乘稀疏矩阵
|
||
* @param ret_mat 乘积稀疏矩阵
|
||
*/
|
||
void multiply_sparse_matrix(const spmat<MatrixValueType> &b, spmat<MatrixValueType> &ret_mat) const;
|
||
|
||
/**
|
||
* @brief 矩阵乘法的母函数
|
||
*
|
||
* @param arr2d 待乘二维数组的指针（需提前开辟好内存空间）
|
||
* @param[in] arr2d_rsize 待乘二维数组的行数
|
||
* @param[in] arr2d_csize 待乘二维数组的列数
|
||
* @param ret_arr2d 乘积二维数组的指针（需提前开辟好内存空间）
|
||
* @param[in] ret_arr2d_rsize 乘积二维数组的行数
|
||
* @param[in] ret_arr2d_csize 乘积二维数组的列数
|
||
*/
|
||
void multiply_matrix(MatrixValueType **arr2d, int arr2d_rsize, int arr2d_csize,
|
||
MatrixValueType **ret_arr2d, int ret_arr2d_rsize, int ret_arr2d_csize) const;
|
||
|
||
/**
|
||
* @brief 矩阵乘法的母函数
|
||
*
|
||
* @note 稀疏矩阵与稠密矩阵的积也会是一个稠密的矩阵
|
||
*
|
||
* @param arr2d 待乘二维数组（需提前开辟好内存空间）
|
||
* @param ret_arr2d 乘积二维数组的指针
|
||
*/
|
||
void multiply_matrix(const matrix<MatrixValueType> &arr2d, matrix<MatrixValueType> &ret_arr2d) const;
|
||
|
||
/**
|
||
* @brief 矩阵与向量相乘
|
||
*
|
||
* @param[in] arr 输入数组指针
|
||
* @param[in] arr_size 输入数组大小
|
||
* @param ret_arr 输出数组指针
|
||
* @param[in] ret_size 输出数组大小
|
||
* @param[in] layout 是否使用稀疏矩阵的转置
|
||
*/
|
||
void multiply_vector(const MatrixValueType *arr, int arr_size, MatrixValueType *ret_arr,
|
||
int ret_size, matrix_layout_e layout = NoTrans) const;
|
||
|
||
/**
|
||
* @brief 稀疏矩阵与向量相乘
|
||
*
|
||
* 此函数会计算稀疏矩阵与一个向量的乘积向量（一维数组）
|
||
*
|
||
* @param[in] arr 待乘的向量数组
|
||
* @param ret_arr 乘积向量的数组
|
||
* @param[in] layout 是否使用稀疏矩阵的转置
|
||
*
|
||
*/
|
||
void multiply_vector(const array<MatrixValueType> &arr, array<MatrixValueType> &ret_arr,
|
||
matrix_layout_e layout = NoTrans) const;
|
||
|
||
/**
|
||
* @brief 稀疏矩阵的某行或列与向量相乘
|
||
*
|
||
* 此函数会计算稀疏矩阵的行或列与一个向量的乘积向量（一维数组）
|
||
*
|
||
* @param[in] arr 待乘的向量数组
|
||
* @param index 与向量相乘的稀疏矩阵的行或列索引
|
||
* @param[in] layout 是否使用稀疏矩阵的转置
|
||
*
|
||
*/
|
||
double multiply_vector(const array<MatrixValueType> &arr, int index,
|
||
matrix_layout_e layout = NoTrans) const;
|
||
|
||
/**
|
||
* @brief 按元素与向量相乘
|
||
*
|
||
* @param[in] oper_id 稀疏矩阵的行（非转置）或者列（转置）序号
|
||
* @param[in] arr 输入数组指针
|
||
* @param[in] arr_size 输入数组大小
|
||
* @param ret_arr 输出数组指针
|
||
* @param[in] layout 是否使用稀疏矩阵的转置
|
||
*/
|
||
void multiply_vector(int oper_id, const MatrixValueType *arr, int arr_size,
|
||
MatrixValueType *ret_arr, matrix_layout_e layout = NoTrans) const;
|
||
|
||
/**
|
||
* @brief 与对角矩阵相乘
|
||
*
|
||
* @param arr 输入的对角矩阵（保存为一个一维数组）
|
||
* @param ret_mat 乘积矩阵
|
||
* @param prepositioned 对角矩阵是否在稀疏矩阵前，为真表示对角矩阵乘稀疏矩阵，为假则表示稀疏矩阵乘对角矩阵
|
||
*/
|
||
void multiply_diagonal_matrix(const array<MatrixValueType> &arr, spmat<MatrixValueType> &ret_mat, bool prepositioned = false) const;
|
||
|
||
/**
|
||
* @brief 返回制定行或列的布尔数组，reverse为真则布尔数组在无元素的位置为真
|
||
*
|
||
* @param[in] oper_id The operator identifier
|
||
* @param ret_arr The ret arr
|
||
* @param[in] ele_mode The ele mode
|
||
* @param[in] reverse The reverse
|
||
*/
|
||
void boolean(int oper_id, array<bool> &ret_arr, matrix_order_e ele_mode = RowMajor,
|
||
bool reverse = false) const;
|
||
|
||
/**
|
||
* @brief 按行或列对矩阵进行归一化
|
||
*
|
||
* @param[in] oper_id 行或列的索引
|
||
* @param[in] norm 归一化后的大小
|
||
* @param[in] line_mode 归一化的方向
|
||
*/
|
||
void normalize(int oper_id, MatrixValueType norm, matrix_order_e line_mode = RowMajor);
|
||
|
||
/**
|
||
* @brief 对矩阵元素进行缩放
|
||
*
|
||
* @param f 缩放倍数
|
||
* @param idx 行或列的索引
|
||
* @param ele_mode 行或列模式
|
||
*/
|
||
void scale(MatrixValueType f, int idx, matrix_order_e ele_mode = RowMajor);
|
||
|
||
/**
|
||
* @brief 返回对角线元素，只对方阵有效
|
||
*
|
||
* @param dia 对角线元素
|
||
*/
|
||
void get_diagonal(array<MatrixValueType> &dia);
|
||
|
||
/**
|
||
* @brief 重载赋值操作符
|
||
*
|
||
* @param[in] b { parameter_description }
|
||
*
|
||
* @return The result of the assignment
|
||
*/
|
||
spmat<MatrixValueType>& operator= (const spmat<MatrixValueType> &b);
|
||
|
||
#ifdef GCTL_EIGEN
|
||
/**
|
||
* @brief 稀疏矩阵与Eigen向量相乘
|
||
*
|
||
* 此函数会计算稀疏矩阵与一个竖或者横向量的乘积向量（一维数组）。
|
||
* 与竖向量相乘为矩阵乘向量，与横向量相乘为向量乘矩阵
|
||
*
|
||
* @param vec 待乘的向量数组
|
||
* @param ret_vec 乘积向量的数组
|
||
* @param row_vector 待乘向量是否为横向量，默认为假（即默认为竖向量）
|
||
*/
|
||
void multiply_vector(const Eigen::VectorXd &vec, Eigen::VectorXd &ret_vec,
|
||
bool row_vector = false) const;
|
||
#endif // GCTL_EIGEN
|
||
|
||
private:
|
||
/**
|
||
* @brief 单个矩阵元素操作
|
||
*
|
||
* @param[in] r 行索引
|
||
* @param[in] c 列索引
|
||
* @param[in] v 操作值
|
||
* @param[in] mode 操作类型（默认为替换，可用的类型请参考枚举类型的声明）
|
||
*
|
||
* @return 若元素操作后为零返回真，否则返回假
|
||
*/
|
||
bool change_element(int r, int c, MatrixValueType v, value_operator_e mode = ReplaceVal);
|
||
|
||
/**
|
||
* @brief 返回行头指针
|
||
* @return 头指针
|
||
*/
|
||
mat_node<MatrixValueType> **row_head() const;
|
||
|
||
/**
|
||
* @brief 返回列头指针
|
||
* @return 头指针
|
||
*/
|
||
mat_node<MatrixValueType> **col_head() const;
|
||
};
|
||
|
||
|
||
template <typename MatrixValueType>
|
||
spmat<MatrixValueType>::spmat()
|
||
{
|
||
r_num = c_num = d_num = 0;
|
||
r_head = c_head = nullptr;
|
||
r_end = c_end = nullptr;
|
||
d_link = nullptr;
|
||
r_count = c_count = nullptr;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
spmat<MatrixValueType>::spmat(int in_rnum, int in_cnum, MatrixValueType null_val)
|
||
{
|
||
malloc(in_rnum, in_cnum, null_val);
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
spmat<MatrixValueType>::spmat(const spmat<MatrixValueType> &b)
|
||
{
|
||
malloc(b.row_size(), b.col_size(), b.zero_value());
|
||
|
||
mat_node<MatrixValueType> **br_head = b.row_head();
|
||
mat_node<MatrixValueType> *bp, *p, *q;
|
||
|
||
int i, j;
|
||
for (int r = 0; r < r_num; r++)
|
||
{
|
||
if (br_head[r] != nullptr)
|
||
{
|
||
bp = br_head[r]->next_right;
|
||
while (bp != nullptr)
|
||
{
|
||
mat_node<MatrixValueType>* new_node =
|
||
new mat_node<MatrixValueType>;
|
||
new_node->r_id = bp->r_id;
|
||
new_node->c_id = bp->c_id;
|
||
new_node->val = bp->val;
|
||
|
||
i = bp->r_id; j = bp->c_id;
|
||
// 我们首先处理行的指针
|
||
if (r_head[i] != nullptr) // 无需检查行尾指针 它们是同步变化的 在尾部插入
|
||
{
|
||
p = r_end[i];
|
||
p->next_right = new_node;
|
||
r_end[i] = new_node;
|
||
r_count[i]++; // 增加行计数
|
||
}
|
||
else
|
||
{
|
||
// 添加行首指针 注意行首指针并不等于第一个元素 而是行首的右向链为第一个元素
|
||
r_head[i] = new mat_node<MatrixValueType>;
|
||
r_head[i]->next_right = new_node;
|
||
// 添加行尾指针 行尾指针直接等于本行的最后一个元素
|
||
r_end[i] = new_node;
|
||
r_count[i]++; // 增加行计数
|
||
}
|
||
|
||
// 接着处理列的指针
|
||
if (c_head[j] != nullptr)
|
||
{
|
||
q = c_end[j];
|
||
q->next_down = new_node;
|
||
c_end[j] = new_node;
|
||
c_count[j]++;
|
||
}
|
||
else
|
||
{
|
||
c_head[j] = new mat_node<MatrixValueType>;
|
||
c_head[j]->next_down = new_node;
|
||
c_end[j] = new_node;
|
||
c_count[j]++;
|
||
}
|
||
|
||
bp = bp->next_right;
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
spmat<MatrixValueType>::~spmat()
|
||
{
|
||
clear();
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::malloc(int in_rnum, int in_cnum, MatrixValueType null_val)
|
||
{
|
||
if (in_rnum <= 0 || in_cnum <= 0)
|
||
{
|
||
throw std::runtime_error("Invalid matrix size. From void spmat<T>::malloc(...)");
|
||
}
|
||
|
||
zero_val = null_val;
|
||
d_num = 0;
|
||
d_link = nullptr;
|
||
r_num = in_rnum; c_num = in_cnum;
|
||
r_head = new mat_node<MatrixValueType>* [in_rnum];
|
||
c_head = new mat_node<MatrixValueType>* [in_cnum];
|
||
r_end = new mat_node<MatrixValueType>* [in_rnum];
|
||
c_end = new mat_node<MatrixValueType>* [in_cnum];
|
||
r_count = new int [in_rnum];
|
||
c_count = new int [in_cnum];
|
||
|
||
for (int i = 0; i < r_num; i++)
|
||
{
|
||
r_head[i] = nullptr; r_end[i] = nullptr; r_count[i] = 0;
|
||
}
|
||
|
||
for (int i = 0; i < c_num; i++)
|
||
{
|
||
c_head[i] = nullptr; c_end[i] = nullptr; c_count[i] = 0;
|
||
}
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::clear(int index, matrix_order_e mode)
|
||
{
|
||
int r, c;
|
||
mat_node<MatrixValueType> *p, *q, *tmp_node;
|
||
|
||
if (index < 0)
|
||
{
|
||
if (r_head != nullptr && c_head != nullptr)
|
||
{
|
||
for (int i = 0; i < r_num; i++)
|
||
{
|
||
p = r_head[i];
|
||
while (p != nullptr)
|
||
{
|
||
q = p; p = p->next_right;
|
||
delete q; q = nullptr;
|
||
}
|
||
}
|
||
|
||
delete[] r_head; r_head = nullptr;
|
||
delete[] c_head; c_head = nullptr;
|
||
delete[] r_end; r_end = nullptr;
|
||
delete[] c_end; c_end = nullptr;
|
||
delete[] r_count; r_count = nullptr;
|
||
delete[] c_count; c_count = nullptr;
|
||
}
|
||
|
||
if (d_num)
|
||
{
|
||
delete[] d_link; d_link = nullptr;
|
||
d_num = 0;
|
||
}
|
||
|
||
r_num = c_num = 0;
|
||
}
|
||
else if (mode == RowMajor)
|
||
{
|
||
if (index >= r_num)
|
||
{
|
||
throw runtime_error("[gctl::spmat<T>::clear] Invalid index.");
|
||
}
|
||
|
||
if (r_head != nullptr && r_head[index] != nullptr)
|
||
{
|
||
p = r_head[index]->next_right;
|
||
while (p != nullptr)
|
||
{
|
||
tmp_node = p;
|
||
p = p->next_right;
|
||
|
||
c = tmp_node->c_id;
|
||
q = c_head[c];
|
||
|
||
while (index != q->next_down->r_id)
|
||
{
|
||
q = q->next_down;
|
||
}
|
||
|
||
if (q == c_head[c] && q->next_down->next_down == nullptr)
|
||
{
|
||
delete c_head[c];
|
||
c_head[c] = nullptr;
|
||
c_end[c] = nullptr;
|
||
}
|
||
else if (q->next_down == c_end[c])
|
||
{
|
||
q->next_down = nullptr;
|
||
c_end[c] = q;
|
||
}
|
||
else
|
||
{
|
||
q->next_down = q->next_down->next_down;
|
||
}
|
||
|
||
c_count[c]--;
|
||
delete tmp_node; tmp_node = nullptr;
|
||
}
|
||
|
||
r_head[index] = nullptr;
|
||
r_end[index] = nullptr;
|
||
r_count[index] = 0;
|
||
|
||
if (d_num)
|
||
{
|
||
delete[] d_link; d_link = nullptr;
|
||
d_num = 0;
|
||
}
|
||
}
|
||
}
|
||
else
|
||
{
|
||
if (index >= c_num)
|
||
{
|
||
throw runtime_error("[gctl::spmat<T>::clear] Invalid index.");
|
||
}
|
||
|
||
if (c_head != nullptr && c_head[index] != nullptr)
|
||
{
|
||
p = c_head[index]->next_down;
|
||
while (p != nullptr)
|
||
{
|
||
tmp_node = p;
|
||
p = p->next_down;
|
||
|
||
r = tmp_node->r_id;
|
||
q = r_head[r];
|
||
|
||
while (index != q->next_right->c_id)
|
||
{
|
||
q = q->next_right;
|
||
}
|
||
|
||
if (q == r_head[r] && q->next_right->next_right == nullptr)
|
||
{
|
||
delete r_head[r];
|
||
r_head[r] = nullptr;
|
||
r_end[r] = nullptr;
|
||
}
|
||
else if (q->next_right == r_end[r])
|
||
{
|
||
q->next_right = nullptr;
|
||
r_end[r] = q;
|
||
}
|
||
else
|
||
{
|
||
q->next_right = q->next_right->next_right;
|
||
}
|
||
|
||
r_count[r]--;
|
||
delete tmp_node; tmp_node = nullptr;
|
||
}
|
||
|
||
c_head[index] = nullptr;
|
||
c_end[index] = nullptr;
|
||
c_count[index] = 0;
|
||
|
||
if (d_num)
|
||
{
|
||
delete[] d_link; d_link = nullptr;
|
||
d_num = 0;
|
||
}
|
||
}
|
||
}
|
||
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::reset()
|
||
{
|
||
mat_node<MatrixValueType> *p, *q;
|
||
for (int i = 0; i < r_num; i++)
|
||
{
|
||
p = r_head[i];
|
||
while (p != nullptr)
|
||
{
|
||
q = p; p = p->next_right;
|
||
delete q; q = nullptr;
|
||
}
|
||
}
|
||
|
||
for (int i = 0; i < r_num; i++)
|
||
{
|
||
r_head[i] = nullptr; r_end[i] = nullptr; r_count[i] = 0;
|
||
}
|
||
|
||
for (int i = 0; i < c_num; i++)
|
||
{
|
||
c_head[i] = nullptr; c_end[i] = nullptr; c_count[i] = 0;
|
||
}
|
||
|
||
if (d_num)
|
||
{
|
||
delete[] d_link; d_link = nullptr;
|
||
d_num = 0;
|
||
}
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::purge()
|
||
{
|
||
mat_node<MatrixValueType> *p, *q, *tmp_node;
|
||
for (int i = 0; i < r_num; i++)
|
||
{
|
||
if (r_head[i] != nullptr)
|
||
{
|
||
p = r_head[i];
|
||
while (p->next_right != nullptr)
|
||
{
|
||
if (p->next_right->val == zero_val)
|
||
{
|
||
// 更改行与列的计数
|
||
r_count[i]--;
|
||
c_count[p->next_right->c_id]--;
|
||
|
||
q = c_head[p->next_right->c_id];
|
||
while (i != q->next_down->r_id)
|
||
{
|
||
q = q->next_down;
|
||
}
|
||
|
||
if (q == c_head[p->next_right->c_id] &&
|
||
q->next_down->next_down == nullptr) // 在列首且仅有列首
|
||
{
|
||
delete c_head[p->next_right->c_id];
|
||
c_head[p->next_right->c_id] = nullptr;
|
||
c_end[p->next_right->c_id] = nullptr;
|
||
}
|
||
else if (q->next_down == c_end[p->next_right->c_id]) // 在列尾
|
||
{
|
||
q->next_down = nullptr;
|
||
c_end[p->next_right->c_id] = q;
|
||
}
|
||
else
|
||
{
|
||
// 此处不删除元素
|
||
q->next_down = q->next_down->next_down;
|
||
}
|
||
|
||
if (p == r_head[i] &&
|
||
p->next_right->next_right == nullptr)
|
||
{
|
||
delete p->next_right;
|
||
p->next_right = nullptr;
|
||
delete r_head[i];
|
||
r_head[i] = nullptr;
|
||
r_end[i] = nullptr;
|
||
break;
|
||
}
|
||
else if (p->next_right == r_end[i])
|
||
{
|
||
tmp_node = p->next_right;
|
||
p->next_right = nullptr;
|
||
delete tmp_node;
|
||
tmp_node = nullptr;
|
||
r_end[i] = p;
|
||
}
|
||
else
|
||
{
|
||
tmp_node = p->next_right;
|
||
p->next_right = p->next_right->next_right;
|
||
delete tmp_node;
|
||
tmp_node = nullptr;
|
||
}
|
||
}
|
||
else p = p->next_right;
|
||
}
|
||
}
|
||
}
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::set_all(MatrixValueType in_val)
|
||
{
|
||
if (in_val == zero_val) // reset the whole matrix to the initial state
|
||
{
|
||
reset();
|
||
return;
|
||
}
|
||
|
||
mat_node<MatrixValueType> *p;
|
||
for (int i = 0; i < r_num; i++)
|
||
{
|
||
if (r_head[i] != nullptr)
|
||
{
|
||
p = r_head[i]->next_right;
|
||
while (p != nullptr)
|
||
{
|
||
p->val = in_val;
|
||
p = p->next_right;
|
||
}
|
||
}
|
||
}
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::set_triplts(const std::vector<mat_node<MatrixValueType> > &triplts, value_operator_e val_mode)
|
||
{
|
||
// 如果存在d_link则只需要更改数值即可
|
||
if (d_num)
|
||
{
|
||
int c = 0;
|
||
int old_r = -1, old_c = -1;
|
||
for (size_t t = 0; t < triplts.size(); t++)
|
||
{
|
||
if(c > d_num) throw std::runtime_error("[gctl::spmat<T>::set_triplts] Invalid triplts' size.");
|
||
|
||
if (triplts[t].r_id == old_r && triplts[t].c_id == old_c) // 重复的元素位置
|
||
{
|
||
c--; // 索引减一 （默认会往前挪一位 所以要往回挪一位）
|
||
if (val_mode == AppendVal) {d_link[c]->val += triplts[t].val; c++;} // 索引加一
|
||
if (val_mode == ReplaceVal) {d_link[c]->val = triplts[t].val; c++;} // 索引加一
|
||
}
|
||
else // 新的元素位置
|
||
{
|
||
old_r = triplts[t].r_id;
|
||
old_c = triplts[t].c_id;
|
||
d_link[c]->val = triplts[t].val;
|
||
c++; // 索引加一
|
||
}
|
||
}
|
||
|
||
return;
|
||
}
|
||
|
||
// 如果矩阵已经初始化 重置稀疏矩阵
|
||
if (!empty()) reset();
|
||
else throw std::runtime_error("The sparse matrix is not initialized. From gctl::spmat<T>::set_triplts(...)");
|
||
|
||
int old_i = -1, old_j = -1;
|
||
int i, j;
|
||
for (size_t t = 0; t < triplts.size(); t++)
|
||
{
|
||
i = triplts[t].r_id;
|
||
j = triplts[t].c_id;
|
||
|
||
if (i < 0 || j < 0 || i >= r_num || j >= c_num)
|
||
{
|
||
throw std::runtime_error("Invalid triplt's index. From gctl::spmat<T>::set_triplts(...)");
|
||
}
|
||
|
||
if (old_i >= i && old_j > j) // 按行按列顺序添加
|
||
{
|
||
throw std::runtime_error("The input triplts are not properly arranged. From gctl::spmat<T>::set_triplts(...)");
|
||
}
|
||
|
||
if (old_i == i && old_j == j) // 重复输入 只需要改变值就行了
|
||
{
|
||
if (val_mode == ReplaceVal) r_end[i]->val = triplts[t].val;
|
||
if (val_mode == AppendVal) r_end[i]->val += triplts[t].val;
|
||
continue;
|
||
}
|
||
|
||
// 新元素
|
||
mat_node<MatrixValueType> *new_node = new mat_node<MatrixValueType>; // 创建新节点 这里必定是需要新建节点的
|
||
|
||
new_node->r_id = triplts[t].r_id;
|
||
new_node->c_id = triplts[t].c_id;
|
||
new_node->val = triplts[t].val;
|
||
|
||
old_i = i;
|
||
old_j = j;
|
||
|
||
// 我们首先处理行的指针
|
||
if (r_head[i] != nullptr) // 无需检查行尾指针 它们是同步变化的
|
||
{
|
||
r_end[i]->next_right = new_node;
|
||
r_end[i] = new_node;
|
||
r_count[i]++;
|
||
}
|
||
else
|
||
{
|
||
// 添加行首指针 注意行首指针并不等于第一个元素 而是行首的右向链为第一个元素
|
||
r_head[i] = new mat_node<MatrixValueType>;
|
||
r_head[i]->next_right = new_node;
|
||
// 添加行尾指针 行尾指针直接等于本行的最后一个元素
|
||
r_end[i] = new_node;
|
||
r_count[i]++;
|
||
}
|
||
|
||
// 接着处理列的指针
|
||
if (c_head[j] != nullptr)
|
||
{
|
||
c_end[j]->next_down = new_node;
|
||
c_end[j] = new_node;
|
||
c_count[j]++;
|
||
}
|
||
else
|
||
{
|
||
c_head[j] = new mat_node<MatrixValueType>;
|
||
c_head[j]->next_down = new_node;
|
||
c_end[j] = new_node;
|
||
c_count[j]++;
|
||
}
|
||
}
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::set_coo(const array<size_t> &row_idx, const array<size_t> &col_idx, const array<MatrixValueType> &data)
|
||
{
|
||
if (row_idx.size() != col_idx.size() || col_idx.size() != data.size() || data.empty())
|
||
{
|
||
throw std::runtime_error("Invalid input COO sparse matrix. From gctl::spmat<T>::set_coo(...)");
|
||
}
|
||
|
||
// 如果矩阵已经初始化 重置稀疏矩阵
|
||
if (!empty()) reset();
|
||
else throw std::runtime_error("The sparse matrix is not initialized. From gctl::spmat<T>::set_coo(...)");
|
||
|
||
size_t i, j;
|
||
for (size_t t = 0; t < data.size(); t++)
|
||
{
|
||
if (row_idx[t] >= r_num || col_idx[t] >= c_num)
|
||
{
|
||
throw std::runtime_error("Invalid COO index. From gctl::spmat<T>::set_coo(...)");
|
||
}
|
||
|
||
// 新元素
|
||
mat_node<MatrixValueType> *new_node = new mat_node<MatrixValueType>; // 创建新节点 这里必定是需要新建节点的
|
||
|
||
i = row_idx[t];
|
||
j = col_idx[t];
|
||
new_node->r_id = row_idx[t];
|
||
new_node->c_id = col_idx[t];
|
||
new_node->val = data[t];
|
||
|
||
// 我们首先处理行的指针
|
||
if (r_head[i] != nullptr) // 无需检查行尾指针 它们是同步变化的
|
||
{
|
||
r_end[i]->next_right = new_node;
|
||
r_end[i] = new_node;
|
||
r_count[i]++;
|
||
}
|
||
else
|
||
{
|
||
// 添加行首指针 注意行首指针并不等于第一个元素 而是行首的右向链为第一个元素
|
||
r_head[i] = new mat_node<MatrixValueType>;
|
||
r_head[i]->next_right = new_node;
|
||
// 添加行尾指针 行尾指针直接等于本行的最后一个元素
|
||
r_end[i] = new_node;
|
||
r_count[i]++;
|
||
}
|
||
|
||
// 接着处理列的指针
|
||
if (c_head[j] != nullptr)
|
||
{
|
||
c_end[j]->next_down = new_node;
|
||
c_end[j] = new_node;
|
||
c_count[j]++;
|
||
}
|
||
else
|
||
{
|
||
c_head[j] = new mat_node<MatrixValueType>;
|
||
c_head[j]->next_down = new_node;
|
||
c_end[j] = new_node;
|
||
c_count[j]++;
|
||
}
|
||
}
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::set_csr(const array<size_t> &row_off, const array<size_t> &col_idx, const array<MatrixValueType> &data)
|
||
{
|
||
if (row_off.size() < 2 || row_off[0] != 0 || row_off.back() != col_idx.size() || col_idx.size() != data.size() || data.empty())
|
||
{
|
||
throw std::runtime_error("Invalid input CSR sparse matrix. From gctl::spmat<T>::set_csr(...)");
|
||
}
|
||
|
||
// 如果矩阵已经初始化 重置稀疏矩阵
|
||
if (!empty()) reset();
|
||
else throw std::runtime_error("The sparse matrix is not initialized. From gctl::spmat<T>::set_csr(...)");
|
||
|
||
size_t i, j;
|
||
for (size_t r = 0; r < row_off.size() - 1; r++)
|
||
{
|
||
for (size_t t = row_off[r]; t < row_off[r+1]; t++)
|
||
{
|
||
if (r >= r_num || col_idx[t] >= c_num)
|
||
{
|
||
throw std::runtime_error("Invalid CSR index. From gctl::spmat<T>::set_csr(...)");
|
||
}
|
||
|
||
// 新元素
|
||
mat_node<MatrixValueType> *new_node = new mat_node<MatrixValueType>; // 创建新节点 这里必定是需要新建节点的
|
||
|
||
i = r;
|
||
j = col_idx[t];
|
||
new_node->r_id = r;
|
||
new_node->c_id = col_idx[t];
|
||
new_node->val = data[t];
|
||
|
||
// 我们首先处理行的指针
|
||
if (r_head[i] != nullptr) // 无需检查行尾指针 它们是同步变化的
|
||
{
|
||
r_end[i]->next_right = new_node;
|
||
r_end[i] = new_node;
|
||
r_count[i]++;
|
||
}
|
||
else
|
||
{
|
||
// 添加行首指针 注意行首指针并不等于第一个元素 而是行首的右向链为第一个元素
|
||
r_head[i] = new mat_node<MatrixValueType>;
|
||
r_head[i]->next_right = new_node;
|
||
// 添加行尾指针 行尾指针直接等于本行的最后一个元素
|
||
r_end[i] = new_node;
|
||
r_count[i]++;
|
||
}
|
||
|
||
// 接着处理列的指针
|
||
if (c_head[j] != nullptr)
|
||
{
|
||
c_end[j]->next_down = new_node;
|
||
c_end[j] = new_node;
|
||
c_count[j]++;
|
||
}
|
||
else
|
||
{
|
||
c_head[j] = new mat_node<MatrixValueType>;
|
||
c_head[j]->next_down = new_node;
|
||
c_end[j] = new_node;
|
||
c_count[j]++;
|
||
}
|
||
}
|
||
}
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::set_2d_vector(const std::vector<std::vector<MatrixValueType> > &den_mat_ptr, MatrixValueType off_set)
|
||
{
|
||
if (den_mat_ptr.empty())
|
||
{
|
||
throw runtime_error("The input 2d vector is empty. From void gctl::spmat<T>::set_2d_vector(...)");
|
||
}
|
||
|
||
// 如果矩阵已经初始化 重置稀疏矩阵
|
||
if (!empty()) reset();
|
||
else throw std::runtime_error("The sparse matrix is not initialized. From gctl::spmat<T>::set_2d_vector(...)");
|
||
|
||
// 确定二维向量的最大列数
|
||
int max_cnum = 0;
|
||
for (int i = 0; i < den_mat_ptr.size(); i++)
|
||
{
|
||
max_cnum = GCTL_MAX(max_cnum, den_mat_ptr[i].size());
|
||
}
|
||
|
||
if (max_cnum > c_num || den_mat_ptr.size() > r_num)
|
||
{
|
||
throw runtime_error("The input 2d vector exceeds the maximal size. From void gctl::spmat<T>::set_2d_vector(...)");
|
||
}
|
||
|
||
mat_node<MatrixValueType> *p, *q;
|
||
for (int i = 0; i < den_mat_ptr.size(); i++)
|
||
{
|
||
for (int j = 0; j < den_mat_ptr[i].size(); j++)
|
||
{
|
||
if (den_mat_ptr[i][j] >= zero_val + off_set ||
|
||
den_mat_ptr[i][j] <= zero_val - off_set) // 输入值不为零则添加
|
||
{
|
||
mat_node<MatrixValueType> *new_node =
|
||
new mat_node<MatrixValueType>; // 创建新节点 这里必定是需要新建节点的
|
||
|
||
new_node->r_id = i;
|
||
new_node->c_id = j;
|
||
new_node->val = den_mat_ptr[i][j];
|
||
|
||
// 我们首先处理行的指针
|
||
if (r_head[i] != nullptr) // 无需检查行尾指针 它们是同步变化的
|
||
{
|
||
p = r_end[i];
|
||
p->next_right = new_node;
|
||
r_end[i] = new_node;
|
||
r_count[i]++;
|
||
}
|
||
else
|
||
{
|
||
// 添加行首指针 注意行首指针并不等于第一个元素 而是行首的右向链为第一个元素
|
||
r_head[i] = new mat_node<MatrixValueType>;
|
||
r_head[i]->next_right = new_node;
|
||
// 添加行尾指针 行尾指针直接等于本行的最后一个元素
|
||
r_end[i] = new_node;
|
||
r_count[i]++;
|
||
}
|
||
|
||
// 接着处理列的指针
|
||
if (c_head[j] != nullptr)
|
||
{
|
||
q = c_end[j];
|
||
q->next_down = new_node;
|
||
c_end[j] = new_node;
|
||
c_count[j]++;
|
||
}
|
||
else
|
||
{
|
||
c_head[j] = new mat_node<MatrixValueType>;
|
||
c_head[j]->next_down = new_node;
|
||
c_end[j] = new_node;
|
||
c_count[j]++;
|
||
}
|
||
}
|
||
}
|
||
}
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::set_2d_array(MatrixValueType **den_mat_ptr, MatrixValueType off_set)
|
||
{
|
||
if (den_mat_ptr == nullptr)
|
||
{
|
||
throw runtime_error("The input pointer is null. From void gctl::spmat<T>::set_2d_array(...)");
|
||
}
|
||
|
||
// 如果矩阵已经初始化 重置稀疏矩阵
|
||
if (!empty()) reset();
|
||
else throw std::runtime_error("The sparse matrix is not initialized. From gctl::spmat<T>::set_2d_array(...)");
|
||
|
||
mat_node<MatrixValueType> *p, *q;
|
||
for (int i = 0; i < r_num; i++)
|
||
{
|
||
for (int j = 0; j < c_num; j++)
|
||
{
|
||
if (den_mat_ptr[i][j] >= zero_val + off_set ||
|
||
den_mat_ptr[i][j] <= zero_val - off_set) // 输入值不为零则添加
|
||
{
|
||
mat_node<MatrixValueType> *new_node =
|
||
new mat_node<MatrixValueType>; // 创建新节点 这里必定是需要新建节点的
|
||
|
||
new_node->r_id = i;
|
||
new_node->c_id = j;
|
||
new_node->val = den_mat_ptr[i][j];
|
||
|
||
// 我们首先处理行的指针
|
||
if (r_head[i] != nullptr) // 无需检查行尾指针 它们是同步变化的
|
||
{
|
||
p = r_end[i];
|
||
p->next_right = new_node;
|
||
r_end[i] = new_node;
|
||
r_count[i]++;
|
||
}
|
||
else
|
||
{
|
||
// 添加行首指针 注意行首指针并不等于第一个元素 而是行首的右向链为第一个元素
|
||
r_head[i] = new mat_node<MatrixValueType>;
|
||
r_head[i]->next_right = new_node;
|
||
// 添加行尾指针 行尾指针直接等于本行的最后一个元素
|
||
r_end[i] = new_node;
|
||
r_count[i]++;
|
||
}
|
||
|
||
// 接着处理列的指针
|
||
if (c_head[j] != nullptr)
|
||
{
|
||
q = c_end[j];
|
||
q->next_down = new_node;
|
||
c_end[j] = new_node;
|
||
c_count[j]++;
|
||
}
|
||
else
|
||
{
|
||
c_head[j] = new mat_node<MatrixValueType>;
|
||
c_head[j]->next_down = new_node;
|
||
c_end[j] = new_node;
|
||
c_count[j]++;
|
||
}
|
||
}
|
||
}
|
||
}
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::set_matrix(const matrix<MatrixValueType> &den_mat_ptr, MatrixValueType off_set)
|
||
{
|
||
if (den_mat_ptr.row_size() > r_num || den_mat_ptr.col_size() > c_num)
|
||
{
|
||
throw runtime_error("The input matrix exceeds the maximal size. From void gctl::spmat<T>::set_matrix(...)");
|
||
}
|
||
|
||
set_2d_array(den_mat_ptr.get(), off_set);
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::init_pattern()
|
||
{
|
||
if (d_num)
|
||
{
|
||
delete[] d_link; d_link = nullptr;
|
||
d_num = 0;
|
||
}
|
||
|
||
d_num = ele_size();
|
||
if (d_num)
|
||
{
|
||
d_link = new mat_node<MatrixValueType>* [d_num];
|
||
for (int i = 0; i < d_num; i++)
|
||
{
|
||
d_link[i] = nullptr;
|
||
}
|
||
}
|
||
else return;
|
||
|
||
int c = 0;
|
||
mat_node<MatrixValueType>* p;
|
||
for (int i = 0; i < r_num; i++)
|
||
{
|
||
if (r_head[i] != nullptr)
|
||
{
|
||
p = r_head[i]->next_right;
|
||
while (p != nullptr)
|
||
{
|
||
d_link[c] = p;
|
||
p = p->next_right;
|
||
c++;
|
||
}
|
||
}
|
||
}
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
bool spmat<MatrixValueType>::empty() const
|
||
{
|
||
if (r_num == 0 && c_num == 0) return true;
|
||
else return false;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
int spmat<MatrixValueType>::row_size() const
|
||
{
|
||
return r_num;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
int spmat<MatrixValueType>::col_size() const
|
||
{
|
||
return c_num;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
int spmat<MatrixValueType>::ele_size(int index, matrix_order_e mode) const
|
||
{
|
||
if (index == -1)
|
||
{
|
||
int sum = 0;
|
||
for (int i = 0; i < r_num; i++)
|
||
{
|
||
sum += r_count[i];
|
||
}
|
||
return sum;
|
||
}
|
||
else if (mode == RowMajor)
|
||
{
|
||
if (index < 0 || index >= r_num)
|
||
{
|
||
throw std::runtime_error("[gctl::spmat<T>::ele_size] Invalid index.");
|
||
}
|
||
|
||
return r_count[index];
|
||
}
|
||
else
|
||
{
|
||
if (index < 0 || index >= c_num)
|
||
{
|
||
throw std::runtime_error("[gctl::spmat<T>::ele_size] Invalid index.");
|
||
}
|
||
|
||
return c_count[index];
|
||
}
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
MatrixValueType spmat<MatrixValueType>::zero_value() const
|
||
{
|
||
return zero_val;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::show_matrix(char delimiter) const
|
||
{
|
||
int st;
|
||
mat_node<MatrixValueType> *p;
|
||
for (int i = 0; i < r_num; i++)
|
||
{
|
||
st = 0;
|
||
if (r_head[i] != nullptr)
|
||
{
|
||
p = r_head[i]->next_right;
|
||
while (p != nullptr)
|
||
{
|
||
for (int s = st; s < p->c_id; s++)
|
||
std::cout << zero_val << delimiter;
|
||
std::cout << p->val << delimiter;
|
||
|
||
st = p->c_id + 1;
|
||
p = p->next_right;
|
||
}
|
||
}
|
||
|
||
if (st <= c_num-1)
|
||
{
|
||
for (int s = st; s < c_num; s++)
|
||
std::cout << zero_val << delimiter;
|
||
}
|
||
std::cout << std::endl;
|
||
}
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::show_list(matrix_order_e mode, char delimiter, bool full) const
|
||
{
|
||
int st;
|
||
if (mode == RowMajor)
|
||
{
|
||
std::cout << "# row" << delimiter << "col" << delimiter << "value" << std::endl;
|
||
|
||
mat_node<MatrixValueType> *p;
|
||
if (full)
|
||
{
|
||
for (int i = 0; i < r_num; i++)
|
||
{
|
||
st = 0;
|
||
if (r_head[i] != nullptr)
|
||
{
|
||
p = r_head[i]->next_right;
|
||
while (p != nullptr)
|
||
{
|
||
for (int s = st; s < p->c_id; s++)
|
||
{
|
||
std::cout << i << delimiter << s << delimiter << zero_val << std::endl;
|
||
}
|
||
std::cout << i << delimiter << p->c_id << delimiter << p->val << std::endl;
|
||
|
||
st = p->c_id + 1;
|
||
p = p->next_right;
|
||
}
|
||
}
|
||
|
||
if (st <= c_num-1)
|
||
{
|
||
for (int s = st; s < c_num; s++)
|
||
std::cout << i << delimiter << s << delimiter << zero_val << std::endl;
|
||
}
|
||
}
|
||
return;
|
||
}
|
||
|
||
for (int i = 0; i < r_num; i++)
|
||
{
|
||
if (r_head[i] != nullptr)
|
||
{
|
||
p = r_head[i]->next_right;
|
||
while (p != nullptr)
|
||
{
|
||
// 输出(row, col, val)三元组
|
||
std::cout << p->r_id << delimiter << p->c_id << delimiter << p->val << std::endl;
|
||
p = p->next_right;
|
||
}
|
||
}
|
||
}
|
||
return;
|
||
}
|
||
|
||
std::cout << "# col" << delimiter << "row" << delimiter << "value" << std::endl;
|
||
|
||
mat_node<MatrixValueType> *q;
|
||
if (full)
|
||
{
|
||
for (int i = 0; i < c_num; i++)
|
||
{
|
||
st = 0;
|
||
if (c_head[i] != nullptr)
|
||
{
|
||
q = c_head[i]->next_down;
|
||
while (q != nullptr)
|
||
{
|
||
for (int s = st; s < q->r_id; s++)
|
||
{
|
||
std::cout << i << delimiter << s << delimiter << zero_val << std::endl;
|
||
}
|
||
std::cout << i << delimiter << q->r_id << delimiter << q->val << std::endl;
|
||
|
||
st = q->r_id + 1;
|
||
q = q->next_down;
|
||
}
|
||
}
|
||
|
||
if (st <= r_num-1)
|
||
{
|
||
for (int s = st; s < r_num; s++)
|
||
std::cout << i << delimiter << s << delimiter << zero_val << std::endl;
|
||
}
|
||
}
|
||
return;
|
||
}
|
||
|
||
for (int i = 0; i < c_num; i++)
|
||
{
|
||
if (c_head[i] != nullptr)
|
||
{
|
||
q = c_head[i]->next_down;
|
||
while (q != nullptr)
|
||
{
|
||
// 输出(row, col, val)三元组
|
||
std::cout << q->c_id << delimiter << q->r_id << delimiter << q->val << std::endl;
|
||
q = q->next_down;
|
||
}
|
||
}
|
||
}
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
MatrixValueType spmat<MatrixValueType>::at(size_t i_index, size_t j_index) const
|
||
{
|
||
#ifdef GCTL_CHECK_BOUNDER
|
||
if (i_index >= r_num || j_index >= c_num)
|
||
{
|
||
throw std::runtime_error("Invalid index. From gctl::spmat<T>::at(...)");
|
||
}
|
||
#endif // GCTL_CHECK_BOUNDER
|
||
|
||
mat_node<MatrixValueType> *p, *q;
|
||
if (r_head[i_index] != nullptr && c_head[j_index] != nullptr)
|
||
{
|
||
p = r_end[i_index];
|
||
if (j_index > p->c_id) return zero_val;
|
||
p = r_head[i_index]->next_right;
|
||
if (j_index < p->c_id) return zero_val;
|
||
|
||
q = c_end[j_index];
|
||
if (i_index > q->r_id) return zero_val;
|
||
q = c_head[j_index]->next_down;
|
||
if (i_index < q->r_id) return zero_val;
|
||
|
||
if (r_count[i_index] <= c_count[j_index])
|
||
{
|
||
while (p != nullptr)
|
||
{
|
||
if (p->c_id == j_index)
|
||
{
|
||
return p->val;
|
||
}
|
||
p = p->next_right;
|
||
}
|
||
return zero_val;
|
||
}
|
||
|
||
while (q != nullptr)
|
||
{
|
||
if (q->r_id == i_index)
|
||
{
|
||
return q->val;
|
||
}
|
||
q = q->next_down;
|
||
}
|
||
return zero_val;
|
||
}
|
||
return zero_val;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::export_dense(matrix<MatrixValueType> &ret_arr2d) const
|
||
{
|
||
ret_arr2d.resize(r_num, c_num);
|
||
|
||
mat_node<MatrixValueType> *p;
|
||
for (int i = 0; i < r_num; i++)
|
||
{
|
||
for (int j = 0; j < c_num; j++)
|
||
{
|
||
ret_arr2d[i][j] = zero_val;
|
||
}
|
||
|
||
if (r_head[i] != nullptr)
|
||
{
|
||
p = r_head[i]->next_right;
|
||
while (p != nullptr)
|
||
{
|
||
ret_arr2d[i][p->c_id] = p->val;
|
||
p = p->next_right;
|
||
}
|
||
}
|
||
}
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::export_dense(int out_id, array<MatrixValueType> &ret_arr,
|
||
double multipler, matrix_order_e line_mode, value_operator_e val_mode) const
|
||
{
|
||
if (line_mode == RowMajor)
|
||
{
|
||
if (out_id < 0 || out_id >= r_num)
|
||
{
|
||
throw std::runtime_error("[gctl::spmat<T>::export_dense] Invalid index.");
|
||
}
|
||
|
||
if (val_mode != AppendVal)
|
||
{
|
||
ret_arr.resize(c_num, zero_val);
|
||
}
|
||
else if (ret_arr.size() != c_num)
|
||
{
|
||
throw std::runtime_error("[gctl::spmat<T>::export_dense] Incompatible array size.");
|
||
}
|
||
|
||
mat_node<MatrixValueType> *p;
|
||
if (r_head[out_id] != nullptr)
|
||
{
|
||
p = r_head[out_id]->next_right;
|
||
while (p != nullptr)
|
||
{
|
||
ret_arr[p->c_id] += p->val * multipler;
|
||
p = p->next_right;
|
||
}
|
||
}
|
||
}
|
||
else
|
||
{
|
||
if (out_id < 0 || out_id >= c_num)
|
||
{
|
||
throw std::runtime_error("[gctl::spmat<T>::export_dense] Invalid index.");
|
||
}
|
||
|
||
if (val_mode != AppendVal)
|
||
{
|
||
ret_arr.resize(r_num, zero_val);
|
||
}
|
||
else if (ret_arr.size() != r_num)
|
||
{
|
||
throw std::runtime_error("[gctl::spmat<T>::export_dense] Incompatible array size.");
|
||
}
|
||
|
||
mat_node<MatrixValueType> *q;
|
||
if (c_head[out_id] != nullptr)
|
||
{
|
||
q = c_head[out_id]->next_down;
|
||
while (q != nullptr)
|
||
{
|
||
ret_arr[q->r_id] += q->val * multipler;
|
||
q = q->next_down;
|
||
}
|
||
}
|
||
}
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::export_coo(array<size_t> &row_idx, array<size_t> &col_idx, array<MatrixValueType> &data) const
|
||
{
|
||
if (empty()) throw std::runtime_error("The sparse matrix is not initialized. From gctl::spmat<T>::export_coo(...)");
|
||
|
||
size_t nz_num = ele_size();
|
||
row_idx.resize(nz_num);
|
||
col_idx.resize(nz_num);
|
||
data.resize(nz_num);
|
||
|
||
size_t c = 0;
|
||
mat_node<MatrixValueType> *p;
|
||
for (size_t i = 0; i < r_num; i++)
|
||
{
|
||
if (r_head[i] != nullptr)
|
||
{
|
||
p = r_head[i]->next_right;
|
||
while (p != nullptr)
|
||
{
|
||
row_idx[c] = i;
|
||
col_idx[c] = p->c_id;
|
||
data[c] = p->val;
|
||
p = p->next_right;
|
||
|
||
c++;
|
||
}
|
||
}
|
||
}
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::export_csr(array<size_t> &row_off, array<size_t> &col_idx, array<MatrixValueType> &data) const
|
||
{
|
||
if (empty()) throw std::runtime_error("The sparse matrix is not initialized. From gctl::spmat<T>::export_coo(...)");
|
||
|
||
size_t nz_num = ele_size();
|
||
col_idx.resize(nz_num);
|
||
data.resize(nz_num);
|
||
|
||
row_off.resize(r_num + 1);
|
||
row_off[0] = 0;
|
||
|
||
size_t r, c = 0;
|
||
mat_node<MatrixValueType> *p;
|
||
for (size_t i = 0; i < r_num; i++)
|
||
{
|
||
r = row_off[i];
|
||
if (r_head[i] != nullptr)
|
||
{
|
||
p = r_head[i]->next_right;
|
||
while (p != nullptr)
|
||
{
|
||
col_idx[c] = p->c_id;
|
||
data[c] = p->val;
|
||
p = p->next_right;
|
||
|
||
c++;
|
||
r++;
|
||
}
|
||
}
|
||
row_off[i+1] = r;
|
||
}
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::insert(int r, int c, MatrixValueType v)
|
||
{
|
||
change_element(r, c, v);
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::remove(int r, int c)
|
||
{
|
||
if (r >= r_num || r < 0 || c >= c_num || c < 0)
|
||
{
|
||
std::string err_str = "Invalid index. From gctl::spmat<T>::remove(...)";
|
||
throw runtime_error(err_str);
|
||
}
|
||
|
||
mat_node<MatrixValueType> *p, *q, *tmp_node;
|
||
p = r_head[r];
|
||
q = c_head[c];
|
||
|
||
if (p == nullptr || q == nullptr) // 元素不在序列中 直接返回
|
||
return;
|
||
|
||
while (p->next_right != nullptr)
|
||
{
|
||
if (c == p->next_right->c_id)
|
||
{
|
||
while (r != q->next_down->r_id)
|
||
{
|
||
q = q->next_down;
|
||
}
|
||
|
||
if (q == c_head[c] &&
|
||
q->next_down->next_down == nullptr)
|
||
{
|
||
delete c_head[c];
|
||
c_head[c] = nullptr;
|
||
c_end[c] = nullptr;
|
||
}
|
||
else if (q->next_down == c_end[c])
|
||
{
|
||
q->next_down = nullptr;
|
||
c_end[c] = q;
|
||
}
|
||
else
|
||
{
|
||
q->next_down = q->next_down->next_down;
|
||
}
|
||
|
||
if (p == r_head[r] &&
|
||
p->next_right->next_right == nullptr)
|
||
{
|
||
delete r_head[r]->next_right;
|
||
r_head[r]->next_right = nullptr;
|
||
delete r_head[r];
|
||
r_head[r] = nullptr;
|
||
r_end[r] = nullptr;
|
||
}
|
||
else if (p->next_right == r_end[r])
|
||
{
|
||
tmp_node = p->next_right;
|
||
p->next_right = nullptr;
|
||
delete tmp_node;
|
||
tmp_node = nullptr;
|
||
r_end[r] = p;
|
||
}
|
||
else
|
||
{
|
||
tmp_node = p->next_right;
|
||
p->next_right = p->next_right->next_right;
|
||
delete tmp_node;
|
||
tmp_node = nullptr;
|
||
}
|
||
|
||
// 更改行与列的计数
|
||
r_count[r]--;
|
||
c_count[c]--;
|
||
break;
|
||
}
|
||
p = p->next_right;
|
||
}
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
MatrixValueType spmat<MatrixValueType>::sum(int index, matrix_order_e mode)
|
||
{
|
||
MatrixValueType sum_val = 0;
|
||
if (index < 0)
|
||
{
|
||
mat_node<MatrixValueType> *p;
|
||
for (int i = 0; i < r_num; i++)
|
||
{
|
||
if (r_head[i] != nullptr)
|
||
{
|
||
p = r_head[i]->next_right;
|
||
while (p != nullptr)
|
||
{
|
||
sum_val = sum_val + p->val;
|
||
p = p->next_right;
|
||
}
|
||
}
|
||
}
|
||
}
|
||
else if (mode == RowMajor)
|
||
{
|
||
if (index < 0 || index >= r_num)
|
||
{
|
||
throw std::runtime_error("[gctl::spmat<T>::sum] Invalid index");
|
||
}
|
||
|
||
mat_node<MatrixValueType> *p;
|
||
if (r_head[index] != nullptr)
|
||
{
|
||
p = r_head[index]->next_right;
|
||
while (p != nullptr)
|
||
{
|
||
sum_val = sum_val + p->val;
|
||
p = p->next_right;
|
||
}
|
||
}
|
||
}
|
||
else
|
||
{
|
||
if (index < 0 || index >= c_num)
|
||
{
|
||
throw std::runtime_error("[gctl::spmat<T>::sum] Invalid index");
|
||
}
|
||
|
||
mat_node<MatrixValueType> *q;
|
||
if (c_head[index] != nullptr)
|
||
{
|
||
q = c_head[index]->next_down;
|
||
while (q != nullptr)
|
||
{
|
||
sum_val = sum_val + q->val;
|
||
q = q->next_down;
|
||
}
|
||
}
|
||
}
|
||
return sum_val;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
MatrixValueType spmat<MatrixValueType>::module(int index, matrix_order_e mode)
|
||
{
|
||
MatrixValueType sum_val = 0;
|
||
if (mode == RowMajor)
|
||
{
|
||
if (index < 0 || index >= r_num)
|
||
{
|
||
throw std::runtime_error("[gctl::spmat<T>::module] Invalid index");
|
||
}
|
||
|
||
mat_node<MatrixValueType> *p;
|
||
if (r_head[index] != nullptr)
|
||
{
|
||
p = r_head[index]->next_right;
|
||
while (p != nullptr)
|
||
{
|
||
sum_val = sum_val + pow(p->val, 2);
|
||
p = p->next_right;
|
||
}
|
||
|
||
return sqrt(sum_val);
|
||
}
|
||
}
|
||
else
|
||
{
|
||
if (index < 0 || index >= c_num)
|
||
{
|
||
throw std::runtime_error("[gctl::spmat<T>::module] Invalid index");
|
||
}
|
||
|
||
mat_node<MatrixValueType> *q;
|
||
if (c_head[index] != nullptr)
|
||
{
|
||
q = c_head[index]->next_down;
|
||
while (q != nullptr)
|
||
{
|
||
sum_val = sum_val + pow(q->val, 2);
|
||
q = q->next_down;
|
||
}
|
||
|
||
return sqrt(sum_val);
|
||
}
|
||
}
|
||
return 0;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::add_scalar(int r, int c, MatrixValueType v)
|
||
{
|
||
if (change_element(r, c, v, AppendVal))
|
||
remove(r, c);
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::minus_scalar(int r, int c, MatrixValueType v)
|
||
{
|
||
if (change_element(r, c, v, RemoveVal))
|
||
remove(r, c);
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::add_sparse_matrix(const spmat<MatrixValueType> &b)
|
||
{
|
||
if (b.row_size() > r_num || b.col_size() > c_num)
|
||
{
|
||
throw runtime_error("Invalid maximal matrix sizes. From gctl::spmat<T>::add_sparse_matrix(...)"); // 大小不匹配 不能操作
|
||
}
|
||
|
||
mat_node<MatrixValueType> **br_head = b.row_head();
|
||
|
||
bool re_flag = false;
|
||
mat_node<MatrixValueType> *pb; // 遍历b的指针
|
||
for (int i = 0; i < r_num; i++)
|
||
{
|
||
pb = br_head[i];
|
||
if (pb != nullptr)
|
||
{
|
||
pb = pb->next_right;
|
||
while (pb != nullptr)
|
||
{
|
||
// 调用逐个加法私有方法
|
||
if (change_element(pb->r_id, pb->c_id, pb->val, AppendVal))
|
||
{
|
||
re_flag = true;
|
||
}
|
||
pb = pb->next_right;
|
||
}
|
||
}
|
||
}
|
||
|
||
if (re_flag) purge();
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::minus_sparse_matrix(const spmat<MatrixValueType> &b)
|
||
{
|
||
if (b.row_size() > r_num || b.col_size() > c_num)
|
||
{
|
||
throw runtime_error("Invalid maximal matrix sizes. From gctl::spmat<T>::minus_sparse_matrix(...)"); // 大小不匹配 不能操作
|
||
}
|
||
|
||
mat_node<MatrixValueType> **br_head = b.row_head();
|
||
|
||
bool re_flag = false;
|
||
mat_node<MatrixValueType> *pb; // 遍历b的指针
|
||
for (int i = 0; i < r_num; i++)
|
||
{
|
||
pb = br_head[i];
|
||
if (pb != nullptr)
|
||
{
|
||
pb = pb->next_right;
|
||
while (pb != nullptr)
|
||
{
|
||
// 调用逐个减法私有方法
|
||
if (change_element(pb->r_id, pb->c_id, pb->val, RemoveVal))
|
||
{
|
||
re_flag = true;
|
||
}
|
||
pb = pb->next_right;
|
||
}
|
||
}
|
||
}
|
||
|
||
if (re_flag) purge();
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::add_matrix(matrix<MatrixValueType> &arr2d) const
|
||
{
|
||
if (r_num != arr2d.row_size() || c_num != arr2d.col_size())
|
||
{
|
||
throw std::runtime_error("Incompatible matrix sizes. From gctl::spmat<T>::add_matrix(...)");
|
||
}
|
||
|
||
mat_node<MatrixValueType> *p;
|
||
for (int i = 0; i < r_num; i++)
|
||
{
|
||
p = r_head[i];
|
||
|
||
if (p != nullptr)
|
||
{
|
||
p = p->next_right;
|
||
while (p != nullptr)
|
||
{
|
||
arr2d[i][p->c_id] += p->val;
|
||
p = p->next_right;
|
||
}
|
||
}
|
||
}
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::minus_matrix(matrix<MatrixValueType> &arr2d, bool preposition) const
|
||
{
|
||
if (r_num != arr2d.row_size() || c_num != arr2d.col_size())
|
||
{
|
||
throw std::runtime_error("Incompatible matrix sizes. From gctl::spmat<T>::minus_matrix(...)");
|
||
}
|
||
|
||
mat_node<MatrixValueType> *p;
|
||
for (int i = 0; i < r_num; i++)
|
||
{
|
||
p = r_head[i];
|
||
|
||
if (p != nullptr)
|
||
{
|
||
p = p->next_right;
|
||
while (p != nullptr)
|
||
{
|
||
if (preposition) arr2d[i][p->c_id] -= p->val;
|
||
else arr2d[i][p->c_id] = p->val - arr2d[i][p->c_id];
|
||
p = p->next_right;
|
||
}
|
||
}
|
||
}
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::map(int from_index, int to_index, MatrixValueType multipler,
|
||
matrix_order_e line_mode, value_operator_e val_mode)
|
||
{
|
||
std::string err_str;
|
||
if (line_mode == RowMajor)
|
||
{
|
||
if (from_index < 0 || from_index >= r_num ||
|
||
to_index < 0 || to_index >= r_num)
|
||
{
|
||
err_str = "Invalid index. From gctl::spmat<T>::map(...)";
|
||
throw runtime_error(err_str);
|
||
}
|
||
|
||
mat_node<MatrixValueType> *p;
|
||
if (r_head[from_index] != nullptr)
|
||
{
|
||
p = r_head[from_index]->next_right;
|
||
while (p != nullptr)
|
||
{
|
||
change_element(to_index, p->c_id, multipler * p->val, val_mode);
|
||
p = p->next_right;
|
||
}
|
||
}
|
||
}
|
||
else
|
||
{
|
||
if (from_index < 0 || from_index >= c_num ||
|
||
to_index < 0 || to_index >= c_num)
|
||
{
|
||
err_str = "Invalid index. From gctl::spmat<T>::map(...)";
|
||
throw runtime_error(err_str);
|
||
}
|
||
|
||
mat_node<MatrixValueType> *q;
|
||
if (c_head[from_index] != nullptr)
|
||
{
|
||
q = c_head[from_index]->next_down;
|
||
while (q != nullptr)
|
||
{
|
||
change_element(q->r_id, to_index, multipler * q->val, val_mode);
|
||
q = q->next_down;
|
||
}
|
||
}
|
||
}
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::transpose(spmat<MatrixValueType> &ret_mat) const
|
||
{
|
||
if (!ret_mat.empty()) ret_mat.clear();
|
||
ret_mat.malloc(c_num, r_num, zero_val);
|
||
|
||
mat_node<MatrixValueType> *p;
|
||
for (int i = 0; i < r_num; i++)
|
||
{
|
||
if (r_head[i] != nullptr)
|
||
{
|
||
p = r_head[i]->next_right;
|
||
while (p != nullptr)
|
||
{
|
||
ret_mat.insert(p->c_id, p->r_id, p->val);
|
||
p = p->next_right;
|
||
}
|
||
}
|
||
}
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::multiply_scalar(int r, int c, MatrixValueType v)
|
||
{
|
||
if (r >= r_num || r < 0 || c >= c_num || c < 0)
|
||
{
|
||
std::string err_str = "Invalid index. From gctl::spmat<T>::multiply(...)";
|
||
throw runtime_error(err_str);
|
||
}
|
||
|
||
if (v == zero_val)
|
||
return;
|
||
|
||
if (r_head != nullptr && r_head[r] != nullptr)
|
||
{
|
||
mat_node<MatrixValueType> *p = r_head[r]->next_right;
|
||
while (p != nullptr)
|
||
{
|
||
if (p->c_id == c)
|
||
{
|
||
p->val *= v;
|
||
break;
|
||
}
|
||
p = p->next_right;
|
||
}
|
||
}
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::multiply_sparse_matrix(const spmat<MatrixValueType> &b,
|
||
spmat<MatrixValueType> &ret_mat) const
|
||
{
|
||
int rb = b.row_size();
|
||
int cb = b.col_size();
|
||
if (c_num != rb)
|
||
{
|
||
throw runtime_error("Incompatible matrix sizes. From spmat<T>::multiply(...)"); // 大小不匹配 不能操作
|
||
}
|
||
|
||
if (!ret_mat.empty()) ret_mat.clear();
|
||
ret_mat.malloc(r_num, cb, zero_val);
|
||
|
||
mat_node<MatrixValueType> **bc_head = b.col_head();
|
||
mat_node<MatrixValueType> *p, *q;
|
||
MatrixValueType compute = zero_val;
|
||
bool has_value = false;
|
||
|
||
for (int i = 0; i < r_num; i++)
|
||
{
|
||
for (int j = 0; j < cb; j++)
|
||
{
|
||
has_value = false;
|
||
p = r_head[i];
|
||
q = bc_head[j];
|
||
compute = zero_val;
|
||
|
||
if (p != nullptr && q != nullptr)
|
||
{
|
||
p = p->next_right;
|
||
q = q->next_down;
|
||
while (p != nullptr && q != nullptr)
|
||
{
|
||
if (p->c_id < q->r_id) // q在p后面
|
||
{
|
||
p = p->next_right; // p往后赶
|
||
}
|
||
else if (p->c_id > q->r_id) // p在q后面
|
||
{
|
||
q = q->next_down; // q往后赶
|
||
}
|
||
else // p->col == q->row
|
||
{
|
||
has_value = true;
|
||
compute += p->val * q->val; // ans[i,j] += A[i,p] * B[p,j];
|
||
p = p->next_right; // p,q一起往后赶
|
||
q = q->next_down; // p,q一起往后赶
|
||
}
|
||
}
|
||
}
|
||
|
||
if (has_value && compute != zero_val) // 如果有非零值
|
||
{
|
||
ret_mat.insert(i, j, compute);
|
||
}
|
||
}
|
||
}
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::multiply_matrix(MatrixValueType **arr2d, int arr2d_rsize,
|
||
int arr2d_csize, MatrixValueType **ret_arr2d, int ret_arr2d_rsize, int ret_arr2d_csize) const
|
||
{
|
||
if (c_num != arr2d_rsize || ret_arr2d_csize != arr2d_csize || ret_arr2d_rsize != r_num)
|
||
{
|
||
std::string err_str = "Incompatible matrix sizes. From gctl::spmat<T>::multiply(...)";
|
||
throw runtime_error(err_str); // 大小不匹配 不能操作
|
||
}
|
||
|
||
mat_node<MatrixValueType> *p;
|
||
MatrixValueType sum;
|
||
|
||
for (int i = 0; i < r_num; i++)
|
||
{
|
||
for (int j = 0; j < arr2d_csize; j++)
|
||
{
|
||
p = r_head[i];
|
||
sum = zero_val;
|
||
|
||
if (p != nullptr)
|
||
{
|
||
p = p->next_right;
|
||
while (p != nullptr)
|
||
{
|
||
sum += p->val * arr2d[p->c_id][j];
|
||
p = p->next_right;
|
||
}
|
||
}
|
||
|
||
ret_arr2d[i][j] = sum;
|
||
}
|
||
}
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::multiply_matrix(const matrix<MatrixValueType> &arr2d,
|
||
matrix<MatrixValueType> &ret_arr2d) const
|
||
{
|
||
if (c_num != arr2d.row_size())
|
||
{
|
||
std::string err_str = "Incompatible matrix sizes. From gctl::spmat<T>::multiply(...)";
|
||
throw runtime_error(err_str); // 大小不匹配 不能操作
|
||
}
|
||
|
||
ret_arr2d.resize(r_num, arr2d.col_size());
|
||
|
||
mat_node<MatrixValueType> *p;
|
||
MatrixValueType sum;
|
||
|
||
for (int i = 0; i < r_num; i++)
|
||
{
|
||
for (int j = 0; j < arr2d.col_size(); j++)
|
||
{
|
||
p = r_head[i];
|
||
sum = zero_val;
|
||
|
||
if (p != nullptr)
|
||
{
|
||
p = p->next_right;
|
||
while (p != nullptr)
|
||
{
|
||
sum += p->val * arr2d[p->c_id][j];
|
||
p = p->next_right;
|
||
}
|
||
}
|
||
|
||
ret_arr2d[i][j] = sum;
|
||
}
|
||
}
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::multiply_vector(const MatrixValueType *arr, int arr_size,
|
||
MatrixValueType *ret_arr, int ret_size, matrix_layout_e layout) const
|
||
{
|
||
size_t i;
|
||
MatrixValueType sum;
|
||
// 待乘向量为横向量 按列与矩阵做乘法与加法
|
||
if (layout == Trans)
|
||
{
|
||
if (arr_size != r_num || ret_size != c_num) // 向量大小必须与矩阵行数相等
|
||
{
|
||
std::string err_str = "Incompatible matrix sizes. From gctl::spmat<T>::multiply(...)";
|
||
throw runtime_error(err_str); // 大小不匹配 不能操作
|
||
}
|
||
|
||
mat_node<MatrixValueType> *q;
|
||
#pragma omp parallel for private (i, sum, q) schedule(guided)
|
||
for (i = 0; i < c_num; i++)
|
||
{
|
||
sum = zero_val;
|
||
|
||
q = c_head[i];
|
||
if (q != nullptr)
|
||
{
|
||
q = q->next_down;
|
||
while (q != nullptr)
|
||
{
|
||
sum += q->val * arr[q->r_id];
|
||
q = q->next_down;
|
||
}
|
||
}
|
||
|
||
ret_arr[i] = sum;
|
||
}
|
||
}
|
||
else // 待乘向量为竖向量 按行与矩阵做乘法与加法
|
||
{
|
||
if (arr_size != c_num || ret_size != r_num)
|
||
{
|
||
std::string err_str = "Incompatible matrix sizes. From gctl::spmat<T>::multiply(...)";
|
||
throw runtime_error(err_str); // 大小不匹配 不能操作
|
||
}
|
||
|
||
mat_node<MatrixValueType> *p;
|
||
#pragma omp parallel for private (i, sum, p) schedule(guided)
|
||
for (i = 0; i < r_num; i++)
|
||
{
|
||
sum = zero_val;
|
||
|
||
p = r_head[i];
|
||
if (p != nullptr)
|
||
{
|
||
p = p->next_right;
|
||
while (p != nullptr)
|
||
{
|
||
sum += p->val * arr[p->c_id];
|
||
p = p->next_right;
|
||
}
|
||
}
|
||
|
||
ret_arr[i] = sum;
|
||
}
|
||
}
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::multiply_vector(const array<MatrixValueType> &arr,
|
||
array<MatrixValueType> &ret_arr, matrix_layout_e layout) const
|
||
{
|
||
size_t i;
|
||
MatrixValueType sum;
|
||
// 待乘向量为横向量 按列与矩阵做乘法与加法
|
||
if (layout == Trans)
|
||
{
|
||
if (arr.size() != r_num) // 向量大小必须与矩阵行数相等
|
||
{
|
||
std::string err_str = "Incompatible matrix sizes. From gctl::spmat<T>::multiply(...)";
|
||
throw runtime_error(err_str); // 大小不匹配 不能操作
|
||
}
|
||
|
||
ret_arr.resize(c_num);
|
||
|
||
mat_node<MatrixValueType> *q;
|
||
#pragma omp parallel for private (i, sum, q) schedule(guided)
|
||
for (i = 0; i < c_num; i++)
|
||
{
|
||
sum = zero_val;
|
||
|
||
q = c_head[i];
|
||
if (q != nullptr)
|
||
{
|
||
q = q->next_down;
|
||
while (q != nullptr)
|
||
{
|
||
sum += q->val * arr[q->r_id];
|
||
q = q->next_down;
|
||
}
|
||
}
|
||
|
||
ret_arr[i] = sum;
|
||
}
|
||
}
|
||
else // 待乘向量为竖向量 按行与矩阵做乘法与加法
|
||
{
|
||
if (arr.size() != c_num)
|
||
{
|
||
std::string err_str = "Incompatible matrix sizes. From gctl::spmat<T>::multiply(...)";
|
||
throw runtime_error(err_str); // 大小不匹配 不能操作
|
||
}
|
||
|
||
ret_arr.resize(r_num);
|
||
|
||
mat_node<MatrixValueType> *p;
|
||
#pragma omp parallel for private (i, sum, p) schedule(guided)
|
||
for (i = 0; i < r_num; i++)
|
||
{
|
||
sum = zero_val;
|
||
|
||
p = r_head[i];
|
||
if (p != nullptr)
|
||
{
|
||
p = p->next_right;
|
||
while (p != nullptr)
|
||
{
|
||
sum += p->val * arr[p->c_id];
|
||
p = p->next_right;
|
||
}
|
||
}
|
||
|
||
ret_arr[i] = sum;
|
||
}
|
||
}
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
double spmat<MatrixValueType>::multiply_vector(const array<MatrixValueType> &arr,
|
||
int index, matrix_layout_e layout) const
|
||
{
|
||
MatrixValueType sum = zero_val;
|
||
|
||
std::string err_str;
|
||
if (layout == Trans)
|
||
{
|
||
if (arr.size() != r_num)
|
||
{
|
||
err_str = "Incompatible matrix sizes. From gctl::spmat<T>::multiply(...)";
|
||
throw runtime_error(err_str); // 大小不匹配 不能操作
|
||
}
|
||
|
||
if (index < 0 || index >= c_num)
|
||
{
|
||
err_str = "Invalid index. From gctl::spmat<T>::multiply(...)";
|
||
throw runtime_error(err_str);
|
||
}
|
||
|
||
mat_node<MatrixValueType> *q = c_head[index];
|
||
if (q != nullptr)
|
||
{
|
||
q = q->next_down;
|
||
while (q != nullptr)
|
||
{
|
||
sum += q->val * arr[q->r_id];
|
||
q = q->next_down;
|
||
}
|
||
}
|
||
}
|
||
else
|
||
{
|
||
if (arr.size() != c_num)
|
||
{
|
||
err_str = "Incompatible matrix sizes. From gctl::spmat<T>::multiply(...)";
|
||
throw runtime_error(err_str); // 大小不匹配 不能操作
|
||
}
|
||
|
||
if (index < 0 || index >= r_num)
|
||
{
|
||
err_str = "Invalid index. From gctl::spmat<T>::multiply(...)";
|
||
throw runtime_error(err_str);
|
||
}
|
||
|
||
mat_node<MatrixValueType> *p = r_head[index];
|
||
if (p != nullptr)
|
||
{
|
||
p = p->next_right;
|
||
while (p != nullptr)
|
||
{
|
||
sum += p->val * arr[p->c_id];
|
||
p = p->next_right;
|
||
}
|
||
}
|
||
}
|
||
|
||
return sum;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::multiply_vector(int oper_id, const MatrixValueType *arr, int arr_size,
|
||
MatrixValueType *ret_arr, matrix_layout_e layout) const
|
||
{
|
||
std::string err_str;
|
||
if (layout == NoTrans)
|
||
{
|
||
if (oper_id < 0 || oper_id >= r_num)
|
||
{
|
||
err_str = "Invalid index. From gctl::spmat<T>::multiply(...)";
|
||
throw runtime_error(err_str);
|
||
}
|
||
|
||
if (arr_size != c_num) // 向量大小必须与矩阵列数相等
|
||
{
|
||
err_str = "Incompatible matrix sizes. From gctl::spmat<T>::multiply(...)";
|
||
throw runtime_error(err_str); // 大小不匹配 不能操作
|
||
}
|
||
|
||
for (int i = 0; i < arr_size; i++)
|
||
ret_arr[i] = zero_val;
|
||
|
||
mat_node<MatrixValueType> *p;
|
||
if (r_head[oper_id] != nullptr)
|
||
{
|
||
p = r_head[oper_id]->next_right;
|
||
while (p != nullptr)
|
||
{
|
||
ret_arr[p->c_id] = arr[p->c_id] * p->val;
|
||
p = p->next_right;
|
||
}
|
||
}
|
||
}
|
||
// 待乘向量为竖向量 按行与矩阵做乘法与加法
|
||
else
|
||
{
|
||
if (oper_id < 0 || oper_id >= c_num)
|
||
{
|
||
err_str = "Invalid index. From gctl::spmat<T>::multiply(...)";
|
||
throw runtime_error(err_str);
|
||
}
|
||
|
||
if (arr_size != r_num) // 向量大小必须与矩阵列数相等
|
||
{
|
||
std::string err_str = "Incompatible matrix sizes. From gctl::spmat<T>::multiply(...)";
|
||
throw runtime_error(err_str); // 大小不匹配 不能操作
|
||
}
|
||
|
||
for (int i = 0; i < arr_size; i++)
|
||
ret_arr[i] = zero_val;
|
||
|
||
mat_node<MatrixValueType> *q;
|
||
if (c_head[oper_id] != nullptr)
|
||
{
|
||
q = c_head[oper_id]->next_down;
|
||
while (q != nullptr)
|
||
{
|
||
ret_arr[q->r_id] = arr[q->r_id] * q->val;
|
||
q = q->next_down;
|
||
}
|
||
}
|
||
}
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::multiply_diagonal_matrix(const array<MatrixValueType> &arr, spmat<MatrixValueType> &ret_mat, bool prepositioned) const
|
||
{
|
||
if (!ret_mat.empty()) ret_mat.clear();
|
||
|
||
double tmp_val;
|
||
if (prepositioned)
|
||
{
|
||
if (arr.size() != r_num)
|
||
{
|
||
throw runtime_error("Incompatible array sizes. From spmat<T>::multiply(...)");
|
||
}
|
||
|
||
ret_mat.malloc(r_num, c_num, zero_val);
|
||
|
||
mat_node<MatrixValueType> *q;
|
||
for (int j = 0; j < c_num; j++)
|
||
{
|
||
q = c_head[j];
|
||
if (q != nullptr)
|
||
{
|
||
q = q->next_down;
|
||
while (q != nullptr)
|
||
{
|
||
tmp_val = q->val * arr[q->r_id];
|
||
if (tmp_val != zero_val)
|
||
{
|
||
ret_mat.insert(q->r_id, j, tmp_val);
|
||
}
|
||
|
||
q = q->next_down;
|
||
}
|
||
}
|
||
}
|
||
}
|
||
else
|
||
{
|
||
if (arr.size() != c_num)
|
||
{
|
||
throw runtime_error("Incompatible array sizes. From spmat<T>::multiply(...)");
|
||
}
|
||
|
||
ret_mat.malloc(r_num, c_num, zero_val);
|
||
|
||
mat_node<MatrixValueType> *p;
|
||
for (int i = 0; i < r_num; i++)
|
||
{
|
||
p = r_head[i];
|
||
if (p != nullptr)
|
||
{
|
||
p = p->next_right;
|
||
while (p != nullptr)
|
||
{
|
||
tmp_val = p->val * arr[p->c_id];
|
||
if (tmp_val != zero_val)
|
||
{
|
||
ret_mat.insert(i, p->c_id, tmp_val);
|
||
}
|
||
|
||
p = p->next_right;
|
||
}
|
||
}
|
||
}
|
||
}
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::boolean(int oper_id, array<bool> &ret_arr, matrix_order_e ele_mode,
|
||
bool reverse) const
|
||
{
|
||
if (ele_mode == RowMajor)
|
||
{
|
||
if (reverse)
|
||
{
|
||
ret_arr.resize(c_num, true);
|
||
|
||
mat_node<MatrixValueType> *p;
|
||
if (r_head[oper_id] != nullptr)
|
||
{
|
||
p = r_head[oper_id]->next_right;
|
||
while (p != nullptr)
|
||
{
|
||
ret_arr.at(p->c_id) = false;
|
||
p = p->next_right;
|
||
}
|
||
}
|
||
}
|
||
else
|
||
{
|
||
ret_arr.resize(c_num, false);
|
||
|
||
mat_node<MatrixValueType> *p;
|
||
if (r_head[oper_id] != nullptr)
|
||
{
|
||
p = r_head[oper_id]->next_right;
|
||
while (p != nullptr)
|
||
{
|
||
ret_arr.at(p->c_id) = true;
|
||
p = p->next_right;
|
||
}
|
||
}
|
||
}
|
||
}
|
||
else
|
||
{
|
||
if (reverse)
|
||
{
|
||
ret_arr.resize(r_num, true);
|
||
|
||
mat_node<MatrixValueType> *q;
|
||
if (c_head[oper_id] != nullptr)
|
||
{
|
||
q = c_head[oper_id]->next_down;
|
||
while (q != nullptr)
|
||
{
|
||
ret_arr.at(q->r_id) = false;
|
||
q = q->next_down;
|
||
}
|
||
}
|
||
}
|
||
else
|
||
{
|
||
ret_arr.resize(r_num, false);
|
||
|
||
mat_node<MatrixValueType> *q;
|
||
if (c_head[oper_id] != nullptr)
|
||
{
|
||
q = c_head[oper_id]->next_down;
|
||
while (q != nullptr)
|
||
{
|
||
ret_arr.at(q->r_id) = true;
|
||
q = q->next_down;
|
||
}
|
||
}
|
||
}
|
||
}
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::normalize(int oper_id, MatrixValueType norm, matrix_order_e line_mode)
|
||
{
|
||
MatrixValueType sum;
|
||
std::string err_str;
|
||
if (line_mode == RowMajor)
|
||
{
|
||
if (oper_id < 0 || oper_id >= r_num)
|
||
{
|
||
err_str = "Invalid index. From gctl::spmat<T>::normalize(...)";
|
||
throw runtime_error(err_str);
|
||
}
|
||
|
||
mat_node<MatrixValueType> *p;
|
||
if (r_head[oper_id] != nullptr)
|
||
{
|
||
p = r_head[oper_id]->next_right;
|
||
sum = p->val;
|
||
p = p->next_right;
|
||
while (p != nullptr)
|
||
{
|
||
sum += p->val;
|
||
p = p->next_right;
|
||
}
|
||
|
||
p = r_head[oper_id]->next_right;
|
||
while (p != nullptr)
|
||
{
|
||
p->val = p->val*norm/sum;
|
||
p = p->next_right;
|
||
}
|
||
}
|
||
}
|
||
else
|
||
{
|
||
if (oper_id < 0 || oper_id >= c_num)
|
||
{
|
||
err_str = "Invalid index. From gctl::spmat<T>::normalize(...)";
|
||
throw runtime_error(err_str);
|
||
}
|
||
|
||
mat_node<MatrixValueType> *q;
|
||
if (c_head[oper_id] != nullptr)
|
||
{
|
||
q = c_head[oper_id]->next_down;
|
||
sum = q->val;
|
||
q = q->next_down;
|
||
while (q != nullptr)
|
||
{
|
||
sum += q->val;
|
||
q = q->next_down;
|
||
}
|
||
|
||
q = c_head[oper_id]->next_down;
|
||
while (q != nullptr)
|
||
{
|
||
q->val = q->val*norm/sum;
|
||
q = q->next_down;
|
||
}
|
||
}
|
||
}
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::scale(MatrixValueType f, int idx, matrix_order_e ele_mode)
|
||
{
|
||
if (ele_mode == RowMajor)
|
||
{
|
||
if (idx < 0 || idx >= r_num)
|
||
{
|
||
throw std::runtime_error("[gctl::spmat<T>::scale] Invalid index.");
|
||
}
|
||
|
||
mat_node<MatrixValueType> *p;
|
||
if (r_head[idx] != nullptr)
|
||
{
|
||
p = r_head[idx]->next_right;
|
||
while (p != nullptr)
|
||
{
|
||
p->val = p->val*f;
|
||
p = p->next_right;
|
||
}
|
||
}
|
||
}
|
||
else
|
||
{
|
||
if (idx < 0 || idx >= c_num)
|
||
{
|
||
throw std::runtime_error("[gctl::spmat<T>::scale] Invalid index.");
|
||
}
|
||
|
||
mat_node<MatrixValueType> *q;
|
||
if (c_head[idx] != nullptr)
|
||
{
|
||
q = c_head[idx]->next_down;
|
||
while (q != nullptr)
|
||
{
|
||
q->val = q->val*f;
|
||
q = q->next_down;
|
||
}
|
||
}
|
||
}
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::get_diagonal(array<MatrixValueType> &dia)
|
||
{
|
||
if (r_num != c_num)
|
||
{
|
||
throw gctl::runtime_error("Not a square matrix. From gctl::spmat<T>::get_diagonal(...)");
|
||
}
|
||
|
||
dia.resize(r_num, zero_val);
|
||
|
||
mat_node<MatrixValueType> *p;
|
||
for (size_t i = 0; i < r_num; i++)
|
||
{
|
||
if (r_head[i] != nullptr)
|
||
{
|
||
p = r_head[i]->next_right;
|
||
while (p != nullptr)
|
||
{
|
||
if (p->c_id == i)
|
||
{
|
||
dia[i] = p->val;
|
||
break;
|
||
}
|
||
|
||
p = p->next_right;
|
||
}
|
||
}
|
||
}
|
||
return;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
spmat<MatrixValueType>& spmat<MatrixValueType>::operator= (const spmat<MatrixValueType> &b)
|
||
{
|
||
if (!empty()) clear();
|
||
|
||
r_num = b.row_size();
|
||
c_num = b.col_size();
|
||
r_head = new mat_node<MatrixValueType>* [r_num];
|
||
c_head = new mat_node<MatrixValueType>* [c_num];
|
||
r_end = new mat_node<MatrixValueType>* [r_num];
|
||
c_end = new mat_node<MatrixValueType>* [c_num];
|
||
r_count = new int [r_num];
|
||
c_count = new int [c_num];
|
||
mat_node<MatrixValueType> **br_head = b.row_head();
|
||
|
||
int i, j;
|
||
for (i = 0; i < r_num; i++)
|
||
{
|
||
r_head[i] = nullptr; r_end[i] = nullptr; r_count[i] = 0;
|
||
}
|
||
for (i = 0; i < c_num; i++)
|
||
{
|
||
c_head[i] = nullptr; c_end[i] = nullptr; c_count[i] = 0;
|
||
}
|
||
|
||
mat_node<MatrixValueType> *bp, *p, *q;
|
||
for (int r = 0; r < r_num; r++)
|
||
{
|
||
if (br_head[r] != nullptr)
|
||
{
|
||
bp = br_head[r]->next_right;
|
||
while (bp != nullptr)
|
||
{
|
||
mat_node<MatrixValueType>* new_node =
|
||
new mat_node<MatrixValueType>;
|
||
new_node->r_id = bp->r_id;
|
||
new_node->c_id = bp->c_id;
|
||
new_node->val = bp->val;
|
||
|
||
i = bp->r_id; j = bp->c_id;
|
||
// 我们首先处理行的指针
|
||
if (r_head[i] != nullptr) // 无需检查行尾指针 它们是同步变化的
|
||
{
|
||
p = r_end[i];
|
||
p->next_right = new_node;
|
||
r_end[i] = new_node;
|
||
r_count[i]++;
|
||
}
|
||
else
|
||
{
|
||
// 添加行首指针 注意行首指针并不等于第一个元素 而是行首的右向链为第一个元素
|
||
r_head[i] = new mat_node<MatrixValueType>;
|
||
r_head[i]->next_right = new_node;
|
||
// 添加行尾指针 行尾指针直接等于本行的最后一个元素
|
||
r_end[i] = new_node;
|
||
r_count[i]++;
|
||
}
|
||
|
||
// 接着处理列的指针
|
||
if (c_head[j] != nullptr)
|
||
{
|
||
q = c_end[j];
|
||
q->next_down = new_node;
|
||
c_end[j] = new_node;
|
||
c_count[j]++;
|
||
}
|
||
else
|
||
{
|
||
c_head[j] = new mat_node<MatrixValueType>;
|
||
c_head[j]->next_down = new_node;
|
||
c_end[j] = new_node;
|
||
c_count[j]++;
|
||
}
|
||
|
||
bp = bp->next_right;
|
||
}
|
||
}
|
||
}
|
||
return *this;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
bool spmat<MatrixValueType>::change_element(int r, int c, MatrixValueType v,
|
||
value_operator_e mode)
|
||
{
|
||
if (r >= r_num || r < 0 || c >= c_num || c < 0)
|
||
{
|
||
throw std::runtime_error("Invalid index. From gctl::spmat<T>::change_element(...)");
|
||
}
|
||
|
||
// 输入的值为0 无效的操作返回假
|
||
if (v == zero_val)
|
||
return false;
|
||
|
||
mat_node<MatrixValueType> *new_node = nullptr;
|
||
mat_node<MatrixValueType> *p, *q, *tmp_node;
|
||
|
||
MatrixValueType in_v;
|
||
if (mode == RemoveVal) in_v = -v;
|
||
else in_v = v;
|
||
|
||
bool d_flag = false;
|
||
|
||
// 首先处理行的指针序列
|
||
// 如果指针不为空 则需要判断插入的位置
|
||
if (r_head[r] != nullptr)
|
||
{
|
||
// 首先检查行尾值
|
||
if (c > r_end[r]->c_id) // 大于行尾的列值 在队尾添加
|
||
{
|
||
// 需要新建节点
|
||
new_node = new mat_node<MatrixValueType>;
|
||
new_node->r_id = r;
|
||
new_node->c_id = c;
|
||
new_node->val = in_v;
|
||
|
||
r_end[r]->next_right = new_node;
|
||
r_end[r] = new_node;
|
||
r_count[r]++;
|
||
}
|
||
else if (c == r_end[r]->c_id)
|
||
{
|
||
if (mode == ReplaceVal)
|
||
{
|
||
r_end[r]->val = in_v;
|
||
}
|
||
else if (mode == AppendVal || mode == RemoveVal)
|
||
{
|
||
r_end[r]->val += in_v;
|
||
if (r_end[r]->val == zero_val)
|
||
d_flag = true;
|
||
}
|
||
else // 操作类型是EraseVal 将元素值置为0 （一般不会用到）
|
||
{
|
||
r_end[r]->val = zero_val;
|
||
d_flag = true;
|
||
}
|
||
}
|
||
else
|
||
{
|
||
p = r_head[r]->next_right; // 该行的首个元素
|
||
if (c < p->c_id) // 输入列值小于行首元素的列值 更改行首元素
|
||
{
|
||
// 需要新建节点
|
||
new_node = new mat_node<MatrixValueType>;
|
||
new_node->r_id = r;
|
||
new_node->c_id = c;
|
||
new_node->val = in_v;
|
||
|
||
r_head[r]->next_right = new_node;
|
||
new_node->next_right = p;
|
||
r_count[r]++;
|
||
}
|
||
else if (c == p->c_id) // 列序列等于行首值，替换节点值就行了
|
||
{
|
||
if (mode == ReplaceVal)
|
||
{
|
||
p->val = in_v;
|
||
}
|
||
else if (mode == AppendVal || mode == RemoveVal)
|
||
{
|
||
p->val += in_v;
|
||
if (p->val == zero_val)
|
||
d_flag = true;
|
||
}
|
||
else // 操作类型是EraseVal 将元素值置为0 （一般不会用到）
|
||
{
|
||
p->val = zero_val;
|
||
d_flag = true;
|
||
}
|
||
}
|
||
else
|
||
{
|
||
// 这个寻址过程是影响效率的主要部分
|
||
while (p->next_right != nullptr)
|
||
{
|
||
if (p->next_right->c_id == c) // 该行的某个节点的列号与输入点一致
|
||
{
|
||
if (mode == ReplaceVal)
|
||
{
|
||
p->next_right->val = in_v;
|
||
}
|
||
else if (mode == AppendVal || mode == RemoveVal)
|
||
{
|
||
p->next_right->val += in_v;
|
||
if (p->next_right->val == zero_val)
|
||
d_flag = true;
|
||
}
|
||
else // 操作类型是EraseVal 将元素值置为0 （一般不会用到）
|
||
{
|
||
p->next_right->val = zero_val;
|
||
d_flag = true;
|
||
}
|
||
break;
|
||
}
|
||
else if (c > p->c_id && c < p->next_right->c_id)
|
||
{
|
||
// 需要新建节点
|
||
new_node = new mat_node<MatrixValueType>;
|
||
new_node->r_id = r;
|
||
new_node->c_id = c;
|
||
new_node->val = in_v;
|
||
|
||
tmp_node = p->next_right;
|
||
p->next_right = new_node;
|
||
new_node->next_right= tmp_node;
|
||
r_count[r]++;
|
||
break;
|
||
}
|
||
p = p->next_right;
|
||
}
|
||
}
|
||
}
|
||
}
|
||
// 若行指针为空 直接添加新元素即可
|
||
else
|
||
{
|
||
// 需要新建节点
|
||
new_node = new mat_node<MatrixValueType>;
|
||
new_node->r_id = r;
|
||
new_node->c_id = c;
|
||
new_node->val = in_v;
|
||
|
||
r_head[r] = new mat_node<MatrixValueType>;
|
||
r_head[r]->next_right = new_node;
|
||
r_end[r] = new_node;
|
||
r_count[r]++;
|
||
}
|
||
|
||
// 下面来处理列序列
|
||
// 如果指针不为空 则需要判断插入的位置
|
||
if (c_head[c] != nullptr)
|
||
{
|
||
// 首先检查列尾值
|
||
if (r > c_end[c]->r_id) // 大于行尾的列值 在队尾添加
|
||
{
|
||
c_end[c]->next_down = new_node;
|
||
c_end[c] = new_node;
|
||
c_count[c]++;
|
||
}
|
||
else
|
||
{
|
||
q = c_head[c]->next_down;
|
||
if (r < q->r_id) // 输入行小于列首行 替换即可
|
||
{
|
||
c_head[c]->next_down = new_node;
|
||
new_node->next_down = q;
|
||
c_count[c]++;
|
||
}
|
||
else
|
||
{
|
||
// 这个寻址过程是影响效率的主要部分
|
||
while (q->next_down != nullptr)
|
||
{
|
||
if (r == q->next_down->r_id)
|
||
break;
|
||
else if (r > q->r_id && r < q->next_down->r_id)
|
||
{
|
||
tmp_node = q->next_down;
|
||
q->next_down = new_node;
|
||
new_node->next_down = tmp_node;
|
||
c_count[c]++;
|
||
break;
|
||
}
|
||
q = q->next_down;
|
||
}
|
||
}
|
||
}
|
||
}
|
||
// 若行指针为空 直接添加新元素即可
|
||
else
|
||
{
|
||
c_head[c] = new mat_node<MatrixValueType>;
|
||
c_head[c]->next_down = new_node;
|
||
c_end[c] = new_node;
|
||
c_count[c]++;
|
||
}
|
||
|
||
// 如果触发删除标签则删除相应元素
|
||
if (d_flag) return true;
|
||
return false;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
mat_node<MatrixValueType> **spmat<MatrixValueType>::row_head() const
|
||
{
|
||
return r_head;
|
||
}
|
||
|
||
template <typename MatrixValueType>
|
||
mat_node<MatrixValueType> **spmat<MatrixValueType>::col_head() const
|
||
{
|
||
return c_head;
|
||
}
|
||
|
||
#ifdef GCTL_EIGEN
|
||
|
||
template <typename MatrixValueType>
|
||
void spmat<MatrixValueType>::multiply_vector(const Eigen::VectorXd &vec, Eigen::VectorXd &ret_vec, bool row_vector) const
|
||
{
|
||
size_t i;
|
||
MatrixValueType sum;
|
||
// 待乘向量为横向量 按列与矩阵做乘法与加法
|
||
if (row_vector)
|
||
{
|
||
if (vec.size() != r_num) // 向量大小必须与矩阵行数相等
|
||
{
|
||
std::string err_str = "Incompatible matrix sizes. From gctl::spmat<T>::multiply(...)";
|
||
throw runtime_error(err_str); // 大小不匹配 不能操作
|
||
}
|
||
|
||
ret_vec.resize(c_num);
|
||
|
||
mat_node<MatrixValueType> *q;
|
||
#pragma omp parallel for private (i, sum, q) schedule(guided)
|
||
for (i = 0; i < c_num; i++)
|
||
{
|
||
sum = zero_val;
|
||
|
||
q = c_head[i];
|
||
if (q != nullptr)
|
||
{
|
||
q = q->next_down;
|
||
while (q != nullptr)
|
||
{
|
||
sum += q->val * vec[q->r_id];
|
||
q = q->next_down;
|
||
}
|
||
}
|
||
|
||
ret_vec[i] = sum;
|
||
}
|
||
}
|
||
else // 待乘向量为竖向量 按行与矩阵做乘法与加法
|
||
{
|
||
if (vec.size() != c_num)
|
||
{
|
||
std::string err_str = "Incompatible matrix sizes. From gctl::spmat<T>::multiply(...)";
|
||
throw runtime_error(err_str); // 大小不匹配 不能操作
|
||
}
|
||
|
||
ret_vec.resize(r_num);
|
||
|
||
mat_node<MatrixValueType> *p;
|
||
#pragma omp parallel for private (i, sum, p) schedule(guided)
|
||
for (i = 0; i < r_num; i++)
|
||
{
|
||
sum = zero_val;
|
||
|
||
p = r_head[i];
|
||
if (p != nullptr)
|
||
{
|
||
p = p->next_right;
|
||
while (p != nullptr)
|
||
{
|
||
sum += p->val * vec[p->c_id];
|
||
p = p->next_right;
|
||
}
|
||
}
|
||
|
||
ret_vec[i] = sum;
|
||
}
|
||
}
|
||
return;
|
||
}
|
||
|
||
#endif // GCTL_EIGEN
|
||
|
||
}
|
||
|
||
#endif //_GCTL_SPMAT_H
|