# Нижне-верхняя (LU) декомпозиция матрицы (алгоритм Краута)

Первым из алгоритмов, посвященным большому разделу решения систем линейных уравнений, представляем алгоритм Краута. Это фактически метод решения систем общего вида, конкурирующий по быстродействию с общеизвестным методом Гаусса-Жордана, но позволяющий более эффективно использовать решение.

Если мы можем разложить матрицу линейной системы A в произведение A=L*U, где L - нижняя, а U - верхняя треугольные матрицы, то решение системы уравнений с произвольной правой частью производится весьма просто, применением двух обратных подстановок. Более того, в отличие от известного метода Гаусса-Жордана, разложенная матрица позволяет быстро решать серии линейных уравнений с различными правыми частями при одной и той же матрице.

Метод Краута позволяет провести LU-декомпозицию матрицы примерно за то же число операций, что и "прямая" часть метода Гаусса-Жордана. Итоговые коэффициенты двух треугольных матриц упаковываются в матрицу того же размера, что и A, и на том же месте в памяти; при этом верхняя матрица U размещается в наддиагональной части и на диагонали; нижняя L в поддиагональной части, а диагональные элементы L считаются все равными 1 (без ограничения общности) и не выводятся.

Алгоритм Краута представляет из себя следующее:

1. Положить Lii=1 i=0,...,N-1 (здесь ничего не производится, а только имеется в виду для дальнейших шагов;
2. Для каждого j=0,1,...,N-1 проделать:
   1. Для всех i=0,...,j вычислить Uij:
      Uij=Aij - SUM(0<=k<i)(Lik*Uik)
      (при i=0 сумму полагать нулем);
   2. Для всех i=j+1,...,N-1 вычислить:
      Lij = Aij - SUM(0<=k<j)(Lik*Uik))/Uii
      Обе процедуры выполняются до перехода к новому j.

Устойчивая работа алгоритма Краута возможно только при условии выбора ведущего элемента для каждого обращения к j из п.2 алгоритма. Выбор производится перед выполнением п.2 для очередного j перестановки необработанных строк матрицы так, чтобы строка, подставляемая на место j, содержала наибольший по модулю элемент в колонке j. Порядок перестановки необходимо запомнить, для чего достаточно использовать дополнительный вектор целых величин. Еще один вектор используется внутри алгоритма для масштабирования.

Алгоритм Краута имеет несколько приложений, одно из которых - решение системы линейных уравнений (обратная подстановка с учетом порядка перестановки строк), другие - вычисление обратной матрицы и вычисление детерминанта. Подробнее вызовы функций, связанных с алгоритмом Краута и реализация деталей алгоритма находятся в комментариях к программе и в ней самой.

Программа самой декомпозиции и некоторых ее приложений приводится ниже

    /*
      LU-decomposition according to Crout's algorithm with pivoting.
            Description:
       int LU_decompos(double **a,int n,int *indx,int *d,double *vv);
            Parameters:
       a - source matrix (n x n) on input, destination on output;
       n - the matrix size;
       indx - integer array (size n) to remember permutations;
       d - on output, contains +1 or -1 for even or odd permutations number.
       vv - temporary array (size n).
               Returns:
       0 - the source matrix is singular (invalid for decomposition),
       1 - if OK.

      Back substitution, using LU decomposed matrix.
            Description:
      void LU_backsub(double **a,int n,int *indx,double *b);
            Parameters:
      a - the matrix decomposed by Crout;
      n - the matrix size;
      indx - permutation order obtained by decomposition algorithm;
      b - the vector (size n) to be substituted on input, the result
          of the substitution on output.
      Note: a and indx are not modified by this routine and could be
      used in multiple calls.

      Invertation of matrix, using LU decomposed matrix.
            Description:
      void LU_invert(double **a,int n,int *indx,double **inv,double *col);
            Parameters:
      a - the matrix decomposed by Crout;
      n - the matrix size;
      indx - permutation order obtained by decomposition algorithm;
      inv - the destination matrix;
      col - temporary array (size n).
      Note: test for singularity has been already obtained on the
      matrix decomposition, a and indx are not modified by this routine,
      the routine uses multiple backsubstitutions (previous algorithm).

      Obtaining the matrix determinant, using LU-decomposed matrix
            Description:
      double LU_determ(double **a,int n,int *indx,int *d);
            Parameters:
      a - the matrix decomposed by Crout;
      n - the matrix size;
      indx - permutation order obtained by decomposition algorithm;
      d - the parity sign (+1 or -1) obtained at decomposition.
            Returns:
      the determinant value. Note: non-zero (the matrix cannot be
      singular, if decomposed properly); a, indx and d are not modified
      by this routine.

    */

    /* for fabs(); inclusion could be removed if fabs() is implemented inline */
    #include <math.h>
    /* "small number" to avoid overflow in some cases */
    #define TINY 1.e-30

    /* the decomposition itself */
    int LU_decompos(double **a, int n, int *indx, int *d, double *vv) {
      register int i,imax,j,k;
      double big,sum,temp;
      /* set 1 for initially even (0) transmutations number */
      *d=1;
      /* search for the largest element in each row; save the scaling in the
         temporary array vv and return zero if the matrix is singular */
      for(i=0;i<n;i++) {
        big=0.;
        for(j=0;j<n;j++) if((temp=fabs(a[i][j]))>big) big=temp;
        if(big==0.) return(0);
        vv[i]=big;
      }
      /* the main loop for the Crout's algorithm */
      for(j=0;j<n;j++) {
        /* this is the part a) of the algorithm except for i==j */
        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;
        }
        /* initialize for the search for the largest pivot element */
        big=0.; imax=j;
        /* this is the part a) for i==j and part b) for i>j
           plus pivot search */
        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;
          /* is the figure of merit for the pivot better than
             the best so far? */
          if((temp=vv[i]*fabs(sum))>=big) {big=temp;imax=i;}
        }
        /* interchange rows, if needed, change parity and the scale factor */
        if(imax!=j) {
          for(k=0;k<n;k++) {
            temp=a[imax][k]; a[imax][k]=a[j][k]; a[j][k]=temp;
          }
          *d=-(*d); vv[imax]=vv[j];
        }
        /* store the index */
        indx[j]=imax;
        /* if the pivot element is zero, the matrix is singular but for some
           applications a tiny number is desirable instead */
        if(a[j][j]==0.) a[j][j]=TINY;
        /* finally, divide by the pivot element */
        if(j<n-1) {
          temp=1./a[j][j];
          for(i=j+1;i<n;i++) a[i][j]*=temp;
        }
      }
      return(1);
    }

    /* the back substitution */
    void LU_backsub(double **a,int n,int *indx,double *b) {
      register int i,j,ip,ii=-1;
      double sum;
      /* First step of backsubstitution; the only wrinkle is to unscramble
         the permutation order. Note: the algorithm is optimized for a
         possibility of large amount of zeroes in b */
      for(i=0;i<n;i++) {
        ip=indx[i];
        sum=b[ip];b[ip]=b[i];
        if(ii>=0) for(j=ii;j<i;j++) sum-=a[i][j]*b[j];
        else if(sum) ii=i;  /* a nonzero element encounted */
        b[i]=sum;
      }
      /* the second step */
      for(i=n-1;i>=0;i--) {
        sum=b[i];
        for(j=i+1;j<n;j++) sum-=a[i][j]*b[j];
        b[i]=sum/a[i][i];
      }
    }

    /* the matrix invertation */
    void LU_invert(double **a,int n,int *indx,double **inv,double *col) {
      register int i,j;
      for(j=0;j<n;j++) {
        for(i=0;i<n;i++) col[i]=0.;
        col[j]=1.;
        LU_backsub(a,n,indx,col);
        for(i=0;i<n;i++) inv[i][j]=col[i];
      }
    }

    /* calculating determinant */
    double LU_determ(double **a,int n,int *indx,int *d) {
      register int j;
      double res=(double)(*d);
      for(j=0;j<n;j++) res*=a[j][j];
      return(res);
    }

Здесь Вы можете скачать текст программы crout.c

Линейные уравнения
Индекс