2

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.

Community
  • 1
  • 1
Michal
  • 43
  • 5

2 Answers2

2

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".

ztik
  • 3,482
  • 16
  • 34
  • Thank you very much for your answer. I want to click an up arrow, but I can't because I have too low reputation score. You are right, when I transpose matrix A function returns right values. Are you sure that LU algorithm always uses pivoting? If yes, maybe I must use LU factorization with partial pivoting. – Michal Feb 02 '15 at 19:19
  • If you have unsymmetric matrices you need to use LU with pivoting. You can not upvote an answer when you have low reputation. You can only accepted as the correct. – ztik Feb 03 '15 at 09:43
  • 2
    Just to provide some more detail, pivoting is required to guarantee that an LU factorization exists for a non-singular matrix A. Consider A = ((0, 1), (1, 0)). There do not exist non-singular triangular L and U such that A = LU. It's necessary to first apply a row-interchange ("pivot") to get a factorization of the form PA = LU, after which we have L = U = identity. – Stephen Canon Feb 11 '15 at 02:04
-1

Isn't MATLAB calculates LU without the use of Pivoting. And from research I found that MATLAB uses LAPACK routines to do the calculation. I wonder how they do it.

HKK
  • 21
  • 7