home *** CD-ROM | disk | FTP | other *** search
/ CP/M / CPM_CDROM.iso / simtel / cpmug / cpmug049.ark / MINV.FOR < prev    next >
Encoding:
Text File  |  1984-04-29  |  2.4 KB  |  119 lines
  1.       programminvt
  2.       integer*1l,n
  3.       dimensiona(5,5),l(5)
  4.       integer*1char
  5.       reallog,log10
  6.       datan/5/
  7.       datakb/3/,kons/3/
  8.       datarpd/.01745329/,dpr/57.29578/
  9.       datapi/3.141593/
  10.       do 23000i=1,n
  11.       do 23002j=1,n
  12.       a(j,i)=1./real(i+j-1)
  13. 23002 continue
  14. 23003 continue
  15. 23000 continue
  16. 23001 continue
  17.       write(kons,20)a
  18. 20    format(1x,5g12.6)
  19.       callminv(a,n,l,det)
  20.       write(kons,20)det
  21.       write(kons,20)a
  22.       callminv(a,n,l,det)
  23.       write(kons,20)det
  24.       write(kons,20)a
  25.       stop
  26.       end
  27.       subroutineminv(a,n,l,det)
  28.       dimensiona(n,n),l(n)
  29.       integer*1l,j,k,n,it,i
  30.       datakb/3/,kons/3/
  31.       datarpd/.01745329/,dpr/57.29578/
  32.       datapi/3.141593/
  33.       nd=n
  34.       do 23004j=1,n
  35.       l(j)=j
  36. 23004 continue
  37. 23005 continue
  38.       det=1.
  39.       do 23006k=1,n
  40.       if(.not.(fmxmgv(a(k,k),n,n-k+1,jr).eq.0.))goto 23008
  41.       write(kons,900)
  42.       goto 23007
  43. 23008 continue
  44. 900   format(1x,15hSingular Matrix)
  45.       if(.not.(jr.ne.1))goto 23010
  46.       j=k-1+jr
  47.       it=l(k)
  48.       l(k)=l(j)
  49.       l(j)=it
  50.       callfvswap(a(1,k),1,a(1,j),1,nd)
  51.       det=-det
  52. 23010 continue
  53.       q=a(k,k)
  54.       det=q*det
  55.       a(k,k)=1.
  56.       do 23012i=1,n
  57.       a(i,k)=a(i,k)/q
  58. 23012 continue
  59. 23013 continue
  60.       continue
  61.       j=1
  62. 23014 if(.not.(j.le.n))goto 23016
  63.       if(.not.(j.eq.k))goto 23017
  64.       goto 23015
  65. 23017 continue
  66.       q=-a(k,j)
  67.       a(k,j)=0.
  68.       do 23019i=1,n
  69.       a(i,j)=a(i,k)*q+a(i,j)
  70. 23019 continue
  71. 23020 continue
  72. 23015 j=j+1
  73.       goto 23014
  74. 23016 continue
  75. 23006 continue
  76. 23007 continue
  77.       n1=n-1
  78.       do 23021k=1,n1
  79.       continue
  80.       i=k
  81. 23023 if(.not.(i.le.n.and.l(i).ne.k))goto 23025
  82. 23024 i=i+1
  83.       goto 23023
  84. 23025 continue
  85.       callfvswap(a(k,1),nd,a(i,1),nd,nd)
  86.       l(i)=l(k)
  87. 23021 continue
  88. 23022 continue
  89.       return
  90.       end
  91.       functionfmxmgv(a,ia,l,k)
  92.       dimensiona(1)
  93.       fmxmgv=0.
  94.       is=1
  95.       do 23026i=1,l
  96.       absa=abs(a(is))
  97.       if(.not.(absa.gt.fmxmgv))goto 23028
  98.       fmxmgv=absa
  99.       k=i
  100. 23028 continue
  101.       is=is+ia
  102. 23026 continue
  103. 23027 continue
  104.       return
  105.       end
  106.       subroutinefvswap(a,ia,b,ib,n)
  107.       dimensiona(1),b(1)
  108.       isa=1
  109.       isb=1
  110.       do 23030i=1,n
  111.       t=a(isa)
  112.       a(isa)=b(isb)
  113.       b(isb)=t
  114.       isa=isa+ia
  115.       isb=isb+ib
  116. 23030 continue
  117. 23031 continue
  118.       return
  119.       end