/* W. Everling, "On the Evaluation of e**AT by power series",   */
/* Proc. IEEE, Vol. 55, 1967, p. 413.                           */

EXPA: PROCEDURE (A, ERROR);
        /* ERROR is relative precision. */
        declare A(*,*) float (18), ERROR FLOAT (18),
                (I, I1, I2, N) fixed binary;
        declare ((A_H, A_I, A_J)(HBOUND(A,1),HBOUND(A,2)),
                  ANORM, ATEST, ETEST ) float (18);

        /* Notation corresponds to that of Liou as follows:     */
        /* A      AT                                            */
        /* I      k                                             */
        /* I1     K                                             */
        /* A_H    auxiliary storage for argument                */
        /* A_I    (AT)**k/k!                                    */
        /* ATEST  10**d * |r(i,j)|                              */

        ON ERROR SNAP SYSTEM;

        N = DIM(A,1);

        ANORM = RSNORM(A);
        I1 = MAX(2*ANORM, 2); I2 = I1/2;
        ANORM = 2*ANORM / ERROR;
        A_H, A_I = A;
        /* Form I + A */
        DO I = 1 TO N;
                A(I,I) = A(I,I) + 1;
        END;

        /* Loop to form A**2/2! + A**3/3! + A**4/4! + ... */
        DO I = 2 BY 1;
                A_J = A_H/I;
                CALL MATMULT (A_I, A_J);
                A = A + A_I;
                IF I = I1 THEN
                        DO;
                                I1 = I1 + I2;
                                ATEST = RSNORM(A_I)*ANORM/(I+1);
                                ETEST = RSNORM(A)*ERROR**2;
                                IF ALL(ATEST < ABS(A) | ABS(A) < ETEST) THEN
                                        GO TO EXPA_END;
                        END;
        END;

EXPA_END:

END EXPA;

/* Forms the maximum_row_sum norm of A. */

RSNORM: PROCEDURE (A) RETURNS (FLOAT (18));
        DECLARE A(*,*) FLOAT (18), ANORM FLOAT (18) INITIAL (0),
                I FIXED BINARY;