ldltdcmp.c

LDL 分解源码。从文件输入数据

#include <stdio.h>

void LDLTDCMP (unsigned int n, double * *a)
/* 刘钦圣等, 数值计算方法教程, 冶金工业出版社, 北京, 2002, page 69.             */
/* LDL' 分解法算法, 可用与求解对称线性方程组. A 的各阶顺序主子式不为零.         */
/* This function implement the LDL' decomposition algorithm.                    */
/* The LDL' decomposition algorithm is used to solve the symmetric linear       */
/* determinantal equation.                                                      */
/* note:                                                                        */
/*    1. The matrix a[][] must be symmetric, and every rank of principal minor  */
/*    determinant of a[][] cannot be zero.                                      */
/*    2. This function must include file "stdio.h".                             */
/* Parameter usage:                                                             */
/*  n    : The dimension of the determinantal equation.                         */
/*  a[][]: Input the coefficient matrix of the determinantal equation.          */
/*         Every rank of principal minor determinant of a[][] cannot be zero.   */
/*         When calculate over, L is saved in the lower triangle of a[][](that  */
/*         is l[i][k] are in a[i][k]), D is saved in the a[i][i](that is d[i]   */
/*         are in a[i][i]), and temporary varible u[i][k] are in a[k][i].       */
/* Revision history:                                                            */
/*    2005.06.29       Zhu Songsheng([e-mail])    initial begin   */
{
    int k;
    int m;
    int i;
    for (k = 0; k < n; k++)
    {
        /* Calculate d[k], and saved in a[k][k] */
        for (m = 0; m < k; m++)
            a[k][k] = a[k][k] - a[m][k] * a[k][m];
        if (a[k][k] == 0)
        {
            printf ("\n\nERROR: LDL\' decompose failed !!\n");
            printf (" Every rank of principal minor determinant of A");
            printf (" cannot be zero.\n\n");
            return;
        }
        for (i = k + 1; i < n; i++)
        {
            /* temporary varible u[i][k] is saved in a[k][i]*/
            for (m = 0; m < k; m++)
                a[k][i] = a[k][i] - a[m][i] * a[k][m];
            /* Calculate l[i][k], and saved in a[i][k]*/
            a[i][k] = a[k][i] / a[k][k];
        }
    }
    printf ("\n\n");
    printf ("** When calculate over, L is saved in the lower triangle of **\n");
    printf ("** a[][](that is l[i][k] are in a[i][k]), D is saved in the **\n");
    printf ("** a[i][i](that is d[i] are in a[i][i]), and temporary      **\n");
    printf ("** varible u[i][k] are in upper triangle of a[k][i].        **\n");
}