矩阵的LU分解(完整版)实用资料.doc

矩阵的LU分解（完整版）实用资料(可以直接使用，可编辑 完整版实用资料，欢迎下载）矩阵的LU分解/* 矩阵的LU分解* To factor the n by n matrix A = (A(I,J) into the product of the* lower triangular matrix L = (L(I,J) and U = (U(I,J), that is* A = LU, where the main diagonal of either L or U consists of all ones:* INPUT: dimension n; the entries A(I,J), 1<=I, J<=n, of A;* the diagonal L(1,1), ., L(N,N) of L or the diagonal* U(1,1), ., U(N,N) of U.* OUTPUT: the entries L(I,J), 1<=J<=I, 1<=I<=n of L and the entries* U(I,J), I<=J<=n, 1<=I<=n of U.*/ 1.C语言版*#include#include#define ZERO 1.0E-20#define true 1#define false 0double absval(double);void INPUT(int *, double 10, int *, int *);void OUTPUT(int, double 10, int);main() double A1010,XL10; double S,SS; int N,M,I,J,ISW,JJ,K,KK,OK; INPUT(&OK, A, &N, &ISW); if (OK) for (I=1; I<=N; I+) XLI-1 = 1.0; /* STEP 1 */ if (absval(A00) <= ZERO) OK = false; else /* the entries of L below the main diagonal will be placed in the corresponding entries of A; the entries of U above the main diagonal will be placed in the corresponding entries of A; the main diagonal which was NOT input will become the main diagonal of A; the input main diagonal of L or U is, of course, placed in XL. */ A00 = A00 / XL0; /* STEP 2 */ for (J=2; J<=N; J+) if (ISW = 0) /* first row of U */ A0J-1 = A0J-1 / XL0; /* first column of L */ AJ-10 = AJ-10 / A00; else /* first row of U */ A0J-1 = A0J-1 / A00; /* first column of L */ AJ-10 = AJ-10 / XL0; /* STEP 3 */ M = N - 1; I = 2; while (I <= M) && OK) /* STEP 4 */ KK = I - 1; S = 0.0; for (K=1; K<=KK; K+) S = S - AI-1K-1 * AK-1I-1; AI-1I-1 = ( AI-1I-1 + S ) / XLI-1; if (absval(AI-1I-1) <= ZERO) OK = false; else /* STEP 5 */ JJ = I + 1; for (J=JJ; J<=N; J+) SS = 0.0; S = 0.0; for (K=1; K<=KK; K+) SS = SS - AI-1K-1 * AK-1J-1; S = S - AJ-1K-1 * AK-1I-1; if (ISW = 0) /* Ith row of U */ AI-1J-1 = (AI-1J-1 + SS) / XLI-1; /* Ith column of L */ AJ-1I-1 = (AJ-1I-1 + S) / AI-1I-1; else /* Ith row of U */ AI-1J-1 = (AI-1J-1 + SS) / AI-1I-1; /* Ith column of L */ AJ-1I-1 = (AJ-1I-1 + S) / XLI-1; I+; if (OK) /* STEP 6 */ S = 0.0; for (K=1; K<=M; K+) S = S - AN-1K-1 * AK-1N-1; AN-1N-1 = (AN-1N-1 + S) / XLN-1; /* If A(N-1,N-1) = 0 then A = LU but the matrix is singular. Process is complete, all entries of A have been determined. */ /* STEP 7 */ OUTPUT(N, A, ISW); if (!OK) printf("System has no unique solutionn"); return 0;void INPUT(int *OK, double A10, int *N, int *ISW) int I, J, FLAG; char AA; char NAME30; FILE *INP; printf("This is the general LU factorization method.n"); printf("The array will be input from a text file in the order:n"); printf("A(1,1), A(1,2), ., A(1,N), A(2,1), A(2,2), ., A(2,N),n"); printf("., A(N,1), A(N,2), ., A(N,N)nn"); printf("Place as many entries as desired on each line, but separate "); printf("entries withn"); printf("at least one blank.nnn"); printf("Has the input file been created? - enter Y or N.n"); scanf("%c",&AA); *OK = false; if (AA = 'Y') | (AA = 'y') printf("Input the file name in the form - drive:name.extn"); printf("for example: A:DATA.DTAn"); scanf("%s", NAME); INP = fopen(NAME, "r"); *OK = false; while (!(*OK) printf("Input the dimension n - an integer.n"); scanf("%d", N); if (*N > 0) for (I=1; I<=*N; I+) for (J=1; J<=*N; J+) fscanf(INP, "%lf", &AI-1J-1); fscanf(INP, "n"); *OK = true; fclose(INP); else printf("The number must be a positive integer.n"); printf("Choice of diagonals:n"); printf("1. Diagonal of L consists of onesn"); printf("2. Diagonal of U consists of onesn"); printf("Please enter 1 or 2.n"); scanf("%d", &FLAG); if (FLAG = 1) *ISW = 0; else *ISW = 1; else printf("The program will end so the input file can be created.n");void OUTPUT(int N, double A10, int ISW) int I, J, FLAG; char NAME30; FILE *OUP; printf("Choice of output method:n"); printf("1. Output to screenn"); printf("2. Output to text filen"); printf("Please enter 1 or 2.n"); scanf("%d", &FLAG); if (FLAG = 2) printf("Input the file name in the form - drive:name.extn"); printf("for example: A:OUTPUT.DTAn"); scanf("%s", NAME); OUP = fopen(NAME, "w"); else OUP = stdout; fprintf(OUP, "GENERAL LU FACTORIZATIONnn"); if (ISW = 0) fprintf(OUP, "The diagonal of L consists of all entries = 1.0n"); else fprintf(OUP, "The diagonal of U consists of all entries = 1.0n"); fprintf(OUP, "nEntries of L below/on diagonal and entries of U above"); fprintf(OUP, "/on diagonaln"); fprintf(OUP, "- output by rows in overwrite format:n"); for (I=1; I<=N; I+) for (J=1; J<=N; J+) fprintf(OUP, " %11.8f", AI-1J-1); fprintf(OUP, "n"); fclose(OUP);/* Absolute Value Function */double absval(double val) if (val >= 0) return val; else return -val;*C语言版END 2.Mathematica语言版*Print"n"Print"A(1,1), A(2,1), ., A(1,n), A(2,1), A(2,2), .n"Print"A(2,n), ., A(n,1), A(n,2), ., A(n,n)n"Print"n"Print"Place as many entries as desired on each line, but n"Print"separate entries with at least one blankn"Print"n"Print"n"AA = InputString"This is General LU factorization methodn The array will be input from a text file in the order:(see screen)n Has the input file been created?n Enter 'yes' or 'no'n"IfAA = "yes" | AA = "y" | AA = "Y", NAME = InputString"Input the file name in the form - n drive: name.ext for examplen A:DATA.DTAn" INP = OpenReadNAME; OK = 0; WhileOK = 0, n = Input"Input the dimension n - an integern" Ifn > 0, Fori = 1,i <= n,i+, ForJ = 1,J <= n,J+, Ai-1,J-1 = ReadINP,Number; ; ; OK = 1; CloseINP, Input"Number must be a positive integern n Press 1 enter to continuen" ; FLAG = Input"Choice of diagonals:n 1. Diagonal of L consists of onesn 2. Diagonal of U consists of onesn Please enter 1 or 2.n" IfFLAG = 1, ISW = 0, ISW = 1; ; , Input"This program will end so the input filen can be created.n n Press 1 enter to continuen" OK = 0;IfOK = 1, Fori = 1,i <= n,i+, XLi-1 = 1; ; (* Step 1 *) IfAbsA0,0 <=10-20, OK = 0, (* The entries below the main diagonal will be placed in the corresponding entries of A; the entries of U above the main diagonal will be placed in the corres- ponding entries of A; the main diagonal which was not inputed will become the main diagonal of A; the input main diagonal of L or U is placed in XL *) A0,0 = A0,0/XL0; (* Step 2 *) ForJ = 2,J <= n,J+, IfISW = 0, (* First row of U *) A0,J-1 = A0,J-1/XL0; (* First column of L *) AJ-1,0 = AJ-1,0/A0,0, (* first row of U *) A0,J-1 = A0,J-1/A0,0; (* First column of L *) AJ-1,0 = AJ-1,0/XL0; ; ; (* Step 3 *) M = n-1; i = 2; Whilei <= M && OK = 1, (* Step 4 *) KK = i-1; S = 0; ForK = 1,K <= KK,K+, S = S-Ai-1,K-1*AK-1,i-1; ; Ai-1,i-1 = (Ai-1,i-1+S)/XLi-1; IfAbsAi-1,i-1 <= 10-20, OK = 0, (* Step 5 *) JJ = i+1; ForJ = JJ,J <= n,J+, SS = 0; S = 0; ForK = 1,K <= KK,K+, SS = SS-Ai-1,K-1*AK-1,J-1; S = S-AJ-1,K-1*AK-1,i-1; ; IfISW = 0, (* Ith row of U *) Ai-1,J-1 = (Ai-1,J-1+SS)/XLi-1; (* Ith column of L *) AJ-1,i-1 = (AJ-1,i-1+S)/Ai-1,i-1, (* Ith row of U *) Ai-1,J-1 = (Ai-1,J-1+SS)/Ai-1,i-1; (* Ith column of L *) AJ-1,i-1 = (AJ-1,i-1+S)/XLi-1; ; ; ; i = i+1; ; IfOK = 1, (* Step 6 *) S = 0; ForK = 1,K <= M,K+, S = S-An-1,K-1*AK-1,n-1 ; An-1,n-1 = (An-1,n-1+S)/XLn-1; (* If A(N-1,N-1)=0 then A=LU but the matrix is singular. Process is complete, all entries of A have been determined. *) (* Step 7 *) FLAG = Input"Select output destinationn 1. Screenn 2. Text filen Enter 1 or 2n" IfFLAG = 2, NAME = InputString"Input the file namen For example: output.dtan" OUP = OpenWriteNAME,FormatType->OutputForm, OUP = "stdout" ; WriteOUP,"GENERAL LU FACTORIZATIONn" WriteOUP,"n" IfISW = 0, WriteOUP,"The diagonal of L consists of all entries = 1.0n", WriteOUP,"The diagonal of U consists of all entries = 1.0n" ; WriteOUP,"n" WriteOUP,"Entries of L below/on diagonal and entriesn of U above/on diagonal n output by rows in overwrite format:n" Fori = 1,i <= n,i+, ForJ = 1,J <= n,J+, WriteOUP,NAi-1,J-1,9; ; WriteOUP,"n" ; WriteOUP,"n" IfOUP = "OutputStream",NAME," 3", Print"Output file: ",NAME," created successfullyn" CloseOUP; ; ; ; IfOK = 0, Print"System has no unique solutionn" ;*Mathematica语言版END 3.Matlab语言版* syms('AA', 'NAME', 'INP', 'OK', 'N', 'I', 'J', 'A'); syms('FLAG', 'ISW', 'XL', 'M', 'KK', 'S', 'K', 'JJ'); syms('SS', 'OUP', 's'); TRUE = 1; FALSE = 0; fprintf(1,'This is the general LU factorization method.n'); fprintf(1,'The array will be input from a text file in the order:n'); fprintf(1,'A(1,1), A(1,2), ., A(1,N), n') fprintf(1,'A(2,1), A(2,2), ., A(2,N),n'); fprintf(1,'., A(N,1), A(N,2), ., A(N,N)nn'); fprintf(1,'Place as many entries as desired on each line, but separaten'); fprintf(1,'entries withn'); fprintf(1,'at least one blank.nnn'); fprintf(1,'Has the input file been created? - enter Y or N.n'); AA = input(' ','s'); if AA = 'Y' | AA = 'y' fprintf(1,'Input the file name in the form - drive:name.extn'); fprintf(1,'for example: A:DATA.DTAn'); NAME = input(' ','s'); INP = fopen(NAME,'rt'); OK = FALSE; while OK = FALSE fprintf(1,'Input the dimension n - an integer.n'); N = input(' '); if N > 0 A = zeros(N,N); XL = zeros(1,N); for I = 1 : N for J = 1 : N A(I,J) = fscanf(INP, '%f',1); end; end; OK = TRUE; fclose(INP); else fprintf(1,'The number must be a positive integer.n'); end; end; fprintf(1,'Choice of diagonals:n'); fprintf(1,'1. Diagonal of L consists of onesn'); fprintf(1,'2. Diagonal of U consists of onesn'); fprintf(1,'Please enter 1 or 2.n'); FLAG = input(' '); if FLAG = 1 ISW = 0; else ISW = 1; end else fprintf(1,'The program will end so the input file can be created.n'); OK = FALSE; end; if OK = TRUE for I = 1 : N XL(I) = 1; end;% STEP 1 if abs(A(1,1) <= 1.0e-20 OK = FALSE; else% the entries of L below the main diagonal will be placed % in the corresponding entries of A; the entries of U % above the main diagonal will be placed in the % corresponding entries of A; the main diagonal which % was not input will become the main diagonal of A; % the input main diagonal of L or U is, % of course, placed in XL A(1,1) = A(1,1)/XL(1);% STEP 2 for J = 2 : N if ISW = 0 % first row of U A(1,J) = A(1,J)/XL(1);% first column of L A(J,1) = A(J,1)/A(1,1); else% first row of U A(1,J) = A(1,J)/A(1,1);% first column of L A(J,1) = A(J,1)/XL(1); end; end;% STEP 3 M = N-1; I = 2; while I <= M & OK = TRUE % STEP 4 KK = I-1; S = 0; for K = 1 : KK S = S-A(I,K)*A(K,I); end; A(I,I) = (A(I,I)+S)/XL(I); if abs(A(I,I) <= 1.0e-20 OK = FALSE; else% STEP 5 JJ = I+1; for J = JJ : N SS = 0; S = 0; for K = 1 : KK SS = SS-A(I,K)*A(K,J); S = S-A(J,K)*A(K,I); end; if ISW = 0 % Ith row of U A(I,J) = (A(I,J)+SS)/XL(I);% Ith column of L A(J,I) = (A(J,I)+S)/A(I,I); else% Ith row of U A(I,J) = (A(I,J)+SS)/A(I,I);% Ith column of L A(J,I) = (A(J,I)+S)/XL(I); end; end; end; I = I+1; end; if OK = TRUE % STEP 6 S = 0; for K = 1 : M S = S-A(N,K)*A(K,N); end; A(N,N) = (A(N,N)+S