LUP分解——矩阵运算

/**
 * n维矩阵LUP分解 PA=LU P置换矩阵 A待分解矩阵 L下三角矩阵 主对角线全为1 U上三角矩阵
 * @param A n维矩阵 A.size() == n
 * @param P 置换矩阵 P.size() == A.size() 置换矩阵以A的行索引表示
 * @return 能否分解
 * 时间复杂度O(n^3) 嵌套三层循环
 */
bool LUP_Decomposition(vector<vector<double> > &A, vector<int> &P) {
    int n = A.size();
    for (int i = 0; i < n; i++)
        P[i] = i;//A置换前每一行的位置
    double temp;
    int flag = 0;
    for (int k = 0; k < n; k++) {
        temp = 0;
        for (int i = k; i < n; i++) {//选定主元
            if (fabs(A[i][k]) > temp) {
                temp = fabs(A[i][k]);
                flag = i;
            }
        }
        if (temp == 0.) return false;//不存在分解
        if (k != flag) {//置换
            P[k] ^= P[flag] ^= P[k] ^= P[flag];//标记行置换
            for (int i = 0; i < n; i++) {//行置换
                temp = A[k][i];
                A[k][i] = A[flag][i];
                A[flag][i] = temp;
            }
        }
        for (int i = k + 1; i < n; i++) {//分解
            A[i][k] = A[i][k] / A[k][k];
            for (int j = k + 1; j < n; j++)
                A[i][j] = A[i][j] - A[i][k] * A[k][j];
        }
    }
    return true;
}

/**
 * LUP分解求n维线性方程组的解
 * @param A 系数矩阵的LUP分解矩阵
 * @param P 置换矩阵
 * @param B 常数列向量
 * @param X 解向量
 * 时间复杂度O(n^3)
 */
void SolvingLinearEquations_LUP(vector<vector<double> > A, vector<int> P, vector<double> B, vector<double> &X) {
    int n = A.size();
    for (int r = 0; r < n; r++) {
        X[r] = B[P[r]];
        for (int c = 0; c < r; c++)
            X[r] -= A[r][c] * X[c];
    }
    for (int r = n - 1; r >= 0; r--) {
        for (int c = n - 1; c > r; c--)
            X[r] -= A[r][c] * X[c];
        X[r] /= A[r][r];
    }
}

/**
 * LUP分解求n维线性方程组的解 AX = B
 * @param A 系数矩阵
 * @param B 常数列向量
 * @param X 解向量
 * @return 是否有解
 * 时间复杂度O(n^3)
 */
bool SolvingLinearEquations(vector<vector<double> > A, vector<double> B, vector<double> &X){
    int n = A.size();
    vector<int> P(n);
    bool flag = LUP_Decomposition(A, P);
    if(!flag) return flag;
    SolvingLinearEquations_LUP(A, P, B, X);
    return true;
}

/**
 * 求逆矩阵 AI = IA = E E为单位矩阵
 * @param A 原矩阵
 * @param I A的逆矩阵
 * @return 是否可逆
 * 时间复杂度O(n^3)
 */
bool InverseMatrix(vector<vector<double> > A, vector<vector<double> > &I) {
    int n = A.size();
    vector<double> i(n);//逆矩阵的某一列i
    vector<double> e(n, 0.);//单位矩阵中与i的对应列e(单位向量) 有Ai = e
    vector<int> P(n);//置换矩阵
    bool flag = LUP_Decomposition(A, P);//LUP分解
    if (!flag) return flag;
    for (int k = 0; k < n; k++) {//求I的每一列
        e[k] = 1.;//赋1表示单位向量e
        SolvingLinearEquations_LUP(A, P, e, i);//Ai=e
        for (int r = 0; r < n; r++)//i赋给I当前列
            I[r][k] = i[r];
        e[k] = 0;//重置为零向量
    }
    return true;
}

/**
 * 矩阵乘法 A * B = C
 * @param A m*o矩阵
 * @param B o*n矩阵
 * @return m*n矩阵
 */
vector<vector<double> > MatrixMultiply(vector<vector<double> >A, vector<vector<double> >B) {
    int m = A.size(), o = B.size(), n = B[0].size();
    vector<vector<double> > C(m, vector<double>(n));
    for(int i = 0; i < m; i++)
        for(int j = 0; j < n; j++){
            C[i][j] = 0;
            for(int k = 0; k < o; k++)
                C[i][j] += A[i][k]*B[k][j];
        }
    return C;
}

发表评论

电子邮件地址不会被公开。 必填项已用*标注