Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
Joost Opschoor
NumCSE
Commits
b9817f75
Commit
b9817f75
authored
Sep 17, 2019
by
ralfh
Browse files
Reformatting of codes
parent
023c6ad2
Changes
4
Hide whitespace changes
Inline
Side-by-side
Assignments/Codes/MatVec/ArrowMatrix/solutions/arrowmatvec.cpp
View file @
b9817f75
#include
<iostream>
#include
<iomanip>
#include
<iostream>
#include
<Eigen/Dense>
...
...
@@ -17,31 +17,27 @@ using namespace Eigen;
/* SAM_LISTING_BEGIN_0 */
void
arrow_matrix_2_times_x
(
const
VectorXd
&
d
,
const
VectorXd
&
a
,
const
VectorXd
&
x
,
VectorXd
&
y
)
{
assert
(
d
.
size
()
==
a
.
size
()
&&
a
.
size
()
==
x
.
size
()
&&
"Vector size must be the same!"
);
int
n
=
d
.
size
();
// In this lines, we extract the blocks used to construct the matrix A.
VectorXd
d_head
=
d
.
head
(
n
-
1
);
VectorXd
a_head
=
a
.
head
(
n
-
1
);
MatrixXd
d_diag
=
d_head
.
asDiagonal
();
MatrixXd
A
(
n
,
n
);
// We build the matrix A using the "comma initialization": each expression separated
// by a comma is a "block" of the matrix we are building.
// d\_diag is the top left (n-1)x(n-1) block
// a\_head is the top right vertical vector
// a\_head.transpose() is the bottom left horizontal vector
// d(n-1) is a single element (a 1x1 matrix), on the bottom right corner
// This is how the matrix looks like:
// A = | D | a |
// |-----+--------|
// | a\^T | d(n-1) |
A
<<
d_diag
,
a_head
,
a_head
.
transpose
(),
d
(
n
-
1
);
y
=
A
*
A
*
x
;
assert
(
d
.
size
()
==
a
.
size
()
&&
a
.
size
()
==
x
.
size
()
&&
"Vector size must be the same!"
);
int
n
=
d
.
size
();
// In this lines, we extract the blocks used to construct the matrix A.
VectorXd
d_head
=
d
.
head
(
n
-
1
);
VectorXd
a_head
=
a
.
head
(
n
-
1
);
MatrixXd
d_diag
=
d_head
.
asDiagonal
();
MatrixXd
A
(
n
,
n
);
// We build the matrix A using the "comma initialization": each expression
// separated by a comma is a "block" of the matrix we are building. d\_diag is
// the top left (n-1)x(n-1) block a\_head is the top right vertical vector
// a\_head.transpose() is the bottom left horizontal vector
// d(n-1) is a single element (a 1x1 matrix), on the bottom right corner
// This is how the matrix looks like:
// A = | D | a |
// |-----+--------|
// | a\^T | d(n-1) |
A
<<
d_diag
,
a_head
,
a_head
.
transpose
(),
d
(
n
-
1
);
y
=
A
*
A
*
x
;
}
/* SAM_LISTING_END_0 */
...
...
@@ -53,52 +49,47 @@ void arrow_matrix_2_times_x(const VectorXd &d, const VectorXd &a,
* @param[out] y The vector y = A*A*x
*/
/* SAM_LISTING_BEGIN_1 */
void
efficient_arrow_matrix_2_times_x
(
const
VectorXd
&
d
,
const
VectorXd
&
a
,
const
VectorXd
&
x
,
VectorXd
&
y
)
{
assert
(
d
.
size
()
==
a
.
size
()
&&
a
.
size
()
==
x
.
size
()
&&
void
efficient_arrow_matrix_2_times_x
(
const
VectorXd
&
d
,
const
VectorXd
&
a
,
const
VectorXd
&
x
,
VectorXd
&
y
)
{
assert
(
d
.
size
()
==
a
.
size
()
&&
a
.
size
()
==
x
.
size
()
&&
"Vector size must be the same!"
);
int
n
=
d
.
size
();
// Notice that we can compute (A*A)*x more efficiently using
// A*(A*x). This is, in fact, performing two matrix vector
// multiplications
// instead of a more expensive matrix-matrix multiplication.
// Therefore, as first step, we need a way to efficiently
// compute A*x
// This function computes A*x. you can use it
// by calling A\_times\_x(x).
// This is the syntax for lambda functions: notice the extra
// [variables] code. Each variable written within [] brackets
// will be captured (i.e. seen) inside the lambda function.
// Without \&, a copy of the variable will be performed.
// Notice that A = D + H + V, s.t. A*x = D*x + H*x + V*x
// D*x can be rewritten as d*x componentwise
// H*x is zero, except at the last component
// V*x is only affected by the last component of x
auto
A_times_x
=
[
&
a
,
&
d
,
n
]
(
const
VectorXd
&
x
)
{
// This takes care of the diagonal (D*x)
// Notice: we use d.array() to tell Eigen to treat
// a vector as an array. As a result: each operation
// is performed componentwise.
VectorXd
Ax
=
(
d
.
array
()
*
x
.
array
()
).
matrix
();
// H*x only affects the last component of A*x
// This is a dot product between a and x with the last
// component removed
Ax
(
n
-
1
)
+=
a
.
head
(
n
-
1
).
dot
(
x
.
head
(
n
-
1
));
// V*x is equal to the vector
// (a(0)*x(n-1), ..., a(n-2)*x(n-1), 0)
Ax
.
head
(
n
-
1
)
+=
x
(
n
-
1
)
*
a
.
head
(
n
-
1
);
return
Ax
;
};
// <=> y = A*(A*x)
y
=
A_times_x
(
A_times_x
(
x
));
int
n
=
d
.
size
();
// Notice that we can compute (A*A)*x more efficiently using
// A*(A*x). This is, in fact, performing two matrix vector
// multiplications
// instead of a more expensive matrix-matrix multiplication.
// Therefore, as first step, we need a way to efficiently
// compute A*x
// This function computes A*x. you can use it
// by calling A\_times\_x(x).
// This is the syntax for lambda functions: notice the extra
// [variables] code. Each variable written within [] brackets
// will be captured (i.e. seen) inside the lambda function.
// Without \&, a copy of the variable will be performed.
// Notice that A = D + H + V, s.t. A*x = D*x + H*x + V*x
// D*x can be rewritten as d*x componentwise
// H*x is zero, except at the last component
// V*x is only affected by the last component of x
auto
A_times_x
=
[
&
a
,
&
d
,
n
](
const
VectorXd
&
x
)
{
// This takes care of the diagonal (D*x)
// Notice: we use d.array() to tell Eigen to treat
// a vector as an array. As a result: each operation
// is performed componentwise.
VectorXd
Ax
=
(
d
.
array
()
*
x
.
array
()).
matrix
();
// H*x only affects the last component of A*x
// This is a dot product between a and x with the last
// component removed
Ax
(
n
-
1
)
+=
a
.
head
(
n
-
1
).
dot
(
x
.
head
(
n
-
1
));
// V*x is equal to the vector
// (a(0)*x(n-1), ..., a(n-2)*x(n-1), 0)
Ax
.
head
(
n
-
1
)
+=
x
(
n
-
1
)
*
a
.
head
(
n
-
1
);
return
Ax
;
};
// <=> y = A*(A*x)
y
=
A_times_x
(
A_times_x
(
x
));
}
/* SAM_LISTING_END_1 */
...
...
@@ -106,48 +97,41 @@ void efficient_arrow_matrix_2_times_x(const VectorXd &d,
* Repeat tests 10 times, and output the minimal runtime
* amongst all times. Test both the inefficient and the efficient
* versions.
*/
*/
void
runtime_arrow_matrix
()
{
/* SAM_LISTING_BEGIN_3 */
std
::
cout
<<
std
::
setw
(
8
)
<<
"n"
<<
std
::
setw
(
15
)
<<
"original"
<<
std
::
setw
(
15
)
<<
"efficient"
<<
std
::
endl
;
for
(
unsigned
n
=
4
;
n
<=
2048
;
n
*=
2
)
{
// Number of repetitions
unsigned
int
repeats
=
10
;
Timer
timer
,
timer_eff
;
// Repeat test many times
for
(
unsigned
int
r
=
0
;
r
<
repeats
;
++
r
)
{
// Create test input using random vectors
Eigen
::
VectorXd
a
=
Eigen
::
VectorXd
::
Random
(
n
),
d
=
Eigen
::
VectorXd
::
Random
(
n
),
x
=
Eigen
::
VectorXd
::
Random
(
n
),
y
;
// Compute times
timer
.
start
();
arrow_matrix_2_times_x
(
d
,
a
,
x
,
y
);
timer
.
stop
();
// Compute times for efficient implementation
timer_eff
.
start
();
efficient_arrow_matrix_2_times_x
(
d
,
a
,
x
,
y
);
timer_eff
.
stop
();
}
// Print results
std
::
cout
<<
std
::
setw
(
8
)
<<
n
<<
std
::
scientific
<<
std
::
setprecision
(
3
)
<<
std
::
setw
(
15
)
<<
timer
.
min
()
<<
std
::
setw
(
15
)
<<
timer_eff
.
min
()
<<
std
::
endl
;
/* SAM_LISTING_BEGIN_3 */
std
::
cout
<<
std
::
setw
(
8
)
<<
"n"
<<
std
::
setw
(
15
)
<<
"original"
<<
std
::
setw
(
15
)
<<
"efficient"
<<
std
::
endl
;
for
(
unsigned
n
=
4
;
n
<=
2048
;
n
*=
2
)
{
// Number of repetitions
unsigned
int
repeats
=
10
;
Timer
timer
,
timer_eff
;
// Repeat test many times
for
(
unsigned
int
r
=
0
;
r
<
repeats
;
++
r
)
{
// Create test input using random vectors
Eigen
::
VectorXd
a
=
Eigen
::
VectorXd
::
Random
(
n
),
d
=
Eigen
::
VectorXd
::
Random
(
n
),
x
=
Eigen
::
VectorXd
::
Random
(
n
),
y
;
// Compute times
timer
.
start
();
arrow_matrix_2_times_x
(
d
,
a
,
x
,
y
);
timer
.
stop
();
// Compute times for efficient implementation
timer_eff
.
start
();
efficient_arrow_matrix_2_times_x
(
d
,
a
,
x
,
y
);
timer_eff
.
stop
();
}
/* SAM_LISTING_END_3 */
// Print results
std
::
cout
<<
std
::
setw
(
8
)
<<
n
<<
std
::
scientific
<<
std
::
setprecision
(
3
)
<<
std
::
setw
(
15
)
<<
timer
.
min
()
<<
std
::
setw
(
15
)
<<
timer_eff
.
min
()
<<
std
::
endl
;
}
/* SAM_LISTING_END_3 */
}
// Rename long variable name to duration_t (easy to change)
...
...
@@ -161,29 +145,28 @@ using duration_t = std::chrono::nanoseconds;
* \param[in] repeats Number of repetitions.
*/
template
<
class
Function
>
duration_t
timing
(
const
Function
&
F
,
int
repeats
=
10
)
{
// Shortcut for time_point
using
time_point_t
=
std
::
chrono
::
high_resolution_clock
::
time_point
;
duration_t
timing
(
const
Function
&
F
,
int
repeats
=
10
)
{
// Shortcut for time_point
using
time_point_t
=
std
::
chrono
::
high_resolution_clock
::
time_point
;
// Loop many times
duration_t
min_elapsed
;
for
(
int
r
=
0
;
r
<
repeats
;
r
++
)
{
// Start clock (MATLAB: tic)
time_point_t
start
=
std
::
chrono
::
high_resolution_clock
::
now
();
// Loop many times
duration_t
min_elapsed
;
for
(
int
r
=
0
;
r
<
repeats
;
r
++
)
{
// Start clock (MATLAB: tic)
time_point_t
start
=
std
::
chrono
::
high_resolution_clock
::
now
();
// Run function
F
();
// Run function
F
();
// Stop clock (MATLAB: toc) and measure difference
duration_t
elapsed
=
std
::
chrono
::
duration_cast
<
duration_t
>
(
std
::
chrono
::
high_resolution_clock
::
now
()
-
start
);
// Stop clock (MATLAB: toc) and measure difference
duration_t
elapsed
=
std
::
chrono
::
duration_cast
<
duration_t
>
(
std
::
chrono
::
high_resolution_clock
::
now
()
-
start
);
// Compute min between all runs
min_elapsed
=
r
==
0
?
elapsed
:
std
::
min
(
elapsed
,
min_elapsed
);
}
return
min_elapsed
;
// Compute min between all runs
min_elapsed
=
r
==
0
?
elapsed
:
std
::
min
(
elapsed
,
min_elapsed
);
}
return
min_elapsed
;
}
/* \brief Compute timing using chrono
...
...
@@ -191,72 +174,67 @@ duration_t timing(const Function & F, int repeats = 10) {
*/
void
runtime_arrow_matrix_with_chrono
()
{
// Table header
std
::
cout
<<
std
::
setw
(
8
)
<<
"n"
<<
std
::
scientific
<<
std
::
setprecision
(
3
)
<<
std
::
setw
(
15
)
<<
"original"
<<
std
::
setw
(
15
)
<<
"efficient"
// Table header
std
::
cout
<<
std
::
setw
(
8
)
<<
"n"
<<
std
::
scientific
<<
std
::
setprecision
(
3
)
<<
std
::
setw
(
15
)
<<
"original"
<<
std
::
setw
(
15
)
<<
"efficient"
<<
std
::
endl
;
// Run from $2^5$ to $2^11$ with powers of two
for
(
unsigned
int
n
=
(
1
<<
5
);
n
<
(
1
<<
12
);
n
=
n
<<
1
)
{
// Create random vectors
VectorXd
d
=
VectorXd
::
Random
(
n
);
VectorXd
a
=
VectorXd
::
Random
(
n
);
VectorXd
x
=
VectorXd
::
Random
(
n
);
VectorXd
y
(
n
);
// Call "timing", using a lambda function for F
// Remember: we cannot pass arrow\_matrix\_2\_times\_x directly to timing
// the timing function expects a n object with operator()(void)
duration_t
elapsed
=
timing
([
&
a
,
&
d
,
&
x
,
&
y
]()
{
arrow_matrix_2_times_x
(
d
,
a
,
x
,
y
);
},
10
);
// Call "timing", using a lambda function for F
duration_t
elapsed_efficient
=
timing
(
[
&
a
,
&
d
,
&
x
,
&
y
]()
{
efficient_arrow_matrix_2_times_x
(
d
,
a
,
x
,
y
);
},
10
);
// Output timings
std
::
cout
<<
std
::
setw
(
8
)
<<
n
<<
std
::
scientific
<<
std
::
setprecision
(
3
)
<<
std
::
setw
(
15
)
<<
elapsed
.
count
()
*
1e-9
// ns to s
<<
std
::
setw
(
15
)
<<
elapsed_efficient
.
count
()
*
1e-9
// ns to s
<<
std
::
endl
;
// Run from $2^5$ to $2^11$ with powers of two
for
(
unsigned
int
n
=
(
1
<<
5
);
n
<
(
1
<<
12
);
n
=
n
<<
1
)
{
// Create random vectors
VectorXd
d
=
VectorXd
::
Random
(
n
);
VectorXd
a
=
VectorXd
::
Random
(
n
);
VectorXd
x
=
VectorXd
::
Random
(
n
);
VectorXd
y
(
n
);
// Call "timing", using a lambda function for F
// Remember: we cannot pass arrow\_matrix\_2\_times\_x directly to timing
// the timing function expects a n object with operator()(void)
duration_t
elapsed
=
timing
([
&
a
,
&
d
,
&
x
,
&
y
]
()
{
arrow_matrix_2_times_x
(
d
,
a
,
x
,
y
);
},
10
);
// Call "timing", using a lambda function for F
duration_t
elapsed_efficient
=
timing
([
&
a
,
&
d
,
&
x
,
&
y
]
()
{
efficient_arrow_matrix_2_times_x
(
d
,
a
,
x
,
y
);
},
10
);
// Output timings
std
::
cout
<<
std
::
setw
(
8
)
<<
n
<<
std
::
scientific
<<
std
::
setprecision
(
3
)
<<
std
::
setw
(
15
)
<<
elapsed
.
count
()
*
1e-9
// ns to s
<<
std
::
setw
(
15
)
<<
elapsed_efficient
.
count
()
*
1e-9
// ns to s
<<
std
::
endl
;
}
}
}
int
main
(
void
)
{
// Test vectors
VectorXd
a
(
5
);
a
<<
1.
,
2.
,
3.
,
4.
,
5.
;
VectorXd
d
(
5
);
d
<<
1.
,
3.
,
4.
,
5.
,
6.
;
VectorXd
x
(
5
);
x
<<
-
5.
,
4.
,
6.
,
-
8.
,
5.
;
VectorXd
yi
;
// Run both functions
arrow_matrix_2_times_x
(
a
,
d
,
x
,
yi
);
VectorXd
ye
(
yi
.
size
());
efficient_arrow_matrix_2_times_x
(
a
,
d
,
x
,
ye
);
// Compute error
double
err
=
(
yi
-
ye
).
norm
();
// Output error
std
::
cout
<<
"--> Correctness test."
<<
std
::
endl
;
std
::
cout
<<
"Error: "
<<
err
<<
std
::
endl
;
// Print out runtime
std
::
cout
<<
"--> Runtime test."
<<
std
::
endl
;
runtime_arrow_matrix
();
// Print out runtime with chrono
std
::
cout
<<
"--> Runtime test. with chrono."
<<
std
::
endl
;
runtime_arrow_matrix_with_chrono
();
// Final test: exit with error if error is too big
double
eps
=
std
::
numeric_limits
<
double
>::
denorm_min
();
exit
(
err
<
eps
);
// Test vectors
VectorXd
a
(
5
);
a
<<
1.
,
2.
,
3.
,
4.
,
5.
;
VectorXd
d
(
5
);
d
<<
1.
,
3.
,
4.
,
5.
,
6.
;
VectorXd
x
(
5
);
x
<<
-
5.
,
4.
,
6.
,
-
8.
,
5.
;
VectorXd
yi
;
// Run both functions
arrow_matrix_2_times_x
(
a
,
d
,
x
,
yi
);
VectorXd
ye
(
yi
.
size
());
efficient_arrow_matrix_2_times_x
(
a
,
d
,
x
,
ye
);
// Compute error
double
err
=
(
yi
-
ye
).
norm
();
// Output error
std
::
cout
<<
"--> Correctness test."
<<
std
::
endl
;
std
::
cout
<<
"Error: "
<<
err
<<
std
::
endl
;
// Print out runtime
std
::
cout
<<
"--> Runtime test."
<<
std
::
endl
;
runtime_arrow_matrix
();
// Print out runtime with chrono
std
::
cout
<<
"--> Runtime test. with chrono."
<<
std
::
endl
;
runtime_arrow_matrix_with_chrono
();
// Final test: exit with error if error is too big
double
eps
=
std
::
numeric_limits
<
double
>::
denorm_min
();
exit
(
err
<
eps
);
}
Assignments/Codes/MatVec/ArrowMatrix/templates/arrowmatvec.cpp
View file @
b9817f75
#include
<iostream>
#include
<iomanip>
#include
<iostream>
#include
<Eigen/Dense>
...
...
@@ -17,20 +17,15 @@ using namespace Eigen;
/* SAM_LISTING_BEGIN_0 */
void
arrow_matrix_2_times_x
(
const
VectorXd
&
d
,
const
VectorXd
&
a
,
const
VectorXd
&
x
,
VectorXd
&
y
)
{
assert
(
d
.
size
()
==
a
.
size
()
&&
a
.
size
()
==
x
.
size
()
&&
"Vector size must be the same!"
);
int
n
=
d
.
size
();
VectorXd
d_head
=
d
.
head
(
n
-
1
);
VectorXd
a_head
=
a
.
head
(
n
-
1
);
MatrixXd
d_diag
=
d_head
.
asDiagonal
();
MatrixXd
A
(
n
,
n
);
A
<<
d_diag
,
a_head
,
a_head
.
transpose
(),
d
(
n
-
1
);
y
=
A
*
A
*
x
;
assert
(
d
.
size
()
==
a
.
size
()
&&
a
.
size
()
==
x
.
size
()
&&
"Vector size must be the same!"
);
int
n
=
d
.
size
();
VectorXd
d_head
=
d
.
head
(
n
-
1
);
VectorXd
a_head
=
a
.
head
(
n
-
1
);
MatrixXd
d_diag
=
d_head
.
asDiagonal
();
MatrixXd
A
(
n
,
n
);
A
<<
d_diag
,
a_head
,
a_head
.
transpose
(),
d
(
n
-
1
);
y
=
A
*
A
*
x
;
}
/* SAM_LISTING_END_0 */
...
...
@@ -42,15 +37,13 @@ void arrow_matrix_2_times_x(const VectorXd &d, const VectorXd &a,
* @param[out] y The vector y = A*A*x
*/
/* SAM_LISTING_BEGIN_1 */
void
efficient_arrow_matrix_2_times_x
(
const
VectorXd
&
d
,
const
VectorXd
&
a
,
const
VectorXd
&
x
,
VectorXd
&
y
)
{
assert
(
d
.
size
()
==
a
.
size
()
&&
a
.
size
()
==
x
.
size
()
&&
void
efficient_arrow_matrix_2_times_x
(
const
VectorXd
&
d
,
const
VectorXd
&
a
,
const
VectorXd
&
x
,
VectorXd
&
y
)
{
assert
(
d
.
size
()
==
a
.
size
()
&&
a
.
size
()
==
x
.
size
()
&&
"Vector size must be the same!"
);
int
n
=
d
.
size
();
int
n
=
d
.
size
();
// TODO: Implement an efficient version of arrow\_matrix\_2\_times\_x
// TODO: Implement an efficient version of arrow\_matrix\_2\_times\_x
}
/* SAM_LISTING_END_1 */
...
...
@@ -58,40 +51,38 @@ void efficient_arrow_matrix_2_times_x(const VectorXd &d,
* Repeat tests 10 times, and output the minimal runtime
* amongst all times. Test both the inefficient and the efficient
* versions.
*/
*/
void
runtime_arrow_matrix
()
{
// TODO: your code here, time the codes
// TODO: your code here, time the codes
}
int
main
(
void
)
{
// Test vectors
VectorXd
a
(
5
);
a
<<
1.
,
2.
,
3.
,
4.
,
5.
;
VectorXd
d
(
5
);
d
<<
1.
,
3.
,
4.
,
5.
,
6.
;
VectorXd
x
(
5
);
x
<<
-
5.
,
4.
,
6.
,
-
8.
,
5.
;
VectorXd
yi
;
// Run both functions
arrow_matrix_2_times_x
(
a
,
d
,
x
,
yi
);
VectorXd
ye
(
yi
.
size
());
efficient_arrow_matrix_2_times_x
(
a
,
d
,
x
,
ye
);
// Compute error
double
err
=
(
yi
-
ye
).
norm
();
// Output error
std
::
cout
<<
"--> Correctness test."
<<
std
::
endl
;
std
::
cout
<<
"Error: "
<<
err
<<
std
::
endl
;
// Print out runtime
std
::
cout
<<
"--> Runtime test."
<<
std
::
endl
;
runtime_arrow_matrix
();
// Final test: exit with error if error is too big
double
eps
=
std
::
numeric_limits
<
double
>::
denorm_min
();
exit
(
err
<
eps
);
// Test vectors
VectorXd
a
(
5
);
a
<<
1.
,
2.
,
3.
,
4.
,
5.
;
VectorXd
d
(
5
);
d
<<
1.
,
3.
,
4.
,
5.
,
6.
;
VectorXd
x
(
5
);
x
<<
-
5.
,
4.
,
6.
,
-
8.
,
5.
;
VectorXd
yi
;
// Run both functions
arrow_matrix_2_times_x
(
a
,
d
,
x
,
yi
);
VectorXd
ye
(
yi
.
size
());
efficient_arrow_matrix_2_times_x
(
a
,
d
,
x
,
ye
);
// Compute error
double
err
=
(
yi
-
ye
).
norm
();
// Output error
std
::
cout
<<
"--> Correctness test."
<<
std
::
endl
;
std
::
cout
<<
"Error: "
<<
err
<<
std
::
endl
;
// Print out runtime
std
::
cout
<<
"--> Runtime test."
<<
std
::
endl
;
runtime_arrow_matrix
();
// Final test: exit with error if error is too big
double
eps
=
std
::
numeric_limits
<
double
>::
denorm_min
();
exit
(
err
<
eps
);
}
Assignments/Codes/MatVec/GramSchmidt/solutions/gramschmidt.cpp
View file @
b9817f75
...
...
@@ -11,55 +11,46 @@ using namespace Eigen;
* \return Matrix with ONB of $span(a_1, \cdots, a_n)$ as columns
*/
/* SAM_LISTING_BEGIN_1 */
MatrixXd
gram_schmidt
(
const
MatrixXd
&
A
)
{
// We create a matrix Q with the same size and data of A
MatrixXd
Q
(
A
);
// The first vector just gets normalized
Q
.
col
(
0
).
normalize
();
// Iterate over all other columns of A
for
(
unsigned
int
j
=
1
;
j
<
A
.
cols
();
++
j
)
{
// See eigen documentation for usage of col and leftCols
Q
.
col
(
j
)
-=
Q
.
leftCols
(
j
)
*
(
Q
.
leftCols
(
j
).
transpose
()
*
A
.
col
(
j
));
// Normalize vector, if possible
// (otherwise it means columns of $\mathbf{A}$ are
// almost linear dependant)
double
eps
=
std
::
numeric_limits
<
double
>::
denorm_min
();
if
(
Q
.
col
(
j
).
norm
()
<=
eps
*
A
.
col
(
j
).
norm
()
)
{
std
::
cerr
<<
"Gram-Schmidt failed because "
<<
"A has (almost) linear dependant "
<<
"columns. Bye."
<<
std
::
endl
;
break
;
}
else
{
Q
.
col
(
j
).
normalize
();
}
MatrixXd
gram_schmidt
(
const
MatrixXd
&
A
)
{
// We create a matrix Q with the same size and data of A
MatrixXd
Q
(
A
);
// The first vector just gets normalized
Q
.
col
(
0
).
normalize
();
// Iterate over all other columns of A
for
(
unsigned
int
j
=
1
;
j
<
A
.
cols
();
++
j
)
{
// See eigen documentation for usage of col and leftCols
Q
.
col
(
j
)
-=
Q
.
leftCols
(
j
)
*
(
Q
.
leftCols
(
j
).
transpose
()
*
A
.
col
(
j
));
// Normalize vector, if possible
// (otherwise it means columns of $\mathbf{A}$ are
// almost linearly dependant)
double
eps
=
std
::
numeric_limits
<
double
>::
denorm_min
();