From: kovarik@mcmail.cis.McMaster.CA (Zdislav V. Kovarik) Subject: Re: periodic tridiagonal solver Date: 17 Jan 2000 20:09:57 -0500 Newsgroups: sci.math.num-analysis Summary: [missing] In article , Tony Keating wrote: : :Could someone please post some pseudo-code for a periodic tridiagonal :solver? I've got a standard tridiagonal solver written, but I'm having :difficulties with implementing periodic boundary conditions. Not quite a pseudocode but explanation of the structure. I presume you are talking about matrices of the form (rigid spacing on the screen assumed) [** *] [*** ] [ *** ] [ *** ] [ *** ] [ *** ] [ *** ] [ ***] [* **] and you want to know how to handle the NE and SW positions. This pattern of non-zeros is a spacial case of the "arrow pattern" [** *] [*** *] [ *** *] [ *** *] [ *** *] [ *** *] [ ****] [ ***] [*********] If it is safe to perform LU-decomposition without pivoting (e.g. if the matrix is positive definite), the factors have sparsity patterns [* ] [** ] [ ** ] [ ** ] [ ** ] [ ** ] [ ** ] [ ** ] [*********] and [** *] [ ** *] [ ** *] [ ** *] [ ** *] [ ** *] [ ***] [ **] [ *] so forward and back substitution are about twice as costly as for a "clean" tridiagonal matrix. The storage can be economized accordingly. The "fill-in" seems to be inevitable, unless someone came up with a neat trick how to get rid of it. Good luck, ZVK(Slavek). ============================================================================== From: Christof Vomel Subject: Re: periodic tridiagonal solver Date: Tue, 18 Jan 2000 11:04:36 +0100 Newsgroups: sci.math.num-analysis I just wanted to add, that for a system [a* -b] [*** ] [ *** ] [ *** ] [ *** ] [ *** ] [ *** ] [ ***] [b *a] you will perform a Givens rotation first. ============================================================================== From: elsner@mathematik.tu-chemnitz.de (Ulrich Elsner) Subject: Re: periodic tridiagonal solver Date: 18 Jan 2000 16:33:28 GMT Newsgroups: sci.math.num-analysis According to Tony Keating : > >Could someone please post some pseudo-code for a periodic tridiagonal >solver? I've got a standard tridiagonal solver written, but I'm having >difficulties with implementing periodic boundary conditions. Your matrix looks like: xx x xxx xxx xxx xxx x xx You don't metion what your 'stadard' solver is. If it is LU-decomposition without pivoting, you will just get one row and column of fill-in, i.e.: xx x xxx * xxx * xxx* xxx x****xx For a more general method (if you use some other method for the triangular system) I'll quote from an old posting of mine: Your matrix is an (easy to solve) tridiagonal matrix disturbed by a low rank matrix. So one general way is the use of the Sherman-Morrison- Woodbury formula (see e.g. Golub, van Loan: Matrix Computations): (A+UV^T)^{-1}=A^{-1} - A^{-1} U (I + V^T A^{-1} U)^{-1} V^T A^{-1} You write your matrix as A+UV^T where A is easy to solve (i.e. tridiagonal) plus a low rank update, here U*V^T = (1 0 0 ... 0) * (0 0 0 .. x) (0 0 0 ... 1)^T (x 0 0 ... 0) a rank 2 modification (U,V are 2 by n matrices). This leaves the tridiagonal part of your matrix undisturbed or U*V^T = (1 0 0 ... 0 1)^T * (x 0 0 ... 0 x) a rank 1 modification (this changes the (1,1) and the (n,n) element of the tridiagonal part which might be a problem if e.g. your tridiagonal part is positive definite). So, instead of solving one linear system with the complicated matrix, you solve 2+rank system with the (easy to solve) matrix A and the small system in the inner bracket (size = rank of U*V^T). (Note: of course you do NOT invert these matrices but only solve the linear systems). Hope that helps, Ulrich -- Ulrich Elsner, Fak. fuer Mathematik, TU Chemnitz, D-09107 Chemnitz, Germany