0716 LU decomposition

[markoo (上行下笑阿..)]

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#define Num 5
void mprint(double A[Num][Num], double B[Num]);
time_t t1;
int main()
{
        int i,j,k;
        int o_row,n_row;
        double A[Num][Num];
        double A2[Num][Num];
        double B[Num];
        double B2[Num];
        double LU[Num][Num];
        double temp;
        double B1st[Num];
        double ans[Num];//final answer
        srand(time(&t1));
        for (i=0;i<Num;i++)
        {
                B[i] = rand()/32767.0;
                B2[i]= B[i];
                for (j=0;j<Num;j++)
                {
                        A[i][j] = rand()/32767.0;
                        A2[i][j] = A[i][j];
                }
        }
        printf("The matrix is now as following:\n");
        mprint(A,B);
        for (i=0;i<Num;i++)     LU[i][0] = A[i][0];
        for (i=1;i<Num;i++)     LU[0][i] = A[0][i]/LU[0][0];
        for (i=1;i<Num;i++)
        {
                for (j=i;j<Num;j++)
                {
                        LU[j][i] = A[j][i];
                        for (k=i-1;k>=0;k--)
                        {
                                LU[j][i] -= LU[j][k]*LU[k][i];
                        }
                        if (j<Num-1)
                        {
                                LU[i][j+1] = A[i][j+1];
                                for (k=i-1;k>=0;k--)
                                {
                                                LU[i][j+1] -= LU[i][k]*LU[k][j+1];
                                }
                                LU[i][j+1] = LU[i][j+1]/LU[i][i];
                        }
                }

        }
        printf("The matrix is now as following:\n");
        mprint(LU,B);
        // LB' = B, find B'
        B1st[0] = B[0]/LU[0][0];
        for (i=1;i<Num;i++)
        {
                B1st[i]=B[i];
                for (j=0;j<i;j++)       B1st[i]-=LU[i][j]*B1st[j];
                B1st[i] = B1st[i]/LU[i][i];
        }
        // Ux = B', find x
        ans[Num-1] = B1st[Num-1];
        for (i=Num-2;i>=0;i--)
        {
                ans[i]=B1st[i];
                for ( j=Num-1;j>i;j--)  ans[i]-=LU[i][j]*ans[j];
        }
        for (i=0;i<Num;i++)     printf("solution x[%d] = %lf\n",i+1,ans[i]);
        for (i=0;i<Num;i++)
        {
                temp = 0;
                for (j=0;j<Num;j++)     temp+=A2[i][j]*ans[j];
                printf("%dth:old value = %lf, new value =%lf\n",i+1,B2[i],temp);
        }
        return 0;
}

void mprint(double A[Num][Num], double B[Num])
{
        int i,j;
        for (i=0;i<Num;i++)
        {
                for (j=0;j<Num-1;j++) printf("%lf ",A[i][j]);
                printf("%lf | ",A[i][Num-1]);
                printf("%lf\n",B[i]);
        }
}

--
※ 發信站: 批踢踢實業坊(ptt.cc) 
◆ From: 140.113.7.248
推 wanting0605 :上行下笑阿..OK~~我笑我笑XDDDDD