

c      Dear Physics 210 students,
c      Here is a fortran code that calls the library routine `invertr'
c      to invert a test matrix.  
c      This routine destroys your input matrix A, so use a copy
c      of it instead of your real matrix, if you want to save it.
c      N is the declared dimension of the array holding your matrix.
c      Np is actual dimension of the matrix you want to invert.
c      idet controls whether the you calculate the determinant only (idet=0)
c      the inverse only (idet=1) or both (idet=2)
c      Richard Scalettar
 
        implicit none
        integer Np
        parameter (Np=64)
        real*8 A(Np,Np),Ainv(Np,Np),det
        real*8 Acopy(Np,Np),ident(Np,Np)
        integer i,j,k,N,idet,iran,ipvt(Np)
 
        N=4
        A(1,1)=3.4d0
        A(1,2)=2.4d0
        A(1,3)=1.4d0
        A(1,4)=0.4d0
        A(2,1)=3.3d0
        A(2,2)=5.3d0
        A(2,3)=3.3d0
        A(2,4)=1.3d0
        A(3,1)=0.7d0
        A(3,2)=0.7d0
        A(3,3)=6.7d0
        A(3,4)=0.7d0
        A(4,1)=2.1d0
        A(4,2)=1.1d0
        A(4,3)=2.1d0
        A(4,4)=4.1d0

        do 110 i=1,N
        do 100 j=1,N
            Acopy(i,j)=A(i,j)
100     continue
110     continue

	     idet=2
        call invertr(Acopy,Ainv,det,Np,N,idet,ipvt)   
 
        write (67,*) '  '
        write (67,*) ' matrix inverse is '
        write (67,*) '  '
        do 210 i=1,N
        do 200 j=1,N
             write (67,990) i,j,Ainv(i,j)
200     continue
210     continue
990     format(2i5,f12.6)

c       TEST

        write (67,*) '  '
        write (67,*) ' A*Ainv is '
        write (67,*) '  '
        do 320 i=1,N
        do 310 j=1,N
            ident(i,j)=0.d0
            do 300 k=1,N
                ident(i,j)=ident(i,j)+A(i,k)*Ainv(k,j)
300         continue
            write (67,990) i,j,ident(i,j)
310     continue
320     continue
 
        end
 

        subroutine invertr(b,c,d,dim,toff,idet,ipvt)   
                   
c       This subroutine finds the inverse of the matrix b
c       using Gaussian elimination with partial pivoting.
c       b is the input matrix.
c       c is the inverse of b.
c       d is the determinent of b.
c       toff is the size of the matrix
c       If idet=0 the inverse of b is returned but not the determinent.
c       If idet=1 the determinent is returned,but not the inverse.
c       If idet=2 both the inverse and the determinent are returned.
c                                     
        implicit none
        integer dim,toff,ipvt(dim) 
        real*8 b(dim,dim),c(dim,dim),d,tb            
        integer iwn,nm1,k,kp1,mb,i,j,idet,kb,km1
 
        ipvt(toff)=1
       iwn=1
       nm1=toff-1
c                                   
c       Finding the pivot
c
       do 35 k=1,nm1
       kp1=k+1
       mb=k
       do 15 i=kp1,toff
       if (dabs(b(i,k)).gt.dabs(b(mb,k))) mb=i
15       continue
       ipvt(k)=mb                                          
       tb=b(mb,k)
       if (mb.ne.k) then
       ipvt(toff)=-ipvt(toff)
       b(mb,k)=b(k,k)
       b(k,k)=tb
       endif
       if (dabs(tb).lt.1.d-10) then
       iwn=0
       go to 35
       endif
c
c       Computing multipliers
c
       do 20 i=kp1,toff
20       b(i,k)=-b(i,k)/tb
c
c       Interchange and elimination by columns
c
       do 30 j=kp1,toff
       tb=b(mb,j)
       if (mb.ne.k) then
       b(mb,j)=b(k,j)
        b(k,j)=tb
       endif
       do 25 i=kp1,toff
25       b(i,j)=b(i,j)+b(i,k)*tb
cray       call saxpy(toff-k,tb,b(kp1,k),1,b(kp1,j),1)
30       continue
35       continue
       if (iwn.eq.0) then
       write (6,*) 'Warning: Pivot=0'
       d=0.d0    
       go to 150
       endif
       if (idet.gt.0) then
c      
c       Compute the determinent
c
       d=ipvt(toff)
       do 40 i=1,toff
40       d=d*b(i,i)
       endif
       if (idet.ne.1) then
c
c       Compute the inverse
c       Ax=b goes to A~x=Lb~
c
       do 50 i=1,toff
       do 45 j=1,toff
45       c(i,j)=0.d0
50       c(i,i)=1.d0
       do 100 j=1,toff
       do 60 k=1,nm1
       kp1=k+1
       mb=ipvt(k)
       tb=c(mb,j)
       if (mb.ne.k) then
       c(mb,j)=c(k,j)
       c(k,j)=tb
       endif
       do 55 i=kp1,toff
55       c(i,j)=c(i,j)+b(i,k)*tb
cray       call saxpy(toff-k,tb,b(kp1,k),1,c(kp1,j),1)
60       continue
c
c       Inverting A~
c
       do 70 kb=1,nm1
       km1=toff-kb
       k=km1+1
       c(k,j)=c(k,j)/b(k,k)
       tb=-c(k,j)
       do 65 i=1,km1    
65       c(i,j)=c(i,j)+b(i,k)*tb
cray       call saxpy(km1,tb,b(1,k),1,c(1,j),1)
70       continue
       c(1,j)=c(1,j)/b(1,1)
100       continue
       endif
150       continue
       return
       end
     
