!*************************************************************************
!*                                                                       *
!*     Computation of function and subgradients values needed to define  *
!*     the clusterwise support vector linear regression problem.         *
!*     (last modified 01.09.2019).                                       *
!*                                                                       *
!*************************************************************************
!*
!*     Modules included:
!*
!*     clrfunctionmod         !

MODULE clrfunctionmod ! Computation of function and (sub)gradients values for
                      ! clusterwise support vector linear regression problem.
    USE omp_lib

    USE r_precision, ONLY : prec      ! Precision for reals.
    IMPLICIT NONE

    ! MODULE clrfunctionmod includes the following subroutines.

    PUBLIC :: &
        hauxfunc, &        ! Computation of h(x^i,y_i,a,b) in aux clusterwise linear regression problem.
        hauxfunc_obj4, &   ! Computation of the value |(x^i)^T a +y_i - b| used in auxiliary clust. prob. objfunc=4
        hfunc, &           ! Computation of h(x^i,y_i,a,b) in clusterwise linear regression problem.
        hauxsub , &        ! Computation of subgradient of h(x^i,y_i,a,b) in aux clusterwise linear regression problem.
        hauxsub_obj4, &    ! Computation of the subgradient of |(x^i)^T a +y_i - b| used in auxiliary clust. prob. objfunc=4
        auxfunc, &         ! Computation of aux clusterwise linear regression problem.
        auxfunc1, &        ! Computation of f1 for aux clusterwise linear regression problem.
        auxfunc2, &        ! Computation of f2 for aux clusterwise linear regression problem.
        clrfunc, &         ! Computation of clusterwise linear regression problem.
        clrfunc1, &        ! Computation of f1 for clusterwise linear regression problem.
        clrfunc2, &        ! Computation of f2 for clusterwise linear regression problem.
        fgrad, &           ! Computation of subgradient.
        fsubgrad1, &       ! Computation of subgradient of the first DC component when objfunc=2.
        fsubgrad2, &       ! Computation of subgradient of the second DC component when objfunc=2.
        fgrad1, &          ! Computation of subgradient of the first DC component.
        fgrad2, &          ! Computation of subgradient of the second DC component.
        fvb, &
        fvb1, &
        fvb2, &
        gradient, &
        gradient1, &
        gradient2


CONTAINS


    !=============================================================
    ! Computation of h(x^i,y_i,a,b) used in auxiliary clust. prob.
    !=============================================================
    SUBROUTINE hauxfunc(x,clust,point,f)

        USE param, ONLY : zero, large
        USE initclr, ONLY : &
            maxdim, &      ! Maximum number of variables in optimization,
                           !   maxdim >= nclust * nft.
            a4, &          ! Data matrix a4(maxrec,maxnft).
            nft, &         ! Number of used features.
            nc, &          ! Current number of clusters, loops from 1 to nclust.
            nrecord1, &    ! Number of records in data.
            nfit, &        ! The number of function evaluations.
            objfunc, &     ! The objective function used.
            objeps         ! The parameter used in objective function objfunc=2, 3 and 4.

        IMPLICIT NONE

        ! Subroutine arguments
        REAL(KIND=prec) ::  &
            x(:), &
            f

        INTEGER :: &
            clust, &   !number of the cluster whose regression function is considered
            point      !data point used in computation   
            
        ! Local variables
        REAL(KIND=prec) ::  &
            ff,d1,d2
        INTEGER :: &
            i,j,k

        INTRINSIC ::  DMIN1
        
        f=zero      
        k=point
        i=clust

        IF (objfunc==1) THEN                            ! objective 1
          ff=x(nft*i) 
            DO j=1,nft-1
              ff=ff+a4(k,j)*x(j+(i-1)*nft)
            END DO
          f=(ff-a4(k,nft))**2           
        
        ELSE IF (objfunc==2 .OR. objfunc==4) THEN       ! objective 2 or 4
          ff=x(nft*i)
            DO j=1,nft-1
              ff=ff+a4(k,j)*x(j+(i-1)*nft)
            END DO
          ff= ABS(ff-a4(k,nft))-objeps
          IF (ff>0.0_prec) THEN 
              f = ff
          ELSE    
              f = 0.0_prec
          END IF          
        
        ELSE IF (objfunc==3) THEN                       ! objective 3
        
          ff=x(nft*i)
            DO j=1,nft-1
              ff=ff+a4(k,j)*x(j+(i-1)*nft)
            END DO
          f= (ff-a4(k,nft))**2-objeps
          IF (f>0.0_prec) THEN 
              f = f
          ELSE    
              f = 0.0_prec
          END IF         
        
        END IF
    
    
        RETURN

    END SUBROUTINE hauxfunc
    
    !=============================================================
    ! Computation of the value |(x^i)^T a +y_i - b| used in auxiliary clust. prob. objfunc=4
    !=============================================================
    SUBROUTINE hauxfunc_obj4(x,clust,point,f)

        USE param, ONLY : zero, large
        USE initclr, ONLY : &
            maxdim, &      ! Maximum number of variables in optimization,
                           !   maxdim >= nclust * nft.
            a4, &          ! Data matrix a4(maxrec,maxnft).
            nft, &         ! Number of used features.
            nc, &          ! Current number of clusters, loops from 1 to nclust.
            nrecord1, &    ! Number of records in data.
            nfit, &        ! The number of function evaluations.
            objfunc, &     ! The objective function used.
            objeps         ! The parameter used in objective function objfunc=2, 3 and 4.

        IMPLICIT NONE

        ! Subroutine arguments
        REAL(KIND=prec) ::  &
            x(:), &
            f

        INTEGER :: &
            clust, &   !number of the cluster whose regression function is considered
            point      !data point used in computation   
            
        ! Local variables
        REAL(KIND=prec) ::  &
            ff,d1,d2
        INTEGER :: &
            i,j,k

        INTRINSIC ::  DMIN1
        
        f=zero      
        k=point
        i=clust

        
          ff=x(nft*i)
            DO j=1,nft-1
              ff=ff+a4(k,j)*x(j+(i-1)*nft)
            END DO
          ff= ABS(ff-a4(k,nft))
          
          f = ff
    
        RETURN

    END SUBROUTINE hauxfunc_obj4    

    !====================================================
    ! Computation of h(x^i,y_i,a,b) used in clust. prob.
    !====================================================
    SUBROUTINE hfunc(x,clust,point,f)

        USE param, ONLY : zero, large
        USE initclr, ONLY : &
            maxdim, &      ! Maximum number of variables in optimization,
                           !   maxdim >= nclust * nft.
            a, &           ! Data matrix a(nrecord,nft), from input file.
            nft, &         ! Number of used features.
            nc, &          ! Current number of clusters, loops from 1 to nclust.
            nrecord1, &    ! Number of records in data.
            nfit, &        ! The number of function evaluations.
            objfunc, &     ! The objective function used.
            objeps         ! The parameter used in objective function objfunc=2, 3 and 4.
            
        IMPLICIT NONE

        ! Subroutine arguments
        REAL(KIND=prec) ::  &
            x(:), &
            f

        INTEGER :: &
            clust, &   !number of the cluster whose regression function is considered
            point      !data point used in computation   
            
        ! Local variables
        REAL(KIND=prec) ::  &
            ff,d1,d2
        INTEGER :: &
            i,j,k

        INTRINSIC ::  DMIN1
        
        f=zero      
        k=point
        i=clust

        IF (objfunc==1) THEN                       ! objective 1
          ff=x(nft*i)
            DO j=1,nft-1
              ff=ff+a(k,j)*x(j+(i-1)*nft)
            END DO
          f=(ff-a(k,nft))**2        
        
        ELSE IF (objfunc==2 .OR. objfunc==4) THEN  ! objective 2 or 4
          ff=x(nft*i)
            DO j=1,nft-1
              ff=ff+a(k,j)*x(j+(i-1)*nft)
            END DO
          ff= ABS(ff-a(k,nft))-objeps
          IF (ff>0.0_prec) THEN 
              f = ff
          ELSE    
              f = 0.0_prec
          END IF          
        
        ELSE IF (objfunc==3) THEN                  ! objective 3
          ff=x(nft*i)
            DO j=1,nft-1
              ff=ff+a(k,j)*x(j+(i-1)*nft)
            END DO
          f= (ff-a(k,nft))**2-objeps
          IF (f>0.0_prec) THEN 
              f = f
          ELSE    
              f = 0.0_prec
          END IF            
        END IF  
    
        RETURN

    END SUBROUTINE hfunc
	
	
    !=============================================================
    ! Computation of the value of the quadratic term in clusterwise linear regression problem
    !=============================================================
    SUBROUTINE lisa(x,f)

        USE param, ONLY : zero, large
        USE initclr, ONLY : &
            maxdim, &      ! Maximum number of variables in optimization,
                           !   maxdim >= nclust * nft.
            objc, &        ! The parameter used in objective function objfunc=2, 3 and 4.
            a4, &          ! Data matrix a4(maxrec,maxnft).
            nft, &         ! Number of used features.
            nc, &          ! Current number of clusters, loops from 1 to nclust.
            nrecord1, &    ! Number of records in data.
            objfunc, &     ! Objective function used.
            objeps, &      ! The parameter used in objective function objfunc=2, 3 and 4.   
            objkvad, &     ! The parameter used in objective function objfunc=2, 3 and 4.           
            nfit           ! The number of function evaluations.

        IMPLICIT NONE

        ! Subroutine arguments
        REAL(KIND=prec) ::  &
            x(maxdim), &
            f

        ! Local variables
        REAL(KIND=prec) ::  &
            fa,fb,d1,d2,smallest
        INTEGER :: &
            i,j,k

        INTRINSIC ::  DMIN1

        !f=zero
       		
		IF (objfunc==2 .OR. objfunc==3 .OR. objfunc==4) THEN  !objectives 2, 3 and 4
              fa = 0.0_prec
              
            IF (objkvad > (10.0_prec)**(-10)) THEN 
			  DO k=1,nc
                DO i=1,nft-1
                 fa = fa +x(i+(k-1)*nft)*x(i+(k-1)*nft)
                END DO
              END DO
              fa = 0.5_prec * objkvad * fa
            END IF
			
	        IF (objc<=(10.0_prec)**(-10)) THEN 
              f = zero + fa
            END IF
			f = fa + objc * f
        ELSE
		    f = f
        END IF 
		
	    end subroutine lisa
    
    !================================================================================
    ! Computation of the subgradient of h(x^i,y_i,a,b) used in auxiliary clust. prob.
    !================================================================================
    SUBROUTINE hauxsub(x,clust,point,grad)

        USE param, ONLY : zero, large, two
        USE initclr, ONLY : &
            maxdim, &      ! Maximum number of variables in optimization,
                           !   maxdim >= nclust * nft.
            a4, &          ! Data matrix a4(maxrec,maxnft).
            nft, &         ! Number of used features.
            nc, &          ! Current number of clusters, loops from 1 to nclust.
            nrecord1, &    ! Number of records in data.
            nfit, &        ! The number of function evaluations.
            objfunc, &     ! The objective function used.
            objeps         ! The parameter used in objective function objfunc=2, 3 and 4.

        IMPLICIT NONE

        ! Subroutine arguments
        REAL(KIND=prec) ::  &
            x(:)
            
        REAL(KIND=prec) ::  &
            grad(nft)

        INTEGER :: &
            clust, &   !number of the cluster whose regression function is considered
            point      !data point used in computation   
            
        ! Local variables
        REAL(KIND=prec) ::  &
            ff,d1,d2,b,e
        INTEGER :: &
            i,j,k

        INTRINSIC ::  DMIN1
        
        DO i =1,nft
          grad(i)=zero  
        END DO      
        k=point
        i=clust
        e= (10.0_prec)**(-9)

        IF (objfunc==2 .OR. objfunc==4) THEN       ! objective 2 or 4
            ff=x(nft*i)
            DO j=1,nft-1
              ff=ff+a4(k,j)*x(j+(i-1)*nft)
            END DO
            ff=(ff-a4(k,nft))
            IF (ABS(ff)-objeps>0.0_prec) THEN 
              IF (ff>0.0_prec) THEN 
                DO j=1,nft-1
                   grad(j) = grad(j)+a4(k,j)
                END DO
                grad(nft) = 1.0_prec
              ELSE  
                DO j=1,nft-1
                   grad(j) =  grad(j)-a4(k,j)
                END DO      
                grad(nft) = -1.0_prec               
              END IF                
            ELSE
              DO i = 1,nft
                grad(i) = 0.0_prec
              END DO
            END IF  
        ELSE IF (objfunc==1) THEN                ! objective 1
            ff=x(nft*i)
            DO j=1,nft-1
              ff=ff+a4(k,j)*x(j+(i-1)*nft)
            END DO
            ff=(ff-a4(k,nft))
            DO i =1, nft-1
                grad(i)=grad(i)+two*ff*a4(k,i)
            END DO          
            grad(nft)=grad(nft)+two*ff
        ELSE IF (objfunc==3) THEN                ! objective 3
            ff=x(nft*i)
            DO j=1,nft-1
              ff=ff+a4(k,j)*x(j+(i-1)*nft)
            END DO
            ff=(ff-a4(k,nft))
            IF (ff**2-objeps > 0.0_prec) THEN
              DO i =1, nft-1
                grad(i)=grad(i)+two*ff*a4(k,i)
              END DO            
              grad(nft)=grad(nft)+two*ff    
            END IF            
        END IF
        
        
    
        RETURN

    END SUBROUTINE hauxsub
    
    !================================================================================
    ! Computation of the subgradient of |(x^i)^T a +y_i - b| used in auxiliary clust. prob. objfunc=4
    !================================================================================
    SUBROUTINE hauxsub_obj4(x,clust,point,grad)

        USE param, ONLY : zero, large, two
        USE initclr, ONLY : &
            maxdim, &      ! Maximum number of variables in optimization,
                           !   maxdim >= nclust * nft.
            a4, &          ! Data matrix a4(maxrec,maxnft).
            nft, &         ! Number of used features.
            nc, &          ! Current number of clusters, loops from 1 to nclust.
            nrecord1, &    ! Number of records in data.
            nfit, &        ! The number of function evaluations.
            objfunc, &     ! The objective function used.
            objeps         ! The parameter used in objective function objfunc=2, 3 and 4.

        IMPLICIT NONE

        ! Subroutine arguments
        REAL(KIND=prec) ::  &
            x(:)
            
        REAL(KIND=prec) ::  &
            grad(nft)

        INTEGER :: &
            clust, &   !number of the cluster whose regression function is considered
            point      !data point used in computation   
            
        ! Local variables
        REAL(KIND=prec) ::  &
            ff,d1,d2,b,e
        INTEGER :: &
            i,j,k

        INTRINSIC ::  DMIN1
        
        DO i =1,nft
          grad(i)=zero  
        END DO      
        k=point
        i=clust

            ff=x(nft*i)
            DO j=1,nft-1
              ff=ff+a4(k,j)*x(j+(i-1)*nft)
            END DO
            ff=(ff-a4(k,nft))
              IF (ff>0.0_prec) THEN 
                DO j=1,nft-1
                   grad(j) = grad(j)+a4(k,j)
                END DO
                grad(nft) = 1.0_prec
              ELSE  
                DO j=1,nft-1
                   grad(j) =  grad(j)-a4(k,j)
                END DO      
                grad(nft) = -1.0_prec   
              END IF                
    
        RETURN

    END SUBROUTINE hauxsub_obj4 
    

    !==========================================================
    ! Computation of aux clusterwise linear regression problem
    !==========================================================
    SUBROUTINE auxfunc(x,fval)

        USE param, ONLY : zero
        USE initclr, ONLY : &
            maxdim, &      ! Maximum number of variables in optimization,
                           !   maxdim >= nclust * nft.
            a4, &          ! Data matrix a4(maxrec,maxnft)
            nft, &         ! Number of used features.
            nrecord1, &    ! Number of records in data.
            objfunc, &     ! Objective function used.
            optmet, &      ! Optimization method.
            dminim, &      ! dminim(nrecord)
            naux           ! The total number of auxiliary problems computed.

        IMPLICIT NONE

        ! Subroutine arguments
        REAL(KIND=prec) ::  &
            x(maxdim), &
            fval

        ! Local variables
        REAL(KIND=prec) ::  &
            f1,f2,f3
        INTEGER :: &
            i,j

        INTRINSIC ::  DMIN1

        fval=zero
        IF (objfunc==1) THEN       ! objective 1 

          DO i=1,nrecord1
            f1=dminim(i)

            CALL hauxfunc(x,1,i,f3)                     

            fval=fval+DMIN1(f1,f3)
          END DO
          naux=naux+1
        ELSE                             ! objectives 2, 3 and 4
           CALL auxfunc1(x,f1)
           CALL auxfunc2(x,f2)
           fval = f1-f2
        END IF 
        RETURN

    END SUBROUTINE auxfunc
    
    !================================================================
    ! Computation of f1 in aux clusterwise linear regression problem
    !================================================================
    SUBROUTINE auxfunc1(x,fval)

        USE param, ONLY : zero
        USE initclr, ONLY : &
            maxdim, &      ! Maximum number of variables in optimization,
                           !   maxdim >= nclust * nft.
            a4, &          ! Data matrix a4(maxrec,maxnft)
            kvad, &        ! Scaling vector for kvadratic terms
            nft, &         ! Number of used features.
            objfunc, &     ! Objective function used.
            nrecord1, &    ! Number of records in data.
            dminim, &      ! dminim(nrecord)
            naux, &        ! The total number of auxiliary problems computed.
            objkvad, &     ! The parameter used in objective function objfunc=2, 3 and 4.
            objeps         ! The parameter used in objective function objfunc=2, 3 and 4.

        IMPLICIT NONE

        ! Subroutine arguments
        REAL(KIND=prec) ::  &
            x(maxdim), &
            fval

        ! Local variables
        REAL(KIND=prec) ::  &
            f1,f2,f3
        INTEGER :: &
            i,j

        INTRINSIC ::  DMIN1

        fval=zero

        IF (objfunc==2 .OR. objfunc==3 .OR. objfunc==4) THEN    ! objective 2, 3 or 4
		  IF (objkvad > (10.0_prec)**(-10)) THEN
            DO j = 1, nft-1
               fval = fval + objkvad*kvad(j)*0.5_prec*x(j)*x(j)
            END DO  
		  END IF	
        END IF
                             
        IF (objfunc ==1 .OR. objfunc==2 .OR. objfunc==3) THEN   ! objective 1 or 2 or 3
          DO i=1,nrecord1
            f2=dminim(i)

            CALL hauxfunc(x,1,i,f3)                     

            fval=fval+f2+f3
          END DO 
        ELSE IF (objfunc==4) THEN                               ! objective 4
          DO i=1,nrecord1
            f2=dminim(i)

            CALL hauxfunc_obj4(x,1,i,f3)                        

            fval=fval+MAX(f2+f3-objeps, MAX(f2,f3-objeps))
          END DO        
        END IF

        naux=naux+1
        RETURN

    END SUBROUTINE auxfunc1 
    
    !================================================================
    ! Computation of f2 in aux clusterwise linear regression problem
    !================================================================
    SUBROUTINE auxfunc2(x,fval)

        USE param, ONLY : zero
        USE initclr, ONLY : &
            maxdim, &      ! Maximum number of variables in optimization,
                           !   maxdim >= nclust * nft.
            a4, &          ! Data matrix a4(maxrec,maxnft)
            nft, &         ! Number of used features.
            objfunc, &     ! Objective function used.
            nrecord1, &    ! Number of records in data.
            dminim, &      ! dminim(nrecord)
            naux, &        ! The total number of auxiliary problems computed.
            objeps         ! The parameter used in objective function objfunc=2, 3 and 4.           

        IMPLICIT NONE

        ! Subroutine arguments
        REAL(KIND=prec) ::  &
            x(maxdim), &
            fval

        ! Local variables
        REAL(KIND=prec) ::  &
            f1,f2,f3
        INTEGER :: &
            i,j

        INTRINSIC ::  DMIN1

        fval=zero
        
        DO i=1,nrecord1
            f2=dminim(i)
          IF (objfunc==1 .OR. objfunc==2 .OR. objfunc==3) THEN       ! objective 1, 2 or 3   
            CALL hauxfunc(x,1,i,f3)                     

            fval=fval+MAX(f2,f3)     
          ELSE IF (objfunc==4) THEN                                  ! objective 4
            CALL hauxfunc_obj4(x,1,i,f3)
            
            fval=fval+MAX(f2,f3-objeps)
          END IF
        END DO

        !naux=naux+1
        RETURN

    END SUBROUTINE auxfunc2     


    !======================================================
    ! Computation of clusterwise linear regression problem
    !======================================================
    SUBROUTINE clrfunc(x,f)

        USE param, ONLY : zero, large
        USE initclr, ONLY : &
            maxdim, &      ! Maximum number of variables in optimization,
                           !   maxdim >= nclust * nft.
            a4, &          ! Data matrix a4(maxrec,maxnft).
            nft, &         ! Number of used features.
            nc, &          ! Current number of clusters, loops from 1 to nclust.
            nrecord1, &    ! Number of records in data.
            objfunc, &     ! Objective function used.  
            optmet, &      ! Optimization method.           
            nfit           ! The number of function evaluations.

        IMPLICIT NONE

        ! Subroutine arguments
        REAL(KIND=prec) ::  &
            x(maxdim), &
            f

        ! Local variables
        REAL(KIND=prec) ::  &
            f1,f2,d1,d2
        INTEGER :: &
            i,j,k

        INTRINSIC ::  DMIN1

        f=zero 
        IF (objfunc==1) THEN        ! objective 1 
        DO i=1,nrecord1
            f2=large
            DO k=1,nc
                CALL hauxfunc(x,k,i,d2)
                f2=DMIN1(f2,d2)
            END DO
            f=f+f2
        END DO
        nfit=nfit+1

        ELSE                                        ! objective 2, 3 or 4 
           CALL clrfunc1(x,f1)
           CALL clrfunc2(x,f2)
           f = f1-f2
        END IF
        RETURN

    END SUBROUTINE clrfunc
    
    !=============================================================
    ! Computation of f1 for clusterwise linear regression problem
    !=============================================================
    SUBROUTINE clrfunc1(x,f)

        USE param, ONLY : zero, large
        USE initclr, ONLY : &
            maxdim, &      ! Maximum number of variables in optimization,
                           !   maxdim >= nclust * nft.
            objc, &        ! The parameter used in objective function objfunc=2, 3 and 4.
            a4, &          ! Data matrix a4(maxrec,maxnft).
            kvad, &        ! Scaling vector for kvadratic terms
            nft, &         ! Number of used features.
            nc, &          ! Current number of clusters, loops from 1 to nclust.
            nrecord1, &    ! Number of records in data.
            objfunc, &     ! Objective function used.
            objeps, &      ! The parameter used in objective function objfunc=2, 3 and 4.   
            objkvad, &     ! The parameter used in objective function objfunc=2, 3 and 4.           
            nfit           ! The number of function evaluations.

        IMPLICIT NONE

        ! Subroutine arguments
        REAL(KIND=prec) ::  &
            x(maxdim), &
            f

        ! Local variables
        REAL(KIND=prec) ::  &
            fa,fb,d1,d2,smallest
        INTEGER :: &
            i,j,k

        INTRINSIC ::  DMIN1

        f=zero
        
        fa = 0.0_prec
		
		IF (objkvad > (10.0_prec)**(-10)) THEN
          DO k=1,nc
            DO i=1,nft-1
              fa = fa +objkvad*kvad(i)*x(i+(k-1)*nft)*x(i+(k-1)*nft)
            END DO
          END DO
		END IF
		
        fa = 0.5_prec * fa

        fb = 0.0_prec   
        
       IF (objfunc==1 .OR. objfunc==2 .OR. objfunc==3) THEN      ! objective 1, 2 or 3 
         DO i=1,nrecord1
            DO k=1,nc
                CALL hauxfunc(x,k,i,d2)
                fb=fb+d2
            END DO
         END DO
         
       ELSE IF (objfunc==4) THEN                                 ! objective 4
         DO i=1,nrecord1
            smallest = large
            DO k=1,nc
                CALL hauxfunc_obj4(x,k,i,d2)
                fb=fb+d2
                IF (smallest > d2) THEN
                   smallest = d2
                END IF
            END DO
            IF (smallest < objeps) THEN 
              fb = fb - smallest + objeps
            END IF
         END DO
       END IF
        
        IF (objfunc==1) THEN           ! objective 1
          f = fb
        ELSE                           ! objective 2, 3 or 4
          f=fa+objc*fb
        END IF 
        nfit=nfit+1
        RETURN

    END SUBROUTINE clrfunc1 
    
    !=============================================================
    ! Computation of f2 for clusterwise linear regression problem
    !=============================================================
    SUBROUTINE clrfunc2(x,f)

        USE param, ONLY : zero, large
        USE initclr, ONLY : &
            maxdim, &      ! Maximum number of variables in optimization,
                           !   maxdim >= nclust * nft.
            a4, &          ! Data matrix a4(maxrec,maxnft).
            objc, &        ! The parameter used in objective function objfunc=2, 3 and 4.
            nft, &         ! Number of used features.
            nc, &          ! Current number of clusters, loops from 1 to nclust.
            nrecord1, &    ! Number of records in data.
            objfunc, &     ! Objective function used.
            objeps, &      ! The parameter used in objective function objfunc=2, 3 and 4.   
            nfit           ! The number of function evaluations.

        IMPLICIT NONE

        ! Subroutine arguments
        REAL(KIND=prec) ::  &
            x(maxdim), &
            f

        ! Local variables
        REAL(KIND=prec) ::  &
            fsum,d1,d2,fmin, smallest
        INTEGER :: &
            i,j,k

        INTRINSIC ::  DMIN1

        f=zero
        
       IF (objfunc==1 .OR. objfunc==2 .OR. objfunc==3) THEN      ! objective 1, 2 or 3
        DO i=1,nrecord1
            fmin = large
            fsum = 0.0_prec
            DO k=1,nc
                CALL hauxfunc(x,k,i,d2)
                fsum = fsum + d2
                fmin = MIN(fmin, d2)
            END DO
            f = f + fsum - fmin
        END DO
       ELSE IF (objfunc==4) THEN                                 ! objective 4
         DO i=1,nrecord1
            smallest = large
            DO k=1,nc
                CALL hauxfunc_obj4(x,k,i,d2)
                f=f+d2
                IF (smallest > d2) THEN
                   smallest = d2
                END IF
            END DO
            f = f - smallest + objeps
         END DO
       END IF       
        
        IF (objfunc==1) THEN       ! objective 1
          f = f
        ELSE                       ! objective 2, 3, or 4
          f = objc * f
        END IF  
        !nfit=nfit+1
        RETURN

    END SUBROUTINE clrfunc2     
    
    
    
    !================================
    ! Computation of the subgradient
    !================================

    SUBROUTINE fgrad(x,grad)

        USE param, ONLY : zero, two, large
        USE initclr, ONLY : &
            maxdim, &      ! Maximum number of variables in optimization,
                           !   maxdim >= nclust * nft.
            a4, &          ! Data matrix a4(maxrec,maxnft).
            dminim, &      !
            nft, &         ! Number of used features.
            nc, &          ! Current number of clusters, loops from 1 to nclust.
            nrecord1, &    ! Number of records in data.
            objfunc, &     ! Objective function used.
            optmet, &      ! Optimization method.           
            ns, &          ! Switch for auxiliary and real clustering problem.
            m              ! Number of variables in optimization:
                           !   m = nft    if ns = 1,
                           !   m = nft*nc if ns = 2.

        IMPLICIT NONE

        ! Subroutine arguments
        REAL(KIND=prec) ::  &
            x(maxdim), &
            grad(maxdim), &
            grad1(maxdim), &
            grad2(maxdim) 

        ! Local variables
        REAL(KIND=prec) ::  &
            d1,d2,d4,f12
        INTEGER :: &
            i,j,k,jmin


        DO i=1,m
            grad(i)=zero
        END DO

      IF (objfunc==1) THEN  ! objective 1 
        IF (ns == 1) THEN   ! aux clusterwise linear regression problem
            DO i=1,nrecord1
                d1 = x(nft)
                DO j=1,nft-1
                    d1=d1+a4(i,j)*x(j)
                END DO
                d1 = d1-a4(i,nft)
                d2=d1*d1
                IF (d2 <= dminim(i)) THEN
                    DO j=1,nft-1
                        grad(j)=grad(j)+two*d1*a4(i,j)
                    END DO
                    grad(nft)=grad(nft)+two*d1
                END IF
            END DO
        ELSE                ! clusterwise linear regression problem

            DO i=1,nrecord1
                f12 = large
                DO k=1,nc
                    d1=x(k*nft)
                    DO j=1,nft-1
                        d1=d1+a4(i,j)*x(j+(k-1)*nft)
                    END DO
                    d2=(d1-a4(i,nft))**2
                    IF (f12 > d2) THEN
                        f12=d2
                        jmin=k
                        d4=d1-a4(i,nft)
                    END IF
                END DO
                DO j=1,nft-1
                    grad(j+(jmin-1)*nft)=grad(j+(jmin-1)*nft) + two*d4*a4(i,j)
                END DO
                grad(jmin*nft)=grad(jmin*nft)+two*d4

            END DO
        END IF
      ELSE   ! objective 2, 3 or 4 
         CALL fsubgrad1(x,grad1)
         CALL fsubgrad2(x,grad2)
         DO i =1, m
            grad(i) = grad1(i)-grad2(i)
         END DO
      END IF      

        RETURN

    END SUBROUTINE fgrad

    
    !===============================================================
    ! Computation of the subgradient (first DC component) 
    !===============================================================
    SUBROUTINE fsubgrad1(x,grad)

        USE param, ONLY : zero, two, large
        USE initclr, ONLY : &
            maxdim, &      ! Maximum number of variables in optimization,
                           !   maxdim >= nclust * nft.
            objc, &        ! The parameter used in objective function objfunc=2, 3 and 4.    
            a4, &          ! Data matrix a4(maxrec,maxnft).
            kvad, &        ! Scaling vector for kvadratic terms
            nft, &         ! Number of used features.
            nc, &          ! Current number of clusters, loops from 1 to nclust.
            nrecord1, &    ! Number of records in data.
            objfunc, &     ! Objective function used.
            dminim, &      ! dminim(nrecord)
            objeps, &      ! The parameter used in objective function objfunc=2, 3 and 4.    
            objkvad, &     ! The parameter used in objective function objfunc=2, 3 and 4.           
            ns, &          ! Switch for auxiliary and real clustering problem.
            m              ! Number of variables in optimization:
                           !   m = nft    if ns = 1,
                           !   m = nft*nc if ns = 2.

        IMPLICIT NONE

        ! Subroutine arguments
        REAL(KIND=prec) ::  &
            x(maxdim), &
            grad(maxdim)

        ! Local variables
        REAL(KIND=prec) ::  &
            d1, sub(nft),ylimsub(nft), smallest, d2,r
        INTEGER :: &
            i,j,k,ind, indapu 

        DO i=1,m
            grad(i)=zero
        END DO

        IF (ns == 1) THEN   ! aux clusterwise linear regression problem
              IF (objfunc==2 .OR. objfunc==3 .OR. objfunc==4) THEN    ! objective 2, 3 or 4
			    IF (objkvad > (10.0_prec)**(-10)) THEN 
                  DO j=1,nft-1
                    grad(j)=grad(j)+objkvad*kvad(j)* x(j)
                  END DO
                END IF    
              END IF
			   
            IF (objfunc==1 .OR. objfunc==2 .OR. objfunc==3) THEN      ! objective 1, 2 or 3  
              DO i=1,nrecord1
                CALL hauxsub(x,1,i,sub)
                DO j = 1, nft               
                  grad(j)=grad(j)+sub(j)
                END DO  
              END DO
            
            ELSE IF (objfunc==4) THEN                                 ! objective 4
               DO i=1,nrecord1         
                   CALL hauxfunc_obj4(x,1,i,d2)
                   CALL hauxsub_obj4(x,1,i,sub)
                   r = dminim(i)
                   IF ( (r+d2-objeps) > MAX(r,d2-objeps)) THEN
                     DO j = 1, nft              
                       grad(j)=grad(j)+sub(j)
                     END DO                       
                   ELSE
                      IF (r > d2 -objeps) THEN 
                         DO j = 1, nft              
                            grad(j)=grad(j)
                         END DO                      
                      ELSE
                        DO j = 1, nft               
                            grad(j)=grad(j)+ sub(j)
                        END DO                    
                      END IF
                   END IF
               END DO   
            END IF
 
        ELSE                ! clusterwise linear regeression problem
          IF (objfunc==2 .OR. objfunc==3 .OR. objfunc==4) THEN      ! objectives 2, 3 or 4
   	       IF (objkvad > (10.0_prec)**(-10)) THEN 
           	 DO k = 1,nc
               DO i = 1, nft-1
                 ind=i+(k-1)*nft
                 grad(ind)=grad(ind)+objkvad*kvad(i)*x(ind)
               END DO
             END DO
           END IF    
          END IF  
           
          IF (objfunc==1 .OR. objfunc==2 .OR. objfunc==3) THEN      ! objectives 1, 2 or 3
            DO i=1,nrecord1
                DO k=1,nc
                    CALL hauxsub(x,k,i,sub)
                    DO j=1,nft
                        ind=j+(k-1)*nft
                        IF (objfunc==1) THEN 
                          grad(ind)=grad(ind)+sub(j)                        
                        ELSE
                          grad(ind)=grad(ind)+objc*sub(j)
                        END IF  
                    END DO
                END DO
            END DO  
            
          ELSE IF (objfunc==4) THEN                                ! objective 4
            DO i=1,nrecord1
                smallest = large
                indapu = 0
                DO k=1,nc
                    CALL hauxfunc_obj4(x,k,i,d2)
                    CALL hauxsub_obj4(x,k,i,sub)
                    DO j=1,nft
                      ind=j+(k-1)*nft
                      grad(ind)=grad(ind)+objc*sub(j)   
                    END DO
                    IF (smallest>d2) THEN 
                       smallest = d2
                       indapu = k
                       ylimsub = sub                          
                    END IF                  
                END DO
                IF (objeps > smallest) THEN 
                    DO j=1,nft
                      ind=j+(indapu-1)*nft
                      grad(ind)=grad(ind)-objc*ylimsub(j)   
                    END DO                 
                END IF
            END DO  
          END IF          
        END IF

        RETURN

    END SUBROUTINE fsubgrad1    
    
    !===============================================================
    ! Computation of the subgradient (second DC component) 
    !===============================================================
    SUBROUTINE fsubgrad2(x,grad)

        USE param, ONLY : zero, large
        USE initclr, ONLY : &
            maxdim, &      ! Maximum number of variables in optimization,
                           !   maxdim >= nclust * nft.
            objc, &        ! The parameter used in objective function objfunc=2, 3 and 4..
            a4, &          ! Data matrix a4(maxrec,maxnft).
            nft, &         ! Number of used features.
            nc, &          ! Current number of clusters, loops from 1 to nclust.
            nrecord1, &    ! Number of records in data.
            objfunc, &     ! the objective function used.
            dminim, &      ! dminim(nrecord)
            objeps, &      ! The parameter used in objective function objfunc=2, 3 and 4.           
            ns, &          ! Switch for auxiliary and real clustering problem.
            m              ! Number of variables in optimization:
                           !   m = nft    if ns = 1,
                           !   m = nft*nc if ns = 2.

        IMPLICIT NONE

        ! Subroutine arguments
        REAL(KIND=prec) ::  &
            x(maxdim), &
            grad(maxdim)

        ! Local variables
        REAL(KIND=prec) ::  &
            d1, sub(nft), ff, r, &
            fmin, d2, smallest, ylim(nft)
        INTEGER :: &
            i,j,k,ind,minind, indapu

        DO i=1,m
            grad(i)=zero
        END DO

        IF (ns == 1) THEN   ! aux clusterwise linear regression problem
        
          IF (objfunc==1 .OR. objfunc==2 .OR. objfunc==3) THEN     ! objective 1, 2 or 3
            DO i=1,nrecord1
                CALL hauxfunc(x,1,i,ff)
                IF (ff>dminim(i)) THEN 
                  CALL hauxsub(x,1,i,sub)               
                  DO j = 1, nft             
                    grad(j)=grad(j)+sub(j)
                  END DO  
                END IF  
            END DO
            
          ELSE IF (objfunc==4) THEN          ! objective 4
               DO i=1,nrecord1
                   CALL hauxfunc_obj4(x,1,i,d2)
                   CALL hauxsub_obj4(x,1,i,sub)
                   r = dminim(i)
                   IF (d2-objeps >= r) THEN
                     DO j = 1, nft              
                       grad(j)=grad(j)+sub(j)
                     END DO                       
                   ELSE
                     DO j = 1, nft              
                       grad(j)=grad(j)
                     END DO                      
                   END IF
                END DO   
          END IF          
            
        ELSE                ! clusterwise linear regression problem
        
          IF(objfunc==1 .OR. objfunc==2 .OR. objfunc==3) THEN      ! objective 1, 2 or 3
            DO i=1,nrecord1
                fmin=large
                DO k=1,nc
                    CALL hauxfunc(x,k,i,ff)
                    CALL hauxsub(x,k,i,sub)
                    DO j=1,nft
                        ind=j+(k-1)*nft
                        IF (objfunc==1) THEN
                          grad(ind)=grad(ind)+sub(j)
                        ELSE 
                          grad(ind)=grad(ind)+objC*sub(j)
                        END IF  
                    END DO
                    IF (ff<fmin) THEN 
                       fmin=ff
                       minind=k
                    END IF 
                END DO
                CALL hauxsub(x,minind,i,sub)                
                DO j=1,nft
                   ind=j+(minind-1)*nft
                   IF (objfunc==1) THEN
                     grad(ind)=grad(ind)-sub(j)
                   ELSE  
                     grad(ind)=grad(ind)-objc*sub(j)
                   END IF                    
                END DO              
            END DO
          
          ELSE IF (objfunc==4) THEN          ! objective 4
            DO i=1,nrecord1
                smallest = large
                DO k=1,nc
                    CALL hauxfunc_obj4(x,k,i,d2)
                    CALL hauxsub_obj4(x,k,i,sub)
                    DO j=1,nft
                      ind=j+(k-1)*nft
                      grad(ind)=grad(ind)+objc*sub(j)   
                    END DO
                    IF (smallest>d2) THEN 
                       smallest = d2
                       indapu = k
                       ylim = sub                         
                    END IF                  
                END DO
                DO j=1,nft
                   ind=j+(indapu-1)*nft
                   grad(ind)=grad(ind)-objc*ylim(j) 
                END DO                 
            END DO      
          END IF            
            
        END IF

        RETURN

    END SUBROUTINE fsubgrad2        
    

    !===============================================================
    ! Computation of the subgradient (first DC component) objfunc=1   (TAITAA OLLA TURHA)
    !===============================================================
    SUBROUTINE fgrad1(x,grad)

        USE param, ONLY : zero, two
        USE initclr, ONLY : &
            maxdim, &      ! Maximum number of variables in optimization,
                           !   maxdim >= nclust * nft.
            a4, &          ! Data matrix a4(maxrec,maxnft).
            nft, &         ! Number of used features.
            nc, &          ! Current number of clusters, loops from 1 to nclust.
            nrecord1, &    ! Number of records in data.
            ns, &          ! Switch for auxiliary and real clustering problem.
            m              ! Number of variables in optimization:
                           !   m = nft    if ns = 1,
                           !   m = nft*nc if ns = 2.

        IMPLICIT NONE

        ! Subroutine arguments
        REAL(KIND=prec) ::  &
            x(maxdim), &
            grad(maxdim)

        ! Local variables
        REAL(KIND=prec) ::  &
            d1
        INTEGER :: &
            i,j,k

        DO i=1,m
            grad(i)=zero
        END DO

        IF (ns == 1) THEN
            DO i=1,nrecord1
                d1=x(nft)
                DO j=1,nft-1
                    d1=d1+a4(i,j)*x(j)
                END DO
                DO j=1,nft-1
                    grad(j)=grad(j)+two*(d1-a4(i,nft))*a4(i,j)
                END DO
                grad(nft)=grad(nft)+two*(d1-a4(i,nft))
            END DO
        ELSE
            DO i=1,nrecord1
                DO k=1,nc
                    d1=x(k*nft)
                    DO j=1,nft-1
                        d1=d1+a4(i,j)*x(j+(k-1)*nft)
                    END DO
                    d1=d1-a4(i,nft)
                    DO j=1,nft-1
                        grad(j+(k-1)*nft)=grad(j+(k-1)*nft)+two*d1*a4(i,j)
                    END DO
                    grad(k*nft)=grad(k*nft)+two*d1
                END DO
            END DO
        END IF

        RETURN

    END SUBROUTINE fgrad1
    
    
    !======================================================
    ! Computation of the subgradient (second DC component) objfunc=2   (TAITAA OLLA TURHA)
    !======================================================
    SUBROUTINE fgrad2(x,grad)

        USE param, ONLY : zero, two
        USE initclr, ONLY : &
            maxdim, &      ! Maximum number of variables in optimization,
                           !   maxdim >= nclust * nft.
            a4, &          ! Data matrix a4(maxrec,maxnft).
            nft, &         ! Number of used features.
            nc, &          ! Current number of clusters, loops from 1 to nclust.
            nrecord1, &    ! Number of records in data.
            dminim, &      ! dminim(nrecord).
            ns, &          ! Switch for auxiliary and real clustering problem.
            m              ! Number of variables in optimization:
                           !   m = mf    if ns = 1,
                           !   m = mf*nc if ns = 2.


        IMPLICIT NONE

        ! Subroutine arguments
        REAL(KIND=prec) ::  &
            x(maxdim), &
            grad(maxdim)

        ! Local variables
        REAL(KIND=prec), DIMENSION(nc) ::  &
            d2,d3          ! Auxiliary arrays.

        REAL(KIND=prec) ::  &
            d0,d1,d4,d5
        INTEGER :: &
            i,j,k,jindex

        DO i=1,m
            grad(i)=zero
        END DO

        IF (ns == 1) THEN
            PRINT*,'ns=1'
            DO i=1,nrecord1
                d1=x(nft)
                DO j=1,nft-1
                    d1=d1+a4(i,j)*x(j)
                END DO
                d0=(d1-a4(i,nft))**2
                IF (d0 >= dminim(i)) THEN
                    DO j=1,nft-1
                        grad(j)=grad(j)+two*(d1-a4(i,nft))*a4(i,j)
                    END DO
                    grad(nft)=grad(nft)+two*(d1-a4(i,nft))
                END IF
            END DO
        ELSE
            PRINT*,'ns=2'
            DO i=1,nrecord1
                DO k=1,nc
                    d3(k)=x(k*nft)
                    DO j=1,nft-1
                        d3(k)=d3(k)+a4(i,j)*x(j+(k-1)*nft)
                    END DO
                    d2(k)=(d3(k)-a4(i,nft))**2
                END DO
                d4=zero
                DO j=1,nc
                    d5=zero
                    DO k=1,nc
                        IF (k /= j) THEN
                            d5=d5+d2(k)
                        END IF
                    END DO
                    IF (d4 < d5) THEN
                        d4=d5
                        jindex=j
                    END IF
                END DO

                DO j=1,nc
                    IF (j /= jindex) THEN
                        DO k=1,nft-1
                            grad(k+(j-1)*nft)=grad(k+(j-1)*nft)+two*(d3(j)-a4(i,nft))*a4(i,k)
                        END DO
                        grad(j*nft)=grad(j*nft)+two*(d3(j)-a4(i,nft))
                    END IF
                END DO
            END DO
        END IF

        RETURN

    END SUBROUTINE fgrad2


    SUBROUTINE fvb(x,f)

        USE param, ONLY : zero, one
        USE initclr, ONLY : &
            maxdim, &      ! Maximum number of variables.
            objfunc, &     ! Objective function used.
            nft, &
            a4, &
            nel1,nob1, &
            nupdate,naux2

        ! Subroutine arguments
        REAL(KIND=prec) ::  &
            x(maxdim), &
            f

        ! Local variables
        REAL(KIND=prec) ::      &
            d1,d2, f1, f2
        INTEGER :: &
            j,k

        d2=zero
        IF (objfunc==1) THEN          ! objective 1  
          do k=1,nel1
            CALL hauxfunc(x,1,nob1(k),d1)
            d2 = d2 + d1
          end do
          f=d2
          nupdate=nupdate+1 
          naux2=naux2+1        
        ELSE                          ! objective 2, 3 or 4               
           CALL fvb1(x,f1) 
           f = f1
        END IF


        RETURN
    END SUBROUTINE fvb
    
    SUBROUTINE fvb1(x,f)

        USE param, ONLY : zero, one
        USE initclr, ONLY : &
            maxdim, &      ! Maximum number of variables.
            nft, &
            a4, &
            kvad, &        ! Scaling vector for kvadratic terms
            nel1,nob1, &
            objfunc, &     ! Objective function used
            objkvad, &     ! The parameter used in objective function objfunc=2, 3 or 4.
            nupdate,naux2

        ! Subroutine arguments
        REAL(KIND=prec) ::  &
            x(maxdim), &
            f

        ! Local variables
        REAL(KIND=prec) ::      &
            d1,d2
        INTEGER :: &
            j,k

        d2=zero
        
        DO k=1,nel1
            CALL hauxfunc(x,1,nob1(k),d1)
            d2 = d2 + d1
        END DO
        f = d2
        
        IF (objfunc==2 .OR. objfunc==3 .OR. objfunc==4) THEN     ! objectives 2, 3 or 4
	      IF (objkvad > (10.0_prec)**(-10)) THEN 
		    DO j = 1, nft-1
             f = f + objkvad*0.5_prec*kvad(j)*x(j)*x(j)
            END DO  
		  END IF	
        END IF
                
        nupdate=nupdate+1
        naux2=naux2+1

        RETURN
    END SUBROUTINE fvb1

    SUBROUTINE fvb2(x,f)

        USE param, ONLY : zero, one
        USE initclr, ONLY : &
            maxdim, &      ! Maximum number of variables.
            nft, &
            a4, &
            nel1,nob1, &
            objfunc, &     ! Objective function used
            nupdate,naux2

        ! Subroutine arguments
        REAL(KIND=prec) ::  &
            x(maxdim), &
            f

        ! Local variables
        REAL(KIND=prec) ::      &
            d1,d2
        INTEGER :: &
            j,k

        d2=zero
        
        f = zero

        RETURN
    END SUBROUTINE fvb2

    SUBROUTINE gradient(x,grad)

        USE param, ONLY : zero, two
        USE initclr, ONLY : &
            maxdim, &      ! Maximum number of variables.
            nft, &
            a4, &
            objfunc, &     ! Objective function used
            nel1,nob1

        ! Subroutine arguments
        REAL(KIND=prec) ::  &
            x(maxdim), &
            grad(maxdim), & 
            grad1(maxdim), & 
            grad2(maxdim) 

        ! Local variables
        REAL(KIND=prec) ::      &
            d1
        INTEGER :: &
            j,k

        DO j=1,nft
            grad(j)=zero
        END DO
        
        IF (objfunc==1) THEN                ! objective 1
        DO k=1,nel1
            d1=x(nft)
            DO j=1,nft-1
                d1=d1+a4(nob1(k),j)*x(j)
            END DO
            d1=d1-a4(nob1(k),nft)
            DO j=1,nft-1
                grad(j)=grad(j)+two*d1*a4(nob1(k),j)
            END DO
            grad(nft)=grad(nft)+two*d1
        END DO 
        
        ELSE                                ! objective 2, 3 or 4
            CALL gradient1(x,grad1)
            DO j =1 ,nft
              grad(j) = grad1(j)
            END DO  
        END IF 

        RETURN
    END SUBROUTINE gradient
    
    SUBROUTINE gradient1(x,grad)

        USE param, ONLY : zero, two
        USE initclr, ONLY : &
            maxdim, &       ! Maximum number of variables.
            objfunc, &      ! Objective function used.
            objkvad, &      ! Scaling vector for kvadratic terms
            nft, &
            a4, &
            kvad, &
            nel1,nob1

        ! Subroutine arguments
        REAL(KIND=prec) ::  &
            x(maxdim), &
            grad(maxdim)

        ! Local variables
        REAL(KIND=prec) ::      &
            d1, sub(nft), grad1(maxdim)
        INTEGER :: &
            j,k

        DO j=1,nft
            grad(j)=zero
        END DO

             IF (objfunc==2 .OR. objfunc==3 .OR. objfunc==4) THEN    ! objective 2, 3 or 4
               IF (objkvad > (10.0_prec)**(-10)) THEN 
 			     DO j=1,nft-1
                    grad(j)=grad(j)+objkvad*kvad(j)* x(j)
                 END DO
			   END IF 
             END IF  
            
            DO k=1,nel1
              
              CALL hauxsub(x,1,nob1(k),sub)
              DO j=1,nft
                grad(j)=grad(j)+sub(j)
              END DO
              
             END DO

        RETURN
    END SUBROUTINE gradient1

    SUBROUTINE gradient2(x,grad)

        USE param, ONLY : zero, two
        USE initclr, ONLY : &
            maxdim, &      ! Maximum number of variables.
            nft, &
            a4, &
            nel1,nob1

        ! Subroutine arguments
        REAL(KIND=prec) ::  &
            x(maxdim), &
            grad(maxdim)

        ! Local variables
        REAL(KIND=prec) ::      &
            d1
        INTEGER :: &
            j,k

        DO j=1,nft
            grad(j)=zero
        END DO


        RETURN
    END SUBROUTINE gradient2    


END MODULE clrfunctionmod
