----------------------------------------------------------------------------------
@MSGID:
<c8ac7613-a6bc-4804-a668-40dbc9de519en@googlegroups.com> 06fed85f
@REPLY:
<8a0d90c2-1912-41f9-8cf1-76c8f12698adn@googlegroups.com> 74627e8d
@REPLYADDR Farzad Tatar
<tatarfarzad1992@gmail.com>
@REPLYTO 2:5075/128 Farzad Tatar
@CHRS: CP866 2
@RFC: 1 0
@RFC-References:
<aa2d38fc-1df2-4703-b1b3-67adf45e2793n@googlegroups.com> <8a0d90c2-1912-41f9-8cf1-76c8f12698adn@googlegroups.com>
@RFC-Message-ID:
<c8ac7613-a6bc-4804-a668-40dbc9de519en@googlegroups.com>
@TZUTC: -0700
@PID: G2/1.0
@TID: FIDOGATE-5.12-ge4e8b94
> I do not quite understand the results you post: the result for
Intel has 10 numbers, for MicroSoft powerstation there are 12 numbers - >
neither set of numbers fit in a square matrix ;). So, I have no idea if
the two are really close as you say.
Sorry for the misleading information. The numbers are the first few
elements of the first column just to be able to compare them. the matrix
shape is 1400by 1400. for this case.
> You should beware that the inverse of a matrix can be difficult
to compute accurately, though these seem to be small. You can check
the results by multiplying the original matrix with the inverse matrix
(or even by multiplying the inverse matrix by the original). The answer
should be the identity matrix or nearly so. Differences may occur because
of rounding effects.
I will try this though I already know that the result is not
accurate. I actually tried REAL128 from the intrinsic kind function in the
modern fortran and the results got closer to the one that I was receiving
from the old code. It means the reason might be due to the precision.
But aren`t Real*8 and Real64 the same?
> It would help if you could show your code.
This subroutine is a very small part of my project. Here is the
part I am dealing with.
SUBROUTINE RBF inv(iflag)
USE COMMON 1
USE, INTRINSIC :: ISO FORTRAN ENV, DP => REAL64
INTEGER, INTENT(IN) :: iflag
INTEGER :: ii, jj, i1, j1, INFO, NRHS, M, N, ising, lwork
INTEGER, DIMENSION(NNODE LS) :: IPIV
REAL (DP):: xi , yi , xj , yj , dist , temp, determ
REAL (DP):: F intp
REAL (DP):: xni, yni, xnj, ynj,xw ls norm
REAL(DP), ALLOCATABLE:: work(:)
! Determine the coeff of the Radial Basis Functions
! This approach is used to interpolate the function everywhere in the domain
! 1. Invert RBF LS only at the 1st iteration (iflag==0) , then
I will have the RBF coeffs
IF (iflag == 0) THEN
WRITE(*,*) ` Solve global RBF iter`, iter, Nnode LS !nnode ls is
the number of node in each row.
! Matrix of the RBF function values
DO ii = 1, Nnode LS
xi = xy lsgrid(1, ii)
yi = xy lsgrid(2, ii)
DO jj = 1, Nnode LS
xj = xy lsgrid(1, jj)
yj = xy lsgrid(2, jj)
dist = SQRT((xj - xi )**2.0 dp + (yj - yi )**2.0 dp)
RBF LS(ii, jj) = EXP(-((dist * epsw)**2.0 dp))
!IF (RBF LS(II,JJ) <= 1.0E-58 DP) RBF LS(II,JJ)=0.0 DP
END DO
END DO
! Solve the linear system
NRHS=1
print*, "1:", rbf ls(1:10,1)
! Perform LU decomposition of RBF LS
CALL DGETRF(NNODE LS, NNODE LS, RBF LS, NNODE LS, ipiv, info)
print*, "2:composed LU", rbf ls(1:10,1)
IF (info == 0) THEN
! RBF LS was successfully factorized, now compute its inverse
ALLOCATE(work(lwork)) !this is work vector
CALL DGETRI(NNODE LS, RBF LS, NNODE LS, ipiv, work, nnode ls, info)
print*, "3:", rbf ls(1:10,1)
IF (info == 0) THEN
! A inv contains the inverse of A
! Use A inv as needed
ELSE
PRINT *, "Error computing the inverse of A."
END IF
DEALLOCATE(work)
ELSE
PRINT *, "Error in LU decomposition of A."
END IF
> One last remark: MicroSoft Powerstation has been out of support
for several decades.
The point is that this old method is giving a more accurate result
in comparison with the new one.
Thank you,
Farzad
--- G2/1.0
* Origin: usenet.network (2:5075/128)
SEEN-BY: 5001/100 5005/49 5015/255 5019/40 5020/715
848 1042 4441 12000
SEEN-BY: 5030/49 1081 5058/104 5075/128
@PATH: 5075/128 5020/1042 4441