********************************************************************************
*   Copyright(C) 2004-2014 Intel Corporation. All Rights Reserved.
*   
*   The source code, information  and  material ("Material") contained herein is
*   owned  by Intel Corporation or its suppliers or licensors, and title to such
*   Material remains  with Intel Corporation  or its suppliers or licensors. The
*   Material  contains proprietary information  of  Intel or  its  suppliers and
*   licensors. The  Material is protected by worldwide copyright laws and treaty
*   provisions. No  part  of  the  Material  may  be  used,  copied, reproduced,
*   modified, published, uploaded, posted, transmitted, distributed or disclosed
*   in any way  without Intel's  prior  express written  permission. No  license
*   under  any patent, copyright  or  other intellectual property rights  in the
*   Material  is  granted  to  or  conferred  upon  you,  either  expressly,  by
*   implication, inducement,  estoppel or  otherwise.  Any  license  under  such
*   intellectual  property  rights must  be express  and  approved  by  Intel in
*   writing.
*   
*   *Third Party trademarks are the property of their respective owners.
*   
*   Unless otherwise  agreed  by Intel  in writing, you may not remove  or alter
*   this  notice or  any other notice embedded  in Materials by Intel or Intel's
*   suppliers or licensors in any way.
*
********************************************************************************
*   Content : TR Solver Fortran-77 example
*
********************************************************************************

!** NONLINEAR LEAST SQUARE PROBLEM WITHOUT BOUNDARY CONSTRAINTS
      PROGRAM EXAMPLE_DTRNLSP_POWELL
        IMPLICIT NONE
!** HEADER-FILE WITH DEFINITIONS (CONSTANTS, EXTERNALS)
        INCLUDE "mkl_rci.fi"
!** USER’S OBJECTIVE FUNCTION
        EXTERNAL            EXTENDED_POWELL
!** N - NUMBER OF FUNCTION VARIABLES
        INTEGER             N
        PARAMETER           (N = 40)
!** M - DIMENSION OF FUNCTION VALUE
        INTEGER             M
        PARAMETER           (M = 40)
!** SOLUTION VECTOR. CONTAINS VALUES X FOR F(X)
        DOUBLE PRECISION    X (N)
!** PRECISIONS FOR STOP-CRITERIA (SEE MANUAL FOR MORE DETAILS)
        DOUBLE PRECISION    EPS (6)
!** JACOBI CALCULATION PRECISION
        DOUBLE PRECISION    JAC_EPS
!** REVERSE COMMUNICATION INTERFACE PARAMETER
        INTEGER             RCI_REQUEST
!** FUNCTION (F(X)) VALUE VECTOR
        DOUBLE PRECISION    FVEC (M)
!** JACOBI MATRIX
        DOUBLE PRECISION    FJAC (M, N)
!** NUMBER OF ITERATIONS
        INTEGER             ITER
!** NUMBER OF STOP-CRITERION
        INTEGER             ST_CR
!** CONTROLS OF RCI CYCLE
        INTEGER             SUCCESSFUL
!** MAXIMUM NUMBER OF ITERATIONS
        INTEGER             ITER1
!** MAXIMUM NUMBER OF ITERATIONS OF CALCULATION OF TRIAL-STEP
        INTEGER             ITER2
!** INITIAL STEP BOUND
        DOUBLE PRECISION    RS
!** INITIAL AND FINAL RESIDUALS
        DOUBLE PRECISION    R1, R2
!** TR SOLVER HANDLE
        INTEGER*8           HANDLE
!** CYCLE’S COUNTERS
        INTEGER             I, J
!** RESULTS OF INPUT PARAMETER CHECKING
        INTEGER             INFO(6)
!** SET PRECISIONS FOR STOP-CRITERIA
        DO I = 1, 6
            EPS (I) = 1.D-5
        END DO
!** SET MAXIMUM NUMBER OF ITERATIONS
        ITER1 = 1000
!** SET MAXIMUM NUMBER OF ITERATIONS OF CALCULATION OF TRIAL-STEP
        ITER2 = 100
!** SET INITIAL STEP BOUND
        RS = 100.D0
!** PRECISIONS FOR JACOBI CALCULATION
        JAC_EPS = 1.D-8
!** SET THE INITIAL GUESS
        DO I = 1, N/4
            X (4*I - 3) =  3.D0
            X (4*I - 2) = -1.D0
            X (4*I - 1) =  0.D0
            X (4*I)     =  1.D0
        END DO
!** SET INITIAL VALUES
        DO I = 1, M
            FVEC (I) = 0.D0
            DO J = 1, N
                FJAC (I, J) = 0.D0
            END DO
        END DO
!** INITIALIZE SOLVER (ALLOCATE MEMORY, SET INITIAL VALUES)
!**   HANDLE    IN/OUT: TR SOLVER HANDLE
!**   N         IN:     NUMBER OF FUNCTION VARIABLES
!**   M         IN:     DIMENSION OF FUNCTION VALUE
!**   X         IN:     SOLUTION VECTOR. CONTAINS VALUES X FOR F(X)
!**   EPS       IN:     PRECISIONS FOR STOP-CRITERIA
!**   ITER1     IN:     MAXIMUM NUMBER OF ITERATIONS
!**   ITER2     IN:     MAXIMUM NUMBER OF ITERATIONS OF CALCULATION OF TRIAL-STEP
!**   RS        IN:     INITIAL STEP BOUND
        IF (DTRNLSP_INIT (HANDLE, N, M, X, EPS, ITER1, ITER2, RS)
     &      .NE. TR_SUCCESS) THEN
!** IF FUNCTION DOES NOT COMPLETE SUCCESSFULLY THEN PRINT ERROR MESSAGE
            PRINT *, '| ERROR IN DTRNLSP_INIT'
!** RELEASE INTERNAL MKL MEMORY THAT MIGHT BE USED FOR COMPUTATIONS.
!** NOTE: IT IS IMPORTANT TO CALL THE ROUTINE BELOW TO AVOID MEMORY LEAKS
!** UNLESS YOU DISABLE MKL MEMORY MANAGER
            CALL MKL_FREE_BUFFERS
!** AND STOP
            STOP 1
        END IF
!** CHECKS THE CORRECTNESS OF HANDLE AND ARRAYS CONTAINING JACOBIAN MATRIX, 
!** OBJECTIVE FUNCTION, LOWER AND UPPER BOUNDS, AND STOPPING CRITERIA.
        IF (DTRNLSP_CHECK (HANDLE, N, M, FJAC, FVEC, EPS, 
     &     INFO) .NE. TR_SUCCESS) THEN
!** IF FUNCTION DOES NOT COMPLETE SUCCESSFULLY THEN PRINT ERROR MESSAGE
            PRINT *, '| ERROR IN DTRNLSPBC_INIT'
!** RELEASE INTERNAL MKL MEMORY THAT MIGHT BE USED FOR COMPUTATIONS.
!** NOTE: IT IS IMPORTANT TO CALL THE ROUTINE BELOW TO AVOID MEMORY LEAKS
!** UNLESS YOU DISABLE MKL MEMORY MANAGER
            CALL MKL_FREE_BUFFERS
!** AND STOP
            STOP 1
        ELSE
!** THE HANDLE IS NOT VALID.
            IF( INFO(1) .NE. 0 .OR. 
!** THE FJAC ARRAY IS NOT VALID.
     &              INFO(2) .NE. 0 .OR. 
!** THE FVEC ARRAY IS NOT VALID.
     &              INFO(3) .NE. 0 .OR. 
!** THE EPS ARRAY IS NOT VALID.
     &              INFO(4) .NE. 0 ) THEN
                PRINT *, '| INPUT PARAMETERS ARE NOT VALID'
!** RELEASE INTERNAL MKL MEMORY THAT MIGHT BE USED FOR COMPUTATIONS.
!** NOTE: IT IS IMPORTANT TO CALL THE ROUTINE BELOW TO AVOID MEMORY LEAKS
!** UNLESS YOU DISABLE MKL MEMORY MANAGER
                CALL MKL_FREE_BUFFERS
!** AND STOP
                STOP 1
            END IF
        END IF
!** SET INITIAL RCI CYCLE VARIABLES
        RCI_REQUEST = 0
        SUCCESSFUL = 0
!** RCI CYCLE
        DO WHILE (SUCCESSFUL == 0)
!** CALL TR SOLVER
!**   HANDLE        IN/OUT: TR SOLVER HANDLE
!**   FVEC          IN:     VECTOR
!**   FJAC          IN:     JACOBI MATRIX
!**   RCI_REQUEST   IN/OUT: RETURN NUMBER WHICH DENOTE NEXT STEP FOR PERFORMING
            IF (DTRNLSP_SOLVE (HANDLE, FVEC, FJAC, RCI_REQUEST)
     &              .NE. TR_SUCCESS) THEN
!** IF FUNCTION DOES NOT COMPLETE SUCCESSFULLY THEN PRINT ERROR MESSAGE
                PRINT *, '| ERROR IN DTRNLSP_SOLVE'
!** RELEASE INTERNAL MKL MEMORY THAT MIGHT BE USED FOR COMPUTATIONS.
!** NOTE: IT IS IMPORTANT TO CALL THE ROUTINE BELOW TO AVOID MEMORY LEAKS
!** UNLESS YOU DISABLE MKL MEMORY MANAGER
                CALL MKL_FREE_BUFFERS
!** AND STOP
                STOP 1
            END IF
!** ACCORDING WITH RCI_REQUEST VALUE WE DO NEXT STEP
            SELECT CASE (RCI_REQUEST)
            CASE (-1, -2, -3, -4, -5, -6)
!**   STOP RCI CYCLE
                SUCCESSFUL = 1
            CASE (1)
!**   RECALCULATE FUNCTION VALUE
!**     M               IN:     DIMENSION OF FUNCTION VALUE
!**     N               IN:     NUMBER OF FUNCTION VARIABLES
!**     X               IN:     SOLUTION VECTOR
!**     FVEC            OUT:    FUNCTION VALUE F(X)
                CALL EXTENDED_POWELL (M, N, X, FVEC)
            CASE (2)
!**   COMPUTE JACOBI MATRIX
!**     EXTENDED_POWELL IN:     EXTERNAL OBJECTIVE FUNCTION
!**     N               IN:     NUMBER OF FUNCTION VARIABLES
!**     M               IN:     DIMENSION OF FUNCTION VALUE
!**     FJAC            OUT:    JACOBI MATRIX
!**     X               IN:     SOLUTION VECTOR
!**     JAC_EPS         IN:     JACOBI CALCULATION PRECISION
                IF (DJACOBI (EXTENDED_POWELL, N, M, FJAC, X, JAC_EPS)
     &                  .NE. TR_SUCCESS) THEN
!** IF FUNCTION DOES NOT COMPLETE SUCCESSFULLY THEN PRINT ERROR MESSAGE
                    PRINT *, '| ERROR IN DJACOBI'
!** RELEASE INTERNAL MKL MEMORY THAT MIGHT BE USED FOR COMPUTATIONS.
!** NOTE: IT IS IMPORTANT TO CALL THE ROUTINE BELOW TO AVOID MEMORY LEAKS
!** UNLESS YOU DISABLE MKL MEMORY MANAGER
                    CALL MKL_FREE_BUFFERS
!** AND STOP
                    STOP 1
                END IF
            ENDSELECT
        END DO
!** GET SOLUTION STATUSES
!**   HANDLE            IN: TR SOLVER HANDLE
!**   ITER              OUT: NUMBER OF ITERATIONS
!**   ST_CR             OUT: NUMBER OF STOP CRITERION
!**   R1                OUT: INITIAL RESIDUALS
!**   R2                OUT: FINAL RESIDUALS
        IF (DTRNLSP_GET (HANDLE, ITER, ST_CR, R1, R2)
     &          .NE. TR_SUCCESS) THEN
!** IF FUNCTION DOES NOT COMPLETE SUCCESSFULLY THEN PRINT ERROR MESSAGE
            PRINT *, '| ERROR IN DTRNLSP_GET'
!** RELEASE INTERNAL MKL MEMORY THAT MIGHT BE USED FOR COMPUTATIONS.
!** NOTE: IT IS IMPORTANT TO CALL THE ROUTINE BELOW TO AVOID MEMORY LEAKS
!** UNLESS YOU DISABLE MKL MEMORY MANAGER
            CALL MKL_FREE_BUFFERS
!** AND STOP
            STOP 1
        END IF
!** FREE HANDLE MEMORY
        IF (DTRNLSP_DELETE (HANDLE) .NE. TR_SUCCESS) THEN
!** IF FUNCTION DOES NOT COMPLETE SUCCESSFULLY THEN PRINT ERROR MESSAGE
            PRINT *, '| ERROR IN DTRNLSP_DELETE'
!** RELEASE INTERNAL MKL MEMORY THAT MIGHT BE USED FOR COMPUTATIONS.
!** NOTE: IT IS IMPORTANT TO CALL THE ROUTINE BELOW TO AVOID MEMORY LEAKS
!** UNLESS YOU DISABLE MKL MEMORY MANAGER
            CALL MKL_FREE_BUFFERS
!** AND STOP
            STOP 1
        END IF

!** RELEASE INTERNAL MKL MEMORY THAT MIGHT BE USED FOR COMPUTATIONS.
!** NOTE: IT IS IMPORTANT TO CALL THE ROUTINE BELOW TO AVOID MEMORY LEAKS
!** UNLESS YOU DISABLE MKL MEMORY MANAGER
        CALL MKL_FREE_BUFFERS
!** IF FINAL RESIDUAL LESS THEN REQUIRED PRECISION THEN PRINT PASS
        IF (R2 .LT. 1.D-5) THEN
            PRINT *, '|         DTRNLSP POWELL............PASS'
            STOP 0
!** ELSE PRINT FAILED
        ELSE
            PRINT *, '|         DTRNLSP POWELL............FAILED'
            STOP 1
        END IF
      END PROGRAM EXAMPLE_DTRNLSP_POWELL

!** ROUTINE FOR EXTENDED POWELL FUNCTION CALCULATION
!**   M     IN:     DIMENSION OF FUNCTION VALUE
!**   N     IN:     NUMBER OF FUNCTION VARIABLES
!**   X     IN:     VECTOR FOR FUNCTION CALCULATING
!**   F     OUT:    FUNCTION VALUE F(X)
      SUBROUTINE EXTENDED_POWELL (M, N, X, F)
        IMPLICIT NONE
        INTEGER M, N
        DOUBLE PRECISION X (*), F (*)
        INTEGER I

        DO I = 1, N/4
            F (4*I-3) = X(4*I - 3) + 10.D0 * X(4*I - 2)
            F (4*I-2) = 2.2360679774998D0 * (X(4*I-1) - X(4*I))
            F (4*I-1) = ( X(4*I-2) - 2.D0*X(4*I-1) )**2
            F (4*I)   = 3.1622776601684D0 * (X(4*I-3) - X(4*I))**2
        END DO

      END SUBROUTINE EXTENDED_POWELL