C     metodo de Cholesky. Retorna TRUE quando matriz e def pos
      subroutine sCholesky(N,pM,pRes)
      integer N
      double precision pM(N,N)
      logical pRes
C     *****
      integer i,j,k
      double precision vL(N,N),vA(N,N),vAux
      logical vTeste

      pRes = .true.
      vTeste = .true.

      call sCopiaMatriz(N,N,pM,vA)
      call sInicializaMatriz(N,N,vL)

      i = 1
      do while(vTeste)

        if(vA(i,i).lt.0) then
          vTeste = .false.
          pRes = .false.
        else
          vL(i,i) = sqrt(vA(i,i))

          do j=i+1,N
            vL(j,i) = vA(j,i)/vL(i,i)

            do k=i+1,j
              vAux = vL(j,i)*vL(k,i)
              vA(j,k) = vA(j,k) - vAux
            end do

          end do

          if(i.ge.N) then
            vTeste = .false.
          else
            i = i+1
          end if

        end if

      end do

      end


C     ******************************************************************
C     ******************************************************************

C     rotina que arredonda um real em n casa decimais
      subroutine sArredondaNum(N,pNum,pRes)
      integer N
      double precision pNum, pRes
C     *****
      integer i,vNumInt, vExp

      vExp = 1
      do i=1,N
        vExp = vExp*10
      end do

      pRes = pNum*vExp
      vNumInt = nint(pRes)
      pRes = vNumInt/vExp

      end

C     ******************************************************************
C     ******************************************************************

C     rotina que arredonda os elem de um vetor em n casas decimais
      subroutine sArredondaVetor(N,M,pV,pRes)
      integer i,N,M
      double precision pV(M),pRes(M)

      do i=1,M
        call sArredondaNum(N,pV(i),pRes(i))
      end do

      end


C     ******************************************************************
C     ******************************************************************

C     rotina que identifica se um vetor e nulo
      subroutine sVetorENulo(N,pVetor,pErr,pRes)
      integer i,N
      double precision pVetor(N),pErr
      logical pRes

      pRes = .true.

      do i=1,N
        if( abs(pVetor(i)).gt.pErr) then
          pRes = .false.
        end if
      end do

      end

C     ******************************************************************
C     ******************************************************************

C      rotina para inicializar uma matriz
       subroutine sInicializaMatriz(N,M,pRes)
       integer i,j,N,M
       double precision pRes(N,M)

       do i=1,N
         do j=1,M
           pRes(i,j) = 0.0d0
         end do
       end do

       end

C     ******************************************************************
C     ******************************************************************

C      rotina para inicializar um vetor
       subroutine sInicializaVetor(N,pRes)
       integer i,N
       double precision pRes(N)

       do i=1,N
         pRes(i) = 0.0d0
       end do

       end

C     ******************************************************************
C     ******************************************************************

C      rotina para copiar uma matriz
       subroutine sCopiaMatriz(N,M,pM,pRes)
       integer i,j,N,M
       double precision pM(N,M),pRes(N,M)

       do i=1,N
         do j=1,M
           pRes(i,j) = pM(i,j)
         end do
       end do

       end

C     ******************************************************************
C     ******************************************************************

       subroutine sCopiaVetor(N,pM,pRes)

C      rotina para copiar um vetor
C
C      Parametros de Entrada:
C
C      N    integer,
C           numero de variaveis,
C
C      pM   double precision pM(N)
C           vetor a ser copiado
C
C      pRes double precision pRes(N)
C           vetor copia

C      SCALAR ARGUMENTS
       integer N

C      ARRAY ARGUMENTS
       double precision pM(N),pRes(N)

C      LOCAL SCALARS
       integer i

       do i=1,N
         pRes(i) = pM(i)
       end do

       end

C     ******************************************************************
C     ******************************************************************

C      rotina para copiar uma coluna de uma matriz em um vetor
       subroutine sCopiaMatVetor(N,M,pCol,pM,pRes)
       integer N,M,pCol,i
       double precision pM(N,M),pRes(N)

       do i=1,N
         pRes(i) = pM(i,pCol)
       end do

       end

C     ******************************************************************
C     ******************************************************************

C      rotina para soma de matrizes
       subroutine sSomaMatriz(N,M,pMatriz1,pMatriz2,pRes)
       integer N,M,i,j
       double precision pMatriz1(N,M),pMatriz2(N,M),pRes(N,M)

       do i=1,N
         do j=1,M
           pRes(i,j) = pMatriz1(i,j) + pMatriz2(i,j)
         end do
       end do

       end

C     ******************************************************************
C     ******************************************************************

C      rotina para soma de vetores
       subroutine sSomaVetor(N,pMatriz1,pMatriz2,pRes)
       integer N,i
       double precision pMatriz1(N),pMatriz2(N),pRes(N)

       do i=1,N
           pRes(i) = pMatriz1(i) + pMatriz2(i)
       end do

       end

C     ******************************************************************
C     ******************************************************************

C      rotina para  subtracao de vetores
       subroutine sDiferencaVetor(N,pMatriz1,pMatriz2,pRes)
       integer N,i
       double precision pMatriz1(N),pMatriz2(N),pRes(N)

       do i=1,N
           pRes(i) = pMatriz1(i) - pMatriz2(i)
       end do

       end

C     ******************************************************************
C     ******************************************************************

      subroutine sum_vv(oper,x,y,n,z)

C     This subroutine makes a sum or a subtraction of two vectors

C     Parameters of the subroutine:
C
C     On Entry:
C
C     oper     integer,
C              indicates if it was a sum or a subtraction:
C              1: indicates a sum
C              other value: indicates a subtraction
C
C     x        double precision x(n),
C              vector 1,
C
C     y        double precision y(n),
C              vector 2,
C
C     n        integer,
C              number of variables,
C
C     On Return
C
C     z        double precision z(n)
C              sum's result
C
C     PARAMETERS

C     SCALAR ARGUMENTS
      integer n,oper

C     ARRAY ARGUMENTS
      double precision x(n),y(n),z(n)

C     LOCAL SCALARS
      integer i
      double precision signal

C     LOCAL ARRAYS

C     SET LOCAL SCALARS
      if ( oper .eq. 1 ) then

          signal = 1.0d0

      else

          signal = -1.0d0

      end if

C     CALCULATE THE OPERATION

      do i = 1,n

          z(i) = x(i) + ( signal * y(i) )

      end do

      end


C     ******************************************************************
C     ******************************************************************

C      rotina para multiplica�o escalar x matriz
       subroutine sMultNumMatriz(pNum,N,M,pMatriz,pRes)
       integer N,M,i,j
       double precision pNum, pMatriz(N,M),pRes(N,M)

       do i=1,N
         do j=1,M
           pRes(i,j) = pNum * pMatriz(i,j)
         end do
       end do

       end

C     ******************************************************************
C     ******************************************************************

       subroutine sMultNumVetor(pNum,N,pVetor,pRes)
C
C      Rotina para calcular a multiplicacao de um scalar com um vetor
C
C      Parametros de Entrada:
C
C      pNum double precision
C           escalar
C      N integer
C        numero de variaveis
C      pVetor double precision pVetor(N)
C             vetor
C
C      Parametros de Saida:
C
C      pRes double precision pRes(N)
C           resultado da multiplicao
C
C      SCALAR ARGUMENTS
       integer N,i
       double precision pNum

C      ARRAY ARGUMENTS
       double precision pVetor(N),pRes(N)

       do i=1,N
           pRes(i) = pNum * pVetor(i)
       end do

       end

C     ******************************************************************
C     ******************************************************************

C      rotina para multiplica�o matriz x matriz
       subroutine sMultMatrizMatriz(N,M,pM1,O,pM2,pRes)
       integer N,M,O,i,j,k
       double precision pM1(N,M),pM2(M,O),pRes(N,O)

       call sInicializaMatriz(N,O,pRes)

       do i=1,N
         do j=1,O
           do k=1,M
             pRes(i,j)=pRes(i,j)+(pM1(i,k)*pM2(k,j))
           end do
         end do
       end do

       end


C     ******************************************************************
C     ******************************************************************

C      rotina para multiplica�o vetor x matriz
       subroutine sMultVetMatriz(N,M,pV,pM,pRes)
       integer N,M,i,j
       double precision pV(N),pM(N,M),pRes(M)

       call sInicializaVetor(M,pRes)

       do i=1,M
         do j=1,N
           pRes(i)=pRes(i)+(pV(i)*pM(j,i))
         end do
       end do

       end

C     ******************************************************************
C     ******************************************************************

C      rotina para multiplica�o matriz x vetor
       subroutine sMultMatVetor(N,M,pM,pV,pRes)
       integer N,M,i,j
       double precision pV(M),pM(N,M),pRes(N)

       call sInicializaVetor(N,pRes)

       do i=1,N
         do j=1,M
           pRes(i)=pRes(i)+(pM(i,j)*pV(j))
         end do
       end do

       end

C     ******************************************************************
C     ******************************************************************

       subroutine sMultTranspVetor(N,pV1,pV2,pRes)

C      rotina para multiplicacao vetor transposto x vetor
C
C      Parametros de Entrada:
C
C      N integer
C        numero de veriaveis
C      pV1 double precision  pV1(N)
C          vetor transposto
C      pV2 double precision pv2(N)
C          segundo vetor
C
C      Parametros de Saida:
C
C      pRes double precision
C           resultado da multiplicacao
C

C      SCALARS ARGUMENTS
       integer N
       double precision pRes

C      ARRAYS ARGUMENTS
       double precision pV1(N),pV2(N)

C      LOCAL SCALARS
       integer i

       pRes = 0.0d0

       do i=1,N
         pRes=pRes+(pV1(i)*pV2(i))
       end do

       end

C     ******************************************************************
C     ******************************************************************

      subroutine mul_vtv(x,y,n,z)

C     This subroutine multiplies two vectors on the form:
C            transpost vector 1 * vector 2 = return scalar
C     Parameters of the subroutine:
C
C     On Entry:
C
C     x        double precision x(n),
C              vector 1,
C
C     y        double precision y(n),
C              vector 2,
C
C     n        integer,
C              number of variables,
C
C     On Return
C
C     z        double precision
C              return scalar
C
C     PARAMETERS

C     SCALAR ARGUMENTS
      integer n
      double precision z

C     ARRAY ARGUMENTS
      double precision x(n),y(n)

C     LOCAL SCALARS
      integer i

C     LOCAL ARRAYS

C     SET LOCAL SCALARS
      z = 0.0d0

C     CALCULATE THE MULTIPLICATION

      do i = 1,n

          z = z + x(i) * y(i)

      end do

      end

C     ******************************************************************
C     ******************************************************************

C      rotina para multiplica�o vetor x vetor transposto
       subroutine sMultVetTransposto(N,pV1,M,pV2,pRes)
       integer N,M,i,j
       double precision pV1(N),pV2(M),pRes(N,M)

       call sInicializaMatriz(N,M,pRes)

       do i=1,N
         do j=1,M
           pRes(i,j)=pV1(i)*pV2(j)
         end do
       end do

       end

C     ******************************************************************
C     ******************************************************************

C      rotina para transpor uma matriz
       subroutine sTranspoeMatriz(N,M,pMatriz,pRes)
       integer N,M,i,j
       double precision pMatriz(N,M),pRes(M,N)

       call sInicializaMatriz(M,N,pRes)

       do i=1,N
         do j=1,M
           pRes(j,i)=pMatriz(i,j)
         end do
       end do

       end

C     ******************************************************************
C     ******************************************************************

C      rotina para calcular a norma de uma matriz
       subroutine sNormaMatriz(N,M,pMatriz,pNorma)
       integer N,M,i,j
       double precision pMatriz(N,M),pNorma,vAux

       pNorma = 0.0d0
       vAux = 0.0d0

       do i=1,N
         do j=1,M
           vAux = pMatriz(i,j)*pMatriz(i,j)
           pNorma = pNorma + vAux
         end do
       end do
       vAux = pNorma
       pNorma = sqrt(vAux)

       end


C     ******************************************************************
C     ******************************************************************

C      rotina para calcular a norma infinita de um vetor
       subroutine sNormaInfVetor(N,pVetor,pNorma)
       integer N,i
       double precision pVetor(N),pNorma

       pNorma = abs(pVetor(1))

       do i=2,N
         if(pNorma.lt.abs(pVetor(i))) then
           pNorma = abs(pVetor(i))
         end if
       end do

       end

C     ******************************************************************
C     ******************************************************************

       subroutine sNormaVetor(N,pVetor,pNorma)
C
C      Rotina para calcular a norma quadr�ica de um vetor
C
C      Parametros de Entrada:
C
C      N integer
C        numero de variaveis
C      pVetor double precision pVetor(N)
C             vetor que sera calculada a norma
C
C      Parametros de Saida:
C
C      pNorma double precision
C             valor da norma quadratica de pVetor

C      SCALAR ARGUMENTS
       integer N
       double precision pNorma

C      ARRAY ARGUMENTS
       double precision pVetor(N)

C      LOCAL SCALARS
       integer i
       double precision vAux

       pNorma = 0.0d0

       do i=1,N
           pNorma = pNorma + pVetor(i)*pVetor(i)
       end do
       vAux = pNorma
       pNorma = sqrt(vAux)

       end


C     ******************************************************************
C     ******************************************************************

C       Calculo da Matriz Inversa
        subroutine matinv(pM,N,pInv)
C       Calcula a inversa de uma matriz NxN
C
        integer N,i,j,p,o,q,r
        double precision pM(N,N),pInv(N,N),vM(N,N)
        double precision vPivo,vAux,vAuxInv

C       Montando a matriz identidade e copiando a matriz pM
        do i = 1,N
          do j = 1,N
            if( i.eq.j )then
              pInv(i,j) = 1.0d0
            else
              pInv(i,j) = 0.0d0
            end if
            vM(i,j) = pM(i,j)
          end do
        end do

C       Calculando a inversa por Eliminicao de Gauss-Jordan
        do p = 1,N
          do o = 1,N

            if(o.ne.p) then
              if(vM(p,p).eq.0.0d0) then
                write(*,*) 'Matriz indicada e singular!'
                exit
              end if

              vPivo = - vM(o,p)/vM(p,p)

              do q = 1,N
                vAux=vPivo*vM(p,q)
                vM(o,q)=vM(o,q)+vAux
                vAuxInv=vPivo*pInv(p,q)
                pInv(o,q)=pInv(o,q)+vAuxInv
              end do

            else
              vPivo = vM(o,p)

              do r = 1,N
                vM(o,r) = vM(o,r)/vPivo
                pInv(o,r) = pInv(o,r)/vPivo
              end do

            end if
          end do
        end do
       end


C     ******************************************************************
C     THIS SECTION CONTAINS DEVELOPMENT OR EXPERIMENTAL CODE ONLY.
c     DO NOT USE THIS SOURCES!!!
C     ******************************************************************

C     metodo de Cholesky Modif. Retorna TRUE quando matriz e def pos
C     retorna false e a matriz modificada quando nao e
      subroutine sCholeskyModif(N,pM,pRes)
      integer N
      double precision pM(N,N)
      logical pRes
C     *****
      integer i,j,k,l
      double precision vL(N,N),vA(N,N),vAux,vAprox,vE(N)
      logical vTeste

      vAprox = 1
      pRes = .true.
      vTeste = .true.

      call sCopiaMatriz(N,N,pM,vA)
      call sInicializaMatriz(N,N,vL)

      do j=1,N
        vE(j) = abs(vA(j,j))
      end do

      i = 1
      do while(vTeste)

        if(vA(i,i).le.0) then
C          vTeste = .false.
           vE(i) = vE(i)+abs(vA(i,i))+vAprox
           vA(i,i) = 0
          pRes = .false.
        else
          vL(i,i) = sqrt(vA(i,i))

          do j=i+1,N
            vL(j,i) = vA(j,i)/vL(i,i)

            do k=i+1,j
              vAux = vL(j,i)*vL(k,i)
              vA(j,k) = vA(j,k) - vAux
            end do

          end do

          if(i.ge.N) then
            vTeste = .false.
          else
            i = i+1
          end if

        end if

      end do

C     soma o valor de E na diagonal principal
      do l=1,N
        pM(l,l) = vE(l)
      end do

C      write(*,*) 'M :'
C      write(*,*) pM(1,1),pM(1,2),pM(1,3)
C      write(*,*) pM(2,1),pM(2,2),pM(2,3)
C      write(*,*) pM(3,1),pM(3,2),pM(3,3)
C      write(*,*) ' '
C      write(*,*) 'A :'
C      write(*,*) vA(1,1),vA(1,2),vA(1,3)
C      write(*,*) vA(2,1),vA(2,2),vA(2,3)
C      write(*,*) vA(3,1),vA(3,2),vA(3,3)
C      write(*,*) ' '
C      write(*,*) 'L :'
C      write(*,*) vL(1,1),vL(1,2),vL(1,3)
C      write(*,*) vL(2,1),vL(2,2),vL(2,3)
C      write(*,*) vL(3,1),vL(3,2),vL(3,3)
C      write(*,*) ' '
C      write(*,*) 'E: '
C      write(*,*) vE
C      write(*,*) ' '
C      write(*,*) 'RES => ',pRes
C      write(*,*) 'i   => ',i
C      call sTranspoeMatriz(N,N,vL,vTL)
C      call sMultMatrizMatriz(N,N,vL,N,vTL,vML)
C      write(*,*) ' '
C      write(*,*) 'ML :'
C      write(*,*) vML(1,1),vML(1,2),vML(1,3)
C      write(*,*) vML(2,1),vML(2,2),vML(2,3)
C      write(*,*) vML(3,1),vML(3,2),vML(3,3)

      end

C     ******************************************************************
C     ******************************************************************