如何使用DLSACG解线性方程AX=B？？

xiaomelody · 发表于 2015-1-10 01:06:12

解线性方程 Ax = b con A(n,n), n=40;  A, b, x 为复数
A(k,j) = 1+(k+j)/(kj+1) + (|k-j|/(k+j+1))i
k,j 属于 1:n

data.txt
b(1:n)的系数为:

b( 1) = (  44.54490425, 250.18068298) b( 2) = (  27.00252589, 201.18580799)
b( 3) = (  24.92829934, 181.70600784) b( 4) = (  26.13000225, 170.32385033)
b( 5) = (  28.11808687, 162.60101547) b( 6) = (  30.29764162, 156.93859099)
b( 7) = (  32.37274813, 152.51444544) b( 8) = (  34.34191794, 148.91242944)
b( 9) = (  36.08469956, 145.99615297) b(10) = (  37.67551424, 143.58667430)

b(11) = (  39.03974183, 141.56238561) b(12) = (  40.22426534, 139.84721623)
b(13) = (  41.21625847, 138.39768230) b(14) = (  42.06822441, 137.19325965)
b(15) = (  42.75210007, 136.18257629) b(16) = (  43.32238389, 135.33706972)
b(17) = (  43.77547980, 134.62809246) b(18) = (  44.09373685, 134.01838944)
b(19) = (  44.31534306, 133.53859998) b(20) = (  44.44533734, 133.13095692)

b(21) = (  44.46167540, 132.79173557) b(22) = (  44.37426847, 132.54243979)
b(23) = (  44.19829310, 132.37563081) b(24) = (  43.94053061, 132.26491425)
b(25) = (  43.63025583, 132.21516615) b(26) = (  43.24852817, 132.20963631)
b(27) = (  42.82035002, 132.25508895) b(28) = (  42.34420273, 132.34102784)
b(29) = (  41.81558928, 132.46138741) b(30) = (  41.22640369, 132.63698669)

b(31) = (  40.58568110, 132.85259849) b(32) = (  39.89902576, 133.11778084)
b(33) = (  39.15822301, 133.42519531) b(34) = (  38.39083704, 133.77522902)
b(35) = (  37.57507213, 134.16202444) b(36) = (  36.71999067, 134.57527665)
b(37) = (  35.83927537, 135.02432074) b(38) = (  34.92706608, 135.50543584)
b(39) = (  33.99396423, 136.00449838) b(40) = (  33.02425789, 136.52679242)

  result.txt：

x( 1) = (          ,          ) x( 2) = (          ,          )
x( 3) = (          ,          ) x( 4) = (          ,          )
x( 5) = (          ,          ) x( 6) = (          ,          )
x( 7) = (          ,          ) x( 8) = (          ,          )
x( 9) = (          ,          ) x(10) = (          ,          )

x(11) = (          ,          ) x(12) = (          ,          )
x(13) = (          ,          ) x(14) = (          ,          )
x(15) = (          ,          ) x(16) = (          ,          )
x(17) = (          ,          ) x(18) = (          ,          )
x(19) = (          ,          ) x(20) = (          ,          )

x(21) = (          ,          ) x(22) = (          ,          )
x(23) = (          ,          ) x(24) = (          ,          )
x(25) = (          ,          ) x(26) = (          ,          )
x(27) = (          ,          ) x(28) = (          ,          )
x(29) = (          ,          ) x(30) = (          ,          )

x(31) = (          ,          ) x(32) = (          ,          )
x(33) = (          ,          ) x(34) = (          ,          )
x(35) = (          ,          ) x(36) = (          ,          )
x(37) = (          ,          ) x(38) = (          ,          )
x(39) = (          ,          ) x(40) = (          ,          )

[Fortran] 纯文本查看 复制代码

PROGRAM sislincom
IMPLICIT NONE

PARAMETER (IPATH=1,LDA=40,n=40)
INTEGER,PARAMETER::ne=40
INTEGER,PARAMETER::nsal=40  
COMPLEX a(LDA,LDA),b(n),x(n)

OPEN(nsal,FILE='RESULT.SAL')
WRITE (nsal,'(A/35("=")/)')
CALL erset (0, 0, 0) 

OPEN(ne,FILE='DATO.TXT')
WRITE (ne,'(A/35("=")/)')

!LECTURA DE DATOS

i=o
DO
 i=i+1
 READ(ne,'(6X,4F8.1)',END=40) datos(i,1:n)
ENDDO


CONTAINS

  SUBROUTINE sistema
    INTEGER code, tipo
    INTEGER, EXTERNAL:: iercd, n1rty

    CALL dlsacg (n, a, LDA, b, IPATH, x)
    code = iercd()
    tipo = n1rty(1)

    IF (code == 0) THEN
      WRITE (nsal,'(/A,3F10.6)') 'Sistema compatible: x = ', x
    ELSE
      WRITE (nsal,'(/A,I2)') 'Error en la resolucion code = ', code
      IF (code == 2) WRITE (nsal,'(A)') 'Matriz singular'
      IF (tipo == 4) WRITE (nsal,'(A)') 'Error fatal'
    ENDIF
  ENDSUBROUTINE sistema

ENDPROGRAM sislincom