Alicja Smoktunowicz

Jesse L. Barlow, Alicja Smoktunowicz

A new reorthogonalized block classical Gram--Schmidt algorithm is proposed that factorizes a full column rank matrix $A$ into $A=QR$ where $Q$ is left orthogonal (has orthonormal columns) and $R$ is upper triangular and nonsingular. With appropriate assumptions on the diagonal blocks of $R$, the algorithm, when implemented in floating point arithmetic with machine unit $\macheps$, produces $Q$ and $R$ such that $\| I- Q^{T} Q \|_2 =O(\macheps)$ and $\| A-QR \|_2 =O(\macheps \| A \|_2)$. The resulting bounds also improve a previous bound by Giraud et al. [Num. Math., 101(1):87-100,\ 2005] on the CGS2 algorithm originally developed by Abdelmalek [BIT, 11(4):354--367,\ 1971]. \medskip Keywords: Block matrices, Q--R factorization, Gram-Schmidt process, Condition numbers, Rounding error analysis.

Alicja Smoktunowicz, Jakub Kierzkowski, Iwona Wrobel

Stability analysis of Wilkinson's iterative refinement with a relaxation IR(omega) for solving linear systems is given. It extends existing results for omega=1, i.e., for Wilkinson's iterative refinement. We assume that all computations are performed in fixed (working) precision arithmetic. Numerical tests were done in MATLAB to illustrate our theoretical results. A particular emphasis is given on convergence of iterative refinement with a relaxation. Our tests confirm that the choice omega=1 is the best choice from the point of numerical stability.

Gianluca Oderda, Alicja Smoktunowicz, Ryszard Kozera

A real quadratic matrix is generalized doubly stochastic (g.d.s.) if all of its row sums and column sums equal one. We propose numerically stable methods for generating such matrices having possibly orthogonality property or/and satisfying Yang-Baxter equation (YBE). Additionally, an inverse eigenvalue problem for finding orthogonal generalized doubly stochastic matrices with prescribed eigenvalues is solved here. The tests performed in \textsl{MATLAB} illustrate our proposed algorithms and demonstrate their useful numerical properties.

Alicja Smoktunowicz, Jesse L. Barlow, Julien Langou

An error analysis result is given for classical Gram--Schmidt factorization of a full rank matrix $A$ into $A=QR$ where $Q$ is left orthogonal (has orthonormal columns) and $R$ is upper triangular. The work presented here shows that the computed $R$ satisfies $\normal{R}=\normal{A}+E$ where $E$ is an appropriately small backward error, but only if the diagonals of $R$ are computed in a manner similar to Cholesky factorization of the normal equations matrix. A similar result is stated in [Giraud at al, Numer. Math. 101(1):87--100,2005]. However, for that result to hold, the diagonals of $R$ must be computed in the manner recommended in this work.

Alicja Smoktunowicz, Iwona Wróbel

It is known that Goertzel's algorithm is much less numerically accurate than the Fast Fourier Transform (FFT)(Cf. \cite{gen:69}). In order to improve accuracy we propose modifications of both Goertzel's and Horner's algorithms based on the divide-and-conquer techniques. The proof of the numerical stability of these two modified algorithms is given. The numerical tests in Matlab demonstrate the computational advantages of the proposed modifications. The appendix contains the proof of numerical stability of Goertzel's algorithm of polynomial evaluation.

Alicja Smoktunowicz, Przemyslaw Kosowski, Iwona Wrobel

The scheme of divided differences is widely used in many approximation and interpolation problems. Computing the Newton coefficients of the interpolating polynomial is the first step of the Björck and Pereyra algorithm for solving Vandermonde systems of equations (Cf. \cite{bjorck: 70}). Very often this algorithm produces very accurate solution. The problem of determining the Newton coefficients is intimately related with the problem of evaluation the Lagrange interpolating polynomial, which can be realized by many algorithms. For these reasons we use the uniform approach and analyze also Aitken's algorithm of the evaluation of an interpolating polynomial. We propose new algorithms that are always numerically stable with respect to perturbation in the function values and more accurate than the Aitken's algorithm and the scheme of divided differences, even for complex data.