@MSGID:
<26454e02-d247-4d86-a721-6166cdf11d56n@googlegroups.com> 0d4e1842
@REPLY:
<c8ac7613-a6bc-4804-a668-40dbc9de519en@googlegroups.com> 06fed85f
@REPLYADDR jfh <harperjf2@gmail.com>
@REPLYTO 2:5075/128 jfh
@CHRS: CP866 2
@RFC: 1 0
@RFC-References:
<aa2d38fc-1df2-4703-b1b3-67adf45e2793n@googlegroups.com> <8a0d90c2-1912-41f9-8cf1-76c8f12698adn@googlegroups.com>
<c8ac7613-a6bc-4804-a668-40dbc9de519en@googlegroups.com>
@RFC-Message-ID:
<26454e02-d247-4d86-a721-6166cdf11d56n@googlegroups.com>
@TZUTC: -0700
@PID: G2/1.0
@TID: FIDOGATE-5.12-ge4e8b94
On Monday, September 4, 2023 at 8:34:45 AM UTC+12, Farzad Tatar wrote:
 > > 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. 
>
 How do you know the old method is more accurate? I have heard of
cases where people thought a result was right because it agreed with a
previous one, but both turned out to be wrong. We have not been shown what
either matmul(A inv ,A) or matmul(A,A inv) is. Nor have we been told how
ill-conditioned A is. I recall years ago that one of my students thought of a
way to do a problem in quantum mechanics that involved solving a linear
system of the form A*x=y. His pure mathematics was right, but he had
never learnt about ill-conditioned matrices. I said that the elements of y
were obtained from physical experiments and were known only to within
about 1 part in 10**6 and so I asked him to try changing some of them
by that amount and see what happened to x. That showed him why his
method was not in general use for his problem.
--- 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