Nп/п : 49 из 100
 От   : Farzad Tatar                        2:5075/128        03 сен 23 13:34:41
 К    : Arjen Markus                                          03 сен 23 23:36:02
 Тема : Re: Inversion problem with minimum accuracy in Dgetri in LAPack
----------------------------------------------------------------------------------
                                                                                 
@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



   GoldED+ VK   │                                                 │   09:55:30    
                                                                                
В этой области больше нет сообщений.

Остаться здесь
Перейти к списку сообщений
Перейти к списку эх