微软公司内部培训程序员资料---操作矩阵类

/** 操作矩阵的类 Matrix* * 周长发编制*/
using System;namespace MSAlgorithm
{public class Matrix{/// <summary>/// 矩阵列数/// </summary>private int numColumns = 0;/// <summary>/// 矩阵行数/// </summary>private int numRows = 0;/// <summary>/// 缺省精度/// </summary>private double eps = 0.0;/// <summary>/// 矩阵数据缓冲区/// </summary>private double[] elements = null;/// <summary>/// 矩阵列数/// </summary>public int Columns{get { return numColumns; }}/// <summary>/// 矩阵行数/// </summary>public int Rows{get { return numRows; }}/// <summary>/// 索引器: 访问矩阵元素/// </summary>/// <param name="row">元素的行</param>/// <param name="col">元素的列</param>/// <returns></returns>public double this[int row, int col]{get { return elements[col + row * numColumns]; }set { elements[col + row * numColumns] = value; }}/// <summary>/// 缺省精度/// </summary>public double Eps{get { return eps; }set { eps = value; }}/// <summary>/// 基本构造函数/// </summary>public Matrix(){numColumns = 1;numRows = 1;Init(numRows, numColumns);}/// <summary>/// 指定行列构造函数/// </summary>/// <param name="nRows">指定的矩阵行数</param>/// <param name="nCols">指定的矩阵列数</param>public Matrix(int nRows, int nCols){numRows = nRows;numColumns = nCols;Init(numRows, numColumns);}/// <summary>/// 指定值构造函数/// </summary>/// <param name="value">二维数组，存储矩阵各元素的值</param>public Matrix(double[,] value){numRows = value.GetLength(0);numColumns = value.GetLength(1);double[] data = new double[numRows * numColumns];int k = 0;for (int i = 0; i < numRows; ++i){for (int j = 0; j < numColumns; ++j){data[k++] = value[i, j];}}Init(numRows, numColumns);SetData(data);}/// <summary>/// 指定值构造函数/// </summary>/// <param name="nRows">指定的矩阵行数</param>/// <param name="nCols">指定的矩阵列数</param>/// <param name="value">一维数组，长度为nRows*nCols，存储矩阵各元素的值</param>public Matrix(int nRows, int nCols, double[] value){numRows = nRows;numColumns = nCols;Init(numRows, numColumns);SetData(value);}/// <summary>/// 方阵构造函数/// </summary>/// <param name="nSize">方阵行列数</param>public Matrix(int nSize){numRows = nSize;numColumns = nSize;Init(nSize, nSize);}/// <summary>/// 方阵构造函数/// </summary>/// <param name="nSize">方阵行列数</param>/// <param name="value">一维数组，长度为nRows*nRows，存储方阵各元素的值</param>public Matrix(int nSize, double[] value){numRows = nSize;numColumns = nSize;Init(nSize, nSize);SetData(value);}/// <summary>/// 拷贝构造函数/// </summary>/// <param name="other">源矩阵</param>public Matrix(Matrix other){numColumns = other.GetNumColumns();numRows = other.GetNumRows();Init(numRows, numColumns);SetData(other.elements);}/// <summary>/// 初始化函数/// </summary>/// <param name="nRows">指定的矩阵行数</param>/// <param name="nCols">指定的矩阵列数</param>/// <returns>bool, 成功返回true, 否则返回false</returns>public bool Init(int nRows, int nCols){numRows = nRows;numColumns = nCols;int nSize = nCols * nRows;if (nSize < 0)return false;//分配内存elements = new double[nSize];return true;}/// <summary>/// 设置矩阵运算的精度/// </summary>/// <param name="newEps">新的精度值</param>public void SetEps(double newEps){eps = newEps;}/// <summary>/// 取矩阵的精度值/// </summary>/// <returns>double型，矩阵的精度值</returns>public double GetEps(){return eps;}/// <summary>/// 重载 + 运算符/// </summary>/// <param name="m1">指定的矩阵1</param>/// <param name="m2">指定的矩阵2</param>/// <returns>Matrix对象</returns>public static Matrix operator +(Matrix m1, Matrix m2){return m1.Add(m2);}/// <summary>/// 重载 - 运算符/// </summary>/// <param name="m1">指定的矩阵1</param>/// <param name="m2">指定的矩阵2</param>/// <returns>Matrix对象</returns>public static Matrix operator -(Matrix m1, Matrix m2){return m1.Substract(m2);}/// <summary>/// 重载 * 运算符/// </summary>/// <param name="m1">指定的矩阵1</param>/// <param name="m2">指定的矩阵2</param>/// <returns>Matrix对象</returns>public static Matrix operator *(Matrix m1, Matrix m2){return m1.Multiply(m2);}/// <summary>/// 重载 double[] 运算符/// </summary>/// <param name="m">Matrix对象</param>/// <returns>double[]对象</returns>public static implicit operator double[](Matrix m){return m.elements;}/// <summary>/// 将方阵初始化为单位矩阵/// </summary>/// <param name="nSize">nSize - 方阵行列数</param>/// <returns>bool 型，初始化是否成功</returns>public bool MakeUnitMatrix(int nSize){if (!Init(nSize, nSize))return false;for (int i = 0; i < nSize; ++i)for (int j = 0; j < nSize; ++j)if (i == j)SetElement(i, j, 1);return true;}/// <summary>/// 将矩阵各元素的值转化为字符串, 元素之间的分隔符为",", 行与行之间有回车换行符/// </summary>/// <returns>string 型，转换得到的字符串</returns>public override string ToString(){return ToString(",", true);}/// <summary>/// 将矩阵各元素的值转化为字符串/// </summary>/// <param name="sDelim">元素之间的分隔符</param>/// <param name="bLineBreak">行与行之间是否有回车换行符</param>/// <returns>string 型，转换得到的字符串</returns>public string ToString(string sDelim, bool bLineBreak){string s = "";for (int i = 0; i < numRows; ++i){for (int j = 0; j < numColumns; ++j){string ss = GetElement(i, j).ToString("F");s += ss;if (bLineBreak){if (j != numColumns - 1)s += sDelim;}else{if (i != numRows - 1 || j != numColumns - 1)s += sDelim;}}if (bLineBreak)if (i != numRows - 1)s += "\r\n";}return s;}/// <summary>/// 将矩阵指定行中各元素的值转化为字符串/// </summary>/// <param name="nRow">指定的矩阵行，nRow = 0表示第一行</param>/// <param name="sDelim">元素之间的分隔符</param>/// <returns>string 型，转换得到的字符串</returns>public string ToStringRow(int nRow, string sDelim){string s = "";if (nRow >= numRows)return s;for (int j = 0; j < numColumns; ++j){string ss = GetElement(nRow, j).ToString("F");s += ss;if (j != numColumns - 1)s += sDelim;}return s;}/// <summary>/// 将矩阵指定列中各元素的值转化为字符串/// </summary>/// <param name="nCol">指定的矩阵行，nCol = 0表示第一列</param>/// <param name="sDelim">元素之间的分隔符</param>/// <returns>string 型，转换得到的字符串</returns>public string ToStringCol(int nCol, string sDelim){string s = "";if (nCol >= numColumns)return s;for (int i = 0; i < numRows; ++i){string ss = GetElement(i, nCol).ToString("F");s += ss;if (i != numRows - 1)s += sDelim;}return s;}/// <summary>/// 设置矩阵各元素的值/// </summary>/// <param name="value">一维数组，长度为numColumns*numRows，存储矩阵各元素的值</param>public void SetData(double[] value){elements = (double[])value.Clone();}/// <summary>/// 设置指定元素的值/// </summary>/// <param name="nRow">元素的行</param>/// <param name="nCol">元素的列</param>/// <param name="value">指定元素的值</param>/// <returns>bool 型，说明设置是否成功</returns>public bool SetElement(int nRow, int nCol, double value){if ((nCol < 0) || (nCol >= numColumns) || (nRow < 0) || (nRow >= numColumns))return false;elements[nCol + nRow * numColumns] = value;return true;}/// <summary>/// 获取指定元素的值/// </summary>/// <param name="nRow">元素的行</param>/// <param name="nCol">元素的列</param>/// <returns>double 型，指定元素的值</returns>public double GetElement(int nRow, int nCol){return elements[nCol + nRow * numColumns];}/// <summary>/// 获取矩阵的列数/// </summary>/// <returns>int 型，矩阵的列数</returns>public int GetNumColumns(){return numColumns;}/// <summary>/// 获取矩阵的行数/// </summary>/// <returns>int 型，矩阵的行数</returns>public int GetNumRows(){return numRows;}/// <summary>/// 获取矩阵的数据/// </summary>/// <returns>double型数组，指向矩阵各元素的数据缓冲区</returns>public double[] GetData(){return elements;}/// <summary>/// 获取指定行的向量/// </summary>/// <param name="nRow">向量所在的行</param>/// <param name="pVector">指向向量中各元素的缓冲区</param>/// <returns>int 型，向量中元素的个数，即矩阵的列数</returns>public int GetRowVector(int nRow, double[] pVector){for (int j = 0; j < numColumns; ++j)pVector[j] = GetElement(nRow, j);return numColumns;}/// <summary>/// 获取指定列的向量/// </summary>/// <param name="nCol">向量所在的列</param>/// <param name="pVector">指向向量中各元素的缓冲区</param>/// <returns>int 型，向量中元素的个数，即矩阵的行数</returns>public int GetColVector(int nCol, double[] pVector){for (int i = 0; i < numRows; ++i)pVector[i] = GetElement(i, nCol);return numRows;}/// <summary>/// 给矩阵赋值/// </summary>/// <param name="other">用于给矩阵赋值的源矩阵</param>/// <returns>Matrix型，阵与other相等</returns>public Matrix SetValue(Matrix other){if (other != this){Init(other.GetNumRows(), other.GetNumColumns());SetData(other.elements);}return this;}/// <summary>/// 判断矩阵否相等/// </summary>/// <param name="obj">用于比较的矩阵</param>/// <returns>bool 型，两个矩阵相等则为true，否则为false</returns>public override bool Equals(object other){Matrix m = other as Matrix;if (m == null)return false;//首先检查行列数是否相等if (numColumns != m.GetNumColumns() || numRows != m.GetNumRows())return false;for (int i = 0; i < numRows; ++i){for (int j = 0; j < numColumns; ++j){if (Math.Abs(GetElement(i, j) - m.GetElement(i, j)) > eps)return false;}}return true;}/// <summary>/// 因为重写了Equals，因此必须重写GetHashCode/// </summary>/// <returns>int型，返回复数对象散列码</returns>public override int GetHashCode(){double sum = 0;for (int i = 0; i < numRows; ++i){for (int j = 0; j < numColumns; ++j){sum += Math.Abs(GetElement(i, j));}}return (int)Math.Sqrt(sum);}/// <summary>/// 实现矩阵的加法/// </summary>/// <param name="other">与指定矩阵相加的矩阵</param>/// <returns>Matrix型，指定矩阵与other相加之和,如果矩阵的行/列数不匹配，则会抛出异常</returns>      public Matrix Add(Matrix other){//首先检查行列数是否相等if ((numColumns != other.GetNumColumns()) ||(numRows != other.GetNumRows()))throw new Exception("矩阵的行/列数不匹配。");//构造结果矩阵,拷贝构造Matrix result = new Matrix(this);//矩阵加法for (int i = 0; i < numRows; ++i){for (int j = 0; j < numColumns; ++j)result.SetElement(i, j, result.GetElement(i, j) + other.GetElement(i, j));}return result;}/// <summary>/// 实现矩阵的减法/// </summary>/// <param name="other">与指定矩阵相减的矩阵</param>/// <returns>Matrix型，指定矩阵与other相减之差,如果矩阵的行/列数不匹配，则会抛出异常</returns>public Matrix Substract(Matrix other){if ((numColumns != other.GetNumColumns()) ||(numRows != other.GetNumRows()))throw new Exception("矩阵的行/列数不匹配。");//构造结果矩阵,拷贝构造Matrix result = new Matrix(this);//进行减法操作for (int i = 0; i < numRows; ++i){for (int j = 0; j < numColumns; ++j)result.SetElement(i, j, result.GetElement(i, j) - other.GetElement(i, j));}return result;}/// <summary>/// 实现矩阵的数乘/// </summary>/// <param name="value">与指定矩阵相乘的实数</param>/// <returns>Matrix型，指定矩阵与value相乘之积</returns>public Matrix Multiply(double value){//构造目标矩阵,拷贝构造Matrix result = new Matrix(this);//进行数乘for (int i = 0; i < numRows; ++i){for (int j = 0; j < numColumns; ++j)result.SetElement(i, j, result.GetElement(i, j) * value);}return result;}/// <summary>/// 实现矩阵的乘法/// </summary>/// <param name="other">与指定矩阵相乘的矩阵</param>/// <returns>Matrix型，指定矩阵与other相乘之积,如果矩阵的行/列数不匹配，则会抛出异常</returns>public Matrix Multiply(Matrix other){//首先检查行列数是否符合要求if (numColumns != other.GetNumColumns())throw new Exception("矩阵的行/列数不匹配。");//构造目标矩阵Matrix result = new Matrix(numRows, other.GetNumColumns());//矩阵乘法，即//// [A][B][C]   [G][H]     [A*G + B*I + C*K][A*H + B*J + C*L]// [D][E][F] * [I][J] =   [D*G + E*I + F*K][D*H + E*J + F*L]//             [K][L]//double value;for (int i = 0; i < result.GetNumRows(); ++i){for (int j = 0; j < other.GetNumColumns(); ++j){value = 0.0;for (int k = 0; k < numColumns; ++k){value += GetElement(i, k) * other.GetElement(k, j);}result.SetElement(i, j, value);}}return result;}/// <summary>/// 复矩阵的乘法/// </summary>/// <param name="AR">左边复矩阵的实部矩阵</param>/// <param name="AI">左边复矩阵的虚部矩阵</param>/// <param name="BR">右边复矩阵的实部矩阵</param>/// <param name="BI">右边复矩阵的虚部矩阵</param>/// <param name="CR">乘积复矩阵的实部矩阵</param>/// <param name="CI">乘积复矩阵的虚部矩阵</param>/// <returns>bool型，复矩阵乘法是否成功</returns>public bool Multiply(Matrix AR, Matrix AI, Matrix BR, Matrix BI, Matrix CR, Matrix CI){//首先检查行列数是否符合要求if (AR.GetNumColumns() != AI.GetNumColumns() ||AR.GetNumRows() != AI.GetNumRows() ||BR.GetNumColumns() != BI.GetNumColumns() ||BR.GetNumRows() != BI.GetNumRows() ||AR.GetNumColumns() != BR.GetNumRows())return false;//构造乘积矩阵实部矩阵和虚部矩阵Matrix mtxCR = new Matrix(AR.GetNumRows(), BR.GetNumColumns());Matrix mtxCI = new Matrix(AR.GetNumRows(), BR.GetNumColumns());//复矩阵相乘for (int i = 0; i < AR.GetNumRows(); ++i){for (int j = 0; j < BR.GetNumColumns(); ++j){double vr = 0;double vi = 0;for (int k = 0; k < AR.GetNumColumns(); ++k){double p = AR.GetElement(i, k) * BR.GetElement(k, j);double q = AI.GetElement(i, k) * BI.GetElement(k, j);double s = (AR.GetElement(i, k) + AI.GetElement(i, k)) * (BR.GetElement(k, j) + BI.GetElement(k, j));vr += p - q;vi += s - p - q;}mtxCR.SetElement(i, j, vr);mtxCI.SetElement(i, j, vi);}}CR = mtxCR;CI = mtxCI;return true;}/// <summary>/// 矩阵的转置/// </summary>/// <returns>Matrix型，指定矩阵转置矩阵</returns>public Matrix Transpose(){//构造目标矩阵Matrix trans = new Matrix(numColumns, numRows);//转置各元素for (int i = 0; i < numRows; ++i){for (int j = 0; j < numColumns; ++j)trans.SetElement(j, i, GetElement(i, j));}return trans;}/// <summary>/// 实矩阵求逆的全选主元高斯 - 约当法/// </summary>/// <returns>bool型，求逆是否成功</returns>public bool InvertGaussJordan(){int i, j, k, l, u, v;double d = 0, p = 0;//分配内存int[] pnRow = new int[numColumns];int[] pnCol = new int[numColumns];//消元for (k = 0; k <= numColumns - 1; k++){d = 0.0;for (i = k; i <= numColumns - 1; i++){for (j = k; j <= numColumns - 1; j++){l = i * numColumns + j;p = Math.Abs(elements[l]);if (p > d){d = p;pnRow[k] = i;pnCol[k] = j;}}}//失败if(d==0.0)return false;if (pnRow[k] != k){for (j = 0; j <= numColumns - 1; j++){u = k * numColumns + j;v = pnRow[k] * numColumns + j;p = elements[u];elements[u] = elements[v];elements[v] = p;}}if (pnCol[k] != k){for (i = 0; i <= numColumns - 1; i++){u = i * numColumns + k;v = i * numColumns + pnCol[k];p = elements[u];elements[u] = elements[v];elements[v] = p;}}l = k * numColumns + k;elements[l] = 1.0 / elements[l];for (j = 0; j <= numColumns - 1; j++){if (j != k){u = k * numColumns + j;elements[u] = elements[u] * elements[l];}}for (i = 0; i <= numColumns - 1; i++){if (i != k){for (j = 0; j <= numColumns - 1; j++){if (j != k){u = i * numColumns + j;elements[u] = elements[u] - elements[i * numColumns + k] * elements[k * numColumns + j];}}}}for (i = 0; i <= numColumns - 1; i++){if (i != k){u = i * numColumns + k;elements[u] = -elements[u] * elements[l];}}}//调整恢复行列次序for (k = numColumns - 1; k >= 0; k--){if (pnCol[k] != k){for (j = 0; j <= numColumns - 1; j++){u = k * numColumns + j;v = pnCol[k] * numColumns + j;p = elements[u];elements[u] = elements[v];elements[v] = p;}}if (pnRow[k] != k){for (i = 0; i <= numColumns - 1; i++){u = i * numColumns + k;v = i * numColumns + pnRow[k];p = elements[u];elements[u] = elements[v];elements[v] = p;}}}return true;}/// <summary>/// 复矩阵求逆的全选主元高斯 - 约当法/// </summary>/// <param name="mtxImag">复矩阵的虚部矩阵，当前矩阵为复矩阵的实部</param>/// <returns>bool型，求逆是否成功</returns>public bool InvertGaussJordan(Matrix mtxImag){int i, j, k, l, u, v, w;double p, q, s, t, d, b;//分配内存int[] pnRow = new int[numColumns];int[] pnCol = new int[numColumns];//消元for (k = 0; k <= numColumns - 1; k++){d = 0.0;for (i = k; i <= numColumns - 1; i++){for (j = k; j <= numColumns - 1; j++){u = i * numColumns + j;p = elements[u] * elements[u] + mtxImag.elements[u] * mtxImag.elements[u];if (p > d){d = p;pnRow[k] = i;pnCol[k] = j;}}}//失败if (d == 0.0)return false;if (pnRow[k] != k){for (j = 0; j <= numColumns - 1; j++){u = k * numColumns + j;v = pnRow[k] * numColumns + j;t = elements[u];elements[u] = elements[v];elements[v] = t;t = mtxImag.elements[u];mtxImag.elements[u] = mtxImag.elements[v];mtxImag.elements[v] = t;}}if (pnCol[k] != k){for (i = 0; i <= numColumns - 1; i++){u = i * numColumns + k;v = i * numColumns + pnCol[k];t = elements[u];elements[u] = elements[v];elements[v] = t;t = mtxImag.elements[u];mtxImag.elements[u] = mtxImag.elements[v];mtxImag.elements[v] = t;}}l = k * numColumns + k;elements[l] = elements[l] / d; mtxImag.elements[l] = -mtxImag.elements[l] / d;for (j = 0; j <= numColumns - 1; j++){if (j != k){u = k * numColumns + j;p = elements[u] * elements[l];q = mtxImag.elements[u] * mtxImag.elements[l];s = (elements[u] + mtxImag.elements[u]) * (elements[l] + mtxImag.elements[l]);elements[u] = p - q;mtxImag.elements[u] = s - p - q;}}for (i = 0; i <= numColumns - 1; i++){if (i != k){v = i * numColumns + k;for (j = 0; j <= numColumns - 1; j++){if (j != k){u = k * numColumns + j;w = i * numColumns + j;p = elements[u] * elements[v];q = mtxImag.elements[u] * mtxImag.elements[v];s = (elements[u] + mtxImag.elements[u]) * (elements[v] + mtxImag.elements[v]);t = p - q;b = s - p - q;elements[w] = elements[w] - t;mtxImag.elements[w] = mtxImag.elements[w] - b;}}}}for (i = 0; i <= numColumns - 1; i++){if (i != k){u = i * numColumns + k;p = elements[u] * elements[l];q = mtxImag.elements[u] * mtxImag.elements[l];s = (elements[u] + mtxImag.elements[u]) * (elements[l] + mtxImag.elements[l]);elements[u] = q - p;mtxImag.elements[u] = p + q - s;}}}//调整恢复行列次序for (k = numColumns - 1; k >= 0; k--){if (pnCol[k] != k){for (j = 0; j <= numColumns - 1; j++){u = k * numColumns + j;v = pnCol[k] * numColumns + j;t = elements[u];elements[u] = elements[v];elements[v] = t;t = mtxImag.elements[u];mtxImag.elements[u] = mtxImag.elements[v];mtxImag.elements[v] = t;}}if (pnRow[k] != k){for (i = 0; i <= numColumns - 1; i++){u = i * numColumns + k;v = i * numColumns + pnRow[k];t = elements[u];elements[u] = elements[v];elements[v] = t;t = mtxImag.elements[u];mtxImag.elements[u] = mtxImag.elements[v];mtxImag.elements[v] = t;}}}return true;}/// <summary>/// 对称正定矩阵的求逆/// </summary>/// <returns>bool型，求逆是否成功</returns>public bool InvertSsgj(){int i, j, k, m;double w, g;//临时内存double[] pTmp = new double[numColumns];//逐列处理for (k = 0; k <= numColumns - 1; k++){w = elements[0];if (w == 0.0){return false;}m = numColumns - k - 1;for (i = 1; i <= numColumns - 1; i++){g = elements[i * numColumns];pTmp[i] = g / w;if (i <= m)pTmp[i] = -pTmp[i];for (j = 1; j <= i; j++)elements[(i - 1) * numColumns + j - 1] = elements[i * numColumns + j] + g * pTmp[j];}elements[numColumns * numColumns - 1] = 1.0 / w;for (i = 1; i <= numColumns - 1; i++)elements[(numColumns - 1) * numColumns + i - 1] = pTmp[i];}//行列调整for (i = 0; i <= numColumns - 2; i++)for (j = i + 1; j <= numColumns - 1; j++)elements[i * numColumns + j] = elements[j * numColumns + i];return true;}/// <summary>/// 托伯利兹矩阵求逆的埃兰特方法/// </summary>/// <returns>bool型，求逆是否成功</returns>public bool InvertTrench(){int i, j, k;double a, s;//上三角元素double[] t = new double[numColumns];//下三角元素double[] tt = new double[numColumns];//上、下三角元素赋值for (i = 0; i < numColumns; ++i){t[i] = GetElement(0, i);tt[i] = GetElement(i, 0);}//临时缓冲区double[] c = new double[numColumns];double[] r = new double[numColumns];double[] p = new double[numColumns];//非Toeplitz矩阵，返回if (t[0] == 0.0)return false;a = t[0];c[0] = tt[1] / t[0];r[0] = t[1] / t[0];for (k = 0; k <= numColumns - 3; k++){s = 0.0;for (j = 1; j <= k + 1; j++)s = s + c[k + 1 - j] * tt[j];s = (s - tt[k + 2]) / a;for (i = 0; i <= k; i++)p[i] = c[i] + s * r[k - i];c[k + 1] = -s;s = 0.0;for (j = 1; j <= k + 1; j++)s = s + r[k + 1 - j] * t[j];s = (s - t[k + 2]) / a;for (i = 0; i <= k; i++){r[i] = r[i] + s * c[k - i];c[k - i] = p[k - i];}r[k + 1] = -s;a = 0.0;for (j = 1; j <= k + 2; j++)a = a + t[j] * c[j - 1];a = t[0] - a;//求解失败if (a == 0.0)               return false;                }elements[0] = 1.0 / a;for (i = 0; i <= numColumns - 2; i++){k = i + 1;j = (i + 1) * numColumns;elements[k] = -r[i] / a;elements[j] = -c[i] / a;}for (i = 0; i <= numColumns - 2; i++){for (j = 0; j <= numColumns - 2; j++){k = (i + 1) * numColumns + j + 1;elements[k] = elements[i * numColumns + j] - c[i] * elements[j + 1];elements[k] = elements[k] + c[numColumns - j - 2] * elements[numColumns - i - 1];}}return true;}/// <summary>/// 求行列式值的全选主元高斯消去法/// </summary>/// <returns>double型，行列式的值</returns>public double ComputeDetGauss(){int i, j, k, nis = 0, js = 0, l, u, v;double f, det, q, d;//初值f = 1.0;det = 1.0;//消元for (k = 0; k <= numColumns - 2; k++){q = 0.0;for (i = k; i <= numColumns - 1; i++){for (j = k; j <= numColumns - 1; j++){l = i * numColumns + j;d = Math.Abs(elements[l]);if (d > q){q = d;nis = i;js = j;}}}if (q == 0.0){det = 0.0;return (det);}if (nis != k){f = -f;for (j = k; j <= numColumns - 1; j++){u = k * numColumns + j;v = nis * numColumns + j;d = elements[u];elements[u] = elements[v];elements[v] = d;}}if (js != k){f = -f;for (i = k; i <= numColumns - 1; i++){u = i * numColumns + js;v = i * numColumns + k;d = elements[u];elements[u] = elements[v];elements[v] = d;}}l = k * numColumns + k;det = det * elements[l];for (i = k + 1; i <= numColumns - 1; i++){d = elements[i * numColumns + k] / elements[l];for (j = k + 1; j <= numColumns - 1; j++){u = i * numColumns + j;elements[u] = elements[u] - d * elements[k * numColumns + j];}}}//求值det = f * det * elements[numColumns * numColumns - 1];return (det);}/// <summary>/// 求矩阵秩的全选主元高斯消去法/// </summary>/// <returns>int型，矩阵的秩</returns>public int ComputeRankGauss(){int i, j, k, nn, nis = 0, js = 0, l, ll, u, v;double q, d;//秩小于等于行列数nn = numRows;if (numRows >= numColumns)nn = numColumns;k = 0;//消元求解for (l = 0; l <= nn - 1; l++){q = 0.0;for (i = l; i <= numRows - 1; i++){for (j = l; j <= numColumns - 1; j++){ll = i * numColumns + j;d = Math.Abs(elements[ll]);if (d > q){q = d;nis = i;js = j;}}}if (q == 0.0)return (k);k = k + 1;if (nis != l){for (j = l; j <= numColumns - 1; j++){u = l * numColumns + j;v = nis * numColumns + j;d = elements[u];elements[u] = elements[v];elements[v] = d;}}if (js != l){for (i = l; i <= numRows - 1; i++){u = i * numColumns + js;v = i * numColumns + l;d = elements[u];elements[u] = elements[v];elements[v] = d;}}ll = l * numColumns + l;for (i = l + 1; i <= numColumns - 1; i++){d = elements[i * numColumns + l] / elements[ll];for (j = l + 1; j <= numColumns - 1; j++){u = i * numColumns + j;elements[u] = elements[u] - d * elements[l * numColumns + j];}}}return (k);}/// <summary>/// 对称正定矩阵的乔里斯基分解与行列式的求值/// </summary>/// <param name="realDetValue">返回行列式的值</param>/// <returns>bool型，求解是否成功</returns>public bool ComputeDetCholesky(ref double realDetValue){int i, j, k, u, l;double d;//不满足求解要求if (elements[0] <= 0.0)return false;//乔里斯基分解elements[0] = Math.Sqrt(elements[0]);d = elements[0];for (i = 1; i <= numColumns - 1; i++){u = i * numColumns;elements[u] = elements[u] / elements[0];}for (j = 1; j <= numColumns - 1; j++){l = j * numColumns + j;for (k = 0; k <= j - 1; k++){u = j * numColumns + k;elements[l] = elements[l] - elements[u] * elements[u];}if (elements[l] <= 0.0)return false;elements[l] = Math.Sqrt(elements[l]);d = d * elements[l];for (i = j + 1; i <= numColumns - 1; i++){u = i * numColumns + j;for (k = 0; k <= j - 1; k++)elements[u] = elements[u] - elements[i * numColumns + k] * elements[j * numColumns + k];elements[u] = elements[u] / elements[l];}}//行列式求值realDetValue = d * d;//下三角矩阵for (i = 0; i <= numColumns - 2; i++)for (j = i + 1; j <= numColumns - 1; j++)elements[i * numColumns + j] = 0.0;return true;}/// <summary>/// 矩阵的三角分解，分解成功后，原矩阵将成为Q矩阵/// </summary>/// <param name="mtxL">返回分解后的L矩阵</param>/// <param name="mtxU">返回分解后的U矩阵</param>/// <returns>bool型，求解是否成功</returns>public bool SplitLU(Matrix mtxL, Matrix mtxU){int i, j, k, w, v, ll;//初始化结果矩阵if (!mtxL.Init(numColumns, numColumns) ||!mtxU.Init(numColumns, numColumns))return false;for (k = 0; k <= numColumns - 2; k++){ll = k * numColumns + k;if (elements[ll] == 0.0)return false;for (i = k + 1; i <= numColumns - 1; i++){w = i * numColumns + k;elements[w] = elements[w] / elements[ll];}for (i = k + 1; i <= numColumns - 1; i++){w = i * numColumns + k;for (j = k + 1; j <= numColumns - 1; j++){v = i * numColumns + j;elements[v] = elements[v] - elements[w] * elements[k * numColumns + j];}}}for (i = 0; i <= numColumns - 1; i++){for (j = 0; j < i; j++){w = i * numColumns + j;mtxL.elements[w] = elements[w];mtxU.elements[w] = 0.0;}w = i * numColumns + i;mtxL.elements[w] = 1.0;mtxU.elements[w] = elements[w];for (j = i + 1; j <= numColumns - 1; j++){w = i * numColumns + j;mtxL.elements[w] = 0.0;mtxU.elements[w] = elements[w];}}return true;}/// <summary>/// 一般实矩阵的QR分解，分解成功后，原矩阵将成为R矩阵/// </summary>/// <param name="mtxQ">返回分解后的Q矩阵</param>/// <returns>bool型，求解是否成功</returns>public bool SplitQR(Matrix mtxQ){int i, j, k, l, nn, p, jj;double u, alpha, w, t;//初始化Q矩阵if (!mtxQ.Init(numRows, numRows))return false;//对角线元素单位化for (i = 0; i <= numRows - 1; i++){for (j = 0; j <= numRows - 1; j++){l = i * numRows + j;mtxQ.elements[l] = 0.0;if (i == j)mtxQ.elements[l] = 1.0;}}//开始分解nn = numColumns;if (numRows == numColumns)nn = numRows - 1;for (k = 0; k <= nn - 1; k++){u = 0.0;l = k * numColumns + k;for (i = k; i <= numRows - 1; i++){w = Math.Abs(elements[i * numColumns + k]);if (w > u)u = w;}alpha = 0.0;for (i = k; i <= numRows - 1; i++){t = elements[i * numColumns + k] / u;alpha = alpha + t * t;}if (elements[l] > 0.0)u = -u;alpha = u * Math.Sqrt(alpha);if (alpha == 0.0)return false;u = Math.Sqrt(2.0 * alpha * (alpha - elements[l]));if ((u + 1.0) != 1.0){elements[l] = (elements[l] - alpha) / u;for (i = k + 1; i <= numRows - 1; i++){p = i * numColumns + k;elements[p] = elements[p] / u;}for (j = 0; j <= numRows - 1; j++){t = 0.0;for (jj = k; jj <= numRows - 1; jj++)t = t + elements[jj * numColumns + k] * mtxQ.elements[jj * numRows + j];for (i = k; i <= numRows - 1; i++){p = i * numRows + j;mtxQ.elements[p] = mtxQ.elements[p] - 2.0 * t * elements[i * numColumns + k];}}for (j = k + 1; j <= numColumns - 1; j++){t = 0.0;for (jj = k; jj <= numRows - 1; jj++)t = t + elements[jj * numColumns + k] * elements[jj * numColumns + j];for (i = k; i <= numRows - 1; i++){p = i * numColumns + j;elements[p] = elements[p] - 2.0 * t * elements[i * numColumns + k];}}elements[l] = alpha;for (i = k + 1; i <= numRows - 1; i++)elements[i * numColumns + k] = 0.0;}}//调整元素for (i = 0; i <= numRows - 2; i++){for (j = i + 1; j <= numRows - 1; j++){p = i * numRows + j;l = j * numRows + i;t = mtxQ.elements[p];mtxQ.elements[p] = mtxQ.elements[l];mtxQ.elements[l] = t;}}return true;}/// <summary>/// 一般实矩阵的奇异值分解，分解成功后，原矩阵对角线元素就是矩阵的奇异值/// </summary>/// <param name="mtxU">返回分解后的U矩阵</param>/// <param name="mtxV">返回分解后的V矩阵</param>/// <param name="eps">计算精度</param>/// <returns>bool型，求解是否成功</returns>public bool SplitUV(Matrix mtxU, Matrix mtxV, double eps){int i, j, k, l, it, ll, kk, ix, iy, mm, nn, iz, m1, ks;double d, dd, t, sm, sm1, em1, sk, ek, b, c, shh;double[] fg = new double[2];double[] cs = new double[2];int m = numRows;int n = numColumns;//初始化U, V矩阵if (!mtxU.Init(m, m) || !mtxV.Init(n, n))return false;//临时缓冲区int ka = Math.Max(m, n) + 1;double[] s = new double[ka];double[] e = new double[ka];double[] w = new double[ka];//指定迭代次数为60it = 60;k = n;if (m - 1 < n)k = m - 1;l = m;if (n - 2 < m)l = n - 2;if (l < 0)l = 0;//循环迭代计算ll = k;if (l > k)ll = l;if (ll >= 1){for (kk = 1; kk <= ll; kk++){if (kk <= k){d = 0.0;for (i = kk; i <= m; i++){ix = (i - 1) * n + kk - 1;d = d + elements[ix] * elements[ix];}s[kk - 1] = Math.Sqrt(d);if (s[kk - 1] != 0.0){ix = (kk - 1) * n + kk - 1;if (elements[ix] != 0.0){s[kk - 1] = Math.Abs(s[kk - 1]);if (elements[ix] < 0.0)s[kk - 1] = -s[kk - 1];}for (i = kk; i <= m; i++){iy = (i - 1) * n + kk - 1;elements[iy] = elements[iy] / s[kk - 1];}elements[ix] = 1.0 + elements[ix];}s[kk - 1] = -s[kk - 1];}if (n >= kk + 1){for (j = kk + 1; j <= n; j++){if ((kk <= k) && (s[kk - 1] != 0.0)){d = 0.0;for (i = kk; i <= m; i++){ix = (i - 1) * n + kk - 1;iy = (i - 1) * n + j - 1;d = d + elements[ix] * elements[iy];}d = -d / elements[(kk - 1) * n + kk - 1];for (i = kk; i <= m; i++){ix = (i - 1) * n + j - 1;iy = (i - 1) * n + kk - 1;elements[ix] = elements[ix] + d * elements[iy];}}e[j - 1] = elements[(kk - 1) * n + j - 1];}}if (kk <= k){for (i = kk; i <= m; i++){ix = (i - 1) * m + kk - 1;iy = (i - 1) * n + kk - 1;mtxU.elements[ix] = elements[iy];}}if (kk <= l){d = 0.0;for (i = kk + 1; i <= n; i++)d = d + e[i - 1] * e[i - 1];e[kk - 1] = Math.Sqrt(d);if (e[kk - 1] != 0.0){if (e[kk] != 0.0){e[kk - 1] = Math.Abs(e[kk - 1]);if (e[kk] < 0.0)e[kk - 1] = -e[kk - 1];}for (i = kk + 1; i <= n; i++)e[i - 1] = e[i - 1] / e[kk - 1];e[kk] = 1.0 + e[kk];}e[kk - 1] = -e[kk - 1];if ((kk + 1 <= m) && (e[kk - 1] != 0.0)){for (i = kk + 1; i <= m; i++)w[i - 1] = 0.0;for (j = kk + 1; j <= n; j++)for (i = kk + 1; i <= m; i++)w[i - 1] = w[i - 1] + e[j - 1] * elements[(i - 1) * n + j - 1];for (j = kk + 1; j <= n; j++){for (i = kk + 1; i <= m; i++){ix = (i - 1) * n + j - 1;elements[ix] = elements[ix] - w[i - 1] * e[j - 1] / e[kk];}}}for (i = kk + 1; i <= n; i++)mtxV.elements[(i - 1) * n + kk - 1] = e[i - 1];}}}mm = n;if (m + 1 < n)mm = m + 1;if (k < n)s[k] = elements[k * n + k];if (m < mm)s[mm - 1] = 0.0;if (l + 1 < mm)e[l] = elements[l * n + mm - 1];e[mm - 1] = 0.0;nn = m;if (m > n)nn = n;if (nn >= k + 1){for (j = k + 1; j <= nn; j++){for (i = 1; i <= m; i++)mtxU.elements[(i - 1) * m + j - 1] = 0.0;mtxU.elements[(j - 1) * m + j - 1] = 1.0;}}if (k >= 1){for (ll = 1; ll <= k; ll++){kk = k - ll + 1;iz = (kk - 1) * m + kk - 1;if (s[kk - 1] != 0.0){if (nn >= kk + 1){for (j = kk + 1; j <= nn; j++){d = 0.0;for (i = kk; i <= m; i++){ix = (i - 1) * m + kk - 1;iy = (i - 1) * m + j - 1;d = d + mtxU.elements[ix] * mtxU.elements[iy] / mtxU.elements[iz];}d = -d;for (i = kk; i <= m; i++){ix = (i - 1) * m + j - 1;iy = (i - 1) * m + kk - 1;mtxU.elements[ix] = mtxU.elements[ix] + d * mtxU.elements[iy];}}}for (i = kk; i <= m; i++){ix = (i - 1) * m + kk - 1;mtxU.elements[ix] = -mtxU.elements[ix];}mtxU.elements[iz] = 1.0 + mtxU.elements[iz];if (kk - 1 >= 1){for (i = 1; i <= kk - 1; i++)mtxU.elements[(i - 1) * m + kk - 1] = 0.0;}}else{for (i = 1; i <= m; i++)mtxU.elements[(i - 1) * m + kk - 1] = 0.0;mtxU.elements[(kk - 1) * m + kk - 1] = 1.0;}}}for (ll = 1; ll <= n; ll++){kk = n - ll + 1;iz = kk * n + kk - 1;if ((kk <= l) && (e[kk - 1] != 0.0)){for (j = kk + 1; j <= n; j++){d = 0.0;for (i = kk + 1; i <= n; i++){ix = (i - 1) * n + kk - 1;iy = (i - 1) * n + j - 1;d = d + mtxV.elements[ix] * mtxV.elements[iy] / mtxV.elements[iz];}d = -d;for (i = kk + 1; i <= n; i++){ix = (i - 1) * n + j - 1;iy = (i - 1) * n + kk - 1;mtxV.elements[ix] = mtxV.elements[ix] + d * mtxV.elements[iy];}}}for (i = 1; i <= n; i++)mtxV.elements[(i - 1) * n + kk - 1] = 0.0;mtxV.elements[iz - n] = 1.0;}for (i = 1; i <= m; i++)for (j = 1; j <= n; j++)elements[(i - 1) * n + j - 1] = 0.0;m1 = mm;it = 60;while (true){if (mm == 0){ppp(elements, e, s, mtxV.elements, m, n);return true;}if (it == 0){ppp(elements, e, s, mtxV.elements, m, n);return false;}kk = mm - 1;while ((kk != 0) && (Math.Abs(e[kk - 1]) != 0.0)){d = Math.Abs(s[kk - 1]) + Math.Abs(s[kk]);dd = Math.Abs(e[kk - 1]);if (dd > eps * d)kk = kk - 1;elsee[kk - 1] = 0.0;}if (kk == mm - 1){kk = kk + 1;if (s[kk - 1] < 0.0){s[kk - 1] = -s[kk - 1];for (i = 1; i <= n; i++){ix = (i - 1) * n + kk - 1;mtxV.elements[ix] = -mtxV.elements[ix];}}while ((kk != m1) && (s[kk - 1] < s[kk])){d = s[kk - 1];s[kk - 1] = s[kk];s[kk] = d;if (kk < n){for (i = 1; i <= n; i++){ix = (i - 1) * n + kk - 1;iy = (i - 1) * n + kk;d = mtxV.elements[ix];mtxV.elements[ix] = mtxV.elements[iy];mtxV.elements[iy] = d;}}if (kk < m){for (i = 1; i <= m; i++){ix = (i - 1) * m + kk - 1;iy = (i - 1) * m + kk;d = mtxU.elements[ix];mtxU.elements[ix] = mtxU.elements[iy];mtxU.elements[iy] = d;}}kk = kk + 1;}it = 60;mm = mm - 1;}else{ks = mm;while ((ks > kk) && (Math.Abs(s[ks - 1]) != 0.0)){d = 0.0;if (ks != mm)d = d + Math.Abs(e[ks - 1]);if (ks != kk + 1)d = d + Math.Abs(e[ks - 2]);dd = Math.Abs(s[ks - 1]);if (dd > eps * d)ks = ks - 1;elses[ks - 1] = 0.0;}if (ks == kk){kk = kk + 1;d = Math.Abs(s[mm - 1]);t = Math.Abs(s[mm - 2]);if (t > d)d = t;t = Math.Abs(e[mm - 2]);if (t > d)d = t;t = Math.Abs(s[kk - 1]);if (t > d)d = t;t = Math.Abs(e[kk - 1]);if (t > d)d = t;sm = s[mm - 1] / d;sm1 = s[mm - 2] / d;em1 = e[mm - 2] / d;sk = s[kk - 1] / d;ek = e[kk - 1] / d;b = ((sm1 + sm) * (sm1 - sm) + em1 * em1) / 2.0;c = sm * em1;c = c * c;shh = 0.0;if ((b != 0.0) || (c != 0.0)){shh = Math.Sqrt(b * b + c);if (b < 0.0)shh = -shh;shh = c / (b + shh);}fg[0] = (sk + sm) * (sk - sm) - shh;fg[1] = sk * ek;for (i = kk; i <= mm - 1; i++){sss(fg, cs);if (i != kk)e[i - 2] = fg[0];fg[0] = cs[0] * s[i - 1] + cs[1] * e[i - 1];e[i - 1] = cs[0] * e[i - 1] - cs[1] * s[i - 1];fg[1] = cs[1] * s[i];s[i] = cs[0] * s[i];if ((cs[0] != 1.0) || (cs[1] != 0.0)){for (j = 1; j <= n; j++){ix = (j - 1) * n + i - 1;iy = (j - 1) * n + i;d = cs[0] * mtxV.elements[ix] + cs[1] * mtxV.elements[iy];mtxV.elements[iy] = -cs[1] * mtxV.elements[ix] + cs[0] * mtxV.elements[iy];mtxV.elements[ix] = d;}}sss(fg, cs);s[i - 1] = fg[0];fg[0] = cs[0] * e[i - 1] + cs[1] * s[i];s[i] = -cs[1] * e[i - 1] + cs[0] * s[i];fg[1] = cs[1] * e[i];e[i] = cs[0] * e[i];if (i < m){if ((cs[0] != 1.0) || (cs[1] != 0.0)){for (j = 1; j <= m; j++){ix = (j - 1) * m + i - 1;iy = (j - 1) * m + i;d = cs[0] * mtxU.elements[ix] + cs[1] * mtxU.elements[iy];mtxU.elements[iy] = -cs[1] * mtxU.elements[ix] + cs[0] * mtxU.elements[iy];mtxU.elements[ix] = d;}}}}e[mm - 2] = fg[0];it = it - 1;}else{if (ks == mm){kk = kk + 1;fg[1] = e[mm - 2];e[mm - 2] = 0.0;for (ll = kk; ll <= mm - 1; ll++){i = mm + kk - ll - 1;fg[0] = s[i - 1];sss(fg, cs);s[i - 1] = fg[0];if (i != kk){fg[1] = -cs[1] * e[i - 2];e[i - 2] = cs[0] * e[i - 2];}if ((cs[0] != 1.0) || (cs[1] != 0.0)){for (j = 1; j <= n; j++){ix = (j - 1) * n + i - 1;iy = (j - 1) * n + mm - 1;d = cs[0] * mtxV.elements[ix] + cs[1] * mtxV.elements[iy];mtxV.elements[iy] = -cs[1] * mtxV.elements[ix] + cs[0] * mtxV.elements[iy];mtxV.elements[ix] = d;}}}}else{kk = ks + 1;fg[1] = e[kk - 2];e[kk - 2] = 0.0;for (i = kk; i <= mm; i++){fg[0] = s[i - 1];sss(fg, cs);s[i - 1] = fg[0];fg[1] = -cs[1] * e[i - 1];e[i - 1] = cs[0] * e[i - 1];if ((cs[0] != 1.0) || (cs[1] != 0.0)){for (j = 1; j <= m; j++){ix = (j - 1) * m + i - 1;iy = (j - 1) * m + kk - 2;d = cs[0] * mtxU.elements[ix] + cs[1] * mtxU.elements[iy];mtxU.elements[iy] = -cs[1] * mtxU.elements[ix] + cs[0] * mtxU.elements[iy];mtxU.elements[ix] = d;}}}}}}}}/// <summary>/// 内部函数，由SplitUV函数调用/// </summary>/// <param name="a"></param>/// <param name="e"></param>/// <param name="s"></param>/// <param name="v"></param>/// <param name="m"></param>/// <param name="n"></param>private void ppp(double[] a, double[] e, double[] s, double[] v, int m, int n){int i, j, p, q;double d;if (m >= n)i = n;elsei = m;for (j = 1; j <= i - 1; j++){a[(j - 1) * n + j - 1] = s[j - 1];a[(j - 1) * n + j] = e[j - 1];}a[(i - 1) * n + i - 1] = s[i - 1];if (m < n)a[(i - 1) * n + i] = e[i - 1];for (i = 1; i <= n - 1; i++){for (j = i + 1; j <= n; j++){p = (i - 1) * n + j - 1;q = (j - 1) * n + i - 1;d = v[p];v[p] = v[q];v[q] = d;}}}/// <summary>/// 内部函数，由SplitUV函数调用/// </summary>/// <param name="fg"></param>/// <param name="cs"></param>private void sss(double[] fg, double[] cs){double r, d;if ((Math.Abs(fg[0]) + Math.Abs(fg[1])) == 0.0){cs[0] = 1.0;cs[1] = 0.0;d = 0.0;}else{d = Math.Sqrt(fg[0] * fg[0] + fg[1] * fg[1]);if (Math.Abs(fg[0]) > Math.Abs(fg[1])){d = Math.Abs(d);if (fg[0] < 0.0)d = -d;}if (Math.Abs(fg[1]) >= Math.Abs(fg[0])){d = Math.Abs(d);if (fg[1] < 0.0)d = -d;}cs[0] = fg[0] / d;cs[1] = fg[1] / d;}r = 1.0;if (Math.Abs(fg[0]) > Math.Abs(fg[1]))r = cs[1];else if (cs[0] != 0.0)r = 1.0 / cs[0];fg[0] = d;fg[1] = r;}/// <summary>/// 求广义逆的奇异值分解法，分解成功后，原矩阵对角线元素就是矩阵的奇异值/// </summary>/// <param name="mtxAP">返回原矩阵的广义逆矩阵</param>/// <param name="mtxU">返回分解后的U矩阵</param>/// <param name="mtxV">返回分解后的V矩阵</param>/// <param name="eps">计算精度</param>/// <returns>bool型，求解是否成功</returns>public bool InvertUV(Matrix mtxAP, Matrix mtxU, Matrix mtxV, double eps){int i, j, k, l, t, p, q, f;//调用奇异值分解if (!SplitUV(mtxU, mtxV, eps))return false;int m = numRows;int n = numColumns;//初始化广义逆矩阵if (!mtxAP.Init(n, m))return false;//计算广义逆矩阵j = n;if (m < n)j = m;j = j - 1;k = 0;while ((k <= j) && (elements[k * n + k] != 0.0))k = k + 1;k = k - 1;for (i = 0; i <= n - 1; i++){for (j = 0; j <= m - 1; j++){t = i * m + j;mtxAP.elements[t] = 0.0;for (l = 0; l <= k; l++){f = l * n + i;p = j * m + l;q = l * n + l;mtxAP.elements[t] = mtxAP.elements[t] + mtxV.elements[f] * mtxU.elements[p] / elements[q];}}}return true;}/// <summary>/// 约化对称矩阵为对称三对角阵的豪斯荷尔德变换法/// </summary>/// <param name="mtxQ">返回豪斯荷尔德变换的乘积矩阵Q</param>/// <param name="mtxT">返回求得的对称三对角阵</param>/// <param name="dblB">一维数组，长度为矩阵的阶数，返回对称三对角阵的主对角线元素</param>/// <param name="dblC">一维数组，长度为矩阵的阶数，前n-1个元素返回对称三对角阵的次对角线元素</param>/// <returns>bool型，求解是否成功</returns>public bool MakeSymTri(Matrix mtxQ, Matrix mtxT, double[] dblB, double[] dblC){int i, j, k, u;double h, f, g, h2;//初始化矩阵Q和Tif (!mtxQ.Init(numColumns, numColumns) ||!mtxT.Init(numColumns, numColumns))return false;if (dblB == null || dblC == null)return false;for (i = 0; i <= numColumns - 1; i++){for (j = 0; j <= numColumns - 1; j++){u = i * numColumns + j;mtxQ.elements[u] = elements[u];}}for (i = numColumns - 1; i >= 1; i--){h = 0.0;if (i > 1){for (k = 0; k <= i - 1; k++){u = i * numColumns + k;h = h + mtxQ.elements[u] * mtxQ.elements[u];}}if (h == 0.0){dblC[i] = 0.0;if (i == 1)dblC[i] = mtxQ.elements[i * numColumns + i - 1];dblB[i] = 0.0;}else{dblC[i] = Math.Sqrt(h);u = i * numColumns + i - 1;if (mtxQ.elements[u] > 0.0)dblC[i] = -dblC[i];h = h - mtxQ.elements[u] * dblC[i];mtxQ.elements[u] = mtxQ.elements[u] - dblC[i];f = 0.0;for (j = 0; j <= i - 1; j++){mtxQ.elements[j * numColumns + i] = mtxQ.elements[i * numColumns + j] / h;g = 0.0;for (k = 0; k <= j; k++)g = g + mtxQ.elements[j * numColumns + k] * mtxQ.elements[i * numColumns + k];if (j + 1 <= i - 1)for (k = j + 1; k <= i - 1; k++)g = g + mtxQ.elements[k * numColumns + j] * mtxQ.elements[i * numColumns + k];dblC[j] = g / h;f = f + g * mtxQ.elements[j * numColumns + i];}h2 = f / (h + h);for (j = 0; j <= i - 1; j++){f = mtxQ.elements[i * numColumns + j];g = dblC[j] - h2 * f;dblC[j] = g;for (k = 0; k <= j; k++){u = j * numColumns + k;mtxQ.elements[u] = mtxQ.elements[u] - f * dblC[k] - g * mtxQ.elements[i * numColumns + k];}}dblB[i] = h;}}for (i = 0; i <= numColumns - 2; i++)dblC[i] = dblC[i + 1];dblC[numColumns - 1] = 0.0;dblB[0] = 0.0;for (i = 0; i <= numColumns - 1; i++){if ((dblB[i] != (double)0.0) && (i - 1 >= 0)){for (j = 0; j <= i - 1; j++){g = 0.0;for (k = 0; k <= i - 1; k++)g = g + mtxQ.elements[i * numColumns + k] * mtxQ.elements[k * numColumns + j];for (k = 0; k <= i - 1; k++){u = k * numColumns + j;mtxQ.elements[u] = mtxQ.elements[u] - g * mtxQ.elements[k * numColumns + i];}}}u = i * numColumns + i;dblB[i] = mtxQ.elements[u]; mtxQ.elements[u] = 1.0;if (i - 1 >= 0){for (j = 0; j <= i - 1; j++){mtxQ.elements[i * numColumns + j] = 0.0;mtxQ.elements[j * numColumns + i] = 0.0;}}}//构造对称三对角矩阵for (i = 0; i < numColumns; ++i){for (j = 0; j < numColumns; ++j){mtxT.SetElement(i, j, 0);k = i - j;if (k == 0)mtxT.SetElement(i, j, dblB[j]);else if (k == 1)mtxT.SetElement(i, j, dblC[j]);else if (k == -1)mtxT.SetElement(i, j, dblC[i]);}}return true;}/// <summary>/// 实对称三对角阵的全部特征值与特征向量的计算/// </summary>/// <param name="dblB">一维数组，长度为矩阵的阶数，传入对称三对角阵的主对角线元素；返回时存放全部特征值。</param>/// <param name="dblC">一维数组，长度为矩阵的阶数，前n-1个元素传入对称三对角阵的次对角线元素</param>/// <param name="mtxQ">如果传入单位矩阵，则返回实对称三对角阵的特征值向量矩阵；/// 如果传入MakeSymTri函数求得的矩阵A的豪斯荷尔德变换的乘积/// 矩阵Q，则返回矩阵A的特征值向量矩阵。其中第i列为与数组dblB/// 中第j个特征值对应的特征向量。</param>/// <param name="nMaxIt">迭代次数</param>/// <param name="eps">计算精度</param>/// <returns>bool型，求解是否成功</returns>public bool ComputeEvSymTri(double[] dblB, double[] dblC, Matrix mtxQ, int nMaxIt, double eps){int i, j, k, m, it, u, v;double d, f, h, g, p, r, e, s;//初值int n = mtxQ.GetNumColumns();dblC[n - 1] = 0.0;d = 0.0;f = 0.0;//迭代计算for (j = 0; j <= n - 1; j++){it = 0;h = eps * (Math.Abs(dblB[j]) + Math.Abs(dblC[j]));if (h > d)d = h;m = j;while ((m <= n - 1) && (Math.Abs(dblC[m]) > d))m = m + 1;if (m != j){do{if (it == nMaxIt)return false;it = it + 1;g = dblB[j];p = (dblB[j + 1] - g) / (2.0 * dblC[j]);r = Math.Sqrt(p * p + 1.0);if (p >= 0.0)dblB[j] = dblC[j] / (p + r);elsedblB[j] = dblC[j] / (p - r);h = g - dblB[j];for (i = j + 1; i <= n - 1; i++)dblB[i] = dblB[i] - h;f = f + h;p = dblB[m];e = 1.0;s = 0.0;for (i = m - 1; i >= j; i--){g = e * dblC[i];h = e * p;if (Math.Abs(p) >= Math.Abs(dblC[i])){e = dblC[i] / p;r = Math.Sqrt(e * e + 1.0);dblC[i + 1] = s * p * r;s = e / r;e = 1.0 / r;}else{e = p / dblC[i];r = Math.Sqrt(e * e + 1.0);dblC[i + 1] = s * dblC[i] * r;s = 1.0 / r;e = e / r;}p = e * dblB[i] - s * g;dblB[i + 1] = h + s * (e * g + s * dblB[i]);for (k = 0; k <= n - 1; k++){u = k * n + i + 1;v = u - 1;h = mtxQ.elements[u];mtxQ.elements[u] = s * mtxQ.elements[v] + e * h;mtxQ.elements[v] = e * mtxQ.elements[v] - s * h;}}dblC[j] = s * p;dblB[j] = e * p;} while (Math.Abs(dblC[j]) > d);}dblB[j] = dblB[j] + f;}for (i = 0; i <= n - 1; i++){k = i;p = dblB[i];if (i + 1 <= n - 1){j = i + 1;while ((j <= n - 1) && (dblB[j] <= p)){k = j;p = dblB[j];j = j + 1;}}if (k != i){dblB[k] = dblB[i];dblB[i] = p;for (j = 0; j <= n - 1; j++){u = j * n + i;v = j * n + k;p = mtxQ.elements[u];mtxQ.elements[u] = mtxQ.elements[v];mtxQ.elements[v] = p;}}}return true;}/// <summary>/// 约化一般实矩阵为赫申伯格矩阵的初等相似变换法/// </summary>public void MakeHberg(){int i = 0, j, k, u, v;double d, t;for (k = 1; k <= numColumns - 2; k++){d = 0.0;for (j = k; j <= numColumns - 1; j++){u = j * numColumns + k - 1;t = elements[u];if (Math.Abs(t) > Math.Abs(d)){d = t;i = j;}}if (d != 0.0){if (i != k){for (j = k - 1; j <= numColumns - 1; j++){u = i * numColumns + j;v = k * numColumns + j;t = elements[u];elements[u] = elements[v];elements[v] = t;}for (j = 0; j <= numColumns - 1; j++){u = j * numColumns + i;v = j * numColumns + k;t = elements[u];elements[u] = elements[v];elements[v] = t;}}for (i = k + 1; i <= numColumns - 1; i++){u = i * numColumns + k - 1;t = elements[u] / d;elements[u] = 0.0;for (j = k; j <= numColumns - 1; j++){v = i * numColumns + j;elements[v] = elements[v] - t * elements[k * numColumns + j];}for (j = 0; j <= numColumns - 1; j++){v = j * numColumns + k;elements[v] = elements[v] + t * elements[j * numColumns + i];}}}}}/// <summary>/// 求赫申伯格矩阵全部特征值的QR方法/// </summary>/// <param name="dblU">一维数组，长度为矩阵的阶数，返回时存放特征值的实部</param>/// <param name="dblV">一维数组，长度为矩阵的阶数，返回时存放特征值的虚部</param>/// <param name="nMaxIt">迭代次数</param>/// <param name="eps">计算精度</param>/// <returns>bool型，求解是否成功</returns>public bool ComputeEvHBerg(double[] dblU, double[] dblV, int nMaxIt, double eps){int m, it, i, j, k, l, ii, jj, kk, ll;double b, c, w, g, xy, p, q, r, x, s, e, f, z, y;int n = numColumns;it = 0;m = n;while (m != 0){l = m - 1;while ((l > 0) && (Math.Abs(elements[l * n + l - 1]) >eps * (Math.Abs(elements[(l - 1) * n + l - 1]) + Math.Abs(elements[l * n + l]))))l = l - 1;ii = (m - 1) * n + m - 1;jj = (m - 1) * n + m - 2;kk = (m - 2) * n + m - 1;ll = (m - 2) * n + m - 2;if (l == m - 1){dblU[m - 1] = elements[(m - 1) * n + m - 1];dblV[m - 1] = 0.0;m = m - 1;it = 0;}else if (l == m - 2){b = -(elements[ii] + elements[ll]);c = elements[ii] * elements[ll] - elements[jj] * elements[kk];w = b * b - 4.0 * c;y = Math.Sqrt(Math.Abs(w));if (w > 0.0){xy = 1.0;if (b < 0.0)xy = -1.0;dblU[m - 1] = (-b - xy * y) / 2.0;dblU[m - 2] = c / dblU[m - 1];dblV[m - 1] = 0.0; dblV[m - 2] = 0.0;}else{dblU[m - 1] = -b / 2.0;dblU[m - 2] = dblU[m - 1];dblV[m - 1] = y / 2.0;dblV[m - 2] = -dblV[m - 1];}m = m - 2;it = 0;}else{if (it >= nMaxIt)return false;it = it + 1;for (j = l + 2; j <= m - 1; j++)elements[j * n + j - 2] = 0.0;for (j = l + 3; j <= m - 1; j++)elements[j * n + j - 3] = 0.0;for (k = l; k <= m - 2; k++){if (k != l){p = elements[k * n + k - 1];q = elements[(k + 1) * n + k - 1];r = 0.0;if (k != m - 2)r = elements[(k + 2) * n + k - 1];}else{x = elements[ii] + elements[ll];y = elements[ll] * elements[ii] - elements[kk] * elements[jj];ii = l * n + l;jj = l * n + l + 1;kk = (l + 1) * n + l;ll = (l + 1) * n + l + 1;p = elements[ii] * (elements[ii] - x) + elements[jj] * elements[kk] + y;q = elements[kk] * (elements[ii] + elements[ll] - x);r = elements[kk] * elements[(l + 2) * n + l + 1];}if ((Math.Abs(p) + Math.Abs(q) + Math.Abs(r)) != 0.0){xy = 1.0;if (p < 0.0)xy = -1.0;s = xy * Math.Sqrt(p * p + q * q + r * r);if (k != l)elements[k * n + k - 1] = -s;e = -q / s;f = -r / s;x = -p / s;y = -x - f * r / (p + s);g = e * r / (p + s);z = -x - e * q / (p + s);for (j = k; j <= m - 1; j++){ii = k * n + j;jj = (k + 1) * n + j;p = x * elements[ii] + e * elements[jj];q = e * elements[ii] + y * elements[jj];r = f * elements[ii] + g * elements[jj];if (k != m - 2){kk = (k + 2) * n + j;p = p + f * elements[kk];q = q + g * elements[kk];r = r + z * elements[kk];elements[kk] = r;}elements[jj] = q; elements[ii] = p;}j = k + 3;if (j >= m - 1)j = m - 1;for (i = l; i <= j; i++){ii = i * n + k;jj = i * n + k + 1;p = x * elements[ii] + e * elements[jj];q = e * elements[ii] + y * elements[jj];r = f * elements[ii] + g * elements[jj];if (k != m - 2){kk = i * n + k + 2;p = p + f * elements[kk];q = q + g * elements[kk];r = r + z * elements[kk];elements[kk] = r;}elements[jj] = q;elements[ii] = p;}}}}}return true;}/// <summary>/// 求实对称矩阵特征值与特征向量的雅可比法/// </summary>/// <param name="dblEigenValue">一维数组，长度为矩阵的阶数，返回时存放特征值</param>/// <param name="mtxEigenVector">返回时存放特征向量矩阵，其中第i列为与数组,/// dblEigenValue中第j个特征值对应的特征向量</param>/// <param name="nMaxIt">迭代次数</param>/// <param name="eps">计算精度</param>/// <returns>bool型，求解是否成功</returns>public bool ComputeEvJacobi(double[] dblEigenValue, Matrix mtxEigenVector, int nMaxIt, double eps){int i, j, p = 0, q = 0, u, w, t, s, l;double fm, cn, sn, omega, x, y, d;if (!mtxEigenVector.Init(numColumns, numColumns))return false;l = 1;for (i = 0; i <= numColumns - 1; i++){mtxEigenVector.elements[i * numColumns + i] = 1.0;for (j = 0; j <= numColumns - 1; j++)if (i != j)mtxEigenVector.elements[i * numColumns + j] = 0.0;}while (true){fm = 0.0;for (i = 1; i <= numColumns - 1; i++){for (j = 0; j <= i - 1; j++){d = Math.Abs(elements[i * numColumns + j]);if ((i != j) && (d > fm)){fm = d;p = i;q = j;}}}if (fm < eps){for (i = 0; i < numColumns; ++i)dblEigenValue[i] = GetElement(i, i);return true;}if (l > nMaxIt)return false;l = l + 1;u = p * numColumns + q;w = p * numColumns + p;t = q * numColumns + p;s = q * numColumns + q;x = -elements[u];y = (elements[s] - elements[w]) / 2.0;omega = x / Math.Sqrt(x * x + y * y);if (y < 0.0)omega = -omega;sn = 1.0 + Math.Sqrt(1.0 - omega * omega);sn = omega / Math.Sqrt(2.0 * sn);cn = Math.Sqrt(1.0 - sn * sn);fm = elements[w];elements[w] = fm * cn * cn + elements[s] * sn * sn + elements[u] * omega;elements[s] = fm * sn * sn + elements[s] * cn * cn - elements[u] * omega;elements[u] = 0.0;elements[t] = 0.0;for (j = 0; j <= numColumns - 1; j++){if ((j != p) && (j != q)){u = p * numColumns + j; w = q * numColumns + j;fm = elements[u];elements[u] = fm * cn + elements[w] * sn;elements[w] = -fm * sn + elements[w] * cn;}}for (i = 0; i <= numColumns - 1; i++){if ((i != p) && (i != q)){u = i * numColumns + p;w = i * numColumns + q;fm = elements[u];elements[u] = fm * cn + elements[w] * sn;elements[w] = -fm * sn + elements[w] * cn;}}for (i = 0; i <= numColumns - 1; i++){u = i * numColumns + p;w = i * numColumns + q;fm = mtxEigenVector.elements[u];mtxEigenVector.elements[u] = fm * cn + mtxEigenVector.elements[w] * sn;mtxEigenVector.elements[w] = -fm * sn + mtxEigenVector.elements[w] * cn;}}}/// <summary>/// 求实对称矩阵特征值与特征向量的雅可比过关法/// </summary>/// <param name="dblEigenValue">一维数组，长度为矩阵的阶数，返回时存放特征值</param>/// <param name="mtxEigenVector">返回时存放特征向量矩阵，其中第i列为与数组,/// dblEigenValue中第j个特征值对应的特征向量</param>/// <param name="eps">计算精度</param>/// <returns>bool型，求解是否成功</returns>public bool ComputeEvJacobi(double[] dblEigenValue, Matrix mtxEigenVector, double eps){int i, j, p, q, u, w, t, s;double ff, fm, cn, sn, omega, x, y, d;if (!mtxEigenVector.Init(numColumns, numColumns))return false;for (i = 0; i <= numColumns - 1; i++){mtxEigenVector.elements[i * numColumns + i] = 1.0;for (j = 0; j <= numColumns - 1; j++)if (i != j)mtxEigenVector.elements[i * numColumns + j] = 0.0;}ff = 0.0;for (i = 1; i <= numColumns - 1; i++){for (j = 0; j <= i - 1; j++){d = elements[i * numColumns + j];ff = ff + d * d;}}ff = Math.Sqrt(2.0 * ff);ff = ff / (1.0 * numColumns);bool nextLoop = false;while (true){for (i = 1; i <= numColumns - 1; i++){for (j = 0; j <= i - 1; j++){d = Math.Abs(elements[i * numColumns + j]);if (d > ff){p = i;q = j;u = p * numColumns + q;w = p * numColumns + p;t = q * numColumns + p;s = q * numColumns + q;x = -elements[u];y = (elements[s] - elements[w]) / 2.0;omega = x / Math.Sqrt(x * x + y * y);if (y < 0.0)omega = -omega;sn = 1.0 + Math.Sqrt(1.0 - omega * omega);sn = omega / Math.Sqrt(2.0 * sn);cn = Math.Sqrt(1.0 - sn * sn);fm = elements[w];elements[w] = fm * cn * cn + elements[s] * sn * sn + elements[u] * omega;elements[s] = fm * sn * sn + elements[s] * cn * cn - elements[u] * omega;elements[u] = 0.0; elements[t] = 0.0;for (j = 0; j <= numColumns - 1; j++){if ((j != p) && (j != q)){u = p * numColumns + j;w = q * numColumns + j;fm = elements[u];elements[u] = fm * cn + elements[w] * sn;elements[w] = -fm * sn + elements[w] * cn;}}for (i = 0; i <= numColumns - 1; i++){if ((i != p) && (i != q)){u = i * numColumns + p;w = i * numColumns + q;fm = elements[u];elements[u] = fm * cn + elements[w] * sn;elements[w] = -fm * sn + elements[w] * cn;}}for (i = 0; i <= numColumns - 1; i++){u = i * numColumns + p;w = i * numColumns + q;fm = mtxEigenVector.elements[u];mtxEigenVector.elements[u] = fm * cn + mtxEigenVector.elements[w] * sn;mtxEigenVector.elements[w] = -fm * sn + mtxEigenVector.elements[w] * cn;}nextLoop = true;break;}}if (nextLoop)break;}if (nextLoop){nextLoop = false;continue;}nextLoop = false;//如果达到精度要求，退出循环，返回结果if (ff < eps){for (i = 0; i < numColumns; ++i)dblEigenValue[i] = GetElement(i, i);return true;}ff = ff / (1.0 * numColumns);}}}
}
注：本代码为原作者所有，本人只是摘录，如有侵犯作者，请与本人联系
Email:359194966@qq.com
微软公司内部培训程序员资料---操作矩阵类相关推荐

微软公司内部培训程序员资料---求解线性方程组的类
/** 求解线性方程组的类 LEquations* * 周长发编制*/using System;namespace MSAlgorithm {public class LEquations{/// & ...
公司停电，程序员去网吧写代码；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上半年在公司内部培训的经历
在这里谈谈我在今年上半年在公司内部培训实习生(今年七月毕业的学生)的经历. 今年三月份我担任公司技术培训的负责人,回顾培训历程,做个总结. 公司目标是通过培训能让实习生快速掌握知识,成为一名优秀的技术 ...
公司KPI考核代码行数，程序员神操作：10行变500行！
"如果你无法度量,就无法管理." 这年头,谁都逃不过被KPI支配的恐惧. KPI,俗称绩效,全称关键绩效指标,也是领导口中常说的"小目标",往往"领导 ...
二：程序员资料大全-各种神奇的资料收集笔记
http://tools.zhaishidan.cn/ 资料篇技术站点必看书籍大牛博客 GitHub篇工具篇平台工具常用工具第三方服务爬虫相关(好玩的工具) 安全相关 Web服务器性能 ...
微软公司内部培训程序员资料---操作矩阵类

微软公司内部培训程序员资料---操作矩阵类相关推荐

最新文章

热门文章
