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

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 have attached the matlab script, fortran code and the make file. Matlab version 2010a, ifort 12.0.0, mkl 10.3.

thanks

Reddy

Link Copied

6 Replies

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

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

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

ax

^{2}+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

^{2})/b => Eqn1

y = 0.75*x+0.25*y

diffX = abs(y-x)

Residue = ay

^{2}+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

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

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

^{2}+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

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

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

- 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

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.

Topic Options

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