       program temp
c      Matrix square root
       integer i,j,k,m,n,isym,niter,ngauss,ierr,index(100,2)
       real*8 pi,anorm,bnorm,r,a(100,100),b(100,100),stiff(100,100),
     & ident(100,100),matrixsqroot(100,100),xi(200),weight(200),
     & matrixsqrootsq(100,100)
       pi=4.0d0*datan(1.0d0)
       open(unit=5,file='matrixsqroot.dat',status='old')
       read(5,*) n,ngauss,r
       do 10 i=1,n
         read(5,*) (a(i,j),j=1,n)
10     continue
        isym=0
        bnorm=0.0d0
        do 4 i=1,n
         do 2 j=1,n
          bnorm=bnorm+dabs(a(i,j)-a(j,i))
2        continue
4       continue
        if(bnorm.ne.0.0d0)isym=1
       anorm=0.0d0
       if(isym.eq.0)then
        do 6 j=1,n
         do 5 i=1,j
          anorm=anorm+dabs(a(i,j))
5        continue
6       continue
       else
        do 8 j=1,n
         do 7 i=1,n
          anorm=anorm+dabs(a(i,j))
7        continue
8       continue
       endif
       ident(1:n,1:n)=0.0d0
       do 20 i=1,n
        ident(i,i)=1.0d0
20     continue
       niter=10
       k=0
       if(isym.eq.0)then
       do 26 j=1,n
        do 24 i=1,j
         k=k+1
         index(k,1)=i
         index(k,2)=j
24      continue
26     continue
       else
       do 260 j=1,n
        do 240 i=1,n
         k=k+1
         index(k,1)=i
         index(k,2)=j
240      continue
260     continue
       endif
       matrixsqroot(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))**(2.0d0*r)*ident(1:n,1:n)+
     & (1.0d0-xi(i))**2*a(1:n,1:n)
       b(1:n,1:n)=2.0d0*(1.0d0+xi(i)+r*(1.0d0-xi(i)))/pi*
     & (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
       matrixsqroot(1:n,1:n)=matrixsqroot(1:n,1:n)+weight(i)*b(1:n,1:n)
28    continue
      if(isym.eq.0)then
      matrixsqroot(1:n,1:n)=(matrixsqroot(1:n,1:n)+
     &  transpose(matrixsqroot(1:n,1:n)))/2.0d0
      endif
      print*, 'Matrix square root: Initial guess'
      print*, ''
      do 29 i=1,n
       print*, (matrixsqroot(i,j),j=1,n)
29    continue
       print*,' '
       if(isym.eq.0)then
       do 500 m=1,niter
       do 100 i=1,n*(n+1)/2
        do 90 j=1,n*(n+1)/2
         if(index(j,1).eq.index(j,2))then
       stiff(i,j)=0.5d0*(matrixsqroot(index(i,1),index(j,1))*
     & ident(index(i,2),index(j,2))+matrixsqroot(index(i,1),index(j,2)) 
     & *ident(index(i,2),index(j,1))+matrixsqroot(index(i,2),index(j,1))
     & *ident(index(i,1),index(j,2))+matrixsqroot(index(i,2),index(j,2))
     & *ident(index(i,1),index(j,1)))
         else
         stiff(i,j)=matrixsqroot(index(i,1),index(j,1))*
     & ident(index(i,2),index(j,2))+matrixsqroot(index(i,1),index(j,2))  
     & *ident(index(i,2),index(j,1))+matrixsqroot(index(i,2),index(j,1))
     & *ident(index(i,1),index(j,2))+matrixsqroot(index(i,2),index(j,2))
     & *ident(index(i,1),index(j,1))
         endif
90      continue
100    continue
       do 140 i=1,n
        do 130 j=1,n
         matrixsqrootsq(i,j)=0.0d0
         do 120 k=1,n
          matrixsqrootsq(i,j)=matrixsqrootsq(i,j)+
     &    matrixsqroot(i,k)*matrixsqroot(k,j)
120      continue
130     continue
140    continue
       do 150 i=1,n*(n+1)/2
        b(i,1)=a(index(i,1),index(i,2))-
     &  matrixsqrootsq(index(i,1),index(i,2))
150    continue
        bnorm=0.0d0
        do 190 i=1,n*(n+1)/2
         bnorm=bnorm+dabs(b(i,1))
190     continue
         bnorm=bnorm/anorm
         print*,'bnorm ',bnorm
         if(bnorm.le.1.0d-14)go to 510
       call gaussj(stiff,n*(n+1)/2,100,b,1,100,ierr)
        if(ierr.eq.1)then
      print*,'Matrix is singular: Newton-Raphson refinement not possible
     & '
         stop
        endif
        do 160 i=1,n*(n+1)/2
         matrixsqroot(index(i,1),index(i,2))=
     &   matrixsqroot(index(i,1),index(i,2))+b(i,1)
160     continue
        do 180 j=1,n
         do 170 i=j+1,n
          matrixsqroot(i,j)=matrixsqroot(j,i)
170      continue
180     continue
500    continue
       print*,'Failed to converge '
       stop
       else
       do 5000 m=1,niter
       do 1000 i=1,n**2
        do 900 j=1,n**2
         stiff(i,j)=matrixsqroot(index(i,1),index(j,1))  
     & *ident(index(i,2),index(j,2))+matrixsqroot(index(j,2),index(i,2))
     & *ident(index(i,1),index(j,1))
900     continue
1000   continue
       do 1400 i=1,n
        do 1300 j=1,n
         matrixsqrootsq(i,j)=0.0d0
         do 1200 k=1,n
          matrixsqrootsq(i,j)=matrixsqrootsq(i,j)+
     &    matrixsqroot(i,k)*matrixsqroot(k,j)
1200     continue
1300    continue
1400   continue
       do 1500 i=1,n**2
        b(i,1)=a(index(i,1),index(i,2))-
     &  matrixsqrootsq(index(i,1),index(i,2))
1500    continue
         bnorm=0.0d0
        do 1900 i=1,n**2
         bnorm=bnorm+dabs(b(i,1))
1900    continue
        bnorm=bnorm/anorm
         print*,'bnorm ',bnorm
         if(bnorm.le.1.0d-14)go to 510
       call gaussj(stiff,n**2,100,b,1,100,ierr)
        if(ierr.eq.1)then
      print*,'Matrix is singular: Newton-Raphson refinement not possible
     & '
         stop
        endif
        do 1600 i=1,n**2
         matrixsqroot(index(i,1),index(i,2))=
     &   matrixsqroot(index(i,1),index(i,2))+b(i,1)
1600     continue
5000   continue
       print*,'Failed to converge '
       stop
       endif
510   print*, 'Matrix square root: Final result'
      print*, ''
      do 200 i=1,n
       print*, (matrixsqroot(i,j),j=1,n)
200   continue
      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
