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

《矩阵的LU分解(完整版)实用资料.doc》由会员分享，可在线阅读，更多相关《矩阵的LU分解(完整版)实用资料.doc（30页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、矩阵的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

2、=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 *, do

3、uble 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

4、 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

5、 (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=

6、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

7、= (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 0) for (I=1; I=*N; I+) for (J=1; J=*N; J+) fscanf(INP, %

8、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

9、; 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

10、(%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); e

11、lse 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= 0) return val; else return -val;*C语言版END 2.Math

12、ematica语言版*Printn;PrintA(1,1), A(2,1), ., A(1,n), A(2,1), A(2,2), .n;PrintA(2,n), ., A(n,1), A(n,2), ., A(n,n)n;Printn;PrintPlace as many entries as desired on each line, but n;Printseparate entries with at least one blankn;Printn;Printn;AA = InputStringThis is General LU factorization methodn The a

13、rray will be input from a text file in the order:(see screen)n Has the input file been created?n Enter yes or non;IfAA = yes | AA = y | AA = Y, NAME = InputStringInput the file name in the form - n drive: name.ext for examplen A:DATA.DTAn; INP = OpenReadNAME; OK = 0; WhileOK = 0, n = InputInput the

14、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, InputNumber must be a positive integern n Press 1 enter to continuen; ; FLAG = InputChoice of diagonals:n 1. Diagonal of L consists of onesn 2. Diagonal of U consists of onesn Ple

15、ase enter 1 or 2.n; IfFLAG = 1, ISW = 0, ISW = 1; ; , InputThis 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 corr

16、esponding 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

17、 = 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

18、*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

19、= (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 OutputForm, OUP = stdout ; WriteOUP,GENERAL LU FACTORIZATIONn; WriteOUP,n; IfISW = 0, WriteOUP,The diagonal of

20、 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 0 A = zeros(N,N); XL = zeros(1,N); for I = 1 :

21、 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 o

22、r 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 plac

23、ed % 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);%

24、 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 : K

25、K 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