For my master's thesis I implemented multiple solvers for structured generalized eigenvalue problems using LAPACK. In this post, I will briefly discuss a method to simplify the memory management and ways to catch programming errors as early as possible when implementing a solver for linear algebra problems that uses LAPACK. The advice in this post only supplements good programming practices like using version control systems and automated tests.

# Type Coercions and Floating-Point Types

Consider a language where we only allow automatic, implicit type conversions (coercion) among numeric types if every value of the source type can be represented as a value of the target type, i.e., there is no truncation and no round-off. Let us call this kind of coercion *value-preserving coercion*. With such a strict coercion rule, an unsigned integer with four bytes can be coerced into a signed integer with eight bytes but the compiler (interpreter) will not coerce signed integers to unsigned integers. In this post, we highlight a reason why coercions from single-precision to double-precision floating-point types may be undesirable although they are value preserving.

# Building NumPy and SciPy with Intel Compilers and Intel MKL

NumPy and SciPy rely on BLAS and LAPACK for basic linear algebra functionality like matrix-vector multiplication, linear system solves, or routines for eigenvalue computation. The Intel Math Kernel Library (MKL) is a mathematics library providing amongst other things fast and multithreaded implementations of BLAS and LAPACK. In this blog post, I describe how to compile NumPy and SciPy with the Intel compilers using Intel MKL on Linux. Note NumPy and SciPy can be linked to MKL without the Intel compilers by providing the proper linker options to, e.g., GCC, and I will briefly explain this as well.

Continue reading Building NumPy and SciPy with Intel Compilers and Intel MKL

# Talk: Master's Thesis

Yesterday I gave a talk on the preliminary results of my master's thesis on *Projection Methods for Generalized Eigenvalue Problems. *Below is the announcement for the talk; it mentions only complex (Hermitian) matrices but all results are immediately applicable to real matrices. The slides are here (PDF).

## Projection Methods for Generalized Eigenvalue Problems

In this talk, we will discuss backward errors, dense solvers, and a new projection method for generalized eigenvalue problems (GEPs) with Hermitian positive semidefinite (HPSD) matrices.

We will present properties and origins of such GEPs. We will also address quickly computable and structure preserving backward error bounds for these kinds of GEPs. There is an abundance of literature on backward error measures possessing one of these features but only recently, the author came across a backward error providing both.

We will elaborate on dense solvers for GEPs with HPSD matrices. The standard solver for GEPs with Hermitian matrices is fast but requires positive definite mass matrices and is only conditionally backward stable; the QZ algorithm for general GEPs is backward stable but it is also magnitudes slower and does not preserve any problem properties. In the talk, we will present two new backward stable and structure preserving solvers, one using deflation, the other one using the generalized singular value decomposition (GSVD). In comparison to the QZ algorithm, both solvers are competitive with the standard solver in our tests.

Finally, we will touch on a new solver for large, sparse GEPs.

# Matlab: Existing Files and mkdir

In Matlab, a directory called `foo`

can be created by calling `mkdir('foo')`

. If `mkdir`

is successful, then the function returns logical 1; otherwise it returns logical 0. Now consider the case where Matlab is run on Linux and there is already a file called `foo`

in the current directory. Linux does not permit the existence of a file and a directory of the same name in the same folder so the call to `mkdir`

should fail. Instead, `mkdir`

will leave the file intact, return logical 1, and print the misleading warning "Directory already exists". That is, there is no *directory* called `foo`

after this supposedly successful function call.

The peculiarity of this behavior has been reported to Mathworks in September 2015. The example above was tested with Matlab R2015b and R2014a on Linux.

# Implicit Type Conversions Gone Wrong, C++ Edition

Consider the following setup:

#include <cstddef> using namespace std; // ... size_t n = 100; size_t* array = new size_t[n];

The goal is to set every element of the array to the fixed value `x`

.

In C++, an idiomatic way to solve this problem are the functions `std::fill`

(link to documentation) and `std::fill_n`

(link to documentation) from the STL header `algorithm`

. In addition to the fixed value `x`

, `std::fill`

requires a pair of output iterators specifying a half-open range as its arguments whereas the additional argument to `std::fill_n`

is an output iterator pointing to the first element and the *length* of the container. That is, we could assign to the elements of the array above as follows:

#include <cstddef> #include <algorithm> using namespace std; // ... size_t n = 100; size_t* array = new size_t[n]; fill(array, array+n, x); fill_n(array, n, x);

What made me write this blog post is the following line of code:

fill_n(array, array+n, x);

Here, I accidentally used `std::fill_n`

in place of `std::fill`

but the code still compiles because of the automatic, implicit type conversion (type coercion) from `std::size_t*`

to `std::size_t`

. With warnings enabled (`-Wextra -Wall -std=c++11 -pedantic`

) neither g++ 4.8.5 nor clang++ 3.5 warn about this line and yet this piece of code causes a segmentation fault on my computer whenever it is executed.

# More Secure Wireless Telephony with VoIP over Wi-Fi

In this post, I will explain why you should start to phase out DECT and how you can do this with VoIP (Voice over IP, also known as Internet telephony) over Wi-Fi. This guide allows you to attach several wireless devices to your landline and it permits calls among all telephones in your household. Moreover, the guide provides encrypted wireless communication, it uses off-the-shelf consumer products, and works worldwide. On the downside, the voice quality may be worse unless your router supports WMM (Wireless Multimedia Extensions) and you may need a WLAN repeater because DECT tends to have higher range than WLAN. (Surprisingly, the electronics markets in my vicinity offer dual-band WLAN repeaters with 600 MBit/s transfer rate fitting into a wall socket costing no more than the cheapest DECT repeater.) If so desired, you can continue to use your existing DECT phones but communication with these devices is not guaranteed to be encrypted.

You can follow this guide if you have

- an account with an Internet telephony provider,
- your Internet router is a Wi-Fi hotspot,
- the Internet router manages your landline,
- the Internet router allows wireless VoIP telephony using SIP, and
- a mobile device with Wi-Fi as well Android, iOS, or Windows.

Continue reading More Secure Wireless Telephony with VoIP over Wi-Fi

# Matlab: Unawareness of Partial Path Names may cause Bugs

In a filesystem, there are two kinds of paths: relative paths and absolute paths. Absolute paths always point to the same target while the target of a relative path depends on the current working directory. In Matlab, paths are called "path names" and there is a third kind of path name: the partial path name. It is defined as follows:

A partial path name is the last portion of a full path name for a location on the MATLAB search path.

(The Matlab search path is similar to the environment variable PATH and is used to locate Matlab files.) In consequence of this definition, a path name may be treated by Matlab as a partial path name even if the user is unaware of them.

In this post, I will demonstrate Matlab behavior that is unexpected if you are not aware of the existence of partial path names. I will also explain how this causes a bug in xUnit 3 and xUnit 4.

Continue reading Matlab: Unawareness of Partial Path Names may cause Bugs

# Evaluating Approximate Solutions using the Backward Error

Given a linear algebra problem and an approximate solution, can we say something about the quality of the given approximate solution without knowing an exact solution?

Continue reading Evaluating Approximate Solutions using the Backward Error

# Windows XP and AHCI

Currently, AHCI is the standard interface between a host and SATA storage devices. Windows XP does not come with a built-in AHCI driver which may prevent it from detecting a computer's storage device. I describe a common solution to this problem and how it worked out for me.