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,std::allocator > >

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

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

}

Hi Alexander,

Thanks for correcting us. I missed the
http://software.intel.com/sites/products/documentation/hpc/mkl/mklman/GU...

(?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