program matrixproot
c      Matrix p'th root
c      *****************************************************************
c      n: Order of the matrix
c      ngauss: order of Gauss quadrature
c      r: User specified real number greater than or equal to 1. 
c      For well-conditioned matrices $r=1$ works well. For 
c      ill-conditioned matrices, try higher values of $r$ to find which 
c      value of $r$ gives the lowest error norm for a given number of 
c      Gauss points.
c      p: Order of the root
c      A(i,j): Input matrix (a(i,j) in the program). 
c      program finds the principal $p$'th root $A^(1/p)$. For principal
c      root to exist, no eigenvalue  of $A$ should lie in the interval 
c      $(-\infty,0]$.
c      Redirect output to file: Example  matrixproot.e >> matrixproot.out
c      *****************************************************************
       integer i,j,k,m,n,p,ngauss,ierr
       real*8 pi,anorm,bnorm,r,sum,a(100,100),b(100,100),stiff(100,100),
     & ident(100,100),matrixproot(100,100),xi(200),weight(200),
     & matrixprootp(100,100),matrixtemp(100,100)
       pi=4.0d0*datan(1.0d0)
       open(unit=5,file='matrixproot.dat',status='old')
       read(5,*) n,ngauss,r,p
       do 10 i=1,n
         read(5,*) (a(i,j),j=1,n)
10     continue
       anorm=0.0d0
        do 8 j=1,n
         do 7 i=1,n
          anorm=anorm+dabs(a(i,j))
7        continue
8       continue
       ident(1:n,1:n)=0.0d0
       do 20 i=1,n
        ident(i,i)=1.0d0
20     continue
       matrixproot(1:n,1:n)=0.0d0
       call gauss(ngauss,xi,weight)
       do 28 i=1,ngauss
        stiff(1:n,1:n)=(1.0d0+xi(i))**(p*r)*ident(1:n,1:n)+
     & (1.0d0-xi(i))**p*a(1:n,1:n)
       b(1:n,1:n)=p*dsin(pi/p)*(1.0d0+xi(i)+r*(1.0d0-xi(i)))/pi*
     & (1-xi(i))**(p-2.0d0)*(1.0d0+xi(i))**(r-1.0d0)*a(1:n,1:n)
      call gaussj(stiff,n,100,b,n,100,ierr)
        if(ierr.eq.1)then
         print*,'Matrix is singular '
         stop
        endif
       matrixproot(1:n,1:n)=matrixproot(1:n,1:n)+weight(i)*b(1:n,1:n)
28    continue
      print*, 'Matrix', p, ' th root'
      print*, ''
      do 29 i=1,n
       print*, (matrixproot(i,j),j=1,n)
29    continue
      matrixprootp(1:n,1:n)=matrixproot(1:n,1:n)
      do 150 m=1,p-1
       do 140 i=1,n
        do 130 j=1,n
         sum=0.0d0
         do 120 k=1,n
          sum=sum+matrixprootp(i,k)*matrixproot(k,j)
120      continue
         matrixtemp(i,j)=sum
130     continue
140    continue
        matrixprootp(1:n,1:n)=matrixtemp(1:n,1:n)
150    continue
       bnorm=0.0d0
        do 80 j=1,n
         do 70 i=1,n
          bnorm=bnorm+dabs(a(i,j)-matrixprootp(i,j))
70       continue
80      continue
       bnorm=bnorm/anorm
       print*,' '
       print*,'Error norm ',bnorm
       stop
       end
      subroutine gaussj(a,n,np,b,m,mp,ierr)

c  Purpose: Solution of the system of linear equations AX = B by
c     Gauss-Jordan elimination, where A is a matrix of order N and B is
c     an N x M matrix.  On output A is replaced by its matrix inverse
c     and B is preplaced by the corresponding set of solution vectors.

c  Source: W.H. Press et al, "Numerical Recipes," 1989, p. 28.

c  Modifications: 
c     1. Double  precision.
c     2. Error parameter IERR included.  0 = no error. 1 = singular 
c        matrix encountered; no inverse is returned.

c  Prepared by J. Applequist, 8/17/91.

      integer i,j,k,l,ll,m,mp,n,np,nmax,ipiv,indxr,indxc,irow,icol,ierr
      real*8 a,b,big,dum,pivinv

c        Set largest anticipated value of N.

      parameter (nmax=500)
      dimension a(np,np),b(np,mp),ipiv(nmax),indxr(nmax),indxc(nmax)
      ierr=0
      irow=0
      icol=0
      do 11 j=1,n
      ipiv(j)=0
 11   continue
      do 22 i=1,n
      big=0.d0
      do 13 j=1,n
      if (ipiv(j).ne.1) then
      do 12 k=1,n
      if (ipiv(k).eq.0) then
      if (dabs(a(j,k)).ge.big) then
      big=dabs(a(j,k))
      irow=j
      icol=k
      endif
      else if (ipiv(k).gt.1) then
      ierr=1
      return
      endif
 12   continue
      endif
 13   continue
      ipiv(icol)=ipiv(icol)+1
      if (irow.ne.icol) then
      do 14 l=1,n
      dum=a(irow,l)
      a(irow,l)=a(icol,l)
      a(icol,l)=dum
 14   continue
      do 15 l=1,m
      dum=b(irow,l)
      b(irow,l)=b(icol,l)
      b(icol,l)=dum
 15   continue
      endif
      indxr(i)=irow
      indxc(i)=icol
      if (a(icol,icol).eq.0.d0) then
      ierr=1
      return
      endif
      pivinv=1.d0/a(icol,icol)
      a(icol,icol)=1.d0
      do 16 l=1,n
      a(icol,l)=a(icol,l)*pivinv
 16   continue
      do 17 l=1,m
      b(icol,l)=b(icol,l)*pivinv
 17   continue
      do 21 ll=1,n
      if (ll.ne.icol) then
      dum=a(ll,icol)
      a(ll,icol)=0.d0
      do 18 l=1,n
      a(ll,l)=a(ll,l)-a(icol,l)*dum
 18   continue
      do 19 l=1,m
      b(ll,l)=b(ll,l)-b(icol,l)*dum
 19   continue
      endif
 21   continue
 22   continue
      do 24 l=n,1,-1
      if (indxr(l).ne.indxc(l)) then
      do 23 k=1,n
      dum=a(k,indxr(l))
      a(k,indxr(l))=a(k,indxc(l))
      a(k,indxc(l))=dum
 23   continue
      endif
 24   continue
      return
      end
        subroutine gauss(n,x,w)
        integer i,j,m,n
        real*8 x(200), w(200), eps,p1,p2,p3,pp,xl,xm,z,z1,w1(200),
     &         x1(200),pi
        pi=datan(1.0d0)*4.0d0
        eps=1.d-15
        m=(n+1)/2
        xm=0.0d0
        xl=1.0d0
        do 12 i=1,m
         z=dcos(pi*(i-0.25d0)/(n+0.5d0))
1        continue
          p1=1.0d0
          p2=0.0d0
          do 11 j=1,n
           p3=p2
           p2=p1
           p1=((2.0d0*j-1.0d0)*z*p2-(j-1.0d0)*p3)/j
11        continue
          pp=float(n)*(z*p1-p2)/(z*z-1.0d0)
          z1=z
          z=z1-p1/pp
          if(abs(z-z1).gt.eps)go to 1
          x(i)=xm-xl*z
          x(n+1-i)=xm+xl*z
          w(i)=2.0d0*xl/((1.0d0-z*z)*pp*pp)
          w(n+1-i)=w(i)
12        continue
         do 13 i=1,n
         x1(i)=x(i)
         w1(i)=w(i)
13       continue
        return
        end