1074 lines
25 KiB
C
1074 lines
25 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_MATRIX_H
|
|||
|
|
#define _GCTL_MATRIX_H
|
|||
|
|
|
|||
|
|
#include "array.h"
|
|||
|
|
|
|||
|
|
#include "random"
|
|||
|
|
#include "chrono"
|
|||
|
|
#include "typeinfo"
|
|||
|
|
|
|||
|
|
namespace gctl
|
|||
|
|
{
|
|||
|
|
template <typename MatValType> class matrix;
|
|||
|
|
|
|||
|
|
/*
|
|||
|
|
* Here we define some commonly used array types.
|
|||
|
|
*/
|
|||
|
|
typedef matrix<int> _2i_matrix; ///< 2-D integer array
|
|||
|
|
typedef matrix<float> _2f_matrix; ///< 2-D single-precision floating-point array
|
|||
|
|
typedef matrix<double> _2d_matrix; ///< 2-D double-precision floating-point array
|
|||
|
|
typedef matrix<std::string> _2s_matrix; ///< 2-D string array
|
|||
|
|
typedef matrix<std::complex<float>> _2cf_matrix; ///< 1-D complex float array
|
|||
|
|
typedef matrix<std::complex<double>> _2cd_matrix; ///< 1-D complex double array
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 二维数组模版
|
|||
|
|
*
|
|||
|
|
* 二维数组模版是一个简单的二维动态数组实现，包含简单的数组操作与IO操作函数。
|
|||
|
|
*
|
|||
|
|
* @tparam MatValType 数组元素模版类型
|
|||
|
|
*/
|
|||
|
|
template <typename MatValType>
|
|||
|
|
class matrix
|
|||
|
|
{
|
|||
|
|
protected:
|
|||
|
|
MatValType **val; ///< 二维数组数据
|
|||
|
|
size_t row_length, col_length; ///< 二维数组的行列数
|
|||
|
|
|
|||
|
|
public:
|
|||
|
|
/**
|
|||
|
|
* @brief 默认构造函数
|
|||
|
|
*/
|
|||
|
|
matrix();
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 构造函数：初始化二维数组空间
|
|||
|
|
*
|
|||
|
|
* @param[in] row_len 二维数组行数
|
|||
|
|
* @param[in] col_len 二维数组列数
|
|||
|
|
*/
|
|||
|
|
matrix(size_t row_len, size_t col_len);
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 构造函数：初始化二维数组空间并赋处值
|
|||
|
|
*
|
|||
|
|
* @param[in] row_len 二维数组行数
|
|||
|
|
* @param[in] col_len 二维数组列数
|
|||
|
|
* @param[in] init_val 数组元素的初始值
|
|||
|
|
*/
|
|||
|
|
matrix(size_t row_len, size_t col_len, MatValType init_val);
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 构造函数：从一维数组初始化二维数组
|
|||
|
|
*
|
|||
|
|
* @param init_mat 输入的一维数组
|
|||
|
|
* @param[in] row_len 二维数组行数
|
|||
|
|
* @param[in] col_len 二维数组列数
|
|||
|
|
*/
|
|||
|
|
matrix(MatValType *init_mat, size_t row_len, size_t col_len);
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 构造函数：从二维数组初始化二维数组
|
|||
|
|
*
|
|||
|
|
* @param init_mat 输入的二维数组
|
|||
|
|
* @param[in] row_len 二维数组行数
|
|||
|
|
* @param[in] col_len 二维数组列数
|
|||
|
|
*/
|
|||
|
|
matrix(MatValType **init_mat, size_t row_len, size_t col_len);
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 构造函数：从一维向量初始化二维数组
|
|||
|
|
*
|
|||
|
|
* @param[in] b 输入向量的引用
|
|||
|
|
* @param[in] row_len 二维数组行数
|
|||
|
|
* @param[in] col_len 二维数组列数
|
|||
|
|
*/
|
|||
|
|
matrix(const std::vector<MatValType> &b, size_t row_len, size_t col_len);
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 构造函数：从二维向量初始化二维数组
|
|||
|
|
*
|
|||
|
|
* @param[in] b 输入向量的引用
|
|||
|
|
*/
|
|||
|
|
matrix(const std::vector<std::vector<MatValType> > &b);
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 构造函数：从一维数组初始化二维数组
|
|||
|
|
*
|
|||
|
|
* @param[in] b 输入数组的引用
|
|||
|
|
* @param[in] row_len 二维数组行数
|
|||
|
|
* @param[in] col_len 二维数组列数
|
|||
|
|
*/
|
|||
|
|
matrix(const array<MatValType> &b, size_t row_len, size_t col_len);
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 拷贝构造函数
|
|||
|
|
*
|
|||
|
|
* @param[in] b 拷贝数组
|
|||
|
|
*/
|
|||
|
|
matrix(const matrix<MatValType> &b);
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 赋值操作符重载
|
|||
|
|
*
|
|||
|
|
* @param[in] b 赋值数组
|
|||
|
|
*
|
|||
|
|
* @return 新建数组
|
|||
|
|
*/
|
|||
|
|
matrix<MatValType>& operator= (const matrix<MatValType> &b);
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 析构函数
|
|||
|
|
*/
|
|||
|
|
virtual ~matrix();
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 初始化二维数组空间
|
|||
|
|
*
|
|||
|
|
* @param[in] row_len 二维数组行数
|
|||
|
|
* @param[in] col_len 二维数组列数
|
|||
|
|
*/
|
|||
|
|
void resize(size_t row_len, size_t col_len);
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 构造函数：初始化二维数组空间并赋处值
|
|||
|
|
*
|
|||
|
|
* @param[in] row_len 二维数组行数
|
|||
|
|
* @param[in] col_len 二维数组列数
|
|||
|
|
* @param[in] init_val 数组元素的初始值
|
|||
|
|
*/
|
|||
|
|
void resize(size_t row_len, size_t col_len, MatValType init_val);
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 构造函数：从一维数组初始化二维数组
|
|||
|
|
*
|
|||
|
|
* @param init_mat 输入的一维数组
|
|||
|
|
* @param[in] row_len 二维数组行数
|
|||
|
|
* @param[in] col_len 二维数组列数
|
|||
|
|
*/
|
|||
|
|
void resize(MatValType *init_mat, size_t row_len, size_t col_len);
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 构造函数：从二维数组初始化二维数组
|
|||
|
|
*
|
|||
|
|
* @param init_mat 输入的二维数组
|
|||
|
|
* @param[in] row_len 二维数组行数
|
|||
|
|
* @param[in] col_len 二维数组列数
|
|||
|
|
*/
|
|||
|
|
void resize(MatValType **init_mat, size_t row_len, size_t col_len);
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 构造函数：从一维向量初始化二维数组
|
|||
|
|
*
|
|||
|
|
* @param[in] b 输入向量的引用
|
|||
|
|
* @param[in] row_len 二维数组行数
|
|||
|
|
* @param[in] col_len 二维数组列数
|
|||
|
|
*/
|
|||
|
|
void resize(const std::vector<MatValType> &b, size_t row_len, size_t col_len);
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 从二维向量初始化二维数组
|
|||
|
|
*
|
|||
|
|
* @param[in] b 输入向量的引用
|
|||
|
|
*/
|
|||
|
|
void resize(const std::vector<std::vector<MatValType> > &b);
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 构造函数：从一维数组初始化二维数组
|
|||
|
|
*
|
|||
|
|
* @param[in] b 输入数组的引用
|
|||
|
|
* @param[in] row_len 二维数组行数
|
|||
|
|
* @param[in] col_len 二维数组列数
|
|||
|
|
*/
|
|||
|
|
void resize(const array<MatValType> &b, size_t row_len, size_t col_len);
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 从二维数组初始化二维数组
|
|||
|
|
*
|
|||
|
|
* @param[in] b 输入数组
|
|||
|
|
*/
|
|||
|
|
void resize(const matrix<MatValType> &b);
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 清除并重置数组空间
|
|||
|
|
*/
|
|||
|
|
void clear();
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief Initialize the array with selected random types.
|
|||
|
|
*
|
|||
|
|
* @param[in] np1 Mean (Gauss) or low bound value (Even)
|
|||
|
|
* @param[in] np2 Standard deviation (Gauss) or hig bound value (Even).
|
|||
|
|
* @param[in] mode Random types. 'RdNormal' for Gaussian distributed numbers and
|
|||
|
|
* 'RdUniform' for even distributed numbers.
|
|||
|
|
*/
|
|||
|
|
void random(MatValType np1, MatValType np2, random_type_e mode = RdNormal, unsigned int seed = 0);
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 对全体元素进行赋值
|
|||
|
|
*
|
|||
|
|
* @param[in] in_val 输入值
|
|||
|
|
*/
|
|||
|
|
void assign_all(MatValType in_val);
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 对全体元素进行缩放
|
|||
|
|
*
|
|||
|
|
* @param s 缩放的比例
|
|||
|
|
*/
|
|||
|
|
void scale(MatValType s);
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 按行或列计算元素平均值
|
|||
|
|
*
|
|||
|
|
* @param m 均值数组
|
|||
|
|
* @param odr 计算行或列的平均值
|
|||
|
|
*/
|
|||
|
|
void mean(array<MatValType> &m, matrix_order_e odr = RowMajor);
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 计算矩阵的模长
|
|||
|
|
*
|
|||
|
|
* @param n_type 模的类型
|
|||
|
|
*/
|
|||
|
|
double norm(norm_type_e n_type = L2) const;
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief Set elements' value as a sequent.
|
|||
|
|
*
|
|||
|
|
* @param st_val Start value.
|
|||
|
|
* @param row_inc Row increasement.
|
|||
|
|
* @param col_inc Column increasement.
|
|||
|
|
*/
|
|||
|
|
void sequent(MatValType st_val, MatValType row_inc, MatValType col_inc);
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 获取索引位置元素值
|
|||
|
|
*
|
|||
|
|
* @note 可用于指针情况
|
|||
|
|
*
|
|||
|
|
* @param[in] row_index 行索引值
|
|||
|
|
* @param[in] col_index 列索引值
|
|||
|
|
*
|
|||
|
|
* @return 索引处元素
|
|||
|
|
*/
|
|||
|
|
MatValType &at(size_t row_index, size_t col_index);
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 获取索引位置元素值（常量重载）
|
|||
|
|
*
|
|||
|
|
* @note 可用于指针情况
|
|||
|
|
*
|
|||
|
|
* @param[in] row_index 行索引值
|
|||
|
|
* @param[in] col_index 列索引值
|
|||
|
|
*
|
|||
|
|
* @return 索引处元素
|
|||
|
|
*/
|
|||
|
|
MatValType &at(size_t row_index, size_t col_index) const;
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 下标符号重载：获取索引位置元素值
|
|||
|
|
*
|
|||
|
|
* @note 不可用于指针情况，对于二维数组只需重载一次下标符号就行。
|
|||
|
|
*
|
|||
|
|
* @param[in] index 索引值
|
|||
|
|
*
|
|||
|
|
* @return 索引处元素
|
|||
|
|
*/
|
|||
|
|
MatValType *operator[](size_t index);
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 下标符号重载：获取索引位置元素值（常量重载）
|
|||
|
|
*
|
|||
|
|
* @note 不可用于指针情况，对于二维数组只需重载一次下标符号就行。
|
|||
|
|
*
|
|||
|
|
* @param[in] index 索引值
|
|||
|
|
*
|
|||
|
|
* @return 索引处元素
|
|||
|
|
*/
|
|||
|
|
MatValType *operator[](size_t index) const;
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 获取成员数组裸指针
|
|||
|
|
*
|
|||
|
|
* @param[in] row_index 行索引值 row_index = 0 时数组指针为全部矩阵数据的头指针
|
|||
|
|
*
|
|||
|
|
* @return 单个元素的指针
|
|||
|
|
*/
|
|||
|
|
MatValType *get(size_t row_index = 0) const;
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 判断数组是否为空
|
|||
|
|
*
|
|||
|
|
* @return 数组为空返回真，否则返回假
|
|||
|
|
*/
|
|||
|
|
bool empty() const;
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 返回矩阵行大小
|
|||
|
|
*
|
|||
|
|
* @return 矩阵行大小
|
|||
|
|
*/
|
|||
|
|
size_t row_size() const;
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 返回矩阵列大小
|
|||
|
|
*
|
|||
|
|
* @return 矩阵列大小
|
|||
|
|
*/
|
|||
|
|
size_t col_size() const;
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 返回矩阵数组大小
|
|||
|
|
*
|
|||
|
|
* @return 数组大小
|
|||
|
|
*/
|
|||
|
|
size_t val_size() const;
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 拷贝数组到新的二维数组
|
|||
|
|
*
|
|||
|
|
* @param b 拷贝对象数组
|
|||
|
|
*/
|
|||
|
|
void copy_to(matrix<MatValType> &b) const;
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 拷贝数组到二维向量
|
|||
|
|
*
|
|||
|
|
* @param b 拷贝对象向量
|
|||
|
|
*/
|
|||
|
|
void copy_to(std::vector<std::vector<MatValType> > &b) const;
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 提取二维数组的行或者列到输出数组
|
|||
|
|
*
|
|||
|
|
* @param b 输出一维数组
|
|||
|
|
* @param[in] index 行或者列的索引值
|
|||
|
|
* @param[in] order 按行或者列输出
|
|||
|
|
*/
|
|||
|
|
void extract(array<MatValType> &b, size_t index, matrix_order_e order = RowMajor) const;
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief 将二维数组转换为一维数组
|
|||
|
|
*
|
|||
|
|
* @warning 此函数会重置输出的一维数组
|
|||
|
|
*
|
|||
|
|
* @param b 转换后数组
|
|||
|
|
* @param[in] order 按行或者列优先排列输出
|
|||
|
|
*/
|
|||
|
|
void reform(array<MatValType> &b, matrix_order_e order = RowMajor) const;
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief element operate function pointer
|
|||
|
|
*
|
|||
|
|
*/
|
|||
|
|
typedef void (*foreach_m_ptr)(MatValType *ele_ptr, size_t r_id, size_t c_id);
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief Operate on each and every element
|
|||
|
|
*
|
|||
|
|
* @param func operation function
|
|||
|
|
*/
|
|||
|
|
void for_each(foreach_m_ptr func);
|
|||
|
|
|
|||
|
|
/**
|
|||
|
|
* @brief Display the elements.
|
|||
|
|
*
|
|||
|
|
*/
|
|||
|
|
void show(std::ostream &os = std::cout, char sep = '\t') const;
|
|||
|
|
};
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
matrix<MatValType>::matrix()
|
|||
|
|
{
|
|||
|
|
row_length = col_length = 0;
|
|||
|
|
val = nullptr;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
matrix<MatValType>::matrix(size_t row_len, size_t col_len)
|
|||
|
|
{
|
|||
|
|
row_length = col_length = 0;
|
|||
|
|
val = nullptr;
|
|||
|
|
resize(row_len, col_len);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
matrix<MatValType>::matrix(size_t row_len, size_t col_len, MatValType init_val)
|
|||
|
|
{
|
|||
|
|
row_length = col_length = 0;
|
|||
|
|
val = nullptr;
|
|||
|
|
resize(row_len, col_len, init_val);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
matrix<MatValType>::matrix(MatValType *init_mat, size_t row_len, size_t col_len)
|
|||
|
|
{
|
|||
|
|
row_length = col_length = 0;
|
|||
|
|
val = nullptr;
|
|||
|
|
resize(init_mat, row_len, col_len);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
matrix<MatValType>::matrix(MatValType **init_mat, size_t row_len, size_t col_len)
|
|||
|
|
{
|
|||
|
|
row_length = col_length = 0;
|
|||
|
|
val = nullptr;
|
|||
|
|
resize(init_mat, row_len, col_len);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
matrix<MatValType>::matrix(const std::vector<MatValType> &b, size_t row_len, size_t col_len)
|
|||
|
|
{
|
|||
|
|
row_length = col_length = 0;
|
|||
|
|
val = nullptr;
|
|||
|
|
resize(b, row_len, col_len);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
matrix<MatValType>::matrix(const std::vector<std::vector<MatValType> > &b)
|
|||
|
|
{
|
|||
|
|
row_length = col_length = 0;
|
|||
|
|
val = nullptr;
|
|||
|
|
resize(b);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
matrix<MatValType>::matrix(const array<MatValType> &b, size_t row_len, size_t col_len)
|
|||
|
|
{
|
|||
|
|
row_length = col_length = 0;
|
|||
|
|
val = nullptr;
|
|||
|
|
resize(b, row_len, col_len);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
matrix<MatValType>::matrix(const matrix<MatValType> &b)
|
|||
|
|
{
|
|||
|
|
row_length = col_length = 0;
|
|||
|
|
val = nullptr;
|
|||
|
|
resize(b);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
matrix<MatValType>& matrix<MatValType>::operator= (const matrix<MatValType> &b)
|
|||
|
|
{
|
|||
|
|
if (b.empty())
|
|||
|
|
{
|
|||
|
|
clear();
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
resize(b.row_size(), b.col_size());
|
|||
|
|
for (size_t i = 0; i < row_length; i++)
|
|||
|
|
{
|
|||
|
|
for (size_t j = 0; j < col_length; j++)
|
|||
|
|
{
|
|||
|
|
val[i][j] = b[i][j];
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
return *this;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
matrix<MatValType>::~matrix()
|
|||
|
|
{
|
|||
|
|
clear();
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
void matrix<MatValType>::resize(size_t row_len, size_t col_len)
|
|||
|
|
{
|
|||
|
|
if (row_len == 0 || col_len == 0)
|
|||
|
|
{
|
|||
|
|
throw std::invalid_argument("Invalid matrix size. gctl::matrix<T>::resize(...)");
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// nothing needs to be done
|
|||
|
|
if (row_len == row_length && col_len == col_length)
|
|||
|
|
{
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// re-indexing
|
|||
|
|
if (row_len*col_len == row_length*col_length)
|
|||
|
|
{
|
|||
|
|
row_length = row_len;
|
|||
|
|
col_length = col_len;
|
|||
|
|
|
|||
|
|
// construct new row-indexing
|
|||
|
|
MatValType *t_ptr = val[0];
|
|||
|
|
delete[] val;
|
|||
|
|
val = new MatValType* [row_len];
|
|||
|
|
val[0] = t_ptr;
|
|||
|
|
|
|||
|
|
for (size_t i = 1; i < row_len; i++)
|
|||
|
|
{
|
|||
|
|
val[i] = val[i-1] + col_len;
|
|||
|
|
}
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
// re-allocate
|
|||
|
|
if (!empty()) clear();
|
|||
|
|
|
|||
|
|
row_length = row_len;
|
|||
|
|
col_length = col_len;
|
|||
|
|
|
|||
|
|
val = new MatValType* [row_len];
|
|||
|
|
val[0] = new MatValType [row_len*col_len];
|
|||
|
|
|
|||
|
|
for (size_t i = 1; i < row_len; i++)
|
|||
|
|
{
|
|||
|
|
val[i] = val[i-1] + col_len;
|
|||
|
|
}
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
void matrix<MatValType>::resize(size_t row_len, size_t col_len, MatValType init_val)
|
|||
|
|
{
|
|||
|
|
resize(row_len, col_len);
|
|||
|
|
|
|||
|
|
for (size_t i = 0; i < row_len; i++)
|
|||
|
|
{
|
|||
|
|
for (size_t j = 0; j < col_len; j++)
|
|||
|
|
{
|
|||
|
|
val[i][j] = init_val;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
void matrix<MatValType>::resize(MatValType *init_mat, size_t row_len, size_t col_len)
|
|||
|
|
{
|
|||
|
|
if (init_mat == nullptr)
|
|||
|
|
{
|
|||
|
|
throw std::domain_error("Invalid pointer. gctl::matrix<T>::resize(...)");
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
resize(row_len, col_len);
|
|||
|
|
|
|||
|
|
for (size_t i = 0; i < row_len; i++)
|
|||
|
|
{
|
|||
|
|
for (size_t j = 0; j < col_len; j++)
|
|||
|
|
{
|
|||
|
|
val[i][j] = init_mat[j + i*col_len];
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
void matrix<MatValType>::resize(MatValType **init_mat, size_t row_len, size_t col_len)
|
|||
|
|
{
|
|||
|
|
if (init_mat == nullptr)
|
|||
|
|
{
|
|||
|
|
throw std::domain_error("Invalid pointer. gctl::matrix<T>::resize(...)");
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
resize(row_len, col_len);
|
|||
|
|
|
|||
|
|
for (size_t i = 0; i < row_len; i++)
|
|||
|
|
{
|
|||
|
|
for (size_t j = 0; j < col_len; j++)
|
|||
|
|
{
|
|||
|
|
val[i][j] = init_mat[i][j];
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
void matrix<MatValType>::resize(const std::vector<MatValType> &b, size_t row_len, size_t col_len)
|
|||
|
|
{
|
|||
|
|
if (!b.empty())
|
|||
|
|
{
|
|||
|
|
resize(row_len, col_len);
|
|||
|
|
|
|||
|
|
for (size_t i = 0; i < row_len; i++)
|
|||
|
|
{
|
|||
|
|
for (size_t j = 0; j < col_len; j++)
|
|||
|
|
{
|
|||
|
|
val[i][j] = b[j + i*col_len];
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
void matrix<MatValType>::resize(const std::vector<std::vector<MatValType> > &b)
|
|||
|
|
{
|
|||
|
|
if (!b.empty())
|
|||
|
|
{
|
|||
|
|
size_t col_s = b[0].size();
|
|||
|
|
for (size_t i = 1; i < b.size(); i++)
|
|||
|
|
{
|
|||
|
|
if (col_s != b[i].size())
|
|||
|
|
{
|
|||
|
|
throw std::runtime_error("The input 2-D vector is not neatly arranged. From matrix<T>::resize(...)");
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
resize(b.size(), col_s);
|
|||
|
|
for (size_t i = 0; i < row_length; i++)
|
|||
|
|
{
|
|||
|
|
for (size_t j = 0; j < col_length; j++)
|
|||
|
|
{
|
|||
|
|
val[i][j] = b[i][j];
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
void matrix<MatValType>::resize(const array<MatValType> &b, size_t row_len, size_t col_len)
|
|||
|
|
{
|
|||
|
|
if (!b.empty())
|
|||
|
|
{
|
|||
|
|
if (b.size() != row_len*col_len)
|
|||
|
|
{
|
|||
|
|
throw std::runtime_error("Invalid initializing parameters. From matrix<T>::resize(...)");
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
resize(row_len, col_len);
|
|||
|
|
|
|||
|
|
for (size_t i = 0; i < row_len; i++)
|
|||
|
|
{
|
|||
|
|
for (size_t j = 0; j < col_len; j++)
|
|||
|
|
{
|
|||
|
|
val[i][j] = b[j + i*col_len];
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
void matrix<MatValType>::resize(const matrix<MatValType> &b)
|
|||
|
|
{
|
|||
|
|
if (!b.empty())
|
|||
|
|
{
|
|||
|
|
resize(b.row_size(), b.col_size());
|
|||
|
|
|
|||
|
|
for (size_t i = 0; i < row_length; i++)
|
|||
|
|
{
|
|||
|
|
for (size_t j = 0; j < col_length; j++)
|
|||
|
|
{
|
|||
|
|
val[i][j] = b[i][j];
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
void matrix<MatValType>::clear()
|
|||
|
|
{
|
|||
|
|
row_length = col_length = 0;
|
|||
|
|
|
|||
|
|
if (val != nullptr)
|
|||
|
|
{
|
|||
|
|
delete[] val[0];
|
|||
|
|
delete[] val;
|
|||
|
|
val = nullptr;
|
|||
|
|
}
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
void matrix<MatValType>::random(MatValType np1, MatValType np2, random_type_e mode, unsigned int seed)
|
|||
|
|
{
|
|||
|
|
if (empty())
|
|||
|
|
{
|
|||
|
|
throw gctl::runtime_error("The matrix is not initialized. From gctl::matrix<T>::random(...)");
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
if (seed == 0) seed = std::chrono::system_clock::now().time_since_epoch().count();
|
|||
|
|
std::default_random_engine generator(seed);
|
|||
|
|
|
|||
|
|
if (mode == RdNormal)
|
|||
|
|
{
|
|||
|
|
std::normal_distribution<MatValType> dist(np1, np2); // 添加高斯分布的随机值
|
|||
|
|
|
|||
|
|
for (size_t i = 0; i < row_length; i++)
|
|||
|
|
{
|
|||
|
|
for (size_t j = 0; j < col_length; j++)
|
|||
|
|
{
|
|||
|
|
val[i][j] = dist(generator);
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
std::uniform_real_distribution<MatValType> dist(np1, np2); // /添加均匀分布的随机值
|
|||
|
|
|
|||
|
|
for (size_t i = 0; i < row_length; i++)
|
|||
|
|
{
|
|||
|
|
for (size_t j = 0; j < col_length; j++)
|
|||
|
|
{
|
|||
|
|
val[i][j] = dist(generator);
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
void matrix<MatValType>::assign_all(MatValType in_val)
|
|||
|
|
{
|
|||
|
|
for (size_t i = 0; i < row_length; i++)
|
|||
|
|
{
|
|||
|
|
for (size_t j = 0; j < col_length; j++)
|
|||
|
|
{
|
|||
|
|
val[i][j] = in_val;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
void matrix<MatValType>::scale(MatValType s)
|
|||
|
|
{
|
|||
|
|
for (size_t i = 0; i < row_length; i++)
|
|||
|
|
{
|
|||
|
|
for (size_t j = 0; j < col_length; j++)
|
|||
|
|
{
|
|||
|
|
val[i][j] *= s;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
void matrix<MatValType>::mean(array<MatValType> &m, matrix_order_e odr)
|
|||
|
|
{
|
|||
|
|
if (odr == RowMajor)
|
|||
|
|
{
|
|||
|
|
m.resize(col_length);
|
|||
|
|
for (size_t j = 0; j < col_length; j++)
|
|||
|
|
{
|
|||
|
|
m[j] = 0.0;
|
|||
|
|
for (size_t i = 0; i < row_length; i++)
|
|||
|
|
{
|
|||
|
|
m[j] += val[i][j];
|
|||
|
|
}
|
|||
|
|
m[j] /= row_length;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
m.resize(row_length);
|
|||
|
|
for (size_t i = 0; i < row_length; i++)
|
|||
|
|
{
|
|||
|
|
m[i] = 0.0;
|
|||
|
|
for (size_t j = 0; j < col_length; j++)
|
|||
|
|
{
|
|||
|
|
m[i] += val[i][j];
|
|||
|
|
}
|
|||
|
|
m[i] /= col_length;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
double matrix<MatValType>::norm(norm_type_e n_type) const
|
|||
|
|
{
|
|||
|
|
MatValType nval = (MatValType) 0;
|
|||
|
|
if (n_type == L0)
|
|||
|
|
{
|
|||
|
|
for (size_t i = 0; i < row_length; i++)
|
|||
|
|
{
|
|||
|
|
for (size_t j = 0; j < col_length; j++)
|
|||
|
|
{
|
|||
|
|
if (val[i][j] != (MatValType) 0)
|
|||
|
|
{
|
|||
|
|
nval += (MatValType) 1;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
return nval;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
if (n_type == L1)
|
|||
|
|
{
|
|||
|
|
for (size_t i = 0; i < row_length; i++)
|
|||
|
|
{
|
|||
|
|
for (size_t j = 0; j < col_length; j++)
|
|||
|
|
{
|
|||
|
|
nval += GCTL_FABS(val[i][j]);
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
return nval;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
if (n_type == L2)
|
|||
|
|
{
|
|||
|
|
for (size_t i = 0; i < row_length; i++)
|
|||
|
|
{
|
|||
|
|
for (size_t j = 0; j < col_length; j++)
|
|||
|
|
{
|
|||
|
|
nval += val[i][j]*val[i][j];
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
return sqrt(nval);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
nval = val[0][0];
|
|||
|
|
for (size_t i = 0; i < row_length; i++)
|
|||
|
|
{
|
|||
|
|
for (size_t j = 0; j < col_length; j++)
|
|||
|
|
{
|
|||
|
|
nval = GCTL_MAX(nval, val[i][j]);
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
return sqrt(nval);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
void matrix<MatValType>::sequent(MatValType st_val, MatValType row_inc, MatValType col_inc)
|
|||
|
|
{
|
|||
|
|
for (size_t i = 0; i < row_length; i++)
|
|||
|
|
{
|
|||
|
|
for (size_t j = 0; j < col_length; j++)
|
|||
|
|
{
|
|||
|
|
val[i][j] = st_val + i*row_inc + j*col_inc;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
MatValType *matrix<MatValType>::get(size_t row_index) const
|
|||
|
|
{
|
|||
|
|
#ifdef GCTL_BOUND_CHECK
|
|||
|
|
if (row_index >= row_length || col_index >= col_length)
|
|||
|
|
throw std::out_of_range("Invalid index. gctl::matrix<T>::get(...)");
|
|||
|
|
#endif // GCTL_BOUND_CHECK
|
|||
|
|
|
|||
|
|
return val[row_index];
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
MatValType &matrix<MatValType>::at(size_t row_index, size_t col_index)
|
|||
|
|
{
|
|||
|
|
#ifdef GCTL_BOUND_CHECK
|
|||
|
|
if (row_index >= row_length || col_index >= col_length)
|
|||
|
|
throw std::out_of_range("Invalid index. gctl::matrix<T>::at(...)");
|
|||
|
|
#endif // GCTL_BOUND_CHECK
|
|||
|
|
|
|||
|
|
return val[row_index][col_index];
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
MatValType &matrix<MatValType>::at(size_t row_index, size_t col_index) const
|
|||
|
|
{
|
|||
|
|
#ifdef GCTL_BOUND_CHECK
|
|||
|
|
if (row_index >= row_length || col_index >= col_length)
|
|||
|
|
throw std::out_of_range("Invalid index. gctl::matrix<T>::at(...)");
|
|||
|
|
#endif // GCTL_BOUND_CHECK
|
|||
|
|
|
|||
|
|
return val[row_index][col_index];
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
MatValType *matrix<MatValType>::operator[](size_t index)
|
|||
|
|
{
|
|||
|
|
return val[index];
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
MatValType *matrix<MatValType>::operator[](size_t index) const
|
|||
|
|
{
|
|||
|
|
return val[index];
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
bool matrix<MatValType>::empty() const
|
|||
|
|
{
|
|||
|
|
if (row_length == 0 && col_length == 0) return true;
|
|||
|
|
else return false;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
size_t matrix<MatValType>::row_size() const
|
|||
|
|
{
|
|||
|
|
return row_length;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
size_t matrix<MatValType>::col_size() const
|
|||
|
|
{
|
|||
|
|
return col_length;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
size_t matrix<MatValType>::val_size() const
|
|||
|
|
{
|
|||
|
|
return row_length*col_length;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
void matrix<MatValType>::copy_to(matrix<MatValType> &b) const
|
|||
|
|
{
|
|||
|
|
b.resize(row_length, col_length);
|
|||
|
|
|
|||
|
|
for (size_t i = 0; i < row_length; i++)
|
|||
|
|
{
|
|||
|
|
for (size_t j = 0; j < col_length; j++)
|
|||
|
|
{
|
|||
|
|
b[i][j] = val[i][j];
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
void matrix<MatValType>::copy_to(std::vector<std::vector<MatValType> > &b) const
|
|||
|
|
{
|
|||
|
|
if (!b.empty())
|
|||
|
|
{
|
|||
|
|
for (size_t i = 0; i < b.size(); i++)
|
|||
|
|
{
|
|||
|
|
b[i].clear();
|
|||
|
|
}
|
|||
|
|
b.clear();
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
b.resize(row_length);
|
|||
|
|
for (size_t i = 0; i < row_length; i++)
|
|||
|
|
{
|
|||
|
|
b[i].resize(col_length);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
for (size_t i = 0; i < row_length; i++)
|
|||
|
|
{
|
|||
|
|
for (size_t j = 0; j < col_length; j++)
|
|||
|
|
{
|
|||
|
|
b[i][j] = val[i][j];
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
void matrix<MatValType>::extract(array<MatValType> &b, size_t index, matrix_order_e order) const
|
|||
|
|
{
|
|||
|
|
if (order == RowMajor)
|
|||
|
|
{
|
|||
|
|
if (index >= row_length)
|
|||
|
|
{
|
|||
|
|
throw std::out_of_range("Invalid index. gctl::matrix<T>::extract(...)");
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
b.resize(col_length);
|
|||
|
|
for (size_t i = 0; i < col_length; i++)
|
|||
|
|
{
|
|||
|
|
b[i] = val[index][i];
|
|||
|
|
}
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
if (index >= col_length)
|
|||
|
|
{
|
|||
|
|
throw std::out_of_range("Invalid index. Thrown by gctl::matrix<T>::extract(...)");
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
b.resize(row_length);
|
|||
|
|
for (size_t i = 0; i < row_length; i++)
|
|||
|
|
{
|
|||
|
|
b[i] = val[i][index];
|
|||
|
|
}
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
void matrix<MatValType>::reform(array<MatValType> &b, matrix_order_e order) const
|
|||
|
|
{
|
|||
|
|
b.resize(row_length*col_length);
|
|||
|
|
|
|||
|
|
if (order == RowMajor)
|
|||
|
|
{
|
|||
|
|
for (size_t i = 0; i < row_length; i++)
|
|||
|
|
{
|
|||
|
|
for (size_t j = 0; j < col_length; j++)
|
|||
|
|
{
|
|||
|
|
b[j+i*col_length] = val[i][j];
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
for (size_t j = 0; j < col_length; j++)
|
|||
|
|
{
|
|||
|
|
for (size_t i = 0; i < row_length; i++)
|
|||
|
|
{
|
|||
|
|
b[i+j*row_length] = val[i][j];
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
void matrix<MatValType>::for_each(foreach_m_ptr func)
|
|||
|
|
{
|
|||
|
|
for (size_t i = 0; i < row_length; i++)
|
|||
|
|
{
|
|||
|
|
for (size_t j = 0; j < col_length; j++)
|
|||
|
|
{
|
|||
|
|
func(&val[i][j], i, j);
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
template <typename MatValType>
|
|||
|
|
void matrix<MatValType>::show(std::ostream &os, char sep) const
|
|||
|
|
{
|
|||
|
|
for (size_t i = 0; i < row_length; i++)
|
|||
|
|
{
|
|||
|
|
for (size_t j = 0; j < col_length; j++)
|
|||
|
|
{
|
|||
|
|
os << val[i][j] << sep;
|
|||
|
|
}
|
|||
|
|
os << std::endl;
|
|||
|
|
}
|
|||
|
|
return;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
#endif // _GCTL_MATRIX_H
|