/*                                                 */

/*      WHITTLE'S ALGORITHM                        */

/*                                                 */

PROC IML;

RESET NOLOG NOCENTER;

/* ******************************* */

/*   MODULE OF WHITTLE ALGORITHM   */

/* ******************************* */

START WHITTLE(C, D, MAXORDER);

 /*  CONSTANTS   */

 IDEN = I(D);

 AROW = D;

 ACOL = D*(MAXORDER+1);

 L = J(AROW, ACOL, 0);

 M = L;

 G = J(D, D, 0);

 LAMBDA = G;

 V = G;

 AROWS = 1:D;

 ACOLS = 1:D;

 /*  INITIAL STEP  K=0  */

 L[AROWS, ACOLS] = -IDEN;

 M[AROWS, ACOLS] = -IDEN;

 ACOLSNXT = (D+1):(2*D);

 G = C[AROWS, ACOLSNXT];

 LAMBDA = C[AROWS, ACOLS];

 V = LAMBDA;

 /*  ITERATIVE STEP  */

 DO K = 1 TO MAXORDER;

    LOLD = L;

    MOLD = M;

    L[AROWS, ACOLS] = -IDEN;

    M[AROWS, ACOLS] = -IDEN;

    ACOLSLST = (K*D+1):(K+1)*D;

    L[AROWS, ACOLSLST] = G*INV(V);

    M[AROWS, ACOLSLST] = G`*INV(LAMBDA);

    DO J = 1 TO K-1;

       ACOLSJ = (J*D+1):(J+1)*D;

       ACOLSKMJ = ((K-J)*D+1):((K-J+1)*D);

       L[AROWS, ACOLSJ] = LOLD[AROWS, ACOLSJ]

         - L[AROWS, ACOLSLST]*MOLD[AROWS, ACOLSKMJ];

       M[AROWS, ACOLSJ] = MOLD[AROWS, ACOLSJ]

         - M[AROWS, ACOLSLST]*LOLD[AROWS, ACOLSKMJ];

       END;

    ACOLSNXT = ((K+1)*D+1):(K+2)*D;

    G = C[AROWS, ACOLSNXT];

    DO J = 1 TO K;

       ACOLSJ = (J*D+1):((J+1)*D);

       ACOLSKMJ = ((K-J+1)*D+1):((K-J+2)*D);

       G = G-L[AROWS, ACOLSKMJ]*C[AROWS, ACOLSJ];

       END;

    LAMBDA = (IDEN -

        L[AROWS, ACOLSLST]*M[AROWS, ACOLSLST])*LAMBDA;

    V = (IDEN -

        M[AROWS, ACOLSLST]*L[AROWS, ACOLSLST])*V;

    /*  PRINT THE SOLUTION  */

    PRINT,, "ORDER =  "  K,,

              LAMBDA  "  " V,;

    DO J = 1 TO K;

       ACOLSJ = (J*D+1):(J+1)*D;

       LKJ = L[AROWS, ACOLSJ];

       MKJ = M[AROWS, ACOLSJ];

       PRINT "LAG = " J,

               LKJ "  "  MKJ, ;

       END;

    END;

 FINISH;