ldlt_pro.c

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

/* This file is used to test the implementation of LDL' decomposition algorithm */
/* Revision history:                                                            */
/*    2005.06.29       Zhu Songsheng([e-mail])    initial begin   */
#include <stdio.h>
#include <stdlib.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");
}

void LDLTBKSB (unsigned int n, double * *a, double *b)
/* This function can solve the symmetric linear determinantal equation after    */
/* the LDT' decomposition(can use different b[] with the same a[][]).           */
/* note:                                                                        */
/*    1. The matrix a[][] must already be decomposed to LDL'. Can use the       */
/*    matrix calculated after LDLTDCMP().                                       */
/* Parameter usage:                                                             */
/*  n    : The dimension of the determinantal equation.                         */
/*  a[][]: Input the LDL' matrix of the determinantal equation.                 */
/*         The form is 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].       */
/*  b[]  : Input the b vector(the dimension of a and b must be accordant),      */
/*         and the result of x[] are also saved in b[] after calculated.        */
/* Revision history:                                                            */
/*    2005.06.29       Zhu Songsheng([e-mail])    initial begin   */
{
    int i;
    int k;
    for (i = 0; i < n; i++)
    {
        /* Calculate y[i], and saved in b[i] */
        for (k = 0; k < i; k++)
            b[i] = b[i] - a[i][k] * b[k];
    }
    for (i = n - 1; i >= 0; i--)
    {
        /* Calculate z[i], and saved in b[i] */
        b[i] = b[i] / a[i][i];
        /* Calculate x[i], and saved in b[i] */
        for (k = i + 1; k < n; k++)
            b[i] = b[i] - a[k][i] * b[k];
    }
    printf ("\n** The result of x[] are saved in b[] after calculated. **\n");
}

void main ()
{
    FILE *fp_ldl_in;    /* Read parameter file pointer. */
    FILE *fp_ldl_out;   /* Save answer file pointer. */
    unsigned int n;
    double * *a;
    double *b;
    double temp;        /* Only used to get real numbers from file. */
    int i, j;

    printf ("\nStart ...\n");
    /* Read parameter from file. */
    if ( (fp_ldl_in = fopen ("../data/LDLTPARA.dat", "r")) == NULL )
    {
        printf ("\n\nERROR: Failed to open ASCII file LDLTPARA.dat !!!\n");
        printf (" Please put the file in folder ../data/.\n");
        exit (1);
    }
    else
        printf ("Open ASCII file LDLTPARA.dat successfully!\n");

    printf ("Reading file LDLTPARA.dat ...\n");

    /* To ignore the line begin with '#' */
    while ((fgetc (fp_ldl_in)) == '#')
        while ((fgetc (fp_ldl_in)) != '\n');
    fseek (fp_ldl_in, -1L, 1);  /* Backword a char */
    /* End of ignore */

    fscanf (fp_ldl_in, "n =%u", &n);

    /* To ignore the line begin with '#' */
    while ((fgetc (fp_ldl_in)) != '\n');
    while ((fgetc (fp_ldl_in)) == '#')
        while ((fgetc (fp_ldl_in)) != '\n');
    fseek (fp_ldl_in, -1L, 1);  /* Backword a char */
    /* End of ignore */

    /* Apply for memory */
    if ((a = (double **) malloc (n * sizeof(double *))) == NULL)
    {
        printf ("\n\nERROR: Not Enough Memory!\n");
        printf (" Please reduce n and try again.\n");
        exit (1);
    }
    for (i = 0; i < n; i++)
    {
        if ((*(a + i) = (double *) malloc (n * sizeof(double))) == NULL)
        {
            printf ("\n\nERROR: Not Enough Memory!\n");
            printf (" Please reduce n and try again.\n");
            exit (1);
        }
    }
    if ((b = (double *) malloc (n * sizeof(double))) == NULL)
    {
        printf ("\nERROR: Not Enough Memory!\n");
        printf (" Please reduce n and try again.\n");
        exit (1);
    }
    /* End of apply */

    for (i = 0; i < n; i++)
    {
        for (j = 0; j < n; j++)
        {
            fscanf (fp_ldl_in, "%lf", &temp);
            a[i][j] = temp;
            if ((i > j) && (a[i][j] != a[j][i]))
            {
                printf ("\n\nERROR: The matrix A is not symmetric");
                printf ("(A not equal to A\').\n");
                printf ("Please change the data of A in the file");
                printf (" ../data/LDLTPARA.DAT\n");
                fclose (fp_ldl_in);
                exit (1);
            }
        }
        fscanf (fp_ldl_in, "%lf", &temp);
        b[i] = temp;

        /* To ignore the line begin with '#' */
        while ((fgetc (fp_ldl_in)) != '\n');
        while ((fgetc (fp_ldl_in)) == '#')
            while ((fgetc (fp_ldl_in)) != '\n');
        fseek (fp_ldl_in, -1L, 1);  /* Backword a char */
        /* End of ignore */
    }
    printf ("Read data sucessfully.\n\n");
    fclose (fp_ldl_in);

    /* Open file to save answers. */
    if ( (fp_ldl_out = fopen ("../Output/LDLT_KEY.dat", "w+")) == NULL )
    {
        printf ("\n\nERROR: Failed to open ASCII file LDLT_KEY.dat !!!\n");
        printf (" You must have the folder ../Output/\n");
        exit (1);
    }
    else
        printf ("Open ASCII file LDLT_KEY.dat successfully!\n");

    /* Show original parameters. */
    printf ("\nOriginal parameters:\n");
    fprintf(fp_ldl_out, "Original parameters:\n");
    printf ("n = %u\n", n);
    fprintf(fp_ldl_out, "n = %u\n", n);
    for (i = 0; i < n; i++)
    {
        for (j = 0; j < n; j++)
        {
            printf ("%8.3g  ", a[i][j]);
            fprintf(fp_ldl_out, "%8.3g  ", a[i][j]);
        }
        printf ("%8.3g\n", b[i]);
        fprintf(fp_ldl_out, "%8.3g\n", b[i]);
    }

    printf ("\n\nCalculate and save result to ASCII file LDLT_KEY.dat ...\n");

    /* Calculate answers. */
    LDLTDCMP (n, a);
    LDLTBKSB (n, a, b);

    /* Show and save answers. */
    printf ("\nAfter calculated:\n");
    fprintf(fp_ldl_out, "\nAfter calculated:\n");
    for (i = 0; i < n; i++)
    {
        for (j = 0; j < n; j++)
        {
            printf ("%8.3g  ", a[i][j]);
            fprintf(fp_ldl_out, "%8.3g  ", a[i][j]);
        }
        printf ("%8.3g\n", b[i]);
        fprintf(fp_ldl_out, "%8.3g\n", b[i]);
    }
    printf ("\nAnswers x[]:\n");
    fprintf(fp_ldl_out, "\nAnswers x[]:\n");
    for (i = 0; i < n; i++)
    {
        printf ("%8.3g", b[i]);
        fprintf(fp_ldl_out, "%8.3g", b[i]);
    }
    fclose (fp_ldl_out);
    printf ("\n\nAll results are saved.\n");
    printf ("Please open ASCII file ../Output/LDLT_KEY.dat to check it.\n");
    for (i = 0; i < n; i++)  
        free (*(a + i));  
    free (a);  
    free (b);  
}