MKL results do not match with MATLAB results

[Dharma]

Hello,
I am trying to solve an equation of this form using mkl.

• This algorithm solves the equation
• betaTrans = transpose(conjg(beta));
• Imat is an Identity matrix
• zplus = 1i*eta;
• g = inv( (E+1i*eta)*Imat - alpha - betaTrans*g*beta);
• The above equation can alos be written as
• A*X*B*X + C*X + I = 0
• where the unknown X = g
• A = betaTrans
• B = beta
• C = - ((E+1i*eta)*Imat - alpha)
• I = Imat
• alpha, beta, and Imat are square matrices of size N=NW*NU

I find that for a speific case of E=0.0, mkl results donot match with matlab results. You will see that the results match upto certain iteration number with in the numerical error. Can anyone please help me fix the issue.
I have attached the matlab script, fortran code and the make file. Matlab version 2010a, ifort 12.0.0, mkl 10.3.

thanks
Reddy

[mecej4]

Your description shows N^2 quadratic equations in N^2 variables. For such problems, there is no unique solution. Even for the case N=2, the solution would be a conic section rather than a point. In such a situation, why do you expect two solutions, one from Matlab and another from a different program, to agree?

Only after this question has been answered satisfactorily can one start talking about algorithms and their behavior.

[Dharma]

True that the system may not necessarily have a unique solution but my understanding is that if i start a iterative procdure (algorithm) with a given intial condiction, the convergance properties should be similar where the algorithm is implimented in fortran or matlab or any othere langauge (assuming the implimentation is correct). Ofcourse if the system (i.e., the equations) itself is chiotic where the solution is extermly sensitive to intial conditon then i would not expect similar results and this system for sure is not chioatic. If you consider the case of N=1, the equations will look like

ax²+bx+c=0 and the algorithm that i use (fixed point iteration)

xinit = -c/b
y = xinit
ic=0
do while( Residue < 1E-5 )
ic=ic+1
x = (-c - ay²)/b => Eqn1
y = 0.75*x+0.25*y
diffX = abs(y-x)
Residue = ay²+by+c
print*,ic,Residue,diffX

end do

Now, in the above algorithm, i would expect the solutions to have a similar convergence path if the intial guess is same even though the quadratic has two possible solutions(which may be same or different).

In the above algorithm the only critical thing which would differe in MATLAB vs fortran implimentation is solving the linear equation (Eqn1). I use a direct solution (LU decompositn and then solve). The iteration path ( values of Residue and diffX) match (upto 3rd or 4th digit) between matlab and fortran for a while and then suddenly fortran iteration takes differen path. Now this happens for a very specific case and for all other cases that I have tested the results seem to match.

Thanks
Reddy

[Dharma]

Hello,
Looks like the problem is due to matrix being ill conditoned. I am doing the following iteration
The system is a N=8. A,B,C and X are all N by N matrices.
do while(fixedpoint)
D = (A*Xold+B)
solve: DX = -C ( i tried both gesv and gesvx)
Residue = clange('F',N,N,DX+C,work)
! inverse Condition number of D goes to 1E-12 as the simulation progresses
! gesvx does imporve the conditon number to 1E-7
! Residue is not going below 1E-2
! where as for the same problem matlab is able to solve the system some how
! and the residue is near 1E-6 as the solution converges.
X = 0.75*X+0.25*Xold
Xold = X
Fnorm = clange('F',N,N,AX²+BX+C,work)
fixedpoint = (Fnorm > 1E-6)
end do

Can you suggest a fix. Also any idea how is matlab able to handle ill conditoned matrix. All of my claculations are being done in single precesion.

Thanks
Reddy

[timintel]

Quoting dharmareddy84@gmail.com

Also any idea how is matlab able to handle ill conditoned matrix. All of my claculations are being done in single precesion.

I thought Matlab supported only double precision.

[Dharma]

Sorry, i should have been clear. All of my calclulations in fortran are being done in single precesion. Yes Matlab seem to do it in double precesion. And i do see the issue of ill conditoned matrix in matrix but only when the inverse of condition number goes less than 1E-16. Even then in matlab the residue of the equation seem to get as less as 1E-6 for the same problem.

[mecej4]

The residual || A x - b || depends on the scaling of the variables, so the value that you quoted for it (1E-6) does not have much significance.

If, indeed, you have a problem + algorithm combination with a condition number of 1E-16, you are liable to lose 16 decimal digits in calculating your results. Neither single-precision (23 bits, or 6 to 7 decimal digits) nor double precision (53 bits, or 15 to 16 decimal digits) will suffice to give meaningful results after a loss of 16 decimal digits.