! PROGRAM BUILDS THE HILBERT MATRIX OF GIVEN DIMENSION WITH 
! ELEMENTS OF GIVEN SIZE (EXPONENT AND MANTISSA LENGTH). 
! COMPUTES SOLUTION OF THE SYSTEM WITH RIGHT HAND SIDE VECTOR
! GIVEN IN TEXT FILE 'HILBERT_RHS.TXT'. 
! FINDS MAX ELEMENT IN EUCLIDIAN NORM OF RESIDUAL VECTOR
! STORES IN TEXT FILE 'HILBERT_SOL.TXT' THE SOLUTION VECTOR
! AND THE MAX ELEMENT OF RESIDUAL VECTOR

      PROGRAM HILBERT

	CHARACTER* 60 TXTDAT, STROUT
      INTEGER IHRTNL, IHVECT, IHVRES, IHLBRT, IHCLU
      INTEGER NEXP, NMNT, NDIM, INDEX, IROW, ICOL, IDENOM, ICODE 

      WRITE( *, '("ENTER DIMENSION OF HILBERT MATRIX: ",$)' )
      READ( *, * ) NDIM
      WRITE( *, '("ENTER EXPONENT LENGTH (IN 32-BIT WORDS) 
     $OF ELEMENT OF MATRIX : ",$)' )
      READ( *, * ) NEXP
      WRITE( *, '("ENTER MANTISSA LENGTH (IN 32-BIT WORDS) 
     $OF ELEMENT OF MATRIX : ",$)' )
      READ( *, * ) NMNT 

! START EXLAF77 1.01 WORKING SESSION
      CALL HSINIT( 'HILBERT', 256 )
! MASK TEXT MESSAGE OF SPECIFIC FOR THE CASE ERROR
      CALL HSEMSK( 503 )
! CREATES NEW RIGHT HAND SIDE EMPTY H-VECTOR WITH HANDLE IHVECT
      CALL HMV( NEXP, NMNT, .TRUE., NDIM, IHVECT, *2000 )	   

! UPDATES H-VECTOR IHVECT WITH DATA FROM EXT FILE 'HILBERT_RHS.TXT'
      OPEN( 8, FILE = 'HILBERT_RHS.TXT' )
      DO INDEX = 1, NDIM
        READ( 8, * ) TXTDAT
     	  CALL HUEVT( TXTDAT, INDEX, IHVECT, *2000 )
      ENDDO
      CLOSE( 8 )
 
! CREATES NEW HILBERT EMPTY H-MATRIX WITH HANDLE IHLBRT (IN PACKED FORM)
      CALL HMMS( NEXP, NMNT, .TRUE., NDIM, 1, IHLBRT, *2000 )

! UPDATES HILBERT H-MATRIX WITH HANDLE IHLBRT (FILLS WITH DATA)
	DO IROW = 1, NDIM
	  DO ICOL = 1, NDIM
	    IF( IROW .LE. ICOL ) THEN
	      IDENOM = IROW + ICOL - 1
	      ! CREATES RATIONAL H-NUMBER WITH HANDLE IHRTNL
	      CALL HAXF( 1, IDENOM, IHRTNL, *2000 ) 
	      ! UPDATE ELEMENT OF H-MATRIX IHLBRT WITH H-NUMBER IHRTNL
	      CALL HUEMN( IHRTNL, IROW, ICOL, IHLBRT, *2000 )
	      ! DELETES H-NUMBER IHRTNL
	      CALL HSDOBJ( IHRTNL, *2000 )
          ENDIF
        ENDDO
      ENDDO

! CREATES COPY OF H-VECTOR IHVECT ( H-VECTOR IHVRES )
	CALL HACPYH( IHVECT, IHVRES, *2000 )
! CREATES COPY OF H-MATRIX IHLBRT ( H-MATRIX IHCLU )
      CALL HACPYH( IHLBRT, IHCLU, *2000 )	

! PERFORMS COMPLETE LU DECOMPOSITION OF H-MATRIX IHCLU 
      CALL HUCLU( IHCLU, *1000 )
! COMPUTES SOLUTION OF THE SYSTEM WITH RESULT IN H-VECTOR IHVRES
      CALL HUMHH( IHCLU, IHVRES, 'R', *2000 )
	    
! DELETES H-MATRIX IHCLU
	CALL HSDOBJ( IHCLU, *2000 )

! OUTPUT RESULT VECTOR TO THE TEXT FILE 'HILBERT_SOL.TXT'
      OPEN( 8, FILE = 'HILBERT_SOL.TXT' )
      WRITE( 8, * ) ' RESULT VECTOR '
      DO INDEX = 1, NDIM
        CALL HETEV( IHVRES, INDEX, 44, 1, 30, 10, STROUT, *2000 )
      WRITE( 8, * ) STROUT
      ENDDO

! FINDING RESIDUAL VECTOR 
! UPDATE MULTIPLICATION OF H-MATRIX IHLBRT BY H-VECTOR IHVRES 
      CALL HUMHH( IHLBRT, IHVRES, 'R', *2000 )
! UPDATE SUBTRACTION H-VECTOR IHVRES FROM RHS H-VECTOR IHVECT
	CALL HUSHH( IHVECT, IHVRES, *2000 )
! FINDS INDEX OF THE GREATEST IN EUCLIDIAN NORM ELEMENT OF RESIDUAL VECTOR
	CALL HGVG2( IHVRES, INDEX, *2000 )

! PRINT THE ELEMENT 
	CALL HETEV( IHVRES, INDEX, 44, 1, 30, 10, STROUT, *2000 )
	PRINT *
	PRINT *,' MAX ELEMENT IN EUCLIDIAN NORM OF RESIDUAL VECTOR IS '
	PRINT *, STROUT

      WRITE( 8, * )'MAX ELEMENT IN EUCLIDIAN NORM OF RESIDUAL VECTOR IS'
      WRITE( 8, * ) STROUT
      CLOSE( 8 )

	GOTO 3000 
! GETS CODE OF ERROR FROM SUBROUTINE HUCLU
 1000 CALL HSERR( ICODE )
      IF( ICODE .EQ. 503 ) THEN
! SPECIFIC ERROR OF INSUFFICIENT SIZE OF ELEMENTS OF HILBERT MATRIX
	  PRINT *, 'THE MATRIX IS ALGORITHMICALLY INDEFINITE'
	  PRINT *, '(INSUFFICIENT WORKING PRECISION)'
        GOTO 3000
	ELSE
	 GOTO 2000 
      ENDIF
	GOTO 3000       
 2000 PRINT*,'THE ERROR JUMP. SEE HILBERT.LOG!'
 3000 CONTINUE

! CLOSE EXLAF77 WORKING SESSION
      CALL HSEXIT()   
	END PROGRAM HILBERT