MKL API for Condition Number "cgecon" Not working correctly

[rohitspandey]

HI,

We are using cgecon api to compute the condition large number of matrices. It seems there is difference in results computed for certain matrices. we have verified results for large set of matrices 10^5 ranging from 2X2 to 10X10 sizes.  The calculation of condition number for about 10%  are wrongly calculated.

We have verifed the results against matlab aswell. Here is an example result compared against matlab.

% input data

N = 2;

H =[ 0.7989 - 0.2170i -2.1232 + 0.8383i;

0.0180 + 1.2296i 0.2920 + 0.8230i; ];

% norm infinite

matlabNormInf = norm(H, Inf)

% = 3.1105

% This is the expected condition number

matlabCond = 1/cond(H,Inf)

% = 0.3582

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Run DLL (complex I/F), with breakpoint before cgecon()

% We use clange() to compute the norm (result in aNorm)

% We use cgetrf() to compute the LU decomposition (result in A)

% When Breakpoint hits on cgecon(), the input is:

% ntran 73 'I' const char

% n 2 int

% +A[0] 0.017999999+i*1.2296000 std::complex

% +A[1] -0.16693313-i*0.65216720 std::complex

% +A[2] 0.29200000+i*0.82300001 std::complex

% +A[3] -2.6111891+i*1.1661187 std::complex

% lda 2 const int

% aNorm 3.1105480 float

% cond -1.0737418e+008 float

% cwork [4](0.0,0.0,0.0,0.0) std::vector<std::complex,std::allocator<std::complex> >

% fwork [4](3.1105480,2.1029973,0.0,0.0) std::vector<float,std::allocator>

% info 0 int

% After stepping over: cgecon(&ntran, &n, A, &lda, &aNorm, &cond, &cwork[0], &fwork[0], &info);

% The only values that changed are:

% cond 0.54947048 float

% cwork [4](-0.041364070+i*0.18391991,0.88580656-i*0.30841386,0.039692443+i*0.23203264,-0.31928882+i*0.14258970)

% fwork [4](0.81910032,0.00000000,0.00000000,1.1150000)

% The condition number is wrong, it should be 0.3582 (computed in matlab as matlabCond)

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Check that the input to cgecon() was correct

 % clange() reported aNorm=3.1105480 which is correct (Matlab matlabNormInf=3.1105)

% cgetrf() reported A (the H=PLU decomposition from cgetrf() ) which is also correct:

A=[ 0.017999999+j*1.2296000;

-0.16693313-j*0.65216720;

0.29200000+j*0.82300001;

-2.6111891+j*1.1661187; ];

P = [0 1; 1 0];

L = [ 1 0;

A(2) 1; ];

U = [ A(1) A(3);

0 A(4); ];

PLU = P*L*U

% PLU = 0.7989 - 0.2170i -2.1232 + 0.8383i

% 0.0180 + 1.2296i 0.2920 + 0.8230i

% A is correct because H=P*L*U

 % The expected condition number using LU instead of H is of course identical

matlabCond = 1/cond(P*L*U,Inf)

% = 0.3582

Kindly verify the observation and let us know whats wrong with the api

Regards
Rohit

[mecej4]

I find your report of error(s) in MKL routine cgecon next to impossible to follow.

Please report the problem using a C/Fortran source for an example case where the routine works incorrectly.

Providing snippets of Matlab scripts and then saying that MKL gave errors is not helpful. Considering that some versions of Matlab call an MKL shared library to perform matrix computations, I suspect that there are problems in how you set up the arguments to cgecon, how you call the routine or how you interpret the results.

[rohitspandey]

I managed to reproduce the problem on a simpler example, where imaginary is always 0:(Now not refering any thing from matlab)

INPUT:

H =
 1 0.9

0.1 -1

EXPECTED RESULT:

NormInf = 1.9

L = 1 0

0.1 1

U = 1 0.9

0 -1.09

ConditionNumber(Inf) = 0.3019

ACTUAL RESULT:

NormInf returned as expected from clange()

L U returned as expected from cgetrf()

condition number is wrong (0.46514937 instead of 0.3019)

 CALL TO cgecon():

input (all correct):

n 2 int

lda 2 const int

aNorm 1.9000000 float

ntran 73 'I' const char

A is decalred as MKL_Complex8*, but we display it as float* :

[0] 1.0000000 float

[1] 0.00000000 float

[2] 0.10000000 float

[3] 0.00000000 float

[4] 0.89999998 float

[5] 0.00000000 float

[6] -1.0900000 float

[7] 0.00000000 float

call:

cgecon(&ntran, &n, A, &lda, &aNorm, &cond, &cwork[0], &fwork[0], &info);

output:

cond 0.46514937 float

cwork [4](0.73394495,2.6605504,0.73394495,2.6605504) std::vector<:COMPLEX>,std::allocator<:COMPLEX> > >

fwork [4](0.10000000,0.00000000,0.00000000,0.89999998) std::vector >

info 0 int

 Can you please check why the condition number is wrong?

#r
Rohit

[Ying_H_Intel]

Hi Rohit, I saw the discussion in http://software.intel.com/en-us/forums/topic/308534, everything looks ok, just one problem, About the matrix order, you may have know, that default lapack routine is FORTRAN interface, the default matrix order is column_major. and if you call the routine in C language, you may need to take care of it. For example, I change the LU matrix (order) A[0] = MKL_Complex8(1.0f, 0.0f); // U(1,1) A[2] = MKL_Complex8(0.1f, 0.0f); // L(2,1) A[1] = MKL_Complex8(0.9f, 0.0f); // U(1,2) A[3] = MKL_Complex8(-1.09f, 0.0f); // U(2,2) const MKL_INT lda = 2; float aNorm = 1.9f; Then i will get expected rcond value: 0.3019. Or in MKL 10.3, we also introduce c interface to LAPACK, which named lapacke_cgecon, you may call LAPACKE_cgecon(LAPACK_ROW_MAJOR,ntran , n, &A[0], lda, aNorm, &rcond ); get same result. Would you please try them and let me know if any problem. Best Regards, Ying H.

[rohitspandey]

Hi Ying, The matrix A was already in column-major order (it was the output of cgetrf). When you tried a different order, you in fact put it in row-major order. Could it be that internally, cgecon transposes the input matrix (like many MKL functions do) based on a non-exposed flag? Important points are: A given input matrix will always give the same output result (As expected) Over many random input matrixes (all with condition number <100), sometimes it works, sometimes it doesn’t. I put the C code with some input that make it FAIL, using clange cgetrf cgecon (but the problem is only with cgecon). Can you advise on why cgecon seem to transpose the input matrix “sometimes”? Or how would you modify the call to cgecon? Thanks PS Here is the C Code void tmpDebug(void) { /*********************************************/ // Set input matrix in COLUMN-major order: // H = 1 0.9 // 0.1 -1 MKL_Complex8 A[4] = { MKL_Complex8(1.0f, 0), MKL_Complex8(0.1f, 0), MKL_Complex8(0.9f, 0), MKL_Complex8(-1.0f, 0)}; /*********************************************/ //Calculate the aNorm Float_t aNorm; { const char ntran = 'I'; //Inf Norm MKL_INT m = 2; MKL_INT n = 2; MKL_INT lda = 2; // COLUMN-major order => lda = num rows float fwork[4]; // temporary buffer (reserved larger than needed) aNorm = clange(&ntran, &m, &n, A, &lda, &fwork[0]); } //Here aNorm = 1.9 (this is correct) /*********************************************/ //Factorise the matrix A=PLU { MKL_INT m = 2; MKL_INT n = 2; MKL_INT lda = 2; // COLUMN-major order => lda = num rows MKL_INT iwork[4]; // temporary buffer (reserved larger than needed) MKL_INT info; cgetrf(&m, &n, A, &lda, &iwork[0], &info); } //Here A[] = {1.0, 0.1, 0.9, -1.09} (this is correct, in COLUMN-major order) //EXPECTED RESULT: // L = 1 0 // 0.1 1 // // U = 1 0.9 // 0 -1.09 // P = identity matrix so we can ignore it // // So L and U "packed" are the matrix // A = 1 0.9 // 0.1 -1.09 // So everything is correct up to here /*********************************************/ //Calcualte the Condition Float_t cond; { const char ntran = 'I'; //Inf Norm MKL_INT n = 2; MKL_INT lda = 2; // COLUMN-major order => lda = num rows float fwork[4]; // temporary buffer (reserved larger than needed) MKL_Complex8 cwork[4]; // temporary buffer (reserved larger than needed) MKL_INT info; cgecon(&ntran, &n, A, &lda, &aNorm, &cond, &cwork[0], &fwork[0], &info); } // Here cond = 0.46514937 (but correct value is 0.3019 // As pointed out by Ying H (Intel), changing A to ROW-major order would make it work }

[Alexander_K_Intel3]

Hi Rohit,

This is normal behavior. The function does estimation of condition number but not exact computation. The estimate is never less than the true value and usually within 10x of true values.

The value of the function to get approximate understanding(order of the condition number) in a few flops and not destroying the original matrix. The function require O(n²) operations whereas exact computation via GETRI requires O(n³).

Details about the functionality could be found in N. J. HIGHAM, FORTRAN codes for estimating the one-norm of a real or complex matrix, with applications to condition estimation, ACM Trans. Math. Softw., 14 (1988), pp. 381-396.

And the exact number is the same which MatLab gives. It could be checked by adding following code in the end of the reproducer: cgetri( &n, A, &lda, ipiv, cwork, &n, &info ); printf("cond A = %f\n",1.0/aNorm/clange("I",&n,&n,A,&lda, fwork));

W.B.R., Alexander

[rohitspandey]

Hi Alexander, Thanks for correcting us. I missed the http://software.intel.com/sites/products/documentation/hpc/mkl/mklman/GUID-63F50607-0569-435D-9FB5-1E3006A91768.htm (?gecon) Application Notes The computed rcond is never less than r (the reciprocal of the true condition number) and in practice is nearly always less than 10r. A call to this routine involves solving a number of systems of linear equations A*x = b or AH*x = b; the number is usually 4 or 5 and never more than 11. Each solution requires approximately 2*n2 floating-point operations for real flavors and 8*n2 for complex flavors. But is there any API that support accurate condtion number calculation. Other wise I we have to write our own. Regards Rohit

[Ying_H_Intel]

Hi Rohit, As Alexander mentioned, the below routines can calculation the accurate condition number. cgetrf( &m, &n, A, &lda, ipiv, &info ); cgetri( &n, A, &lda, ipiv, cwork, &n, &info ); printf("cond A = %f\n",1.0/aNorm/clange("I",&n,&n,A,&lda, fwork)); Best Regards, Ying