- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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

Link Copied

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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>

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

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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^{2}) operations whereas exact computation via GETRI requires O(n^{3}).

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

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page