Chapter 15 C-implementation PA = LU

1
Chapter 15 C-implementationPA LU

• Speaker Lung-Sheng Chien

2
OutLine

• Data structure of full matrix
• Implementation of PALU
• Visual Studio 2005 How To
• Linux machine How to compile

3
row-major versus col-major
Logical index 2D
col-major based
6
-2
2
4
12
-8
6
10
3
-13
9
3
6
-2
2
4
-6
4
1
-18
12
-8
6
10
3
-13
9
3
row-major based
-6
4
1
-18
6
-2
2
4
12
-8
6
10
We choose col-major representation
3
-13
9
3
-6
4
1
-18
4
col-major C starts index from 0
0
6
1
matrix.h
2
12
3
3
-6
4
-2
5
-8
6
-13
7
col-major based
4
8
legal
2
9
6
10
9
11
1
12
is useless
1
6
-2
2
4
13
4
is forbidden
2
12
-8
6
10
10
14
3
-13
9
3
3
3
15
-6
4
1
-18
-18
16
5
Relaxation for high precision package
global.h
Large 1D array, int is 4-byte, only supports up
to 2G elements
We use col-major, however you can implement both.
http//crd.lbl.gov/dhbailey/mpdist/
6
Data structure of matrix and its method
matrix.h
constructor (???) zeros
Arithmetic matvec, matsub, norm
1
4
destructor (???) dealloc
I/O disp
2
5
Index rule we do it explicitly
3
7
matrix.cpp
constructor 1
empty handler
1
1
set parameter
5
2
2
3
contiguous memory block
4
3
6
1
2
12
3
3
5
8
matrix.cpp
constructor 2
0
6
1
4
2
12
3
3
-6
4
low
-2
5
-8
6
1
-13
7
high
4
8
2
2
9
3
6
10
9
11
4
4
1
12
pointer arithmetic
13
4
10
14
3
15
5
-18
16
9
destructor
0
matrix.cpp
6
1
2
12
3
3
1
-6
4
-2
5
2
-8
6
3
-13
7
4
1
free pointer array
8
2
free physic storage
2
9
6
10
9
11
free matrix handler
1
12
3
13
4
10
14
3
15
-18
16
10
Input/Output
matrix.cpp
6
-2
2
4
12
-8
6
10
Show in standard form human readable
3
-13
9
3
-6
4
1
-18
11
Simple driver
main.cpp
Execution result
12
duplicate of a matrix / vector
matrix.cpp
verification
We require matrix handler to verify
configuration, this is very important for
debugging.
13
y Ax
matrix.cpp
Outer-product formulation
14
Matrix / vector norm
matrix.cpp
Maximum of absolute column sum
Maximum of absolute row sum
2-norm is difficult to compute, we will discuss
this later
15
OutLine

• Data structure of full matrix
• Implementation of PALU
• Visual Studio 2005 How To
• Linux machine How to compile

16
Algorithm ( PA LU ) 1
Given full matrix
, construct initial lower triangle matrix
use permutation vector
to record permutation matrix
verification
(1) we store lower triangle matrix L into storage
of matrix A
(2) The reason to construct type int_matrix is
for permutation matrix P,
since we can check dimension between A and P.
17
Algorithm ( PA LU ) 2
let
,
and
Note that P is a column vector
for
we have compute
update original matrix
stores in lower triangle matrix
18
Algorithm ( PA LU ) 3
such that
find a
1
, NOT efficient
swap row
and row
2
define permutation matrix
, we have
then after swapping rows
19
Algorithm ( PA LU ) 4
3
compute
for
if
then
by swapping
and
compute
4
then
endif
20
Algorithm ( PA LU ) 5
we store lower triangle matrix into storage of
matrix A, this code is different from code in
MATLAB.
decompose
5
How to compute
efficiently?
21
Algorithm ( PA LU ) 6
component-wise update
Observation we implement with column-major form,
column-wise update is good
22
Algorithm ( PA LU ) 7
where
by updating
by updating matrix
(recursion is done)
then
endfor
23
Algorithm ( PA LU ) 8
6
-2
2
4
1
12
-8
6
10
12
-8
6
10
0.25
1
-11
7.5
0.5
3
-13
9
3
-0.5
0
1
4
-13
-6
4
1
-18
0.5
-2/11
1/11
1
3/11
Store L and U into original A
Question 1 how to write down a test driver
Question 2 how to do backward and forward
24
Test driver 1
main.cpp
6
-2
2
4
12
-8
6
10
3
-13
9
3
-6
4
1
-18
1
Remember to duplicate matrix A since A is
modified when execute LU decomposition.
25
Test driver 2
2
5
6
3
7
8
4
Linear solver do backward and forward
5
-6.9306
17.9583
26.5833
6
7.3333
26
Exercise

• Implement PA LU.
• Implement forward and backward.
• Write a driver to test your linear solver.
• Any other implementation for full matrix?

27
OutLine

• Data structure of full matrix
• Implementation of PALU
• Visual Studio 2005 How To
• Linux machine How to compile

28
Visual Studio 2005 How To 1
Create a new project File ? Project
29
Visual Studio 2005 How To 2
Choose Win32 console application, this is the
same as what we do in VC 6.0
Project name vc 2005 will create a directory,
the same name as project name
30
Visual Studio 2005 How To 3
We choose next bottom to set an empty project,
DO NOT press Finish bottom.
31
Visual Studio 2005 How To 4
Empty project, this is very important
32
Visual Studio 2005 How To 5
3
Add main.cpp to this project
Empty project
1
Add new item (source file or header file) to this
project
2
33
Visual Studio 2005 How To 6
Write source code hello world
Compile Build ? Build matrix
34
Visual Studio 2005 How To 7
Result of compilation, including linking success
Additional directory debug appears
Executable file matrix.exe in this debug
directory
35
Visual Studio 2005 How To 8
Execute matrix.exe Debug ? Start Without
Debugging
Execution in VC 6.0
36
Visual Studio 2005 How To 9
Suppose we want to add several source files into
this project.
1
Copy files into directory /matrix/matrix
Add source files into the project by Project ?
Add Existing item
2
37
Visual Studio 2005 How To 10
Select the files you want to add
All souce files in this project
38
Visual Studio 2005 How To 11
Remember that executable file is in directory
/matrix/debug, you can use command window to
execute matrix.exe
39
Visual Studio 2005 How To 12
There is another debug directory in
/matrix/matrix, this directory contains object
files, no executable file
40
OutLine

• Data structure of full matrix
• Implementation of PALU
• Visual Studio 2005 How To
• Linux machine How to compile

41
Linux machine How to compile 1
Suppose we put all source file into directory
matrix_constructor and upload to workstation.
1
Remove all VC related files.
2
use command qmake -project to generate project
file
3
42
Linux machine How to compile 2
use command qmake matrix_constructor.pro to
generate Makefile Later on, we will introduce
Makefile
4
use command make to compile your codes and
generate executable file
5
43
Linux machine How to compile 3
6
Execute executable file by ./matrix_constructor