c++编写的 矩阵 matrix 类源码

#ifndef matrix_cpp #define matrix_cpp #include <iostream.h> #include <matrix.h> #include <math.h> #include "stdlib.h" template <class type> inline void swap(type &a, type &b); matrix::matrix(void) { line_n = 0; column_n = 0; p = NULL; } matrix::matrix(float *mtrx,int ln, int cn) { line_n = ln; column_n = cn; p = new float[line_n * column_n]; if(!p) { cout << "Allocation error!" << endl; exit(1); } for(int i = 0; i < line_n * column_n; i++) p[i] = mtrx[i]; } matrix::matrix(const matrix &obj) { line_n = obj.line_n; column_n = obj.column_n; p = new float[line_n * column_n]; if(!p) { cout << "Allocation error!" << endl; exit(1); } for(int i = 0; i < line_n * column_n; i++) p[i] = obj.p[i]; } ostream &operator<<(ostream &stream, const matrix &obj) { float n; for(int i = 0; i < obj.line_n * obj.column_n; i++) { n = obj.p[i]; if(floor(obj.p[i]) - obj.p[i] < mat::precision && floor(obj.p[i]) - obj.p[i] > -mat::precision) n = floor(obj.p[i]); else if(ceil(obj.p[i]) - obj.p[i] < mat::precision && ceil(obj.p[i]) - obj.p[i] > -mat::precision) n = ceil(obj.p[i]); stream.precision(5); stream.width(12); if(i % obj.column_n == 0) stream << endl; stream << n << ' '; } stream << endl; return stream; } istream &operator>>(istream &stream, const matrix &obj) { } matrix &matrix::operator= (const matrix &obj) { if(this == &obj) return *this; delete []p; p = new float[obj.line_n * obj.column_n]; line_n = obj.line_n; column_n = obj.column_n; if(!p) { cout << "Allocation error!" << endl; exit(1); } for(int i = 0; i < obj.line_n * obj.column_n; i++) { p[i] = obj.p[i]; // p[i] = abs(p[i] - floor(p[i]) < precision) ? floor(p[i]):p[i]; // p[i] = abs(p[i] - ceil(p[i]) < precision) ? ceil(p[i]):p[i]; } return *this; } matrix &matrix::operator= (const float &n) { for(int i =0; i < line_n; i++) for(int j = 0; j < column_n; j++) { p[i * column_n + j] = n; } return *this; } //matrix matrix::operator= (float *mtrx,int ln, int cn) //{ // delete p; // p = new float[ln * cn]; // line_n = ln; // column_n = cn; // if(!p) // { // cout << "Allocation error!" << endl; // exit(1); // } // for(int i = 0; i < line_n * column_n; i++) // p[i] = mtrx[i]; // return *this; // } const matrix matrix::operator+ (const matrix &obj) { matrix temp; if(isempty(obj)) return *this; if(line_n != obj.line_n || column_n != obj.column_n) return temp; temp.p = new float[line_n * column_n]; temp.line_n = line_n ; temp.column_n = column_n; for(int i = 0; i < line_n * column_n; i++) temp.p[i] = this->p[i] + obj.p[i]; return temp; } const matrix matrix::operator- (const matrix &obj) { if(isempty(obj)) return *this; matrix temp; if(line_n != obj.line_n || column_n != obj.column_n) return temp; temp.p = new float[line_n * column_n]; temp.line_n = line_n ; temp.column_n = column_n; for(int i = 0; i < line_n * column_n; i++) temp.p[i] = this->p[i] - obj.p[i]; return temp; } bool matrix::operator==(const matrix &obj) const { if(line_n != obj.line_n || column_n != obj.column_n) return 0; for(int i = 0; i < line_n * column_n;i++) if(p[i] != obj.p[i]) return 0; return 1; } bool matrix::isempty(const matrix &obj) { if(!obj.p) return 1; return 0; } void matrix::print(void) //向屏幕打印矩阵 { for(int i = 0; i < line_n * column_n; i++) { cout.width(6); if(i % column_n == 0) cout << endl; cout << p[i] << ' '; } cout << endl; } void matrix::set(float * mtrx, int ln, int cn) { line_n = ln; column_n = cn; p = new float[line_n * column_n]; if(!p) { cout << "Allocation error!" << endl; exit(1); } for(int i = 0; i < line_n * column_n; i++) p[i] = mtrx[i]; } void matrix::set(const matrix &obj) { line_n = obj.line_n; column_n = obj.column_n; p = new float[line_n * column_n]; if(!p) { cout << "Allocation error!" << endl; exit(1); } for(int i = 0; i < line_n * column_n; i++) p[i] = obj.p[i]; } const matrix matrix::operator* (const matrix &obj) const { matrix temp; int i, j, k; if(column_n != obj.line_n) return temp; temp.p = new float[line_n * obj.column_n]; temp.line_n = line_n ; temp.column_n = obj.column_n; for(i = 0; i < line_n; i++) for(j = 0; j < obj.column_n; j++) { temp.p[i * temp.column_n + j] = 0; for(k = 0; k < column_n; k++) temp.p[i * temp.column_n + j] += p[i * column_n + k] * obj.p[k * obj.column_n + j]; } return temp; } bool matrix::operator %=(const matrix &obj) const { if(line_n == obj.line_n && column_n == obj.column_n) return 1; return 0; } bool matrix::lup_decomposition(matrix &L, matrix &U, matrix &P) const//四个矩阵规模一致 { if(!(*this %= L) || !(*this %= U) || !(*this %= P))//规模一致 return 0; if(line_n != column_n) //需是方阵 return 0; L = U = P = 0; matrix temp; temp = *this; int n = line_n; int max_i, i; float max; int *pi = new int[n]; for(i = 0; i < n; i++) pi[i] = i; for(int k = 0; k < n; k++) { max = 0; for(i = k; i < n; i++)//检索一列中最小的元素 { float c = temp.p[i * column_n + k]; if(c > max || c < -max) { max = c > 0? c : -c; max_i = i; } } if((max - 0) < mat::precision && max - 0 > -mat::precision) { cout << "这是一个奇异矩阵！" << endl; return 0; } swap(pi[k], pi[max_i]); for(i = 0; i < n; i++) { swap(temp.p[k * column_n + i], temp.p[max_i * column_n + i]); swap(L.p[k * column_n + i], L.p[max_i * column_n + i]); //要跟着交换 swap(U.p[k * column_n + i], U.p[max_i * column_n + i]); } U.p[k * column_n + k] = temp.p[k * column_n + k]; L.p[k * column_n + k] = 1; for(i = k + 1; i < n; i++) { L.p[i * column_n + k] = temp.p[i * column_n + k] / temp.p[k * column_n + k]; temp.p[i * column_n + k] = temp.p[i * column_n + k] / temp.p[k * column_n + k]; U.p[k * column_n + i] = temp.p[k * column_n + i]; for(int j = k + 1; j < n; j++) temp.p[i * column_n + j] = temp.p[i * column_n + j] - temp.p[i * column_n + k] * temp.p[k * column_n + j]; } } for(int i = 0; i < n; i++) P.p[i * column_n + pi[i]] = 1; delete []pi; return 1; } const matrix matrix::lup_solve(const matrix &L, const matrix &U, const matrix &P, const matrix &b) const { //cout <<"///////////////////////////////////" << endl; matrix b1; b1 = P * b; matrix x; if(!(L %= U) || !(L %= P)) return x; if(L.line_n != L.column_n || b.line_n != L.line_n || b.column_n != 1) //b必须为n*1的矩阵 ，L U P必须是方阵 return x; // cout << "U:" << endl << U; matrix y; x = b; y = b; int n = L.line_n; float sum = 0; int j; for(int i = 0; i < n; i++) { sum = 0; for(j = 0,sum = 0; j < i; j++) sum += L.p[i* L.column_n + j] * y.p[j]; y.p[i] = b1.p[i] - sum; } for(int i = n - 1;