微软公司内部培训程序员资料---求解线性方程组的类

/** 求解线性方程组的类 LEquations* * 周长发编制*/using System;namespace MSAlgorithm
{public class LEquations{/// <summary>/// 系数矩阵/// </summary>private Matrix mtxLECoef;/// <summary>/// 常数矩阵/// </summary>private Matrix mtxLEConst;/// <summary>/// 基本构造函数/// </summary>public LEquations(){ }/// <summary>/// 指定系数和常数构造函数/// </summary>/// <param name="mtxCoef">指定的系数矩阵</param>/// <param name="mtxConst">指定的常数矩阵</param>public LEquations(Matrix mtxCoef, Matrix mtxConst){Init(mtxCoef, mtxConst);}/// <summary>/// 初始化函数/// </summary>/// <param name="mtxCoef">指定的系数矩阵</param>/// <param name="mtxConst">指定的常数矩阵</param>/// <returns>bool 型，初始化是否成功</returns>public bool Init(Matrix mtxCoef, Matrix mtxConst){if (mtxCoef.GetNumRows() != mtxConst.GetNumRows())return false;mtxLECoef = new Matrix(mtxLECoef);mtxConst = new Matrix(mtxConst);return true;}/// <summary>/// 获取系数矩阵/// </summary>/// <returns>Matrix 型，返回系数矩阵</returns>public Matrix GetCoefMatrix(){return mtxLECoef;}/// <summary>/// 获取常数矩阵/// </summary>/// <returns>Matrix 型，返回常数矩阵</returns>public Matrix GetConstMatrix(){return mtxLEConst;}/// <summary>/// 获取方程个数/// </summary>/// <returns>int 型，返回方程组方程的个数</returns>public int GetNumEquations(){return GetCoefMatrix().GetNumRows();}/// <summary>/// 获取未知数个数/// </summary>/// <returns>int 型，返回方程组未知数的个数</returns>public int GetNumUnknowns(){return GetConstMatrix().GetNumColumns();}/// <summary>/// 全选主元高斯消去法/// </summary>/// <param name="mtxResult">Matrix对象，返回方程组的解</param>/// <returns>bool 型，方程组求解是否成功</returns>public bool GetRootsetGauss(Matrix mtxResult){int l, k, i, j, nIs = 0, p, q;double d, t;//方程组的属性，将常数矩阵赋给解矩阵mtxResult.SetValue(mtxLEConst);double[] pDataCoef = mtxLECoef.GetData();double[] pDataConst = mtxResult.GetData();int n = GetNumUnknowns();//临时缓冲区，存放列数int[] pnJs = new int[n];//消元l = 1;for (k = 0; k <= n - 2; k++){d = 0.0;for (i = k; i <= n - 1; i++){for (j = k; j <= n - 1; j++){t = Math.Abs(pDataCoef[i * n + j]);if (t > d){d = t;pnJs[k] = j;nIs = i;}}}if (d == 0.0)l = 0;else{if (pnJs[k] != k){for (i = 0; i <= n - 1; i++){p = i * n + k;q = i * n + pnJs[k];t = pDataCoef[p];pDataCoef[p] = pDataCoef[q];pDataCoef[q] = t;}}if (nIs != k){for (j = k; j <= n - 1; j++){p = k * n + j;q = nIs * n + j;t = pDataCoef[p];pDataCoef[p] = pDataCoef[q];pDataCoef[q] = t;}t = pDataConst[k];pDataConst[k] = pDataConst[nIs];pDataConst[nIs] = t;}}//求解失败if (l == 0)return false;d = pDataCoef[k * n + k];for (j = k + 1; j <= n - 1; j++){p = k * n + j;pDataCoef[p] = pDataCoef[p] / d;}pDataConst[k] = pDataConst[k] / d;for (i = k + 1; i <= n - 1; i++){for (j = k + 1; j <= n - 1; j++){p = i * n + j;pDataCoef[p] = pDataCoef[p] - pDataCoef[i * n + k] * pDataCoef[k * n + j];}pDataConst[i] = pDataConst[i] - pDataCoef[i * n + k] * pDataConst[k];}}//求解失败d = pDataCoef[(n - 1) * n + n - 1];if (d == 0.0){return false;}//求解pDataConst[n - 1] = pDataConst[n - 1] / d;for (i = n - 2; i >= 0; i--){t = 0.0;for (j = i + 1; j <= n - 1; j++)t = t + pDataCoef[i * n + j] * pDataConst[j];pDataConst[i] = pDataConst[i] - t;}//调整解的位置pnJs[n - 1] = n - 1;for (k = n - 1; k >= 0; k--){if (pnJs[k] != k){t = pDataConst[k];pDataConst[k] = pDataConst[pnJs[k]];pDataConst[pnJs[k]] = t;}}return true;}/// <summary>/// 全选主元高斯－约当消去法/// </summary>/// <param name="mtxResult">Matrix对象，返回方程组的解</param>/// <returns>bool 型，方程组求解是否成功</returns>public bool GetRootsetGaussJordan(Matrix mtxResult){int l, k, i, j, nIs = 0, p, q;double d, t;//方程组的属性，将常数矩阵赋给解矩阵mtxResult.SetValue(mtxLEConst);double[] pDataCoef = mtxLECoef.GetData();double[] pDataConst = mtxResult.GetData();int n = GetNumUnknowns();int m = mtxLEConst.GetNumColumns();//临时缓冲区，存放变换的列数int[] pnJs = new int[n];//消元l = 1;for (k = 0; k <= n - 1; k++){d = 0.0;for (i = k; i <= n - 1; i++){for (j = k; j <= n - 1; j++){t = Math.Abs(pDataCoef[i * n + j]);if (t > d){d = t;pnJs[k] = j;nIs = i;}}}if (d + 1.0 == 1.0)l = 0;else{if (pnJs[k] != k){for (i = 0; i <= n - 1; i++){p = i * n + k;q = i * n + pnJs[k];t = pDataCoef[p];pDataCoef[p] = pDataCoef[q];pDataCoef[q] = t;}}if (nIs != k){for (j = k; j <= n - 1; j++){p = k * n + j;q = nIs * n + j;t = pDataCoef[p];pDataCoef[p] = pDataCoef[q];pDataCoef[q] = t;}for (j = 0; j <= m - 1; j++){p = k * m + j;q = nIs * m + j;t = pDataConst[p];pDataConst[p] = pDataConst[q];pDataConst[q] = t;}}}//求解失败if (l == 0)return false;d = pDataCoef[k * n + k];for (j = k + 1; j <= n - 1; j++){p = k * n + j;pDataCoef[p] = pDataCoef[p] / d;}for (j = 0; j <= m - 1; j++){p = k * m + j;pDataConst[p] = pDataConst[p] / d;}for (j = k + 1; j <= n - 1; j++){for (i = 0; i <= n - 1; i++){p = i * n + j;if (i != k)pDataCoef[p] = pDataCoef[p] - pDataCoef[i * n + k] * pDataCoef[k * n + j];}}for (j = 0; j <= m - 1; j++){for (i = 0; i <= n - 1; i++){p = i * m + j;if (i != k)pDataConst[p] = pDataConst[p] - pDataCoef[i * n + k] * pDataConst[k * m + j];}}}//调整for (k = n - 1; k >= 0; k--){if (pnJs[k] != k){for (j = 0; j <= m - 1; j++){p = k * m + j;q = pnJs[k] * m + j;t = pDataConst[p];pDataConst[p] = pDataConst[q];pDataConst[q] = t;}}}return true;}/// <summary>/// 复系数方程组的全选主元高斯消去法/// </summary>/// <param name="mtxCoefImag">系数矩阵的虚部矩阵</param>/// <param name="mtxConstImag">常数矩阵的虚部矩阵</param>/// <param name="mtxResult">Matrix对象，返回方程组解矩阵的实部矩阵</param>/// <param name="mtxResultImag">Matrix对象，返回方程组解矩阵的虚部矩阵</param>/// <returns>bool 型，方程组求解是否成功</returns>public bool GetRootsetGauss(Matrix mtxCoefImag, Matrix mtxConstImag, Matrix mtxResult, Matrix mtxResultImag){int l, k, i, j, nIs = 0, u, v;double p, q, s, d;//方程组的属性，将常数矩阵赋给解矩阵mtxResult.SetValue(mtxLEConst);mtxResultImag.SetValue(mtxConstImag);double[] pDataCoef = mtxLECoef.GetData();double[] pDataConst = mtxResult.GetData();double[] pDataCoefImag = mtxCoefImag.GetData();double[] pDataConstImag = mtxResultImag.GetData();int n = GetNumUnknowns();int m = mtxLEConst.GetNumColumns();//临时缓冲区，存放变换的列数int[] pnJs = new int[n];//消元for (k = 0; k <= n - 2; k++){d = 0.0;for (i = k; i <= n - 1; i++){for (j = k; j <= n - 1; j++){u = i * n + j;p = pDataCoef[u] * pDataCoef[u] + pDataCoefImag[u] * pDataCoefImag[u];if (p > d){d = p;pnJs[k] = j;nIs = i;}}}//求解失败if (d == 0.0)return false;if (nIs != k){for (j = k; j <= n - 1; j++){u = k * n + j;v = nIs * n + j;p = pDataCoef[u];pDataCoef[u] = pDataCoef[v];pDataCoef[v] = p;p = pDataCoefImag[u];pDataCoefImag[u] = pDataCoefImag[v];pDataCoefImag[v] = p;}p = pDataConst[k];pDataConst[k] = pDataConst[nIs];pDataConst[nIs] = p;p = pDataConstImag[k];pDataConstImag[k] = pDataConstImag[nIs];pDataConstImag[nIs] = p;}if (pnJs[k] != k){for (i = 0; i <= n - 1; i++){u = i * n + k;v = i * n + pnJs[k];p = pDataCoef[u];pDataCoef[u] = pDataCoef[v];pDataCoef[v] = p;p = pDataCoefImag[u];pDataCoefImag[u] = pDataCoefImag[v];pDataCoefImag[v] = p;}}v = k * n + k;for (j = k + 1; j <= n - 1; j++){u = k * n + j;p = pDataCoef[u] * pDataCoef[v];q = -pDataCoefImag[u] * pDataCoefImag[v];s = (pDataCoef[v] - pDataCoefImag[v]) * (pDataCoef[u] + pDataCoefImag[u]);pDataCoef[u] = (p - q) / d;pDataCoefImag[u] = (s - p - q) / d;}p = pDataConst[k] * pDataCoef[v];q = -pDataConstImag[k] * pDataCoefImag[v];s = (pDataCoef[v] - pDataCoefImag[v]) * (pDataConst[k] + pDataConstImag[k]);pDataConst[k] = (p - q) / d;pDataConstImag[k] = (s - p - q) / d;for (i = k + 1; i <= n - 1; i++){u = i * n + k;for (j = k + 1; j <= n - 1; j++){v = k * n + j;l = i * n + j;p = pDataCoef[u] * pDataCoef[v];q = pDataCoefImag[u] * pDataCoefImag[v];s = (pDataCoef[u] + pDataCoefImag[u]) * (pDataCoef[v] + pDataCoefImag[v]);pDataCoef[l] = pDataCoef[l] - p + q;pDataCoefImag[l] = pDataCoefImag[l] - s + p + q;}p = pDataCoef[u] * pDataConst[k];q = pDataCoefImag[u] * pDataConstImag[k];s = (pDataCoef[u] + pDataCoefImag[u]) * (pDataConst[k] + pDataConstImag[k]);pDataConst[i] = pDataConst[i] - p + q;pDataConstImag[i] = pDataConstImag[i] - s + p + q;}}u = (n - 1) * n + n - 1;d = pDataCoef[u] * pDataCoef[u] + pDataCoefImag[u] * pDataCoefImag[u];//求解失败if (d == 0.0)return false;//求解p = pDataCoef[u] * pDataConst[n - 1]; q = -pDataCoefImag[u] * pDataConstImag[n - 1];s = (pDataCoef[u] - pDataCoefImag[u]) * (pDataConst[n - 1] + pDataConstImag[n - 1]);pDataConst[n - 1] = (p - q) / d; pDataConstImag[n - 1] = (s - p - q) / d;for (i = n - 2; i >= 0; i--){for (j = i + 1; j <= n - 1; j++){u = i * n + j;p = pDataCoef[u] * pDataConst[j];q = pDataCoefImag[u] * pDataConstImag[j];s = (pDataCoef[u] + pDataCoefImag[u]) * (pDataConst[j] + pDataConstImag[j]);pDataConst[i] = pDataConst[i] - p + q;pDataConstImag[i] = pDataConstImag[i] - s + p + q;}}//调整位置pnJs[n - 1] = n - 1;for (k = n - 1; k >= 0; k--){if (pnJs[k] != k){p = pDataConst[k];pDataConst[k] = pDataConst[pnJs[k]];pDataConst[pnJs[k]] = p;p = pDataConstImag[k];pDataConstImag[k] = pDataConstImag[pnJs[k]];pDataConstImag[pnJs[k]] = p;}}return true;}/// <summary>/// 复系数方程组的全选主元高斯－约当消去法/// </summary>/// <param name="mtxCoefImag">系数矩阵的虚部矩阵</param>/// <param name="mtxConstImag">常数矩阵的虚部矩阵</param>/// <param name="mtxResult">Matrix对象，返回方程组解矩阵的实部矩阵</param>/// <param name="mtxResultImag">Matrix对象，返回方程组解矩阵的虚部矩阵</param>/// <returns>bool 型，方程组求解是否成功</returns>public bool GetRootsetGaussJordan(Matrix mtxCoefImag, Matrix mtxConstImag, Matrix mtxResult, Matrix mtxResultImag){int l, k, i, j, nIs = 0, u, v;double p, q, s, d;//方程组的属性，将常数矩阵赋给解矩阵mtxResult.SetValue(mtxLEConst);mtxResultImag.SetValue(mtxConstImag);double[] pDataCoef = mtxLECoef.GetData();double[] pDataConst = mtxResult.GetData();double[] pDataCoefImag = mtxCoefImag.GetData();double[] pDataConstImag = mtxResultImag.GetData();int n = GetNumUnknowns();int m = mtxLEConst.GetNumColumns();//临时缓冲区，存放变换的列数int[] pnJs = new int[n];//消元for (k = 0; k <= n - 1; k++){d = 0.0;for (i = k; i <= n - 1; i++){for (j = k; j <= n - 1; j++){u = i * n + j;p = pDataCoef[u] * pDataCoef[u] + pDataCoefImag[u] * pDataCoefImag[u];if (p > d){d = p;pnJs[k] = j;nIs = i;}}}//求解失败if (d == 0.0)return false;if (nIs != k){for (j = k; j <= n - 1; j++){u = k * n + j;v = nIs * n + j;p = pDataCoef[u];pDataCoef[u] = pDataCoef[v];pDataCoef[v] = p;p = pDataCoefImag[u];pDataCoefImag[u] = pDataCoefImag[v];pDataCoefImag[v] = p;}for (j = 0; j <= m - 1; j++){u = k * m + j;v = nIs * m + j;p = pDataConst[u];pDataConst[u] = pDataConst[v];pDataConst[v] = p;p = pDataConstImag[u];pDataConstImag[u] = pDataConstImag[v];pDataConstImag[v] = p;}}if (pnJs[k] != k){for (i = 0; i <= n - 1; i++){u = i * n + k;v = i * n + pnJs[k];p = pDataCoef[u];pDataCoef[u] = pDataCoef[v];pDataCoef[v] = p;p = pDataCoefImag[u];pDataCoefImag[u] = pDataCoefImag[v];pDataCoefImag[v] = p;}}v = k * n + k;for (j = k + 1; j <= n - 1; j++){u = k * n + j;p = pDataCoef[u] * pDataCoef[v];q = -pDataCoefImag[u] * pDataCoefImag[v];s = (pDataCoef[v] - pDataCoefImag[v]) * (pDataCoef[u] + pDataCoefImag[u]);pDataCoef[u] = (p - q) / d;pDataCoefImag[u] = (s - p - q) / d;}for (j = 0; j <= m - 1; j++){u = k * m + j;p = pDataConst[u] * pDataCoef[v];q = -pDataConstImag[u] * pDataCoefImag[v];s = (pDataCoef[v] - pDataCoefImag[v]) * (pDataConst[u] + pDataConstImag[u]);pDataConst[u] = (p - q) / d;pDataConstImag[u] = (s - p - q) / d;}for (i = 0; i <= n - 1; i++){if (i != k){u = i * n + k;for (j = k + 1; j <= n - 1; j++){v = k * n + j;l = i * n + j;p = pDataCoef[u] * pDataCoef[v];q = pDataCoefImag[u] * pDataCoefImag[v];s = (pDataCoef[u] + pDataCoefImag[u]) * (pDataCoef[v] + pDataCoefImag[v]);pDataCoef[l] = pDataCoef[l] - p + q;pDataCoefImag[l] = pDataCoefImag[l] - s + p + q;}for (j = 0; j <= m - 1; j++){l = i * m + j;v = k * m + j;p = pDataCoef[u] * pDataConst[v]; q = pDataCoefImag[u] * pDataConstImag[v];s = (pDataCoef[u] + pDataCoefImag[u]) * (pDataConst[v] + pDataConstImag[v]);pDataConst[l] = pDataConst[l] - p + q;pDataConstImag[l] = pDataConstImag[l] - s + p + q;}}}}//求解调整for (k = n - 1; k >= 0; k--){if (pnJs[k] != k){for (j = 0; j <= m - 1; j++){u = k * m + j;v = pnJs[k] * m + j;p = pDataConst[u];pDataConst[u] = pDataConst[v];pDataConst[v] = p;p = pDataConstImag[u];pDataConstImag[u] = pDataConstImag[v];pDataConstImag[v] = p;}}}return true;}/// <summary>/// 求解三对角线方程组的追赶法/// </summary>/// <param name="mtxResult">Matrix对象，返回方程组解矩阵</param>/// <returns>bool 型，方程组求解是否成功</returns>public bool GetRootsetTriDiagonal(Matrix mtxResult){int k, j;double s;//将常数矩阵赋给解矩阵mtxResult.SetValue(mtxLEConst);double[] pDataConst = mtxResult.GetData();int n = GetNumUnknowns();if (mtxLECoef.GetNumRows() != n)return false;//为系数矩阵对角线数组分配内存double[] pDiagData = new double[3 * n - 2];//构造系数矩阵对角线元素数组k = j = 0;if (k == 0){pDiagData[j++] = mtxLECoef.GetElement(k, k);pDiagData[j++] = mtxLECoef.GetElement(k, k + 1);}for (k = 1; k < n - 1; ++k){pDiagData[j++] = mtxLECoef.GetElement(k, k - 1);pDiagData[j++] = mtxLECoef.GetElement(k, k);pDiagData[j++] = mtxLECoef.GetElement(k, k + 1);}if (k == n - 1){pDiagData[j++] = mtxLECoef.GetElement(k, k - 1);pDiagData[j++] = mtxLECoef.GetElement(k, k);}//求解for (k = 0; k <= n - 2; k++){j = 3 * k;s = pDiagData[j];// 求解失败if (Math.Abs(s) + 1.0 == 1.0){return false;}pDiagData[j + 1] = pDiagData[j + 1] / s;pDataConst[k] = pDataConst[k] / s;pDiagData[j + 3] = pDiagData[j + 3] - pDiagData[j + 2] * pDiagData[j + 1];pDataConst[k + 1] = pDataConst[k + 1] - pDiagData[j + 2] * pDataConst[k];}s = pDiagData[3 * n - 3];if (s == 0.0){return false;}//调整pDataConst[n - 1] = pDataConst[n - 1] / s;for (k = n - 2; k >= 0; k--)pDataConst[k] = pDataConst[k] - pDiagData[3 * k + 1] * pDataConst[k + 1];return true;}/// <summary>/// 一般带型方程组的求解/// </summary>/// <param name="nBandWidth">带宽</param>/// <param name="mtxResult">Matrix对象，返回方程组解矩阵</param>/// <returns>bool 型，方程组求解是否成功</returns>public bool GetRootsetBand(int nBandWidth, Matrix mtxResult){int ls, k, i, j, nis = 0, u, v;double p, t;//带宽必须为奇数if ((nBandWidth - 1) % 2 != 0)return false;//将常数矩阵赋给解矩阵mtxResult.SetValue(mtxLEConst);double[] pDataConst = mtxResult.GetData();//方程组属性int m = mtxLEConst.GetNumColumns();int n = GetNumUnknowns();if (mtxLECoef.GetNumRows() != n)return false;//带宽数组：带型矩阵的有效数据double[] pBandData = new double[nBandWidth * n];//半带宽ls = (nBandWidth - 1) / 2;//构造带宽数组for (i = 0; i < n; ++i){j = 0;for (k = Math.Max(0, i - ls); k < Math.Max(0, i - ls) + nBandWidth; ++k){if (k < n)pBandData[i * nBandWidth + j++] = mtxLECoef.GetElement(i, k);elsepBandData[i * nBandWidth + j++] = 0;}}//求解for (k = 0; k <= n - 2; k++){p = 0.0;for (i = k; i <= ls; i++){t = Math.Abs(pBandData[i * nBandWidth]);if (t > p){p = t;nis = i;}}if (p == 0.0)return false;for (j = 0; j <= m - 1; j++){u = k * m + j;v = nis * m + j;t = pDataConst[u];pDataConst[u] = pDataConst[v];pDataConst[v] = t;}for (j = 0; j <= nBandWidth - 1; j++){u = k * nBandWidth + j;v = nis * nBandWidth + j;t = pBandData[u];pBandData[u] = pBandData[v];pBandData[v] = t;}for (j = 0; j <= m - 1; j++){u = k * m + j;pDataConst[u] = pDataConst[u] / pBandData[k * nBandWidth];}for (j = 1; j <= nBandWidth - 1; j++){u = k * nBandWidth + j;pBandData[u] = pBandData[u] / pBandData[k * nBandWidth];}for (i = k + 1; i <= ls; i++){t = pBandData[i * nBandWidth];for (j = 0; j <= m - 1; j++){u = i * m + j;v = k * m + j;pDataConst[u] = pDataConst[u] - t * pDataConst[v];}for (j = 1; j <= nBandWidth - 1; j++){u = i * nBandWidth + j;v = k * nBandWidth + j;pBandData[u - 1] = pBandData[u] - t * pBandData[v];}u = i * nBandWidth + nBandWidth - 1; pBandData[u] = 0.0;}if (ls != (n - 1))ls = ls + 1;}p = pBandData[(n - 1) * nBandWidth];if (p == 0.0){return false;}for (j = 0; j <= m - 1; j++){u = (n - 1) * m + j;pDataConst[u] = pDataConst[u] / p;}ls = 1;for (i = n - 2; i >= 0; i--){for (k = 0; k <= m - 1; k++){u = i * m + k;for (j = 1; j <= ls; j++){v = i * nBandWidth + j;nis = (i + j) * m + k;pDataConst[u] = pDataConst[u] - pBandData[v] * pDataConst[nis];}}if (ls != (nBandWidth - 1))ls = ls + 1;}return true;}/// <summary>/// 求解对称方程组的分解法/// </summary>/// <param name="mtxResult">CMatrix引用对象，返回方程组解矩阵</param>/// <returns>bool 型，方程组求解是否成功</returns>public bool GetRootsetDjn(Matrix mtxResult){int i, j, l, k, u, v, w, k1, k2, k3;double p;//方程组属性，将常数矩阵赋给解矩阵Matrix mtxCoef = new Matrix(mtxLECoef);mtxResult.SetValue(mtxLEConst);int n = mtxCoef.GetNumColumns();int m = mtxResult.GetNumColumns();double[] pDataCoef = mtxCoef.GetData();double[] pDataConst = mtxResult.GetData();//非对称系数矩阵，不能用本方法求解if (pDataCoef[0] == 0.0)return false;for (i = 1; i <= n - 1; i++){u = i * n;pDataCoef[u] = pDataCoef[u] / pDataCoef[0];}for (i = 1; i <= n - 2; i++){u = i * n + i;for (j = 1; j <= i; j++){v = i * n + j - 1;l = (j - 1) * n + j - 1;pDataCoef[u] = pDataCoef[u] - pDataCoef[v] * pDataCoef[v] * pDataCoef[l];}p = pDataCoef[u];if (p == 0.0)return false;for (k = i + 1; k <= n - 1; k++){u = k * n + i;for (j = 1; j <= i; j++){v = k * n + j - 1;l = i * n + j - 1;w = (j - 1) * n + j - 1;pDataCoef[u] = pDataCoef[u] - pDataCoef[v] * pDataCoef[l] * pDataCoef[w];}pDataCoef[u] = pDataCoef[u] / p;}}u = n * n - 1;for (j = 1; j <= n - 1; j++){v = (n - 1) * n + j - 1;w = (j - 1) * n + j - 1;pDataCoef[u] = pDataCoef[u] - pDataCoef[v] * pDataCoef[v] * pDataCoef[w];}p = pDataCoef[u];if (p == 0.0)return false;for (j = 0; j <= m - 1; j++){for (i = 1; i <= n - 1; i++){u = i * m + j;for (k = 1; k <= i; k++){v = i * n + k - 1;w = (k - 1) * m + j;pDataConst[u] = pDataConst[u] - pDataCoef[v] * pDataConst[w];}}}for (i = 1; i <= n - 1; i++){u = (i - 1) * n + i - 1;for (j = i; j <= n - 1; j++){v = (i - 1) * n + j;w = j * n + i - 1;pDataCoef[v] = pDataCoef[u] * pDataCoef[w];}}for (j = 0; j <= m - 1; j++){u = (n - 1) * m + j;pDataConst[u] = pDataConst[u] / p;for (k = 1; k <= n - 1; k++){k1 = n - k;k3 = k1 - 1;u = k3 * m + j;for (k2 = k1; k2 <= n - 1; k2++){v = k3 * n + k2;w = k2 * m + j;pDataConst[u] = pDataConst[u] - pDataCoef[v] * pDataConst[w];}pDataConst[u] = pDataConst[u] / pDataCoef[k3 * n + k3];}}return true;}/// <summary>/// 求解对称正定方程组的平方根法/// </summary>/// <param name="mtxResult">CMatrix引用对象，返回方程组解矩阵</param>/// <returns>bool 型，方程组求解是否成功</returns>public bool GetRootsetCholesky(Matrix mtxResult){int i, j, k, u, v;//方程组属性，将常数矩阵赋给解矩阵Matrix mtxCoef = new Matrix(mtxLECoef);mtxResult.SetValue(mtxLEConst);int n = mtxCoef.GetNumColumns();int m = mtxResult.GetNumColumns();double[] pDataCoef = mtxCoef.GetData();double[] pDataConst = mtxResult.GetData();//非对称正定系数矩阵，不能用本方法求解if (pDataCoef[0] <= 0.0)return false;pDataCoef[0] = Math.Sqrt(pDataCoef[0]);for (j = 1; j <= n - 1; j++)pDataCoef[j] = pDataCoef[j] / pDataCoef[0];for (i = 1; i <= n - 1; i++){u = i * n + i;for (j = 1; j <= i; j++){v = (j - 1) * n + i;pDataCoef[u] = pDataCoef[u] - pDataCoef[v] * pDataCoef[v];}if (pDataCoef[u] <= 0.0)return false;pDataCoef[u] = Math.Sqrt(pDataCoef[u]);if (i != (n - 1)){for (j = i + 1; j <= n - 1; j++){v = i * n + j;for (k = 1; k <= i; k++)pDataCoef[v] = pDataCoef[v] - pDataCoef[(k - 1) * n + i] * pDataCoef[(k - 1) * n + j];pDataCoef[v] = pDataCoef[v] / pDataCoef[u];}}}for (j = 0; j <= m - 1; j++){pDataConst[j] = pDataConst[j] / pDataCoef[0];for (i = 1; i <= n - 1; i++){u = i * n + i;v = i * m + j;for (k = 1; k <= i; k++)pDataConst[v] = pDataConst[v] - pDataCoef[(k - 1) * n + i] * pDataConst[(k - 1) * m + j];pDataConst[v] = pDataConst[v] / pDataCoef[u];}}for (j = 0; j <= m - 1; j++){u = (n - 1) * m + j;pDataConst[u] = pDataConst[u] / pDataCoef[n * n - 1];for (k = n - 1; k >= 1; k--){u = (k - 1) * m + j;for (i = k; i <= n - 1; i++){v = (k - 1) * n + i;pDataConst[u] = pDataConst[u] - pDataCoef[v] * pDataConst[i * m + j];}v = (k - 1) * n + k - 1;pDataConst[u] = pDataConst[u] / pDataCoef[v];}}return true;}/// <summary>/// 求解大型稀疏方程组的全选主元高斯－约去消去法/// </summary>/// <param name="mtxResult">CMatrix引用对象，返回方程组解矩阵</param>/// <returns>bool 型，方程组求解是否成功</returns>public bool GetRootsetGgje(Matrix mtxResult){int i, j, k, nIs = 0, u, v;double d, t;//方程组属性，将常数矩阵赋给解矩阵Matrix mtxCoef = new Matrix(mtxLECoef);mtxResult.SetValue(mtxLEConst);int n = mtxCoef.GetNumColumns();double[] pDataCoef = mtxCoef.GetData();double[] pDataConst = mtxResult.GetData();//临时缓冲区，存放变换的列数int[] pnJs = new int[n];//消元for (k = 0; k <= n - 1; k++){d = 0.0;for (i = k; i <= n - 1; i++){for (j = k; j <= n - 1; j++){t = Math.Abs(pDataCoef[i * n + j]);if (t > d){d = t;pnJs[k] = j;nIs = i;}}}if (d == 0.0)return false;if (nIs != k){for (j = k; j <= n - 1; j++){u = k * n + j;v = nIs * n + j;t = pDataCoef[u];pDataCoef[u] = pDataCoef[v];pDataCoef[v] = t;}t = pDataConst[k];pDataConst[k] = pDataConst[nIs];pDataConst[nIs] = t;}if (pnJs[k] != k){for (i = 0; i <= n - 1; i++){u = i * n + k;v = i * n + pnJs[k];t = pDataCoef[u];pDataCoef[u] = pDataCoef[v];pDataCoef[v] = t;}}t = pDataCoef[k * n + k];for (j = k + 1; j <= n - 1; j++){u = k * n + j;if (pDataCoef[u] != 0.0)pDataCoef[u] = pDataCoef[u] / t;}pDataConst[k] = pDataConst[k] / t;for (j = k + 1; j <= n - 1; j++){u = k * n + j;if (pDataCoef[u] != 0.0){for (i = 0; i <= n - 1; i++){v = i * n + k;if ((i != k) && (pDataCoef[v] != 0.0)){nIs = i * n + j;pDataCoef[nIs] = pDataCoef[nIs] - pDataCoef[v] * pDataCoef[u];}}}}for (i = 0; i <= n - 1; i++){u = i * n + k;if ((i != k) && (pDataCoef[u] != 0.0))pDataConst[i] = pDataConst[i] - pDataCoef[u] * pDataConst[k];}}//调整for (k = n - 1; k >= 0; k--){if (k != pnJs[k]){t = pDataConst[k];pDataConst[k] = pDataConst[pnJs[k]];pDataConst[pnJs[k]] = t;}}return true;}/// <summary>/// 求解托伯利兹方程组的列文逊方法/// </summary>/// <param name="mtxResult">CMatrix引用对象，返回方程组解矩阵</param>/// <returns>bool 型，方程组求解是否成功</returns>public bool GetRootsetTlvs(Matrix mtxResult){int i, j, k;double a, beta, q, c, h;//未知数个数int n = mtxLECoef.GetNumColumns();//初始化解解向量mtxResult.Init(n, 1);double[] x = mtxResult.GetData();//常数数组double[] pDataConst = mtxLEConst.GetData();//建立T数组double[] t = new double[n];//构造T数组for (i = 0; i < n; ++i)t[i] = mtxLECoef.GetElement(0, i);//临时数组double[] s = new double[n];double[] y = new double[n];//非托伯利兹方程组，不能用本方法求解a = t[0];if (a == 0.0)return false;//列文逊方法求解y[0] = 1.0;x[0] = pDataConst[0] / a;for (k = 1; k <= n - 1; k++){beta = 0.0;q = 0.0;for (j = 0; j <= k - 1; j++){beta = beta + y[j] * t[j + 1];q = q + x[j] * t[k - j];}if (a == 0.0)return false;c = -beta / a;s[0] = c * y[k - 1];y[k] = y[k - 1];if (k != 1){for (i = 1; i <= k - 1; i++)s[i] = y[i - 1] + c * y[k - i - 1];}a = a + c * beta;if (a == 0.0)return false;h = (pDataConst[k] - q) / a;for (i = 0; i <= k - 1; i++){x[i] = x[i] + h * s[i];y[i] = s[i];}x[k] = h * y[k];}return true;}/// <summary>/// 高斯－赛德尔迭代法/// </summary>/// <param name="mtxResult">CMatrix引用对象，返回方程组解矩阵</param>/// <param name="eps">控制精度</param>/// <returns>bool 型，方程组求解是否成功</returns>public bool GetRootsetGaussSeidel(Matrix mtxResult, double eps){int i, j, u, v;double p, t, s, q;//未知数个数int n = mtxLECoef.GetNumColumns();//初始化解向量mtxResult.Init(n, 1);double[] x = mtxResult.GetData();//系数与常数double[] pDataCoef = mtxLECoef.GetData();double[] pDataConst = mtxLEConst.GetData();//求解for (i = 0; i <= n - 1; i++){u = i * n + i;p = 0.0;x[i] = 0.0;for (j = 0; j <= n - 1; j++){if (i != j){v = i * n + j;p = p + Math.Abs(pDataCoef[v]);}}if (p >= Math.Abs(pDataCoef[u]))return false;}//精度控制p = eps + 1.0;while (p >= eps){p = 0.0;for (i = 0; i <= n - 1; i++){t = x[i];s = 0.0;for (j = 0; j <= n - 1; j++)if (j != i)s = s + pDataCoef[i * n + j] * x[j];x[i] = (pDataConst[i] - s) / pDataCoef[i * n + i];q = Math.Abs(x[i] - t) / (1.0 + Math.Abs(x[i]));if (q > p)p = q;}}return true;}/// <summary>/// 求解对称正定方程组的共轭梯度法/// </summary>/// <param name="mtxResult">CMatrix引用对象，返回方程组解矩阵</param>/// <param name="eps">控制精度</param>public void GetRootsetGrad(Matrix mtxResult, double eps){int i, k;double alpha, beta, d, e;//未知数个数int n = GetNumUnknowns();//初始化解向量mtxResult.Init(n, 1);double[] x = mtxResult.GetData();//构造临时矩阵Matrix mtxP = new Matrix(n, 1);double[] p = mtxP.GetData();double[] pDataCoef = mtxLECoef.GetData();double[] pDataConst = mtxLEConst.GetData();double[] r = new double[n];for (i = 0; i <= n - 1; i++){x[i] = 0.0;p[i] = pDataConst[i];r[i] = pDataConst[i];}i=0;while (i <= n - 1){Matrix mtxS = mtxLECoef.Multiply(mtxP);double[] s = mtxS.GetData();d = 0.0;e = 0.0;for (k = 0; k <= n - 1; k++){d = d + p[k] * pDataConst[k];e = e + p[k] * s[k];}alpha = d / e;for (k = 0; k <= n - 1; k++)x[k] = x[k] + alpha * p[k];Matrix mtxQ = mtxLECoef.Multiply(mtxResult);double[] q = mtxQ.GetData();d = 0.0;for (k = 0; k <= n - 1; k++){r[k] = pDataConst[k] - q[k];d = d + r[k] * s[k];}beta = d / e; d = 0.0;for (k = 0; k <= n - 1; k++)d = d + r[k] * r[k];//满足精度，求解结束d = Math.Sqrt(d);if (d < eps)break;for (k = 0; k <= n - 1; k++)p[k] = r[k] - beta * p[k];i = i + 1;}}/// <summary>/// 求解线性最小二乘问题的豪斯荷尔德变换法/// </summary>/// <param name="mtxResult">Matrix对象，返回方程组解矩阵</param>/// <param name="mtxQ">Matrix对象，返回豪斯荷尔德变换的Q矩阵</param>/// <param name="mtxR">Matrix对象，返回豪斯荷尔德变换的R矩阵</param>/// <returns>bool型，方程组求解是否成功</returns>public bool GetRootsetMqr(Matrix mtxResult, Matrix mtxQ, Matrix mtxR){int i, j;double d;//方程组的方程数和未知数个数int m = mtxLECoef.GetNumRows();int n = mtxLECoef.GetNumColumns();//奇异方程组if (m < n)return false;//将解向量初始化为常数向量mtxResult.SetValue(mtxLEConst);double[] pDataConst = mtxResult.GetData();//构造临时矩阵，用于QR分解mtxR.SetValue(mtxLECoef);double[] pDataCoef = mtxR.GetData();//QR分解if (!mtxR.SplitQR(mtxQ))return false;//临时缓冲区double[] c = new double[n];double[] q = mtxQ.GetData();//求解for (i=0; i<=n-1; i++){ d=0.0;for (j=0; j<=m-1; j++)d=d+q[j*m+i]*pDataConst[j];c[i]=d;}pDataConst[n - 1] = c[n - 1] / pDataCoef[n * n - 1];for (i = n - 2; i >= 0; i--){d = 0.0;for (j = i + 1; j <= n - 1; j++)d = d + pDataCoef[i * n + j] * pDataConst[j];pDataConst[i] = (c[i] - d) / pDataCoef[i * n + i];}return true;}/// <summary>/// 求解线性最小二乘问题的广义逆法/// </summary>/// <param name="mtxResult">Matrix对象，返回方程组解矩阵</param>/// <param name="mtxAP">Matrix对象，返回系数矩阵的广义逆矩阵</param>/// <param name="mtxU">Matrix对象，返回U矩阵</param>/// <param name="mtxV">Matrix对象，返回V矩阵</param>/// <param name="eps">控制精度</param>/// <returns>bool型，方程组求解是否成功</returns>public bool GetRootsetGinv(Matrix mtxResult, Matrix mtxAP, Matrix mtxU, Matrix mtxV, double eps){int i, j;//方程个数和未知数个数int m = mtxLECoef.GetNumRows();int n = mtxLECoef.GetNumColumns();//初始化解向量mtxResult.Init(n, 1);double[] pDataConst = mtxLEConst.GetData();double[] x = mtxResult.GetData();//临时矩阵Matrix mtxA = new Matrix(mtxLECoef);//求广义逆矩阵if (!mtxA.InvertUV(mtxAP, mtxU, mtxV, eps))return false;double[] pAPData = mtxAP.GetData();//求解for (i = 0; i <= n - 1; i++){x[i] = 0.0;for (j = 0; j <= m - 1; j++)x[i] = x[i] + pAPData[i * m + j] * pDataConst[j];}return true;}/// <summary>/// Matrix对象，返回方程组解矩阵/// </summary>/// <param name="mtxResult">Matrix对象，返回方程组解矩阵</param>/// <param name="nMaxIt">叠加次数</param>/// <param name="eps">控制精度</param>/// <returns>bool型，方程组求解是否成功</returns>public bool GetRootsetMorbid(Matrix mtxResult, int nMaxIt, double eps){int i, k;double q, qq;//方程的阶数int n = GetNumUnknowns();//设定迭代次数, 缺省为60i = nMaxIt;//用全选主元高斯消元法求解LEquations leqs = new LEquations(mtxLECoef, mtxLEConst);if (!leqs.GetRootsetGauss(mtxResult))return false;double[] x = mtxResult.GetData();q = 1.0 + eps;while (q >= eps){// 迭代次数已达最大值，仍为求得结果，求解失败if (i == 0)return false;// 迭代次数减1i = i - 1;// 矩阵运算Matrix mtxE = mtxLECoef.Multiply(mtxResult);Matrix mtxR = mtxLEConst.Substract(mtxE);// 用全选主元高斯消元法求解leqs = new LEquations(mtxLECoef, mtxR);Matrix mtxRR = new Matrix();if (!leqs.GetRootsetGauss(mtxRR))return false;double[] r = mtxRR.GetData();q = 0.0;for (k = 0; k <= n - 1; k++){qq = Math.Abs(r[k]) / (1.0 + Math.Abs(x[k] + r[k]));if (qq > q)q = qq;}for (k = 0; k <= n - 1; k++)x[k] = x[k] + r[k];}return true;}}
}

注：本代码为原作者所有，本人只是摘录，如有侵犯作者，请与本人联系

Email:359194966@qq.com

微软公司内部培训程序员资料---求解线性方程组的类相关推荐

微软公司内部培训程序员资料---操作矩阵类
/** 操作矩阵的类 Matrix* * 周长发编制*/ using System;namespace MSAlgorithm {public class Matrix{/// <summary ...
公司停电，程序员去网吧写代码；iPhone 14将于北京时间9月8日发布；GitLab修复一个关键远程代码执行漏洞|极客头条
「极客头条」-- 技术人员的新闻圈! CSDN 的读者朋友们早上好哇,「极客头条」来啦,快来看今天都有哪些值得我们技术人关注的重要新闻吧. 整理 | 梦依丹出品 | CSDN(ID:CSDNnews ...
公司停电，程序员去网吧写代码；iPhone 14将于北京时间9月8日发布；GitLab修复一个关键远程代码执行漏洞|极客头条...
「极客头条」-- 技术人员的新闻圈! CSDN 的读者朋友们早上好哇,「极客头条」来啦,快来看今天都有哪些值得我们技术人关注的重要新闻吧. 整理 | 梦依丹出品 | CSDN(ID:CSDNnews ...
程序员职场须知：公司如何衡量程序员的价值？别以为是经常加班！
1.如何评价程序员的薪酬? 我认为程序员是这个时代很好的职业,从最现实的薪酬角度来看,跟很多行业相比,互联网行业依然是高薪行业.比如前不久猎聘发布的2019年一季度中高端人才报告中的数据就显示,互联网 ...
一个我自己建的程序员资料分享站
自己建的一个小站,目的是分享收集起来的程序员资料(大部分是pdf和视频). 资料都经过我分类和统计,每个资料都会有简介或者目录,可以看看是否适合再下载. 具体可以到我的小站程序员资料站看看. 资源 ...
公司招个程序员，34岁以上两年一跳的不要，开出工资以为看错了
往期好文推荐 0基础不用怕,从0到1轻松教你入门Python python系统学习流线图,教你一步一步学会python 导读:对于公司来说,肯定是希望花最少的钱招到最优秀的员工,但事实上这个想法是不太 ...
谈谈我2013上半年在公司内部培训的经历
在这里谈谈我在今年上半年在公司内部培训实习生(今年七月毕业的学生)的经历. 今年三月份我担任公司技术培训的负责人,回顾培训历程,做个总结. 公司目标是通过培训能让实习生快速掌握知识,成为一名优秀的技术 ...
java中类图概念,程序员眼中的UML(4)--类图释疑之一，Attribute和Property之区别
程序员眼中的UML(4) --类图释疑之一,Attribute和Property之区别上一篇中提出了很多问题,其中最令人费解的可能就是Attribute和Property之区别了吧.我在网络上寻找良 ...
身为程序员还看不懂UML类图？一文带你零基础学会看UML类图！
身为程序员还看不懂UML类图? 一文带你零基础学会看UML类图! 一,UML类图示例图二,UML类图图例三,分步解析说明 3.1 类图: 3.2 接口: 3.3 实现继承 3.4 实现接口: 3. ...

微软公司内部培训程序员资料---求解线性方程组的类

微软公司内部培训程序员资料---求解线性方程组的类相关推荐

最新文章

热门文章