LUP module

download report

Transcript LUP module

LU Decomposition
Greg Beckham, Michael Sedivy
Overview
Step 1
This is handled implicitly in the code by only
calculating the diagonal for β
Step 2
Calculating βij
for(j = 0; j < n; j++) // This is the loop over columns of Crout's method
{
for(i = 0; i < j; i++) // Equation (2.3.12) except for i = j
{
sum = a[i][j];
for(k = 0; k < i; k++) sum -= a[i][k] * a[k][j];
a[i][j] = sum;
}
…
}
Calculating αij
for(j = 0; j < n; j++) // This is the loop over columns of Crout's method
{
…
for(i = j; i < n; i++)// This is i=j of equation (2.3.12) and i=j+1
{
// N-1 of equation (2.3.13)
sum = a[i][j];
for(k = 0; k < j; k++) sum -= a[i][k] * a[k][j];
a[i][j] = sum;
…
}
…
if(j != n - 1) // Divide by the pivot element
{
dum = 1.0/(a[j][j]);
for(i = j + 1; i < n; i++) a[i][j] *= dum;
}
…
}
Pivoting
• Initially finds largest element in each row
• Used as a “scaling factor”, not sure of use other than
to rollover
for(i = 0; i < n; i++) // Loop over the rows to get implicit scaling
{
// information
big = 0.0;
for(j = 0; j < n; j++)
{
if((temp = fabs(a[i][j])) > big) big = temp;
}
if (big == 0.0)
{
printf("ERROR: Singular matrix\n");
}
// non-zero largest element.
vv[i] = 1.0/big;
// Save the scaling
}
Pivoting
• Finds maximum
if((dum = vv[i] * fabs(sum)) >= big)
{
// Is the figure of merit for the pivot better than the best so far?
big = dum;
imax = i;
}
Pivoting
• Performs row interchanges
if(j != imax)
// Do we need to interchange
rows?
{
for(k = 0; k < n; k++) // Interchange rows
{
dum = a[imax][k];
a[imax][k] = a[j][k];
a[j][k] = dum;
}
d = -d;
// change the parity of d
vv[imax] = vv[j]; // interchange scale
factor
}
indx[j] = imax;
Related Questions
• What is the advantage of LU(P) solver over GJ(P) solver?
(Complexity)
– Both are O(N3)
– After LU(P) is solved, more solutions supposed to be found in O (N2)
• Are you keeping L and U in the same matrix, or separate?
Advantage/disadvantage?
– LU are being created in place in the same matrix.
– The advantage to this strategy is lower memory usage
– The disadvantage is that the original matrix is lost
• I am somewhat confused with extraction of P in decomposition,
and how it is then used in eq solving. Can you elaborate more?
Related Questions
• Cormen et al., p 824, used a single array
instead of P. Needs careful explanation.
– From Cormen et al. 825. “we dynamically
maintain the permutation matrix P as an array π,
where π[i] = j means that the ith row of P contains
a 1 in column j
|2|
|0 1 0|
π = |3| => P | 0 0 1 |
|1|
|1 0 0|
Related Questions
• How complex equations are solved? (in Text)
– If only the right hand side vector is complex then
the operation can be performed by solving for the
real part, then the imaginary
– If the matrix itself is complex then
• Rewrite the algorithm for complex values
• Split the real and imaginary parts into separate real
number and solve using existing algorithm
– A*x–C*y=b
– C*x+A*y=b
GJ vs. LUP: we found lup faster than gj, but…
Variable Array Size
12000
10000
8000
Time (ms) 6000
LUP
Gauss-Jordan
4000
2000
0
Array Size
GJ vs. LUP: lup not faster amortized
Variable Solution Size
900
800
700
600
500
Time (ms)
LUP
400
Gauss-Jordan
300
200
100
0
Number of Solutions
GJ vs. LUP
GJ vs. LUP Difference
Average Difference is 2.765471
4
3.5
3
2.5
2
1.5
GJ vs. LUP
Difference
1
0.5
940
860
780
700
620
540
460
380
300
220
140
60
0
GJ vs. LUP Difference
Average Difference is 1.255924
2.5
2
1.5
GJ vs. LUP
Difference
1
0.5
20
100
180
260
340
420
500
580
660
740
820
900
980
0