| > | restart; LU := proc() local filename,fp,A,opt,m,n,i,j,k,p_r,p_c,rk,rtmp,E,F,G,itmp,PAQ; Digits := 20; ### Input: printf("LU decomposition P A Q = L U.\n"); printf("Enter your filename:\n"); filename := scanf("%s")[1]; fp := fopen(filename,READ); m := fscanf(fp,"%d")[1]; n := fscanf(fp,"%d")[1]; A := array(1..m,1..n); for i from 1 to m do for j from 1 to n do A[i,j] := fscanf(fp,"%f")[1]; od; od; fclose(fp); ### Print A: printf("\nBefore: A =\n"); for i from 1 to m do for j from 1 to n do printf("%8.3f",A[i,j]); od; printf("\n"); od; ### Initial row and column exchanges: for i from 1 to n do E[i] := i; F[i] := i; od; rk := 0; for j from 1 to n do ### Find a pivot: p_r := 0; p_c := 0; rtmp := 1.0*10^(-13); for i from rk+1 to m do for k from j to n do if abs(A[i,k]) > rtmp then rtmp := abs(A[i,k]); p_r := i; p_c := k; fi; od; od; # Skip if no pivot found, i.e. all zeros if p_c < j then next; fi; rk := rk +1; if p_r > rk then # switch rows E[rk] := p_r; for k from j to n do rtmp := A[rk,k]; A[rk,k] := A[p_r,k]; A[p_r,k] := rtmp; od; fi; if p_c > j then # switch columns F[rk] := p_c; for i from 1 to m do rtmp := A[i,j]; A[i,j] := A[i,p_c]; A[i,p_c] := rtmp; od; fi; ### Replace 0 entries by Gauss vectors: for i from rk+1 to m do A[i,j] := A[i,j]/A[rk,j]; for k from j+1 to n do A[i,k] := A[i,k] -A[rk,k] * A[i,j]; od; od; od; ### Permute entries of Gauss vectors: for j from 1 to n-1 do for i from j+1 to n-1 do if E[i] > i then rtmp := A[i,j]; A[i,j] := A[E[i],j]; A[E[i],j] := rtmp; fi; od; od; ### Print [L\U]: printf("\nAfter: A stores [L\\U] =\n"; for i from 1 to m do for j from 1 to n do printf("%8.3f",A[i,j]); od; printf("\n"); od; printf("rank of A = %d\n", rk); printf("\nRow exchanges in order: "); for i from 1 to rk do printf("(%d %d)",i,E[i]); od; printf(",\n"); for i from 1 to m do G[i] := i; od; for i from 1 to m do itmp := G[i]; G[i]:=G[E[i]]; G[E[i]]:= itmp; od; printf(" i.e. permutation P = "); for i from 1 to m do printf("%3d",i); od; printf("\n"); printf(" "); for i from 1 to m do printf("%3d",G[i]); od; printf("\n"); printf("\nColumn exchanges in order: "); for i from 1 to rk do printf("(%d %d)", i,F[i]); od; printf(",\n"); for i from 1 to n do G[i] := i; od; for i from 1 to n do if i<F[i] then itmp := G[i]; G[i]:=G[F[i]]; G[F[i]]:= itmp; fi; od; printf(" i.e. permutation Q = "); for i from 1 to n do printf("%3d",i); od; printf("\n"); printf(" "); for i from 1 to n do printf("%3d",G[i]); od; printf("\n"); # printf("L * U:\n"); for i from 1 to m do for j from 1 to i-1 do PAQ[i,j] := 0; for k from 1 to j do PAQ[i,j] := PAQ[i,j] + A[i,k]*A[k,j]; od; # printf("%8.3f",PAQ[i,j]); od; for j from i to n do PAQ[i,j] := A[i,j]; for k from 1 to i-1 do PAQ[i,j] := PAQ[i,j] + A[i,k]*A[k,j]; od; # printf("%8.3f",PAQ[i,j]); od; # printf("\n"); od; for k from rk to 1 by -1 do if k<E[k] then for j from 1 to n do rtmp:=PAQ[k,j]; PAQ[k,j]:=PAQ[E[k],j]; PAQ[E[k],j]:=rtmp; od; fi; if k<F[k] then for i from 1 to m do rtmp:=PAQ[i,k]; PAQ[i,k]:=PAQ[i,F[k]]; PAQ[i,F[k]]:=rtmp; od; fi; od; printf("\nP^(-1) * L * U * Q^(-1) =\n"); for i from 1 to m do for j from 1 to n do printf("%8.3f",PAQ[i,j]); od; printf("\n"); od; printf("\n"); end: |
| > | LU(); |
LU decomposition P A Q = L U.
Enter your filename:
| > | D:/aaa.txt |
Before: A =
1.000 2.000 3.000 4.000
4.000 5.000 6.000 7.000
7.000 8.000 9.000 10.000
10.000 11.000 12.000 13.000
After: A stores [L\U] =
13.000 10.000 12.000 11.000
.308 -2.077 -.692 -1.385
.769 .333 .000 -.000
.538 .667 -.000 .000
rank of A = 2
Row exchanges in order: (1 4)(2 4),
i.e. permutation P = 1 2 3 4
4 1 3 2
Column exchanges in order: (1 4)(2 4),
i.e. permutation Q = 1 2 3 4
4 1 3 2
P^(-1) * L * U * Q^(-1) =
1.000 2.000 3.000 4.000
4.000 5.000 6.000 7.000
7.000 8.000 9.000 10.000
10.000 11.000 12.000 13.000
| > | LU(); |
LU decomposition P A Q = L U.
Enter your filename:
| > | D:/bbb.txt |
Before: A =
1.000 1.000 0.000 0.000 1.000
1.000 2.000 1.000 0.000 2.000
1.000 2.000 2.000 1.000 3.000
1.000 2.000 2.000 2.000 4.000
After: A stores [L\U] =
4.000 2.000 2.000 1.000 2.000
.500 1.000 0.000 .500 -1.000
.750 .500 .500 0.000 0.000
.250 .500 -1.000 .500 0.000
rank of A = 4
Row exchanges in order: (1 4)(2 2)(3 3)(4 4),
i.e. permutation P = 1 2 3 4
4 2 3 1
Column exchanges in order: (1 5)(2 2)(3 3)(4 5),
i.e. permutation Q = 1 2 3 4 5
5 2 3 1 4
P^(-1) * L * U * Q^(-1) =
1.000 1.000 0.000 0.000 1.000
1.000 2.000 1.000 0.000 2.000
1.000 2.000 2.000 1.000 3.000
1.000 2.000 2.000 2.000 4.000
| > | LU(); |
LU decomposition P A Q = L U.
Enter your filename:
| > | D:/ccc.txt |
Before: A =
3.000 17.000 10.000
2.000 4.000 -2.000
6.000 18.000 -12.000
After: A stores [L\U] =
18.000 -12.000 6.000
.944 21.333 -2.667
.222 .031 .750
rank of A = 3
Row exchanges in order: (1 3)(2 3)(3 3),
i.e. permutation P = 1 2 3
3 1 2
Column exchanges in order: (1 2)(2 3)(3 3),
i.e. permutation Q = 1 2 3
2 3 1
P^(-1) * L * U * Q^(-1) =
3.000 17.000 10.000
2.000 4.000 -2.000
6.000 18.000 -12.000
| > | LU(); |
LU decomposition P A Q = L U.
Enter your filename:
| > | D:/ddd.txt |
Before: A =
-3.000 2.000 2.000 -1.000
-12.000 7.000 9.000 -6.000
-3.000 1.000 2.000 -1.000
24.000 -24.000 -12.000 8.000
After: A stores [L\U] =
24.000 -24.000 8.000 -12.000
-.500 -5.000 -2.000 3.000
-.125 .400 .800 -.700
-.125 .200 .500 .250
rank of A = 4
Row exchanges in order: (1 4)(2 2)(3 3)(4 4),
i.e. permutation P = 1 2 3 4
4 2 3 1
Column exchanges in order: (1 1)(2 2)(3 4)(4 4),
i.e. permutation Q = 1 2 3 4
1 2 4 3
P^(-1) * L * U * Q^(-1) =
-3.000 2.000 2.000 -1.000
-12.000 7.000 9.000 -6.000
-3.000 1.000 2.000 -1.000
24.000 -24.000 -12.000 8.000
| > | LU(); |
LU decomposition P A Q = L U.
Enter your filename:
| > | D:/eee.txt |
Before: A =
-2.000 4.000 -2.000 10.000
-2.000 4.000 -2.000 8.000
3.000 3.000 -3.000 6.000
After: A stores [L\U] =
10.000 -2.000 4.000 -2.000
.600 4.200 .600 -1.800
.800 -.095 .857 -.571
rank of A = 3
Row exchanges in order: (1 1)(2 3)(3 3),
i.e. permutation P = 1 2 3
1 3 2
Column exchanges in order: (1 4)(2 4)(3 4),
i.e. permutation Q = 1 2 3 4
4 1 2 3
P^(-1) * L * U * Q^(-1) =
-2.000 4.000 -2.000 10.000
-2.000 4.000 -2.000 8.000
3.000 3.000 -3.000 6.000
| > |