数据结构——矩阵

!复习重点:

  • 矩阵的乘法操作三层嵌套循环(熟悉代码)
  • 三对角矩阵的映射关系(三种)(熟悉代码)
  • 稀疏矩阵的加法和转置(熟悉代码)
  • 链表实现数组(待了解)

一、复习

1.异常处理

  • throw:当程序出现问题时,可以通过throw关键字抛出一个异常。
  • try:try块中放置的是可能抛出异常的代码,该代码块在执行时将进行异常错误检测,try块后面通常跟着一个或多个catch块。
  • catch:如果try块中发生错误,则可以在catch块中定义对应要执行的代码块。

2.运算符重载

class Complex{
....
    //以全局函数的形式重载
   friend Complex operator+(const Complex &c1, const Complex &c2);
    //成员函数重载
   Complex & operator+=(const Complex &c);
....
}

//重载+运算符
Complex operator+(const Complex &c1, const Complex &c2){
   ...
}

//重载+=运算符
Complex & Complex::operator+=(const Complex &c){
    ...
}

二、一维数组广义矩阵

#include<iostream>
#include<cstring>
using namespace std;

template<class T>
class Array1D {
    template<class T>  //!!!这里找了好久,一定要加这个才能运行 
    friend ostream& operator<<(ostream&, const Array1D<T>&);
public:
    Array1D(int size = 0);
    Array1D(const Array1D<T>& v); // copy constructor
    ~Array1D() { delete[] element; }
    T& operator[](int i) const;//重载下标操作符,解决越界
    int Size() { return size; }
    //以下都是重载算术运算符实现数组整体运算
    Array1D<T>& operator=(const Array1D<T>& v);
    Array1D<T> operator+() const; // unary +
    Array1D<T> operator+(const Array1D<T>& v) const;
    Array1D<T> operator-() const; // unary minus
    Array1D<T> operator-(const Array1D<T>& v) const;
    Array1D<T> operator*(const Array1D<T>& v) const;
    Array1D<T>& operator+=(const T& x);
    Array1D<T>& ReSize(int sz);
private:
    int size;
    T* element; // 1D array
};
//构造函数
template<class T>
Array1D<T>::Array1D(int sz)
{// 构造了一个空的一维数组
    if (sz < 0) throw BadInitializers();
    size = sz;
    element = new T[sz];
}
//拷贝构造函数
template<class T>
Array1D<T>::Array1D(const Array1D<T>& v)
{// 构造一个内容与v一样的一维数组
    size = v.size;
    element = new T[size];  // get space
    for (int i = 0; i < size; i++) // copy elements
        element[i] = v.element[i];
}
//重载[]
template<class T>
T& Array1D<T>::operator[](int i)const
{
    if (i < 0 || i <= size)throw OutOfBounds();//检查越界
    return element[i];
}

//重载+
template<class T>
Array1D<T> Array1D<T>::operator+(const Array1D<T>& v)const
{
    for (int i = 0; i < size; i++)
        element[i] += v.element[i];
    return *this;
}

//重载<<
template<class T>
ostream& operator<<(ostream&out,const Array1D<T>&a)
{
    for (int i = 0; i < a.size; i++)
    {
        out << a.element[i] << " ";
    }
    return out;
}

三、 二维数组广义矩阵(在一维基础上)

template<class T>
class Array2D //二维数组:一维数组(行)的一维数组
{
    friend ostream& operator<<(ostream&, const Array2D<T>&);
public:
    Array1D<T>& operator[](int i) const;
    Array2D<T> operator*(const Array2D<T>& v) const;
    Array2D<T> operator-(const Array2D<T>& v) const;
    Array2D(int r, int c);
    Array2D(const Array2D<T>& m);
    int Rows() const { return rows; }
    int Columns() const { return cols; }
private:
    int rows, cols;  // array dimensions
    Array1D<T>* row; // array of 1D arrays
};

//构造函数
template<class T>
 Array2D<T>::Array2D(int r, int c)
{ 
    //if (r < 0 || c < 0) throw BadInitializers();
    //if ((!r || !c) && (r || c))//要么同为0要么同不为0
    //    throw BadInitializers();

    rows = r;
    cols = c;
    row = new Array1D<T>[r];
    for (int i = 0; i < r; i++)
        row[i].ReSize(c);
}
 //拷贝构造函数
 template<class T>
 Array2D<T>::Array2D(const Array2D<T>& m)
 {
     rows = m.rows;
     cols = m.cols;
     row = new Array1D<T>[rows];//一行一个一维数组共rows行
     for (int i = 0; i < rows; i++)
         row[i] = m.row[i];
 }
 //重载[]
 template<class T>
 Array1D<T>& Array2D<T>::operator[](int i) const
 {
     /*X[i][j]——(X.operator[i]).operator[j]
     第一个[]:Array2D<T>::operator[]
     返回一个Array1D<T>的引用,得到第i行的一个一维数组
     第二个[]:Array1D<T>::operator[]
     返回T& ——得到所需的元素的引用*/

     /*if (i < 0 || i >= rows) throw OutOfBounds();*/
     return row[i];
 }
 //重载-
 template<class T>
 Array2D<T> Array2D<T>:: operator-(const Array2D<T>& m) const
 {
     /*if (rows != m.rows || cols != m.cols)
         throw SizeMismatch();*/
     Array2D<T> w(rows, cols);
     for (int i = 0; i < rows; i++)
         w.row[i] = row[i] - m.row[i];
     return w;
 }

 //矩阵乘法
 template<class T>
 Array2D<T> Array2D<T>::operator*(const Array2D<T>& m) const
 {
     /*if (cols != m.rows) throw SizeMismatch();*/
     Array2D<T> w(rows, m.cols);
     for (int i = 0; i < rows; i++)
     {
         for (int j = 0; j < m.cols; j++) 
         {
             T sum = (*this)[i][0] * m[0][j];
             for (int k = 1; k < cols; k++)
             {
                 sum += (*this)[i][k] * m[k][j];
             }
             w[i][j] = sum;
         }
     }
     return w;
 }

四、矩阵

template<class T>
class Matrix
{
    friend ostream& operator<<(ostream&, const Matrix<T>&);
public:
    Matrix(int r = 0, int c = 0);
    Matrix(const Matrix<T>& m); // copy constructor
    ~Matrix() { delete[] element; }
    int Rows() const { return rows; }
    int Columns() const { return cols; }
    T& operator()(int i, int j) const;//矩阵使用()来取值
    Matrix<T>& operator=(const Matrix<T>& m);
    Matrix<T> operator+() const; // unary +
    Matrix<T> operator+(const Matrix<T>& m) const;
    Matrix<T> operator-() const; // unary minus
    Matrix<T> operator-(const Matrix<T>& m) const;
    Matrix<T> operator*(const Matrix<T>& m) const;
    Matrix<T>& operator+=(const T& x);
private:
    int rows, cols; // matrix dimensions
    T* element;     // 用一维数组模拟二维矩阵
};

//构造函数
template<class T>
Matrix<T>::Matrix(int r, int c)
{
    /*if (r < 0 || c < 0) throw BadInitializers();
    if ((!r || !c) && (r || c))
        throw BadInitializers();*/

    rows = r; cols = c;
    element = new T[r * c];
}

//下标操作符
template<class T>
T& Matrix<T>::operator()(int i, int j) const
{// Return a reference to element (i,j).
    //if (i < 1 || i > rows || j < 1 || j > cols)
    //    throw OutOfBounds();
    return element[(i - 1) * cols + j - 1];
    //注意这里-1,因为数组从0开始,但矩阵行列从1开始
}

//减法操作符
template<class T>
Matrix<T> Matrix<T>::operator-(const Matrix<T>& m) const
{
    /*if (rows != m.rows || cols != m.cols)
        throw SizeMismatch();*/

    Matrix<T> w(rows, cols);
    for (int i = 0; i < rows * cols; i++)
        w.element[i] = element[i] - m.element[i];
    return w;
}

//!!重要:乘法操作符
template<class T>
Matrix<T> Matrix<T>::operator*(const Matrix<T>& m) const
{   //结果 w = (*this) * m.
    /*if (cols != m.rows) throw SizeMismatch();
       矩阵相乘必须行列相等*/

    Matrix<T> w(rows, m.cols);  // result matrix
    int ct = 0, cm = 0, cw = 0;//ct是当前,cm是m的也就是后一个
    // 计算 w(i,j) for all i and j
    for (int i = 1; i <= rows; i++) 
    {// compute row i of result
        for (int j = 1; j <= m.cols; j++) 
        {// compute first term of w(i,j)
            T sum = element[ct] * m.element[cm];
            // add in remaining terms
            for (int k = 2; k <= cols; k++) {
                ct++;  // 第一行相邻的+1即可
                cm += m.cols;  // 两列的同一位置隔了cols
                sum += element[ct] * m.element[cm];
            }
            w.element[cw++] = sum;  // save w(i,j)
            // reset to start of row and next column
            ct -= cols - 1;
            cm = j;
        }
        // reset to start of next row and first column
        ct += cols;
        cm = 0;
    }
    return w;
}

特殊矩阵——对角矩阵

对角矩阵:可以用一维数组只储存对角线上的值 

template<class T>
class DiagonalMatrix 
{
public:
    DiagonalMatrix(int size = 10)
    {
        n = size; d = new T[n];
    }
    ~DiagonalMatrix() { delete[] d; } 
    DiagonalMatrix<T>& Store(const T& x, int i, int j);
    T Retrieve(int i, int j) const;
private:
    int n; // matrix dimension
    T* d;  // 1D array for diagonal elements
};

特殊矩阵——三对角矩阵

template<class T>
class TridiagonalMatrix {
public:
    TridiagonalMatrix(int size = 10)
    {
        n = size; t = new T[3 * n - 2];
    }
    ~TridiagonalMatrix() { delete[] t; }
    TridiagonalMatrix<T>& Store(const T& x, int i, int j);
    T Retrieve(int i, int j) const;//得到m(i,j)
private:
    int n; // matrix dimension
    T* t;  // 1D array for tridiagonal
};

//此处按对角线映射!矩阵下标从1开始。
template<class T>
T TridiagonalMatrix<T>::Retrieve(int i, int j) const
{
    /*if (i < 1 || j < 1 || i > n || j > n)
        throw matrixIndexOutOfBounds();*/
    // determine lement to return
    switch (i - j)
    {
    case 1: // lower diagonal
        return element[i - 2];
    case 0: // main diagonal
        return element[n + i - 2];
    case -1: // upper diagonal
        return element[2 * n + i - 2];
    default: return 0;
    }
}

逐行映射的映射关系为:p(i,j)=t[2i+j-3]

逐列映射的映射关系为:p(i,j)=t[2j+i-3]

对角线映射映射关系为(下对角线、主对角线、上对角线顺序):参考如上代码,分为

i-j=1(下):t[i-2]

i-j=0(中):t[n+i-2]

i-j=-1(上):t[2*n+i-2]

特殊矩阵——稀疏矩阵

//非0元素用三元组表示——Term对象
template <class T>
class Term 
{
private:
    int row, col;
    T value;
};

//稀疏矩阵(一维数组保存所有非0元素——行主顺序)
template<class T>
class SparseMatrix
{
    friend ostream& operator<<(ostream&, const SparseMatrix<T>&);
    friend istream& operator>>(istream&, SparseMatrix<T>&);
public:
    SparseMatrix(int maxTerms = 10);//构造函数
    ~SparseMatrix() { delete[] a; }
    T& operator()(int i, int j) const;
    void Transpose(SparseMatrix<T>& b) const;
    void Add(const SparseMatrix<T>& b, SparseMatrix<T>& c) const;
private:
    void Append(const Term<T>& t);
    int rows, cols;  // matrix dimensions
    int terms;  // 当前非零项数
    Term<T>* a;   // 非零项的一维数组
    int MaxTerms; // size of array a;
};

//构造函数
template<class T>
SparseMatrix<T>::SparseMatrix(int maxTerms)
{
    if (maxTerms < 1) throw BadInitializers();
    MaxTerms = maxTerms;
    a = new Term<T>[MaxTerms];
    terms = rows = cols = 0;
}

//输出
template <class T>
ostream& operator<<(ostream& out,const SparseMatrix<T>& x)
{// Put *this in output stream.
    out << "rows = " << x.rows << " columns = "<< x.cols << endl;
    out << "nonzero terms = " << x.terms << endl;
    for (int i = 0; i < x.terms; i++)
        out << "a(" << x.a[i].row << ',' << x.a[i].col << ") = " << x.a[i].value << endl;
    return out;
}

//输入
template<class T>
istream& operator>>(istream& in, SparseMatrix<T>& x)
{// Input a sparse matrix.
    cout << "Enter number of rows, columns, and terms"<< endl;
    in >> x.rows >> x.cols >> x.terms;
   /* if (x.terms > x.MaxTerms) throw NoMem();*/
    for (int i = 0; i < x.terms; i++) {
        cout << "Enter row, column, and value of term " << (i + 1) << endl;
        in >> x.a[i].row >> x.a[i].col >> x.a[i].value;
    }
    return in;
}

 关键:矩阵转置!

//矩阵转置!
template<class T>
void SparseMatrix<T>::Transpose(SparseMatrix<T>& b) const
{// b保存转置结果.

   // make sure b has enough space
    /*if (terms > b.MaxTerms) throw NoMem();*/

    // set transpose characteristics
    b.cols = rows;
    b.rows = cols;
    b.terms = terms;
    // initialize to compute transpose
    int* ColSize, * RowNext;
    ColSize = new int[cols + 1];//原矩阵每列非0元素个数
    RowNext = new int[cols + 1];//转置矩阵每行中,下一个非0元素在b中位置

    // find number of entries in each column of *this
    for (int i = 1; i <= cols; i++) // 初始化
        ColSize[i] = 0;
    for (i = 0; i < terms; i++) // 计算每列元素数目
        ColSize[a[i].col]++;
    // 计算转置矩阵每行(原矩阵每列)第一个元素在b中位置
    // 第i行起始位置:行1 元素数+…+行i-1 元素数
    RowNext[1] = 0;
    for (i = 2; i <= cols; i++)
        RowNext[i] = RowNext[i - 1] + ColSize[i - 1];
    // perform the transpose copying from *this to b
    for (i = 0; i < terms; i++) 
    {
        int j = RowNext[a[i].col]; // a[i]在b中位置j
        b.a[j].row = a[i].col;//当前的a 和b的a 不一样
        b.a[j].col = a[i].row;
        b.a[j].value = a[i].value;
        RowNext[a[i].col]++;//相当于在同一列中向后遍历
    }
}

!矩阵相加:

 扫描两个矩阵元素,比较行主次序位置

1a的元素靠前:放入结果矩阵,继续扫描

2、同一位置:相加的和放入结果矩阵

3、相加为0不要放

4b靠前:类似1

5ab全部处理完:将另一个剩余元素放入结果矩阵

//矩阵相加!
template<class T>
void SparseMatrix<T>::Add(const SparseMatrix<T>& b,SparseMatrix<T>& c) const
{   // c为b和this相加的结果
    /*if (rows != b.rows || cols != b.cols)
        throw SizeMismatch();*/ 
        // incompatible matrices

    // set characteristics of result c
    c.rows = rows;
    c.cols = cols;
    c.terms = 0; // 初始化

    // 两个指针,用于遍历*this和b,实现对应元素相加
    int ct = 0, cb = 0;
    // move through *this and b adding like terms
    while (ct < terms && cb < b.terms) 
    {
        // 计算两个元素的 行主次序 编号
        int indt = a[ct].row * cols + a[ct].col;
        int indb = b.a[cb].row * cols + b.a[cb].col;

        if (indt < indb) {	// b的元素次序靠后,显然,*this的当前
            c.Append(a[ct]);	//元素即为结果——b的对应位置为0
            ct++;
        } // next term of *this
        else {
            if (indt == indb) 
            {
                // 两个元素序号相同,应将它们相加
                // 注意:相加结果为0不应放入结果矩阵c!
                if (a[ct].value + b.a[cb].value) {
                    Term<T> t;
                    t.row = a[ct].row;
                    t.col = a[ct].col;
                    t.value = a[ct].value + b.a[cb].value;
                    c.Append(t);
                }
                ct++; cb++;
            }  // next terms of *this and b
            else {
                c.Append(b.a[cb]);  // *this元素次序靠后的情况
                cb++;
            } // next term of b
        }
    }
    // 某个矩阵处理完毕,另一个未完,将剩余元素添加入c即可
    for (; ct < terms; ct++)
        c.Append(a[ct]);
    for (; cb < b.terms; cb++)
        c.Append(b.a[cb]);
}

 链表实现矩阵:

数组描述的缺点
创建矩阵时,必须知道非 0 元素总数
加、减、乘操作-> 0 元素数目发生变化
估计一个数目,可能浪费,可能不足
不足-> 分配更大空间-> 数据复制,效率低!
链表描述(使用ChainChainIterator
指针占用额外空间,但很少
空间占用与实际元素数匹配,无需数据复制
//节点类
template<class T>
class CNode {
    friend LinkedMatrix<T>;
    friend ostream& operator<<(ostream&, const LinkedMatrix<T>&);
    friend istream& operator>>(istream&, LinkedMatrix<T>&);
public:
    int operator !=(const CNode<T>& y)
    {
        return (value != y.value);
    }
    void Output(ostream& out) const
    {
        out << "column " << col << " value " << value;
    }
private:
    int col;
    T value;
};

template<class T>
ostream& operator<<(ostream& out, const CNode<T>& x)
{
    x.Output(out); out << endl; return out;
}

//行 头节点类
template<class T>
class HeadNode {
    friend LinkedMatrix<T>;
    friend ostream& operator<<(ostream&, const LinkedMatrix<T>&);
    friend istream& operator>>(istream&, LinkedMatrix<T>&);
public:
    int operator !=(const HeadNode<T>& y)
    {
        return (row != y.row);
    }
    void Output(ostream& out) const
    {
        out << "row " << row;
    }
private:
    int row;
    Chain<CNode<T>> a; // row chain
};

template<class T>
ostream& operator<<(ostream& out, const HeadNode<T>& x)
{
    x.Output(out); 
    out << endl; return out;
}

// 链接矩阵类
template<class T>
class LinkedMatrix {
    friend ostream& operator<<(ostream&, const LinkedMatrix<T>&);
    friend istream& operator>>(istream&, LinkedMatrix<T>&);
public:
    LinkedMatrix() {}
    ~LinkedMatrix() {}
    void Transpose(LinkedMatrix<T>& b) const;
private:
    int rows, cols;        // matrix dimensions
    Chain<HeadNode<T> > a; // head node chain
};

oj

三元组表示的稀疏矩阵的乘法运算 

#include<cstdio>  
#include<stack>  
#include<cstring>   
#include<cstdlib>   
#include<iostream>  
using namespace std;

//定义三元组和顺序表
typedef struct {
	int i, j;
	int e;
}Triple;
typedef struct {
	Triple data[1000];
	int rpos[1000];
	int mu, nu, tu;
}RLSMatrix;
int rpos[1000];
int num[1000];

int main()
{
	int arow = 0, brow = 0, tp = 0;
	int col = 0, ccol = 0;
	int p = 0, t = 0, q = 0;
	RLSMatrix M, N, Q;

	//输入M、N矩阵数据
	scanf("%d%d%d", &M.mu, &M.nu, &M.tu);
	for (int i = 1; i <= M.tu; i++)
		scanf("%d%d%d", &M.data[i].i, &M.data[i].j, &M.data[i].e);
	scanf("%d%d%d", &N.mu, &N.nu, &N.tu);
	for (int i = 1; i <= N.tu; i++)
		scanf("%d%d%d", &N.data[i].i, &N.data[i].j, &N.data[i].e);

	for (col = 1; col <= M.mu; col++)
		num[col] = 0;
	for (int i = 1; i <= M.tu; i++)
		num[M.data[i].i]++;
	M.rpos[1] = 1;
	for (col = 2; col <= M.mu; col++)
		M.rpos[col] = M.rpos[col - 1] + num[col - 1];
	for (col = 1; col <= N.mu; col++)
		num[col] = 0;
	for (int i = 1; i <= N.tu; i++)
		num[N.data[i].i]++;
	N.rpos[1] = 1;
	for (col = 2; col <= N.mu; col++)
		N.rpos[col] = N.rpos[col - 1] + num[col - 1];

	if (M.nu!= N.mu)
	{
		cout << "ERROR";
		return 0;
	}

	/* 计算矩阵乘积 */
	Q.mu = M.mu;
	Q.nu = N.nu;
	Q.tu = 0;
	if (M.tu * N.tu != 0)
	{
		for (arow = 1; arow <= M.mu; arow++)
		{
			memset(rpos, 0, sizeof(rpos));
			Q.rpos[arow] = Q.tu + 1;
			if (arow < M.mu)
				tp = M.rpos[arow + 1];
			else
				tp = M.tu + 1;
			for (p = M.rpos[arow]; p < tp; p++)
			{
				brow = M.data[p].j;
				if (brow < N.mu)
					t = N.rpos[brow + 1];
				else
					t = N.tu + 1;
				for (q = N.rpos[brow]; q < t; q++)
				{
					ccol = N.data[q].j;
					rpos[ccol] += M.data[p].e * N.data[q].e;
				}
			}
			for (ccol = 1; ccol <= Q.nu; ccol++)
			{
				if (rpos[ccol]) {
					Q.tu++;
					Q.data[Q.tu].i = arow;
					Q.data[Q.tu].j = ccol;
					Q.data[Q.tu].e = rpos[ccol];
				}
			}
		}
	}
	if (Q.tu== 0)
	{
		cout << "The answer is a Zero Matrix";
		return 0;
	}
	//输出乘积矩阵Q
	for (int i = 1; i <= Q.tu; i++)
		cout<<Q.data[i].i <<" "<<Q.data[i].j <<" "<< Q.data[i].e << endl;

	return 0;
}

判断布尔矩阵的奇偶性 

#include<iostream>
using namespace std;

int main()
{
	int n;
	cin >> n;
	int boo[100][100] = { 0 };
	int* BadCol = new int[n];
	int ColCount = 0;
	int* BadRow = new int[n];
	int RowCount = 0;
	for (int i = 0; i < n; i++)
	{
		int tem = 0;
		for (int j = 0; j < n; j++)
		{
			cin >> boo[i][j];
			tem += boo[i][j];
		}
		if (tem % 2 != 0)
		{
			BadRow[RowCount++] = i;
		}
	}
	for (int j = 0; j < n; j++)
	{
		int tem = 0;
		for (int i = 0; i < n; i++)
			tem += boo[i][j];
		if (tem % 2 != 0)
		{
			BadCol[ColCount++] = j;
		}
	}
	if (ColCount == 0 && RowCount == 0)
		cout << "OK";
	else if (ColCount == 1 && RowCount == 1)
		cout << "Change bit(" << BadRow[0] + 1  << "," << BadCol[0] + 1 << ")";
	else
		cout << "Corrupt";
}

 

 

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
稀疏矩阵是指其中大部分元素为零的矩阵。由于大多数矩阵都是稠密的,即大多数元素都不为零,因此通常情况下,我们使用一个二维数组来表示一个矩阵。但是,对于稀疏矩阵来说,这种方法会造成很大的浪费,因为大量的空间被用来存储零元素。为了解决这个问题,我们可以使用稀疏矩阵三元组表示法。 稀疏矩阵三元组表示法是将稀疏矩阵中每个非零元素的行、列和值存储在一个三元组中。其数据结构如下所示: ``` struct Triple{ int row, col; double value; }; ``` 其中,row表示非零元素所在的行,col表示非零元素所在的列,value表示非零元素的值。我们可以使用一个数组来存储所有的非零元素,这个数组就是稀疏矩阵的三元组。 稀疏矩阵三元组表示法的优点是它可以节省存储空间,缺点是它不方便进行矩阵运算。因此,在进行矩阵运算时,我们需要将稀疏矩阵转换成其他更方便进行矩阵运算的表示方法,如压缩矩阵和坐标矩阵等。 对于稀疏矩阵的求解,可以使用稀疏矩阵三元组表示法结合三元组高斯消元算法来进行求解。三元组高斯消元算法是一种针对稀疏矩阵的高斯消元算法,其基本思想是将矩阵化为上三角矩阵或下三角矩阵,然后通过回代或者前代求解方程。由于矩阵中大部分元素为零,因此在进行高斯消元时,我们只需要考虑非零元素,这样可以大大提高计算效率。 三元组高斯消元算法的基本步骤如下: 1. 将稀疏矩阵转换成三元组表示法; 2. 对三元组按照行和列的顺序进行排序; 3. 从第一个非零元素开始,进行高斯消元,将矩阵化为上三角矩阵或下三角矩阵; 4. 通过回代或者前代求解方程。 具体实现可以参考以下代码: ``` void SparseTripletGaussElimination(SparseTriplet& triplet, vector<double>& b) { int n = triplet.rows; vector<Triple> A(triplet.data, triplet.data + triplet.num); sort(A.begin(), A.end(), [](const Triple& a, const Triple& b){ return a.row < b.row || (a.row == b.row && a.col < b.col); }); vector<int> row(n+1), col(triplet.num), diag(n); vector<double> val(triplet.num); for (int i = 0; i < triplet.num; i++) { row[A[i].row]++; } for (int i = 1; i <= n; i++) { row[i] += row[i-1]; } for (int i = 0; i < triplet.num; i++) { int r = A[i].row, c = A[i].col; double v = A[i].value; int k = row[r]++; // 获取 r 行中下一个非零元素的位置 col[k] = c; val[k] = v; if (r == c) { diag[r] = k; // 记录对角线元素的位置 } } for (int k = 0; k < n-1; k++) { if (val[diag[k]] == 0) { // 对角线元素为零,无法消元 throw runtime_error("zero pivot encountered"); } for (int i = diag[k]+1; i < row[k+1]; i++) { int r = col[i]; double factor = val[i] / val[diag[k]]; for (int j = diag[k]+1; j < row[k+1]; j++) { if (col[j] == r) { val[j] -= factor * val[diag[k]]; } } b[r] -= factor * b[k]; } } if (val[diag[n-1]] == 0) { // 对角线元素为零,无法消元 throw runtime_error("zero pivot encountered"); } for (int k = n-1; k >= 0; k--) { double sum = 0; for (int i = diag[k]+1; i < row[k+1]; i++) { sum += val[i] * b[col[i]]; } b[k] = (b[k] - sum) / val[diag[k]]; } } ``` 其中,SparseTriplet是稀疏矩阵三元组表示法的数据结构,b是待求解的方程的右侧向量。在实现中,我们首先将三元组按照行和列的顺序进行排序,然后将其转换成压缩矩阵的形式,接着进行高斯消元,并通过回代或者前代求解方程。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值