Search code examples
c++lapackequation-solving

Solving equation using LU factorization without pivoting (Lapack library)


At the beginning I would like to apologise for my English. Now, let's go to my problem.

I try to write a simple code which will find a solution of a system of linear equations:

Ax = b

where A is a square matrix nxn. In this program I use Lapack library (I MUST use LU factorization WITHOUT pivoting).

I found few examples, e. g.: Understanding LAPACK calls in C++ with a simple example where we can see how to use functions: dgetrf_ and dgetrs_. But even if I copy this code (from the best answer) to my program, it return sometimes correct results (e. g. A and b are the same as in the best answer) and sometimes wrong results (e. g. A = {1, -3, 1, -1}, b = {3, 5}, right answer is: {6, 1} and function return {-4, 7}). For greater matrices it return wrong results. Can anyone say why?

In this site: https://software.intel.com/sites/products/documentation/doclib/mkl_sa/11/mkl_lapack_examples/lapackro1.htm is written:

LAPACK routines assume input matrices do not contain IEEE 754 special values such as INF or NaN. Using these special values may cause LAPACK to return unexpected results or become unstable.

I guess that INF means "infinity" and NaN means "Not A Number", right?

And the second problem is that if even above example would work properly, it uses LU factorization WITH partial pivoting. I need functions from Lapack library, which do LU factorization WITHOUT pivoting. I did research for this function, but I didn't find anything. Is there anyone who know what is (or maybe are) this (these) function (functions)? I lose hope for the solution to this problem.


Solution

  • The LAPACK routines are written in FORTRAN and data are stored column major.. You are solving the transposed A matrix system A^T x = b. Try using A = {1, 1, -3, -1}.

    You are correct, INF means "infinity" and NaN means "Not A Number".

    The LU algorithm always uses pivoting. The Cholesky decomposition does not use pivoting (dpotrf,dpotrs). But then your matrix must be "symmetric positive definite".