单晶塑性abaqus-umat 自学—续

本文是针对黄老师单晶umat程序中的一些细节问题。都是本人一下没看明白，通过思考和查资料弄明白后，写在这里。

1.刚度矩阵的旋转

C-----  Rotation matrix: ROTD to transform local elastic matrix DLOCAL 
C     to global elastic matrix D
C
      DO J=1,3
         J1=1+J/3
         J2=2+J/2

         DO I=1,3
            I1=1+I/3
            I2=2+I/2

            ROTD(I,J)=ROTATE(I,J)**2
            ROTD(I,J+3)=2.*ROTATE(I,J1)*ROTATE(I,J2)
            ROTD(I+3,J)=ROTATE(I1,J)*ROTATE(I2,J)
            ROTD(I+3,J+3)=ROTATE(I1,J1)*ROTATE(I2,J2)+
     2                    ROTATE(I1,J2)*ROTATE(I2,J1)

         END DO
      END DO

C-----  Elastic matrix in global system: D
C     {D} = {ROTD} * {DLOCAL} * {ROTD}transpose

旋转矩阵设为 $\mathbf{R}$ ，向量旋转：
$${\vec{u}}^{New}=\mathbf{R}\cdot {\vec{u}}^{Old}{\vec{u}}_i^{New}=R_{ij}{\vec{u}}_j^{Old}$$

二阶张量旋转：
$${\mathbf{\sigma }}^{New}=\mathbf{R}\cdot {\mathbf{\sigma }}^{Old}\cdot {\mathbf{R}}^T{\mathbf{\sigma }}_{ij}^{New}=R_{ik}{\mathbf{\sigma }}_{kl}^{Old}{R}_{jl}$$

改写为：$${\mathbf{\sigma }}^{New}=\mathbb{R}:\mathbf{\sigma }{\mathbf{\sigma }}_{ij}^{New}={\mathbb{R}}_{ijkl}{\mathbf{\sigma }}_{kl}^{Old}$$
其中：$${\mathbb{R}}_{ijkl}=R_{ik}{R}_{jl}$$

对于四阶刚度张量：
$${\mathbf{D}}^{New}=\mathbb{R}:{\mathbf{D}}^{Old}:{\mathbb{R}}^T{\mathbf{D}}_{ijkl}^{New}={\mathbb{R}}_{ijnm}{\mathbf{D}}_{nm\alpha \beta }^{Old}{\mathbb{R}}_{\alpha \beta kl}$$

程序中ROTATE为$R$，ROTD为 $\mathbb{R}$的voigt展开：

VUMAT的应力分量顺序为 { σ 11   σ 22   σ 33   σ 12   σ 23   σ 31 } \{\sigma_{11}\,\sigma_{22}\,\sigma_{33}\, \sigma_{12}\,\sigma_{23}\,\sigma_{31}\} {σ11​σ22​σ33​σ12​σ23​σ31​}
$${\mathbb{R}}_{ijkl}=\left[\begin{array}{ccccccccccc}& kl& 11& 22& 33& 12& 23& 31& 21& 32& 13\\ ij& & & & & & & & & & \\ 11& & R_{11}{R}_{11}& R_{12}{R}_{12}& R_{13}{R}_{13}& R_{11}{R}_{12}& R_{12}{R}_{13}& R_{13}{R}_{11}& R_{12}{R}_{11}& R_{13}{R}_{12}& R_{11}{R}_{13}\\ 22& & R_{21}{R}_{21}& R_{22}{R}_{22}& R_{23}{R}_{23}& R_{21}{R}_{22}& R_{22}{R}_{23}& R_{23}{R}_{21}& R_{22}{R}_{21}& R_{23}{R}_{22}& R_{21}{R}_{23}\\ 33& & R_{31}{R}_{31}& R_{32}{R}_{32}& R_{33}{R}_{33}& R_{31}{R}_{32}& R_{32}{R}_{33}& R_{33}{R}_{31}& R_{32}{R}_{31}& R_{33}{R}_{32}& R_{31}{R}_{33}\\ 12& & R_{11}{R}_{21}& R_{12}{R}_{22}& R_{13}{R}_{23}& R_{11}{R}_{22}& R_{12}{R}_{23}& R_{13}{R}_{21}& R_{12}{R}_{21}& R_{13}{R}_{21}& R_{11}{R}_{23}\\ 23& & R_{21}{R}_{31}& R_{22}{R}_{32}& R_{23}{R}_{33}& R_{21}{R}_{32}& R_{22}{R}_{33}& R_{23}{R}_{31}& R_{22}{R}_{31}& R_{23}{R}_{32}& R_{21}{R}_{33}\\ 31& & R_{31}{R}_{11}& R_{32}{R}_{12}& R_{33}{R}_{13}& R_{31}{R}_{12}& R_{32}{R}_{13}& R_{33}{R}_{11}& R_{32}{R}_{11}& R_{33}{R}_{12}& R_{31}{R}_{13}\\ 21& & R_{21}{R}_{11}& R_{22}{R}_{12}& R_{23}{R}_{13}& R_{21}{R}_{12}& R_{22}{R}_{13}& R_{23}{R}_{11}& R_{22}{R}_{11}& R_{23}{R}_{12}& R_{21}{R}_{13}\\ 32& & R_{31}{R}_{21}& R_{32}{R}_{22}& R_{33}{R}_{23}& R_{31}{R}_{22}& R_{32}{R}_{23}& R_{33}{R}_{21}& R_{32}{R}_{21}& R_{33}{R}_{22}& R_{31}{R}_{23}\\ 13& & R_{11}{R}_{31}& R_{12}{R}_{32}& R_{13}{R}_{33}& R_{11}{R}_{32}& R_{12}{R}_{33}& R_{13}{R}_{31}& R_{12}{R}_{31}& R_{13}{R}_{32}& R_{11}{R}_{33}\end{array}\right]$$

$${\mathbb{R}}_{ijkl}^{Voigt}=\left[\begin{array}{cccccccc}& kl& 11& 22& 33& 12& 23& 31\\ ij& & & & & & & \\ 11& & R_{11}{R}_{11}& R_{12}{R}_{12}& R_{13}{R}_{13}& 2R_{11}{R}_{12}& 2R_{12}{R}_{13}& 2R_{13}{R}_{11}\\ 22& & R_{21}{R}_{21}& R_{22}{R}_{22}& R_{23}{R}_{23}& 2R_{21}{R}_{22}& 2R_{22}{R}_{23}& 2R_{23}{R}_{21}\\ 33& & R_{31}{R}_{31}& R_{32}{R}_{32}& R_{33}{R}_{33}& 2R_{31}{R}_{32}& 2R_{32}{R}_{33}& 2R_{33}{R}_{31}\\ 12& & R_{11}{R}_{21}& R_{12}{R}_{22}& R_{13}{R}_{23}& R_{11}{R}_{22}+R_{12}{R}_{21}& R_{12}{R}_{23}+R_{13}{R}_{21}& R_{13}{R}_{21}+R_{11}{R}_{23}\\ 23& & R_{21}{R}_{31}& R_{22}{R}_{32}& R_{23}{R}_{33}& R_{21}{R}_{32}+R_{22}{R}_{31}& R_{22}{R}_{33}+R_{23}{R}_{32}& R_{23}{R}_{31}+R_{21}{R}_{33}\\ 31& & R_{31}{R}_{11}& R_{32}{R}_{12}& R_{33}{R}_{13}& R_{31}{R}_{12}+R_{32}{R}_{11}& R_{32}{R}_{13}+R_{33}{R}_{12}& R_{33}{R}_{11}+R_{31}{R}_{13}\end{array}\right]$$

UMAT的应力分量顺序为 { σ 11   σ 22   σ 33   σ 12   σ 13   σ 23 } \{\sigma_{11}\,\sigma_{22}\,\sigma_{33}\, \sigma_{12}\,\sigma_{13}\,\sigma_{23}\} {σ11​σ22​σ33​σ12​σ13​σ23​}

$${\mathbb{R}}_{ijkl}^{Voigt}=\left[\begin{array}{cccccccc}& kl& 11& 22& 33& 12& 13& 23\\ ij& & & & & & & \\ 11& & R_{11}{R}_{11}& R_{12}{R}_{12}& R_{13}{R}_{13}& 2R_{11}{R}_{12}& 2R_{11}{R}_{13}& 2R_{12}{R}_{13}\\ 22& & R_{21}{R}_{21}& R_{22}{R}_{22}& R_{23}{R}_{23}& 2R_{21}{R}_{22}& 2R_{21}{R}_{23}& 2R_{22}{R}_{23}\\ 33& & R_{31}{R}_{31}& R_{32}{R}_{32}& R_{33}{R}_{33}& 2R_{31}{R}_{32}& 2R_{31}{R}_{33}& 2R_{32}{R}_{33}\\ 12& & R_{11}{R}_{21}& R_{12}{R}_{22}& R_{13}{R}_{23}& R_{11}{R}_{22}+R_{12}{R}_{21}& R_{11}{R}_{23}+R_{13}{R}_{21}& R_{12}{R}_{23}+R_{13}{R}_{22}\\ 13& & R_{11}{R}_{31}& R_{12}{R}_{32}& R_{13}{R}_{33}& R_{11}{R}_{32}+R_{12}{R}_{31}& R_{11}{R}_{33}+R_{13}{R}_{31}& R_{12}{R}_{33}+R_{13}{R}_{32}\\ 23& & R_{21}{R}_{31}& R_{22}{R}_{32}& R_{23}{R}_{33}& R_{21}{R}_{32}+R_{22}{R}_{31}& R_{21}{R}_{33}+R_{23}{R}_{31}& R_{22}{R}_{33}+R_{23}{R}_{32}\end{array}\right]$$

2. LU分解

      SUBROUTINE LUDCMP (A, N, NP, INDX, D)

C-----  LU decomposition

C-----  Use single precision on cray
C
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (NMAX=200, TINY=1.0E-20)
      DIMENSION A(NP,NP), INDX(N), VV(NMAX)

      D=1.
      DO I=1,N
         AAMAX=0.

         DO J=1,N
            IF (ABS(A(I,J)).GT.AAMAX) AAMAX=ABS(A(I,J))
         END DO

         IF (AAMAX.EQ.0.) PAUSE 'Singular matrix.'
         VV(I)=1./AAMAX
      END DO

      DO J=1,N
         DO I=1,J-1
            SUM=A(I,J)

            DO K=1,I-1
               SUM=SUM-A(I,K)*A(K,J)
            END DO

            A(I,J)=SUM
         END DO
         AAMAX=0.

         DO I=J,N
            SUM=A(I,J)

            DO K=1,J-1
               SUM=SUM-A(I,K)*A(K,J)
            END DO

            A(I,J)=SUM
            DUM=VV(I)*ABS(SUM)
            IF (DUM.GE.AAMAX) THEN
               IMAX=I
               AAMAX=DUM
            END IF
         END DO

         IF (J.NE.IMAX) THEN
            DO K=1,N
               DUM=A(IMAX,K)
               A(IMAX,K)=A(J,K)
               A(J,K)=DUM
            END DO

            D=-D
            VV(IMAX)=VV(J)
         END IF

         INDX(J)=IMAX
         IF (A(J,J).EQ.0.) A(J,J)=TINY
         IF (J.NE.N) THEN
            DUM=1./A(J,J)
            DO I=J+1,N
               A(I,J)=A(I,J)*DUM
            END DO
         END IF

      END DO

      RETURN
      END

黄老师子程序中的LU分解，实际是改进后的LUP分解，LUP分解中每一步是找的是这一列的最大值并通过P矩阵交换的对角线上，黄老师程序并不是找最大值，而是相对最大值，当前值要除以此行的最大值，再比较大小DUM=VV(I)*ABS(SUM) IF (DUM.GE.AAMAX) THEN，找出相对最大值，且下三角矩阵 $L$ 存在变化后矩阵的下三角部分（下三角矩阵主对角线上的值均为1），上三角矩阵 $U$ 存在变化后矩阵的上三角部分。置换关系P，通过INDX(N)记录
以一个4×4的矩阵举例说明（数值任意取的）：
$$A=\left[\begin{array}{cccc}{A}_{11}& A_{12}& A_{13}& A_{14}\\ A_{21}& A_{22}& A_{23}& A_{24}\\ A_{31}& A_{32}& A_{33}& A_{34}\\ A_{41}& A_{42}& A_{43}& A_{44}\end{array}\right]=\left[\begin{array}{cccc}1& 2& 3& 4\\ 2& 3& 4& 1\\ 3& 3& 6& 7\\ 5& 3& 2& 1\end{array}\right]$$
首先，计算每一行的最大值，并求倒数记为 vv(I)，具体值vv(1)=1/4,vv(2)=1/4,vv(3)=1/7,vv(1)=1/5
其次，列循环：
第一列，当J=1时 I=1,J-1不存在，I=J,N不存在，
找出一列的相对最大值，IMAX=4
$\mathrm{J}\ne \mathrm{I}\mathrm{M}\mathrm{A}\mathrm{X}$（J.NE.IMAX） 将第IMAX行和第J行互换，交换后第IMAX行的最大值倒数为 vv(IMAX)=vv(J)，于此同时INDX(1)=4，即第1行和第4行互换
$$A=P_1\cdot A=\left[\begin{array}{cccc}0& 0& 0& 1\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 1& 0& 0& 0\end{array}\right]\left[\begin{array}{cccc}1& 2& 3& 4\\ 2& 3& 4& 1\\ 3& 3& 6& 7\\ 5& 3& 2& 1\end{array}\right]=\left[\begin{array}{cccc}5& 3& 2& 1\\ 2& 3& 4& 1\\ 3& 3& 6& 7\\ 1& 2& 3& 4\end{array}\right]$$
求第1列下对角阵
$$A=L_1\cdot A=\left[\begin{array}{cccc}5& 3& 2& 1\\ 2& 3& 4& 1\\ 3& 3& 6& 7\\ 1& 2& 3& 4\end{array}\right]=\left[\begin{array}{cccc}1& 0& 0& 0\\ 2/5& 1& 0& 0\\ 3/5& 0& 1& 0\\ 1/5& 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}5& 3& 2& 1\\ 0& 3-3(2/5)& 4-2(2/5)& 1-(2/5)\\ 0& 3-3(3/5)& 6-2(3/5)& 7-(3/5)\\ 0& 2-3(1/5)& 3-2(1/5)& 4-(1/5)\end{array}\right]$$
第1列循环完
$$A=\left[\begin{array}{cccc}5& 3& 2& 1\\ 2/5& 3& 4& 1\\ 3/5& 5& 6& 7\\ 1/5& 2& 3& 4\end{array}\right]$$
循环继续，第2列J=2
I=1,J-1$A_{12}=A_{12}$
I=J,N
DO K=1,J-1 SUM=SUM-A(I,K)*A(K,J) END DO
$$A=\left[\begin{array}{cccc}5& 3& 2& 1\\ 2/5& 3-3(2/5)& 4& 1\\ 3/5& 5-3(3/5)& 6& 7\\ 1/5& 2-3(1/5)& 3& 4\end{array}\right]=\left[\begin{array}{cccc}5& 3& 2& 1\\ 2/5& 9/5& 4& 1\\ 3/5& 6/5& 6& 7\\ 1/5& 7/5& 3& 4\end{array}\right]$$
IMAX=2 \quad IMAX=J \quad 于此同时INDX(2)=2 \quad 求下三角：
$$A=\left[\begin{array}{cccc}5& 3& 2& 1\\ 2/5& 3-3(2/5)& 4& 1\\ 3/5& 5-3(3/5)& 6& 7\\ 1/5& 2-3(1/5)& 3& 4\end{array}\right]=\left[\begin{array}{cccc}5& 3& 2& 1\\ 2/5& 9/5& 4& 1\\ 3/5& 6/9& 6& 7\\ 1/5& 7/9& 3& 4\end{array}\right]$$
循环继续，第3列J=3
I=1,J-1
DO K=1,I-1 SUM=SUM-A(I,K)*A(K,J) END DO$A_{23}=A_{23}-A_{21}{A}_{13}$
$$A=\left[\begin{array}{cccc}5& 3& 2& 1\\ 2/5& 9/5& 4-2(2/5)& 1\\ 3/5& 6/9& 6& 7\\ 1/5& 7/9& 3& 4\end{array}\right]$$

I=J,N
DO K=1,J-1 SUM=SUM-A(I,K)*A(K,J) END DO
$$A=\left[\begin{array}{cccc}5& 3& 2& 1\\ 2/5& 9/5& 16/5& 1\\ 3/5& 6/9& 6-2(3/5)-(16/5)(6/9)& 7\\ 1/5& 7/9& 3-2(1/5)-(16/5)(7/9)& 4\end{array}\right]=\left[\begin{array}{cccc}5& 3& 2& 1\\ 0.4& 1.8& 3.2& 1\\ 0.6& 6/9& 2.667& 7\\ 0.2& 7/9& 0.1111& 4\end{array}\right]$$
IMAX=3 \quad IMAX=J \quad \quad 于此同时INDX(3)=3 \quad 求下三角：
$$A=\left[\begin{array}{cccc}5& 3& 2& 1\\ 0.4& 1.8& 3.2& 1\\ 0.6& 6/9& 2.667& 7\\ 0.2& 7/9& 1/24& 4\end{array}\right]$$

循环继续，第4列J=4
I=1,J-1
DO K=1,I-1 SUM=SUM-A(I,K)*A(K,J) END DO
$$A=\left[\begin{array}{cccc}5& 3& 2& 1\\ 0.4& 1.8& 3.2& 1-0.4\\ 0.6& 6/9& 2.667& 7-0.6-(6/9)(1-0.4)\\ 0.2& 7/9& 1/24& 4\end{array}\right]=\left[\begin{array}{cccc}5& 3& 2& 1\\ 0.4& 1.8& 3.2& 0.6\\ 0.6& 6/9& 2.667& 6\\ 0.2& 7/9& 1/24& 4\end{array}\right]$$
I=J,N
DO K=1,J-1 SUM=SUM-A(I,K)*A(K,J) END DO
$$A=\left[\begin{array}{cccc}5& 3& 2& 1\\ 0.4& 1.8& 3.2& 0.6\\ 0.6& 6/9& 2.667& 6\\ 0.2& 7/9& 1/24& 4-0.2-0.6(7/9)-6(1/24)\end{array}\right]=\left[\begin{array}{cccc}5& 3& 2& 1\\ 0.4& 1.8& 3.2& 0.6\\ 0.6& 6/9& 2.667& 6\\ 0.2& 7/9& 1/24& 3.083\end{array}\right]$$
IMAX=4 \quad IMAX=J \quad \quad 于此同时INDX(4)=4 \quad 且J=N \quad 不用求下三角：
最终：$$A=\left[\begin{array}{cccc}5& 3& 2& 1\\ 0.4& 1.8& 3.2& 0.6\\ 0.6& 6/9& 2.667& 6\\ 0.2& 7/9& 1/24& 3.083\end{array}\right]$$

3. 反向代入法

      SUBROUTINE LUBKSB (A, N, NP, INDX, B)

C-----  Linear equation solver based on LU decomposition

C-----  Use single precision on cray
C
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(NP,NP), INDX(N), B(N)

      II=0
      DO I=1,N
         LL=INDX(I)
         SUM=B(LL)
         B(LL)=B(I)

         IF (II.NE.0) THEN
            DO J=II,I-1
               SUM=SUM-A(I,J)*B(J)
            END DO
         ELSE IF (SUM.NE.0.) THEN
            II=I
         END IF

         B(I)=SUM
      END DO

      DO I=N,1,-1
         SUM=B(I)

         IF (I.LT.N) THEN
            DO J=I+1,N
               SUM=SUM-A(I,J)*B(J)
            END DO
         END IF

         B(I)=SUM/A(I,I)
      END DO

      RETURN
      END

继续用上节4×4矩阵演示：
求解方程组：
$$Ax=y$$

$$\left[\begin{array}{cccc}1& 2& 3& 4\\ 2& 3& 4& 1\\ 3& 3& 6& 7\\ 5& 3& 2& 1\end{array}\right]\left\{\begin{array}{c}{x}_1\\ x_2\\ x_3\\ x_4\end{array}\right\}=\left\{\begin{array}{c}{y}_1\\ y_2\\ y_3\\ y_4\end{array}\right\}$$

LUP分解后的矩阵$A^{LUP}$
$$A^{LUP}=\left[\begin{array}{cccc}5& 3& 2& 1\\ 0.4& 1.8& 3.2& 0.6\\ 0.6& 6/9& 2.667& 6\\ 0.2& 7/9& 1/24& 3.083\end{array}\right]$$

置换关系INDX=4,2,3,4
$y$已知求解$x$

进入程序，现根据置换关系和下三角矩阵调整$y$:
I=1$y_1$ 和 $y_4$互换 且II=1$\left\{\begin{array}{cccc}{y}_4& y_2& y_3& y_1\end{array}\right\}$
I=2利用下三角矩阵调整 $y(2)$ 的值$\left\{\begin{array}{cccc}{y}_4& y_2-0.4\ast y(1)& y_3& y_1\end{array}\right\}$
即：$\left\{\begin{array}{cccc}{y}_4& y_2-0.4\ast y_4& y_3& y_1\end{array}\right\}$
I=3利用下三角矩阵调整 $y(2)$ 的值$\left\{\begin{array}{cccc}y(1)& y(2)& y(3)-0.4y(1)-0.6y(2)& y(4)\end{array}\right\}$
I=4利用下三角矩阵调整 $y(2)$ 的值$\left\{\begin{array}{cccc}y(1)& y(2)& y(3)& y(4)-0.4y(1)-0.6y(2)-0.2y(3)\end{array}\right\}$
此时已经转换成求解：
$$A^{U}x=\left[\begin{array}{cccc}5& 3& 2& 1\\ 0& 1.8& 3.2& 0.6\\ 0& 0& 2.667& 6\\ 0& 0& 0& 3.083\end{array}\right]\left\{\begin{array}{c}{x}_1\\ x_2\\ x_3\\ x_4\end{array}\right\}=\left\{\begin{array}{c}y(1)\\ y(2)\\ y(3)\\ y(4)\end{array}\right\}$$
反向代入：
$x(4)=[y(4)]/A(4,4)$
$x(3)=[y(3)-A(3,4)\ast x(4)]/A(3,3)$
$x(2)=[y(2)-A(2,4)\ast x(4)-A(2,3)\ast x(3)]/A(2,2)$
$x(1)=[y(1)-A(1,4)\ast x(4)-A(1,3)\ast x(3)-A(1,2)\ast x(2)]/A(1,1)$

4.利用LUP分解和反向代入法求逆矩阵

如果 $y$ 的初始值为$y_1={\left\{\begin{array}{cccc}1& 0& 0& 0\end{array}\right\}}^T$…$y_4={\left\{\begin{array}{cccc}0& 0& 0& 1\end{array}\right\}}^T$，求出的$\left[x_1{x}_2{x}_3{x}_4\right]$ 为$A$的逆矩阵。
$$L_1{L}_2...L_{N-1}{U}_{N\times N}=P_{N-1}...P_2{P}_1{A}_{N\times N}$$

求出 $U$之后，通过$U$的主对角线乘积算出 $\det (A)=(-1{)}^k\det (U)$（ k为行交换的次数），判断逆矩阵是否存在。
由于$P_i^{-1}=P_i$ 且$L_i^{-1}=2I-L_i$ ，因此：
$$U^{-1}{L}_{N-1}^{-1}...L_2^{-1}{L}_1^{-1}=A_{N\times N}^{-1}{P}_1{P}_2...P_{N-1}$$

$$A_{N\times N}^{-1}=U^{-1}{L}_{N-1}^{-1}...L_2^{-1}{L}_1^{-1}{P}_{N-1}...P_2{P}_1$$
$U^{-1}$通过反向代入法求出。

程序：

      Program Main
      Implicit none
      Integer,parameter:: N=4
      Integer I,J
      Real A(N,N),B(N,N),C(N,N),Det
      DATA A/1,2,3,4,34,15,2,1,3,22,14,7,7,81,2,13/ !任意输入矩阵
      Do I=1,N   !备份一个相同的矩阵，方便后续验证
          B(I,:)=A(I,:)
          write(*,*) B(I,:)
      End Do
      write(*,*) "*****"
      Call matInv(A,N,Det)  !求逆
      write(*,*) "Det=",Det
      write(*,*) "*****"
      Do I=1,N
          write(*,*) A(I,:)
      End Do
      write(*,*) "*****"
      C=matmul(A,B)  !验证
      Do I=1,N
          write(*,*) C(I,:)
      End Do
      write(*,*) "*****"
      End Program Main
      
C *********************************

      Subroutine matInv(A,N,Det)
      Integer I,J,K,P(N)
      Real A(N,N),B(N,N),U(N,N),L(N,N),Tempt(N),Det
     
      
      Call LUP(A,N,P,Det)  !LUP分解

      Do I=1,N
          Det=Det*A(I,I)
      End Do

      Call Back(A,N)  !求上对角矩阵的逆
      Do I=1,N
          Do J=1,N
              If (I .GT. J) Then
                  U(I,J)=0.0
              Else
                  U(I,J)=A(I,J)
              End If
          End Do
      End Do
    
      Do J=N-1,1,-1
                    
          Do K=1,N
              Do I=1,N
                  If (I .EQ. K) Then
                      L(K,K)=1.0
                  Else
                      L(I,K)=0.0
                  End If
              End Do
          End Do
          Do K=J+1,N 
              L(K,J)=-A(K,J)
          End Do
          
          U=matmul(U,L)
      End Do
      Do J=N-1,1,-1
          If(P(J).NE.J) Then
              Do K=1,N
                  Tempt(K)=U(K,J)
                  U(K,J)=U(K,P(J))
                  U(K,P(J))=Tempt(K)
              End Do          
          End If
      End Do
      Do I=1,N
      Do J=1,N
          A(I,J)=U(I,J)
      End Do
      End Do
      
      End Subroutine matInv
      
C *********************************
      Subroutine LUP(A,N,P,Det)
      Integer I,J,K,N,P(N)
      Real A(N,N),Temp(N),Det,Jmax
      
      Det=1.0
      Do J=1,N
          P(J)=J
          Jmax=Abs(A(J,J))  !先假设每一行的最大值是其对角线值，这样就不用了换行了
          Do I=J,N
              If (Abs(A(I,J)).GT.Jmax) Then  !寻找此列最大值
                  Jmax=Abs(A(I,J))
                  P(J)=I
              End If
          End Do
          If (J .NE. P(J)) Then !换行
              Det=-Det
              Do K=1,N
                  Temp(K)=A(J,K)
                  A(J,K)=A(P(J),K)
                  A(P(J),K)=Temp(K)
              End Do
          End If
          
          Do I=J+1,N  !计算下三角矩阵，并更新相应的上三角矩阵
              A(I,J)=A(I,J)/A(J,J)
              Do K=J+1,N
                  A(I,K)=A(I,K)-A(I,J)*A(J,K)
              End Do
          End Do
      End Do
      Return
      End Subroutine LUP
      
C *********************************
      Subroutine Back(A,N)
      Integer I,J,K,N
      Real A(N,N),Temp
      
      Do J=N,1,-1
          Do I=J,1,-1
              If (I.EQ.J) Then
                  A(I,J)=1.0/A(I,J)
              Else
                  Temp=0.0
                  Do K=J,I+1,-1
                      Temp=Temp+A(K,J)*A(I,K)
                  End Do
                  A(I,J)=-Temp/A(I,I)
              End If          
          End Do
      End Do
     
      Return
      End Subroutine Back

5. 行列式知识点补充

5.1 初等变换

定义：如果B可以由A经过一系列初等变换得到，则称矩阵A与B称等价。
初等行（列）变换：
定义：所谓数域P上矩阵的初等行变换是指下列3种变换：
1）以P中一个非零的数乘矩阵的某一行
2）把矩阵的某一行的c倍加到另一行，这里c是P中的任意一个数
3）互换矩阵中两行的位置
一般来说，一个矩阵经过初等行变换后就变成了另一个矩阵，当矩阵A经过初等行变换变成矩阵B时，一般写作
可以证明：任意一个矩阵经过一系列初等行变换总能变成阶梯型矩阵。

5.2 行列式初等变换的相关性质

性质1：行列互换，行列式不变
性质2：一数乘行列式的一行就相当于这个数乘此行列式
性质3：如果行列式中有两行相同，那么行列式为0，所谓两行相同，即两行对应的元素都相等
性质4：如果行列式中，两行成比例，那么该行列式为0
性质5：把一行的倍数加到另一行，行列式不变
性质6：对换行列式中两行的位置，行列式反号

5.3 矩阵乘积的行列式

引理：一个n阶矩阵A总可以通过第二种行和列的初等变换化为一个对角矩阵
$$A=\left[\begin{array}{cccc}{A}_{11}& A_{12}& ...& A_{1N}\\ A_{21}& A_{22}& ...& A_{2N}\\ ...& ...& ...& ...\\ A_{N1}& A_{N2}& ..& A_{NN}\end{array}\right]\to \bar{A}\left[\begin{array}{cccc}{\bar{A}}_{11}& 0& ...& 0\\ 0& {\bar{A}}_{22}& ...& 0\\ ...& ...& ...& ...\\ 0& 0& ..& {\bar{A}}_{NN}\end{array}\right]$$
第二种行和列的初等变换化，不改变行列式大小，因此$\det (A)=\det (\bar{A})$

定理：设A，B是任意两个n阶矩阵，那么$\det (AB)=\det (A)\cdot \det (B)$

由此可知：$\det (A)=(-1{)}^k\det (U)$ k为行交换的次数。