      implicit real*8(a-h,o-z)
      implicit integer*8(i-n)
!
!  Definimos las variables para aplicar dimensionamiento dinamico
!
      allocatable :: lp1(:),lp2(:),rgmat(:),
     &               gload(:),xdisp(:),resdl(:),
     &               ngfix(:),presc(:),treac(:),
     &               wrksp(:)
!
!  Establecemos los parametros del problema
!
      ntotv=500     ! Numero total de variables
      nvfix=2        ! Numero de variables con desplazamientos prescritos
      niter=500      ! Numero de iteraciones por restart (dimension del subespacio de Krylov)
                     ! Como maximo puede ser igual a la dimension ntotv de la matriz
      nrestart=5  ! Numero de veces que se reinicia el GMRES. nrestart=0 --> Sin reiniciar
      nactu=10        ! Separacion entre iteraciones en las que se imprime el residuo

      epslr=1.d-10   !  Error absoluto permitido para el residuo
      rltvr=1.d-07   !  Error relativo permitido para el residuo

      ntotr=3*ntotv-2   ! Numero de componentes a almacenar para la matriz del sistema.
!
!  Se asigna de forma dinamica la memoria necesaria para la definicion de la matriz del sistema
!
      allocate(lp1(ntotr))
      allocate(lp2(ntotr))
      allocate(rgmat(ntotr))
      allocate(ngfix(nvfix))
      allocate(presc(nvfix))
      allocate(treac(nvfix))
!
!  Se asigna la memoria para los siguientes vectores del sistema y para el vector de trabajo
!  del metodo GMRES (wrksp)
!
      allocate(gload(ntotv))
      allocate(xdisp(ntotv))
      allocate(resdl(ntotv))
      allocate(wrksp((ntotv+4)*(niter+3)))
!
!  Definicion de la matriz de rigidez, del vector de terminos independientes y de la solución
!  inicial de partida
!
      do itotv=1,ntotv
        xdisp(itotv)=0.d+00
      enddo
!
!  Matriz de rigidez
!
      lp1(1)= 1
      lp2(1)= 1
      rgmat(1)=2.d+00
      ind=1
      do i=2,ntotv
        lp1(ind+1)= i-1
        lp1(ind+2)= i
        lp1(ind+3)= i
        lp2(ind+1)=i
        lp2(ind+2)=i-1
        lp2(ind+3)=i
        rgmat(ind+1)= -1.d+00
        rgmat(ind+2)= -1.d+00
        rgmat(ind+3)=  2.d+00
        ind=ind+3
      enddo

!
!  Vector de terminos independientes
!
      do i=1,ntotv
        gload(i)= 1.d+00/((dfloat(ntotv+1))**2)
      enddo
!
!  Variables prescritas
!
      ngfix(1)=1
      ngfix(2)=ntotv
      presc(1)=0.d+00
      presc(2)=0.d+00
      treac(1)=0.d+00
      treac(2)=0.d+00
!
!  Llamada a la subrutina Solver_GMRES que aplica el Metodo GMRES al sistema definido anteriormente
!  con precondicionamiento diagonal
!
      call SLE_Solver_GMRES_SP(ntotv,ntotr,nvfix,
     &                      niter,nrestart,nactu,epslr,rltvr,
     &                      lp1,lp2,rgmat,gload,xdisp,resdl,
     &                      ngfix,presc,treac,
     &                      wrksp,
     &                      rnorm,kiter,irestart,ierr,cputime)
!
!  Escritura de resultados 
!
      write(6,'(/a,e15.6)') ' RNORM   = ',rnorm
      write(6,'(a,i8)') ' RESTART = ',irestart
      write(6,'(a,i8,a,i6,a,i4,a,i6,a)') ' ITER    = ',niter*(irestart)
     &                                   +kiter, '  (',
     &                      niter,' * ',irestart,' + ',kiter,' )'
      write(6,'(a,f9.2,a)') ' CPUTIME =',cputime,' segundos'
      if (ierr.eq.0)then
        write(6,'(a,i8,a/)')' FLAG    = ',ierr,'   --->  Convergencia'
      elseif(ierr.eq.1)then
        write(6,'(a,i8,2a/)')' FLAG    = ',ierr,'   --->  No.',
     &                       ' maximo de iteraciones alcanzado'
      elseif(ierr.eq.2)then
        write(6,'(a,i8,2a/)')' FLAG    = ',ierr,'   --->  GMRES',
     &                       ' estancado'
      endif

      read(5,*)

      if (ierr.eq.0)then

      do itotv=1,ntotv
        write(6,'(a,i8,a,e20.10)') ' ITOTV   = ',itotv,'   XDISP=',
     &                              xdisp(itotv)
      enddo
      
      write(6,'()')

      do ivfix=1,nvfix
        write(6,'(a,i8,a,e20.10)') ' ITOTV   = ',ngfix(ivfix),
     &                            '   TREAC=',treac(ivfix)
      enddo
      
      endif

      end

