Formulation of Matrix Pade Approximation in Rectangular Full Packed Storage

The Extended Euclidean algorithm for matrix Pade approximants is applied to compute matrix Pade approximants in rectangular full packed format (RFP) if the coefficient matrices of the input matrix polynomial are triangular. The procedure given by Gustavson et al for packing a triangular matrix in rectangular full packed format is applied to pack sequence of lower triangular matrices of a matrix polynomial in Rectangular Full Packed format. This RFP format of a matrix polynomial is applied to compute matrix Pade approximants of the matrix polynomial using Matrix Pade Extended Euclidean Algorithm. Algorithms for the multiplication of two triangular matrices and inverse of a triangular matrix in RFP format are also presented. The CPU time and memory comparison in computing the matrix Pade approximants of a matrix polynomial between RFP format case and non packed case are elucidated in detail.


Introduction
Pade approximants have been the subject of much recent interest in many fields of applications (Baker, G.A. Jr., 1975).Some of the applications of this efficient technique is summation of infinite series, solution of integral equations, numerical inversion of Laplace transforms etc. Work on Pade or Pade-type approximants ultimately involves the explicit determination of the polynomials forming the numerator and denominator of rational functions can be found in ( Gragg, W.B., 1972;Baker, G.A. Jr.,1975;Robert J. McEliece, James B. Shearer, 1978;Baker, G.A. Jr. and Graves -Morris, P.R.,1981).There exists in literature the use of extended Euclidean algorithm for matrix Pade approximants.Several researchers have worked on the problem of finding simple algorithm for computing matrix Pade approximants from a formal power series of matrix coefficients ( Rissanen, J., 1972 ;Bultheel, A.,1980;Achuthan, P. and Sundar, S., 1988).So it is natural for us to apply these techniques for computing matrix Pade approximants in a special case where the matrices are of very large order lower triangular type.For the sake of completeness, the definition of matrix Pade approximants is presented in Section 2. Section 3 deals with the procedure for a lower triangular matrix in rectangular full packed format and packing of a sequence of lower triangular matrices in rectangular full packed format.Section 4 deals with the lower triangular inverse in rectangular full packed format.The multiplication of two lower triangular matrices in rectangular full packed format is presented in Section 5. Section 6 is concerned with the algorithm for computing matrix Pade approximants in rectangular full packed format.Section 7 brings out CPU time and memory comparison in computing the matrix Pade approximants.The subroutines used in C language are presented in Section 8.The system configuration is given in appendix.The work presented in this paper is an advanced one that of presented in (Achuthan, P. and Sundar, S.,1988 , P. 287-296).
Throughout this paper the following symbols are used : [ f (x)] → f (x) is the polynomial in which the coefficients are square matrices of the same order n; I n → Identity matrix of order n; ∅ n → Null matrix of order n; [RFP − g(x)] → g (x) is the polynomial in which the coefficients are triangular matrices of the same order n and each coefficient matrices are in rectangular full packed format and also all the rectangular full packed matrices are packed according to their degree.

The matrix Pade approximants
The [M/N] matrix Pade approximant of a formal power series where (s n ) i 's are coefficients of the power series which are square matrices of the same order n, is defined as the rational function The numerator and denominator of this matrix Pade approximant are where [P M (x)] and [Q N (x)] are polynomials of degree M and N at most respectively.[P M (x)]/[Q N (x)] is called the normal matrix Pade approximant of order (M, N) and is symbolically denoted by (M/N) [S ] (x).We may define two types of Pade approximants, since the (s n ) i 's need no longer commute, by the equations However we can show that these two Pade approximants are actually identical.
For different chosen values of M(≥ 0) and N(≥ 0) we can construct the Pade table in which the distinct approximants of the function [S (x)] would be the elements.There arise mainly two problems: one is the coefficient problem, in which it is required to find the numerator and denominator coefficient matrices, (p n ) i and (q n ) i , of equations ( 3) and (4) for the input (s n ) i of equation (1); the other is the value problem, which is to determine the value of the desired approximant for an explicit input chosen for x.The algorithm to find matrix Pade approximants by using Extended Euclidean algorithm was given in (Achuthan, P. and Sundar, S.,1988 , P. 287-296).

Lower triangular matrix in rectangular full packed format
Symmetric or triangular matrix may be stored in rectangular full packed format.The advantage of storing symmetric or triangular matrix in Rectangular full packed format is explained in (Gustavson, G. and Jerzy Wasniewsky., 2007, P. 570 -579).We presented here the method for storing a triangular matrix in Rectangular full packed format.Break the lower triangular matrix A into A 11 , A 21 and A 22 as three sub matrices as shown in figure 1.
When n is even , the lower triangular sub matrix of A 11 and the upper triangular sub matrix of A T 22 can be concatenated along their main diagonals into a (k + 1) × k dense matrix.Then the square sub matrix A 21 of order k × k is appended below the dense matrix (k + 1) × k , where k = n/2 .Thus the lower triangular matrix A can be stored as a (n + 1) × k rectangular matrix.When n is odd, the lower triangular sub matrix A 11 is concatenated with the upper triangular matrix A T 22 along their main diagonals into a k × k dense matrix.Then the rectangular sub matrix A 21 of order (k − 1) × k is appended below the k × k dense matrix.Thus A is stored as a n × k rectangular matrix.The above procedure is presented in figure 2 As in Section 3 a triangular matrix of order n can be stored as a (n + 1) × k rectangular matrix if 'n' is even and n × k rectangular matrix if 'n' is odd in rectangular full packed format.If coefficient matrices of a matrix polynomial [S (x)] are triangular then each coefficient matrices of [S (x)] packed as Rectangular full packed format one by one according to their degree.Thus obtained Rectangular full packed matrices are sequentially packed according to the degree of the matrix polynomial [S (x)] as shown in figure 3. Now the resulting matrix is also rectangular.If the number of coefficient matrices of [S (x)] is m, the resulting rectangular matrix must be of order (n + 1) × (mk), if n is even and n × (mk) if n is odd.

Triangular inverse in rectangular full packed format
There are many algorithms in literature to find the inverse of a lower triangular matrix.Let the inverse of lower triangular matrix A be D.
We find inverse of A 11 , A T 22 using the algorithm given below separately and developed an algorithm to compute the inverse of A in rectangular full packed format.for j = k : −1 : 1 In order to find the inverse of A T 22 , the above algorithm modified to compute row by row in the reverse order.

Multiplications of two triangular matrixes in rectangular full packed format
As the lower triangular matrix A and B are in 2− by −2 block form as shown in figure 5, from C = AB, we obtain three block equations Using the three block equations we presented an algorithm which multiplies A and B and returns the resulting matrix C in rectangular full packed format., 1980).The extended Euclidean Algorithm to matrix Pade was presented by P.Achuthan and S.Sundar in (Achuthan, P. and Sundar, S.,1988, P. 287-296) for square matrix coefficients.
We packed the triangular matrix coefficients of the matrix polynomials as explained in Section 3 and applied this packed matrix polynomials in MAtrix Pade Extended Euclidean Algorithm(MAPEA) in (Achuthan, P. and Sundar, S.,1988, P. 287-296).This algorithm gives a sequences of anti-diagonal approximants starting from (M + N, 0) to (M, N) in the Pade OUT PUT: Sequence of anti -diagonal approximants starting from (M + N, 0) to (M, N) if all the approximants exists.
Function MAPEARFP ([S(x)], degree S, N, M, n) EXIT 1: In sufficient degree for the input matrix polynomial [S (x)], so we cannot find approximants.EXIT 2: As the Pade denominator degree is zero , [b(x)] itself is the required Pade approximant.EXIT 3: As the inverse of b n does not exist, Pade approximant cannot be found.

CPU time
We have computed the matrix Pade approximants for the Semi -normal power Series ¢ www.ccsenet.org/jmr in (Leighton, W. and Scott,W.T., 1939) using the algorithms MAPEA in (Achuthan, P. and Sundar, S.,1988 , P.287-296) and MAPEARFP for various matrix orders and various Pade orders when x = 1.We have considered identity matrices in (7) as triangular matrices while using MAPEARFP algorithm.
CPU time to compute (7/7) matrix Pade approximants for the semi-normal power series (7) at x = 1 for various order of square matrices are presented in table 1 also the CPU time to compute (7/7) matrix Pade approximants for the semi-normal power series (7) at x = 1 for various order of triangular matrices in RFP format are presented in table 2.
While comparing the CPU time to compute matrix Pade approximants for the square matrix case using the algorithm MAPEA and for the triangular matrix case using the algorithm MAPEARFP, it is clear that MAPEARFP algorithm is very fast while using higher order triangular matrix coefficients of the matrix polynomial.

Memory comparison
The memory needed to store a matrix polynomial is depends on its matrix coefficients and degree.Hence the memory needed to store a matrix polynomial is (Number of coefficients matrices of the polynomial) (Memory for storing a coefficient matrix) + (Memory for storing exponents of the polynomial).The difference in memory size to compute (M / N) matrix Pade approximants for square matrix case and triangular matrix case in rectangular full packed storage are presented below.

7.2.1
The memory required to compute (M/N) matrix Pade approximants for square matrix case.
To get the (M/N) matrix Pade approximants, degree of the input series [S (x)] must be M + N + 1, hence the number of coefficient matrices of [S (x)] are at most M + N + 2. So the memory needed to store [S (x)] in square matrix case is Similarly, the memory needed to store various matrix polynomials used to compute (M/N) matrix Pade approximants are presented in the table 3.In addition various subroutines are used to compute the matrix Pade approximants, the memory required to execute the subroutines are presented in table 4.
Total memory we required to compute (M/N) matrix Pade approximants for square matrix case is

7.2.2
The memory required to compute the matrix Pade approximants for the triangular matrix case in rectangular full packed storage.
The memory required to store a triangular matrix of order 'n' in rectangular full packed format is n(n + 1)/2.To get the (M/N) matrix Pade approximants , the degree of the input the series [S (x)] must be M + N + 1 , hence the number of coefficient matrices of [S (x)] are at most M + N + 2 .So the memory needed to store [S (x)] in Rectangular full packed format case is (M + N + 2)(n(n + 1)/2) + M + N + 2. Similarly, the memory needed to store various matrix polynomials used to compute (M/N) matrix Pade approximants in rectangular full packed format are presented in the table 5.In addition we allocated 4(M+N +2)n 2 +4(M+N +2) storage space to store [S (x)], [a(x)], [dummy1(x)] and [dummy2(x)] before packing.
To compute matrix Pade approximants various subroutines are used, the memory required to execute the subroutines are presented in the following table 6.
The total memory we required to compute (M/N) matrix Pade approximants in rectangular full packed storage Hence the difference in memory size to compute (M/N) Matrix Pade approximants between square matrix and triangular matrix cases is

Subroutines used in C
The following subroutines were used to calculate matrix Pade approximants in rectangular full packed format: 1 Rect pack Matrix multiply : Returns the result of multiplication of two triangular matrices in rectangular full packed format. 2 Rect pack Matrix inversion : Returns the inverse of a triangular matrix in rectangular full packed format.3 Rect pack : Returns the rectangular full packed format of a triangular matrix 4 Rect polynomial packing : Returns the packed form of the coefficient matrices of polynomial matrix 5 Rect packpoly multiply : Returns the result of multiplication of two matrix polynomials in rectangular full packed format.6 Rect packpoly addition : Returns the result of addition of two matrix polynomials in rectangular full packed format.7 Rect packpoly subtraction : Returns the result of subtraction of two matrix polynomials in rectangular full packed format.8 Rect packpoly division : Returns the quotient and remainder of the two matrix polynomials in rectangular full packed format.

CONCLUDING REMARKS
We have shown here how to construct normal matrix Pade approximants if the coefficients matrices of the input matrix polynomial are triangular matrix by using MAtrix Pade Extended Euclidean Algorithm in Rectangular Full Packed format (MAPEARFP).If the coefficients of the input matrix polynomial are triangular, MAPEARFP algorithm is very useful in getting matrix Pade approximants instead of using square matrix case algorithm MAPEA.This packed algorithm is very fast while using higher order triangular matrix coefficients of the matrix polynomial.Its memory usage is also reduced nearly half compared to the square matrix case.

APPENDIX
The configuration of the system used to find the CPU time to compute matrix Pade approximants for the semi normal power series ( 7) is presented below.

System configuration:
Operating
Since A and D are in 2− by −2 blocking as shown in figure 4, from the identity AD = DA = I, we get three block equations A 11 D 11 = I, A 21 D 11 + A 22 D 21 = 0, A 22 D 22 = I.The above three block equations implies that

Figure 2 .Figure 3 .
Figure 2. Rectangular full packed format of a lower triangular matrix A when n = 6, 7.

Figure 4
Figure 4. 2− by −2 Block form of A and D.

Figure 8 .
Figure 8. Graphical representation for the comparison of Table1 and Table 2.
Algorithm to find inverse of a lower triangular matrix in RFP form Name of the Algorithm: RFPTI (Rectangular Full Packed Triangular Inverse) INPUT:Lower triangular matrix A in Rectangular full packed format.OUT PUT: Inverse of lower triangular matrix A in Rectangular full packed format. 1 D 11 ← In(A 11 ) % find inverse of A 11 .D 21 ← mul(A 21 , D 11 ) % multiplication of A 21 and D 11 4 D 21 ← (−1)mul(D 22 , D 21 ) % multiplication of D 22 and D 21 .
To find the inverse of matrix A 11 of order k, we used the following algorithm in (Jeremy J Du Croz and Nicholas, Higham, J., 1992, P. 1-19) which computes the inverse of A 11 , column by column in the reverse order.Let the inverse of A 11 be X.Algorithm to find inverse of A 11 INPUT:Lower triangular matrix A 11 of order k.OUT PUT: Inverse of lower triangular matrix A 11 of order k.
Name of the algorithm: RFPTRTRM (Rectangular Full Packed Triangular matrix and Triangular matrix Multiplication) INPUT:Two lower triangular matrices A and B in Rectangular full packed format.OUT PUT: Matrix C in Rectangular full packed format 1 C 11 ←LTLTM(A 11 , B 11 ) % multiplication of the lower triangular matrices A 11 and B 11 C 21 ←mul(A 21 , B 21 ) % multiplication of the matrices A 21 and B 11 4 C 21 ← C 21 +mul(A 22 , B 21 ) % multiplication of the matrices A 22 and B 21 and addition with C 21 Note 1: LTLTM = Lower Triangular matrix Lower Triangular matrix Multiplication.UTUTM = Upper Triangular matrix Upper Triangular matrix Multiplication 6. Algorithm for computing matrix Pade approximants in Rectangular Full Packed format (Aho, A.V., Hopcroft, J.E., and Ullman, J.D., 1974)own algorithm, originally suggested by Aho et al(Aho, A.V., Hopcroft, J.E., and Ullman, J.D., 1974).An application of EEA to ordinary non-matrix Pade was presented byMcEliece  and Sherarer (Robert J.McEliece, James B. Shearer., 1978, P. 611-615)and Brent et al(Brent, R.P., Gustavson,F.G., and  Yun, D.Y.Y. table if all the approximants exists.We presented here a new algorithm MAPEARFP (MAtrix Pade Extended Euclidean Algorithm in Rectangular Full Packed format) to compute matrix Pade approximants if coefficient matrices of a matrix polynomial are triangular.Name of the algorithm: MAPEARFP (MAtrix Pade Extended Euclidean Algorithm in Rectangular Full Packed format) INPUT: n, N, M, [S (x)], n is order of the matrix, M is Pade numerator, N is Pade denominator, [S (x)] is the Series whose coefficients are triangular matrices.

Table 6 .
Memory allocation of various subroutines that are required to compute (M/N) matrix Pade approximants in RFP format.