> restart; > # CHOLESKI'S ALGORITHM 6.6 > # > # To factor the positive definite n by n matrix A into LL**T, > # where L is lower triangular. > # > # INPUT: the dimension n; entries A(I,J), 1<=I, J<=n of A. > # > # OUTPUT: the entries L(I,J), 1<=J<=I, 1<=I<=n of L. > # > alg066 := proc() local AA, NAME, INP, OK, N, I, J, A, NN, KK, S, K, JJ, FLAG, OUP; > printf(`This is Choleski 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)\n\n`); > printf(`Place as many entries as desired on each line, but separate `); > printf(`entries with\n`); > printf(`at least one blank.\n\n\n`); > printf(`Has the input file been created? - enter Y or N.\n`); > AA := scanf(`%c`)[1]; > if AA = "Y" or AA = "y" then > printf(`Input the file name in the form - drive:\\name.ext\n`); > printf(`for example: A:\\DATA.DTA\n`); > NAME := scanf(`%s`)[1]; > INP := fopen(NAME,READ,TEXT); > OK := FALSE; > while OK = FALSE do > printf(`Input the dimension n - an integer.\n`); > N := scanf(`%d`)[1]; > if N > 0 then > for I from 1 to N do > for J from 1 to N do > A[I-1,J-1] := fscanf(INP, `%f`)[1]; > od; > od; > OK := TRUE; > fclose(INP); > else printf(`The number must be a positive integer.\n`); > fi; > od; > else > printf(`The program will end so the input file can be created.\n`); > OK := FALSE; > fi; > if OK = TRUE then > # Step 1 > A[0,0] := sqrt(A[0,0]); > # Step 2 > for J from 2 to N do > A[J-1,0] := A[J-1,0]/A[0,0]; > od; > # Step 3 > NN := N-1; > for I from 2 to NN do > # Step 4 > KK := I-1; > S := 0; > for K from 1 to KK do > S := S-A[I-1,K-1]*A[I-1,K-1]; > od; > A[I-1,I-1] := sqrt(A[I-1,I-1]+S); > # Step 5 > JJ := I+1; > for J from JJ to N do > S := 0; > KK := I-1; > for K from 1 to KK do > S := S - A[J-1,K-1]*A[I-1,K-1]; > od; > A[J-1,I-1] := (A[J-1,I-1]+S)/A[I-1,I-1]; > od; > od; > # Step 6 > S := 0; > for K from 1 to NN do > S := S-A[N-1,K-1]*A[N-1,K-1]; > od; > A[N-1,N-1] := sqrt(A[N-1,N-1]+S); > # Step 7 > printf(`Choice of output method:\n`); > printf(`1. Output to screen\n`); > printf(`2. Output to text file\n`); > printf(`Please enter 1 or 2.\n`); > FLAG := scanf(`%d`)[1]; > if FLAG = 2 then > printf(`Input the file name in the form - drive:\\name.ext\n`); > printf(`for example: A:\\OUTPUT.DTA\n`); > NAME := scanf(`%s`)[1]; > OUP := fopen(NAME,WRITE,TEXT); > else > OUP := default; > fi; > fprintf(OUP, `CHOLESKI FACTORIZATION\n\n`); > fprintf(OUP, `The matrix L output by rows:\n`); > for I from 1 to N do > for J from 1 to I do > fprintf(OUP, ` %12.8f`, A[I-1,J-1]); > od; > fprintf(OUP, `\n`); > od; > if OUP <> default then > fclose(OUP): > printf(`Output file %s created successfully`,NAME); > fi; > fi; > RETURN(0); > end; > alg066();