home *** CD-ROM | disk | FTP | other *** search
/ CP/M / CPM_CDROM.iso / simtel / cpmug / cpmug049.ark / MINV.RAT < prev    next >
Encoding:
Text File  |  1984-04-29  |  2.2 KB  |  97 lines
  1. #Test driver for MINV
  2.     PROGRAM MINVT
  3.     INTEGER*1 L,N
  4.     DIMENSION A(5,5),L(5)
  5.     include F77DEF
  6.     DATA N/5/
  7.     include CONSTS
  8.     DO I=1,N
  9.     $(DO J=1,N
  10.     A(J,I)=1./REAL(I+J-1)$)#Set up arbitrary matrix
  11.     write(kons,20) A
  12. 20    format(1x,5G12.6)
  13.     call MINV(A,N,L,DET)
  14.     write(kons,20)DET
  15.     write(kons,20)A
  16. #Now try to get original back
  17.     call MINV(A,N,L,DET)
  18.     WRITE(KONS,20)DET
  19.     WRITE(KONS,20)A
  20.     STOP
  21.     END
  22. #Gauss-Jordan inversion & determinant calculation
  23. SUBROUTINE MINV(A,N,L,DET)
  24. DIMENSION A(N,N),L(N)
  25. #Use FORTRAN (col,row) subscript convention
  26. INTEGER*1 L,J,K,N,IT,I
  27. include CONSTS
  28. ND=N#Convert to *2 for calls
  29. #Tag original column
  30. DO J=1,N
  31.     L(J)=J
  32. DET=1.
  33. DO K=1,N;$(
  34. #Choose pivot row for maximum value and find if 0.
  35. #Start at diagonal and look downwards
  36.   IF(FMXMGV(A(K,K),N,N-K+1,JR)==0.)
  37.     $(WRITE(KONS,900);BREAK$)
  38.     900 FORMAT(1x,'Singular Matrix')
  39.   IF(JR^=1)
  40.     $(J=K-1+JR#Convert from J counting from diagonal
  41. #to J subscript
  42. #Swap rows K & J and tags
  43.     IT=L(K)
  44.     L(K)=L(J)
  45.     L(J)=IT
  46.     CALL FVSWAP(A(1,K),1,A(1,J),1,ND)
  47. #Swap changes determinant sign
  48.     DET=-DET$)#End of swapping block
  49. #Determinant calculation
  50.   Q=A(K,K)
  51.   DET=Q*DET
  52.   A(K,K)=1.#Start changing this part of I matrix to inverse
  53. # Divide row by diagonal, including both modified
  54. # original and developing inverse
  55.   DO I=1,N
  56.     A(I,K)=A(I,K)/Q
  57. #Eliminate off-diagonal elements of col K
  58.   FOR(J=1;J<=N;J=J+1)$(IF(J==K)NEXT
  59.     Q=-A(K,J)
  60. #Bring in off-diagonal element of I matrix
  61.     A(K,J)=0.
  62.     DO I=1,N;A(I,J)=A(I,K)*Q+A(I,J)$)#End FOR
  63. $)#End DO K=1,N
  64. #Undo column swaps
  65. N1=N-1
  66. DO K=1,N1;$(
  67. #Search for K
  68.     FOR(I=K;I<=N&L(I)^=K;I=I+1);#No body of FOR
  69.     CALL FVSWAP(A(K,1),ND,A(I,1),ND,ND)
  70.     L(I)=L(K)$)#Keep track of columns yet to go
  71. RETURN
  72. END
  73. FUNCTION FMXMGV(A,IA,L,K)
  74. #Start at A, at intervals IA find largest of L numbers
  75. #and return position in K
  76. DIMENSION A(1)#Actually (IA,L)
  77. FMXMGV=0.
  78. IS=1
  79. DO I=1,L
  80.   $(ABSA=ABS(A(IS))
  81.   IF(ABSA>FMXMGV)
  82.     $(FMXMGV=ABSA;K=I$)
  83.   IS=IS+IA$)
  84. RETURN
  85. END
  86. SUBROUTINE FVSWAP(A,IA,B,IB,N)
  87. DIMENSION A(1),B(1)#Actually (IA,N),IB,N)
  88. #Swap REAL vector starting at A interval IA
  89. #with B interval IB, length N
  90. ISA=1
  91. ISB=1
  92. DO I=1,N
  93.     $(T=A(ISA);A(ISA)=B(ISB);B(ISB)=T
  94.     ISA=ISA+IA;ISB=ISB+IB$)
  95. RETURN
  96. END
  97.