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.
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.
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
ax2+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 - ay2)/b => Eqn1 y = 0.75*x+0.25*y diffX = abs(y-x) Residue = ay2+by+c print*,ic,Residue,diffX
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.
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,AX2+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.
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.
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.