Задачи по C/С++ - [2295] :: Прикладное программирование

#include <stdio.h>
#include <alloc.h>
#include <math.h>

int cholesky_decomp(float *mat_a[], float *mat_b[], int sz)
{
    float sum;
    int ii, jj, kk;

    for(ii = 0;ii < sz;ii++)
    {
        for(jj = ii;jj < sz;jj++)
        {
            sum = mat_a[ii][jj];

            for(kk = ii - 1;kk >= 0;kk--)
                sum = sum - mat_a[ii][kk]*mat_a[jj][kk];

            if(ii == jj)
            {
                if(sum <= .0F)
                {
                    printf("Cholesky's Method failed!\n");

                    return 0;
                }

                mat_b[ii][ii] = sqrt(sum);
            }
            else
                mat_b[jj][ii] = mat_a[jj][ii] = sum/mat_b[ii][ii];
        }
    }

return 1;
}

void cholesky_solver(float *mat_b[], float right[], float xx[], int sz)
{
    float *p_temp,
     sum;
    int ii, jj;

    p_temp = (float*)malloc(sz*sizeof(right[0]));

    for(ii = 0; ii < 4; ii++)
    {
        sum = 0.0f;

        for(jj = 0; jj < ii; jj++)
            sum += mat_b[ii][jj]*p_temp[jj];

        p_temp[ii] = (right[ii] - sum)/mat_b[ii][ii];
    }

    for(ii = 3; ii >= 0; ii--)
    {
        sum = 0.0f;

        for(jj = 3; jj > ii; jj--)
            sum += mat_b[jj][ii]*xx[jj];

        xx[ii] = (p_temp[ii] - sum)/mat_b[ii][ii];
    }

free(p_temp);
}

float left[4][4] = {{1.0f, 1.0f, 1.0f, 1.0f}, // This matrix generate by pascal(4) from Matlab
                         {1.0f, 2.0f, 3.0f, 4.0f},
                         {1.0f, 3.0f, 6.0f, 10.0f},
                         {1.0f, 4.0f, 10.0f, 20.0f}};
float bb[4][4];

float *p_left[4] = {(float*)&left[0], (float*)&left[1],
                         (float*)&left[2], (float*)&left[3]};

float *p_bb[4] = {(float*)&bb[0], (float*)&bb[1],
                        (float*)&bb[2], (float*)&bb[3]};

float right[4] = {6.0f, 16.0f, 33.0f, 59.0f};

float xx[4];

int main()
{
    int ii;

    if(!cholesky_decomp(p_left, p_bb, 4))
        return 0;

    cholesky_solver(p_bb, right, xx, 4);

    printf("Cholesky's Method\n\n");

    for(ii = 0; ii < 4; ii++)
     printf(" X%d = %f\n", ii, xx[ii]);

    return 0;
}

Модерирует : ShIvADeSt
Версия для печати • Подписаться • Добавить в закладки
На первую страницу • к этому сообщению • к последнему сообщению