#include <iostream>
#include <stdlib.h>
#include <math.h>
#define MaxD 10

using namespace std;

float Matrix[MaxD][MaxD];
int N;

float Sumof(int cyc, int i, int j)
{
    float fReturnVal = 0;
    for(int k=0; k<cyc; k++)
    {
        fReturnVal += Matrix[i][k]*Matrix[k][j];
    }
    return fReturnVal;
}

void Doolittle()
{
    int i = 0,j=0;
    if(Matrix[0][0] == 0.0)
    {
        printf("不能進(jìn)行LU分解");
        return; 
    }
    
    for(i=1; i<N; i++)
    {
        Matrix[i][0] = Matrix[i][0]/Matrix[0][0];
    }
    
    for(i=1; i<N; i++)
    {
        Matrix[i][i] = Matrix[i][i] - Sumof(i,i,i);
        if(Matrix[i][i] == 0.0)
        {
            printf("不能進(jìn)行LU分解");
            return;
        }
        for(j=i+1; j<N; j++)
        {
            Matrix[i][j] = Matrix[i][j] - Sumof(i,i,j);
            Matrix[j][i] = (Matrix[j][i] - Sumof(i,j,i))/Matrix[i][i];
        }
    }
}

int main(int argc, char *argv[])
{
    int i=0,j=0;
    N = 3;
    do
    {
        printf("請輸入矩陣維數(shù) 不大于10\n");
        scanf("%d",&N);
    }while(N > 10); 
  
    printf("請輸入矩陣元素:\n");
    for(i=0; i<N; i++)
    {
        for(j=0; j<N; j++)
        {
            printf("Matrix[%d][%d]=",i,j);
            scanf("%f",&Matrix[i][j]);
        }
    }
  
    for(i=0; i<N; i++)
    {
        for(j=0; j<N; j++)
        {
            printf("Matrix[%d][%d]=%5f    ",i,j,Matrix[i][j]);
        }
        printf("\n");
    }
    
    printf("Doolittle 分解后\n");
    Doolittle();
    for(i=0; i<N; i++)
    {
        for(j=0; j<N; j++)
        {
            printf("Matrix[%d][%d]=%5f    ",i,j,Matrix[i][j]);
        }
        printf("\n");
    }
    system("PAUSE");	
    return 0;
}