nr.cs

C#矩阵运算类

using System;
using System.Collections.Generic;
using System.Text;
using System.Drawing;

namespace NumericRcipes
{
    public class NR
    {
        public static double SQRT(double a, double b)////计算a，b的平方根
        {
            double absa, absb;
            absa = Math.Abs(a);
            absb = Math.Abs(b);
            if (absa > absb) return absa * Math.Sqrt(1.0 + (absb / absa) * (absb / absa));
            else return (absb == 0.0 ? 0.0 : absb * Math.Sqrt(1.0 + (absa / absb) * (absa / absb)));
        }
        public static double SIGN(double a, double b)
        {
            return b >= 0.0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a);
        }
        public static void SWAP<T>(ref T a, ref T b)//交换两个数据对象
        {
            T dum = a; a = b; b = dum;
        }
        public static bool Gauss_Jordan(out NRMatrix inverse, out NRMatrix result,
            NRMatrix left, NRMatrix right) //Gauss_Jordan全主元消元法求解矩阵方程，
        { //求解成功，inverse是逆矩阵,result是解向量。
            if (left.Rows != left.Cols) throw new NRException("系数矩阵不是方阵！");
            if (left.Rows != right.Rows) throw new NRException("左右两端矩阵的行数不同！");
            NRMatrix A = new NRMatrix(left);
            NRMatrix B = new NRMatrix(right);
            int icol = 0, irow = 0;
            double big, dum, pivinv;
            int n = A.Rows;
            int m = B.Cols;
            int[] indexc = new int[n]; int[] indexr = new int[n]; int[] ipiv = new int[n];//这些整型数组用于记录主元。

            for (int s = 0; s < n; s++)
            {//约化列的主循环。
                big = 0.0;
                for (int i = 0; i < n; i++)//寻找主元的外层循环。
                    if (ipiv[i] != 1)
                        for (int j = 0; j < n; j++)
                        {
                            if (ipiv[j] == 0)
                            {
                                if (Math.Abs(A[i, j]) >= big)
                                {
                                    big = Math.Abs(A[i, j]);
                                    irow = i;
                                    icol = j;
                                }
                            }
                        }
                ++(ipiv[icol]);
                //找到主元，交换到对角线上。列不进行实际的交换，只进行重新标注：
                if (irow != icol)
                {
                    A = NRMatrix.ExchangeRows(A, irow, icol);
                    B = NRMatrix.ExchangeRows(B, irow, icol);
                }
                indexr[s] = irow;
                indexc[s] = icol;
                if (A[icol, icol] == 0.0)
                {
                    Console.WriteLine("gaussj:奇异矩阵!");
                    inverse = null;
                    result = null;
                    return false;
                }
                pivinv = 1.0 / A[icol, icol];
                A[icol, icol] = 1.0;//归一化对角线上的行
                for (int k = 0; k < A.Cols; k++)
                    A[icol, k] *= pivinv;
                for (int k = 0; k < B.Cols; k++)
                    B[icol, k] *= pivinv;
                for (int i = 0; i < n; i++)//约化各列
                    if (i != icol)
                    {
                        dum = A[i, icol];
                        A[i, icol] = 0.0;
                        for (int k = 0; k < A.Cols; k++)
                            A[i, k] += A[icol, k] * (-dum);
                        for (int k = 0; k < B.Cols; k++)
                            B[i, k] += B[icol, k] * (-dum);
                    }
            }//这是列约化主循环的结尾
            //考虑已进行列交换，为了使解向量保持原来的次序，
            //再根据其交换的相反顺序交换各列，以整理得到原始方程的解。
            for (int j = (int)n - 1; j >= 0; j--)
            {
                if (indexr[j] != indexc[j])
                    A = NRMatrix.ExchangeCols(A, indexr[j], indexc[j]);
            }
            inverse = A;
            result = B;
            return true;
        }
        public static bool LUdcmp(out NRMatrix result, out int[] index, out int d, NRMatrix A) //矩阵的LU分解
        { //结果保存在该result中，index[0..n-1]为输出向量，用来记录因部分主元法而改变了的行排列次序。
            //d为输出变量，值为＋1或－1，分别表示进行行交换次数为偶数还是奇数
            const double TINY = 1.0e-20;//一个小数
            int i, imax = 0, j, k;
            double big, dum, sum, temp;
            result = new NRMatrix(A);
            int n = result.Cols;
            double[] vv = new double[n];//用于保存每行的内含比例因子
            index = new int[n];
            d = 1;
            for (i = 0; i < n; i++)
            {//按行循环，求内含比例因子。
                big = 0.0;
                for (j = 0; j < n; j++)
                    if ((temp = Math.Abs(result[i, j])) > big) big = temp;
                if (big == 0.0) throw new NRException("LU分解中遇到奇异矩阵！");
                vv[i] = 1.0 / big;//保存比例因子
            }
            for (j = 0; j < n; j++)
            {//Grout方法的列循环
                for (i = 0; i < j; i++)
                {
                    sum = result[i, j];
                    for (k = 0; k < i; k++) sum -= result[i, k] * result[k, j];
                    result[i, j] = sum;
                }
                big = 0.0;//初始化最大主元素变量
                for (i = j; i < n; i++)
                {
                    sum = result[i, j];
                    for (k = 0; k < j; k++) sum -= result[i, k] * result[k, j];
                    result[i, j] = sum;
                    if ((dum = vv[i] * Math.Abs(sum)) >= big)
                    {//新的候选主元是否比已选中的更合适？
                        big = dum;
                        imax = i;
                    }
                }
                if (j != imax)
                {//是否需要进行交换
                    for (k = 0; k < n; k++)
                    {
                        dum = result[imax, k];
                        result[imax, k] = result[j, k];
                        result[j, k] = dum;
                    }
                    d = -d;
                    vv[imax] = vv[j];
                }
                index[j] = imax;
                if (result[j, j] == 0.0) result[j, j] = TINY;
                if (j != n - 1)
                {
                    dum = 1.0 / (result[j, j]);
                    for (i = j + 1; i < n; i++) result[i, j] *= dum;
                }
            }
            return true;
        }
        public static bool LUdcmp(NRMatrix left, NRMatrix right, out NRMatrix result) //LU分解求解方程组解
        { //解矩阵保存在矩阵result中
            const double TINY = 1.0e-20;//一个小数
            int i, imax = 0, j, k, ip, ii = 0;
            double big, dum, sum, temp;
            NRMatrix A = new NRMatrix(left);
            NRMatrix b = new NRMatrix(right);
            int n = A.Cols;
            double[] vv = new double[n];//用于保存每行的内含比例因子
            int[] index = new int[n];
            for (i = 0; i < n; i++)
            {//按行循环，求内含比例因子。
                big = 0.0;
                for (j = 0; j < n; j++)
                    if ((temp = Math.Abs(A[i, j])) > big) big = temp;
                if (big == 0.0) throw new NRException("LU分解中遇到奇异矩阵！");
                vv[i] = 1.0 / big;//保存比例因子
            }
            for (j = 0; j < n; j++)
            {//Grout方法的列循环
                for (i = 0; i < j; i++)
                {
                    sum = A[i, j];
                    for (k = 0; k < i; k++) sum -= A[i, k] * A[k, j];
                    A[i, j] = sum;
                }
                big = 0.0;//初始化最大主元素变量
                for (i = j; i < n; i++)
                {
                    sum = A[i, j];
                    for (k = 0; k < j; k++) sum -= A[i, k] * A[k, j];
                    A[i, j] = sum;
                    if ((dum = vv[i] * Math.Abs(sum)) >= big)
                    {//新的候选主元是否比已选中的更合适？
                        big = dum;
                        imax = i;
                    }
                }
                if (j != imax)
                {//是否需要进行交换
                    A = NRMatrix.ExchangeRows(A, imax, j);
                    vv[imax] = vv[j];
                }
                index[j] = imax;
                if (A[j, j] == 0.0) A[j, j] = TINY;
                if (j != n - 1)
                {
                    dum = 1.0 / (A[j, j]);
                    for (i = j + 1; i < n; i++) A[i, j] *= dum;
                }
            }

            for (k = 0; k < b.Cols; k++)
            {
                for (i = 0; i < n; i++)
                {
                    ip = index[i];
                    sum = b[ip, k];
                    b[ip, k] = b[i, k];
                    if (ii != 0)
                        for (j = ii - 1; j < i; j++) sum -= A[i, j] * b[j, k];
                    else if (sum != 0.0) ii = i + 1;
                    b[i, k] = sum;
                }
                for (i = n - 1; i >= 0; i--)
                {
                    sum = b[i, k];
                    for (j = i + 1; j < n; j++) sum -= A[i, j] * b[j, k];
                    b[i, k] = sum / A[i, i];
                }
                ii = 0;
            }
            result = new NRMatrix(b);
            return true;
        }
        public static void SVdcmp(NRMatrix a, out NRVec W, out NRMatrix V)//SVD分解A=UWV',
        //矩阵a[m,n]返回U,向量W[n]返回奇异值，矩阵V返回V[n，n]
        {
            bool flag;
            int i = 0, its = 0, j = 0, jj = 0, k = 0, l = 0, nm = 0;
            double anorm, c, f, g, h, s, scale, x, y, z;

            int m = a.Rows;
            int n = a.Cols;
            NRVec w = new NRVec(a.Cols);
            NRMatrix v = new NRMatrix(a.Cols, a.Cols);
            NRVec rv1 = new NRVec(n);
            g = scale = anorm = 0.0;
            for (i = 0; i < n; i++)
            {
                l = i + 2;
                rv1[i] = scale * g;
                g = s = scale = 0.0;
                if (i < m)
                {
                    for (k = i; k < m; k++) scale += Math.Abs(a[k, i]);
                    if (scale != 0.0)
                    {
                        for (k = i; k < m; k++)
                        {
                            a[k, i] /= scale;
                            s += a[k, i] * a[k, i];
                        }
                        f = a[i, i];
                        g = -SIGN(Math.Sqrt(s), f);
                        h = f * g - s;
                        a[i, i] = f - g;
                        for (j = l - 1; j < n; j++)
                        {
                            for (s = 0.0, k = i; k < m; k++) s += a[k, i] * a[k, j];
                            f = s / h;
                            for (k = i; k < m; k++) a[k, j] += f * a[k, i];
                        }
                        for (k = i; k < m; k++) a[k, i] *= scale;
                    }
                }
                w[i] = scale * g;
                g = s = scale = 0.0;
                if (i + 1 <= m && i != n)
                {
                    for (k = l - 1; k < n; k++) scale += Math.Abs(a[i, k]);
                    if (scale != 0.0)
                    {
                        for (k = l - 1; k < n; k++)
                        {
                            a[i, k] /= scale;
                            s += a[i, k] * a[i, k];
                        }
                        f = a[i, l - 1];
                        g = -SIGN(Math.Sqrt(s), f);
                        h = f * g - s;
                        a[i, l - 1] = f - g;
                        for (k = l - 1; k < n; k++) rv1[k] = a[i, k] / h;
                        for (j = l - 1; j < m; j++)
                        {
                            for (s = 0.0, k = l - 1; k < n; k++) s += a[j, k] * a[i, k];
                            for (k = l - 1; k < n; k++) a[j, k] += s * rv1[k];
                        }
                        for (k = l - 1; k < n; k++) a[i, k] *= scale;
                    }
                }