tridiagonal matrix inversion

[nurulhudaismail]

May I know anybody here kindly can share tridiagonal matrix inversion?
Thank you.

[salem c]

Try "tridiagonal matrix inversion c++" in your favourite search engine.

[lastchance]

#include <iostream>
#include <vector>
#include <cmath>     // EDITED, to define abs() properly
using namespace std;

const double SMALL = 1.0E-30;          // used to stop divide-by-zero

//======================================================================

bool tdma( const vector<double> &a, const vector<double> &b, const vector<double> &c, const vector<double> &d, 
           vector<double> &X )
//*********************************************************************************
// Solve, using the Thomas Algorithm (TDMA), the tri-diagonal system              *
//     a_i X_i-1 + b_i X_i + c_i X_i+1 = d_i,     i = 0, n - 1                    *
//                                                                                *
// Effectively, this is the n x n matrix equation.                                *
// a[i], b[i], c[i] are the non-zero diagonals of the matrix and d[i] is the rhs. *
// a[0] and c[n-1] aren't used.                                                   *
//*********************************************************************************
{
   int n = d.size();
   vector<double> P( n, 0 );
   vector<double> Q( n, 0 );
   X = P;

   // Forward pass
   int i = 0;
   double denominator = b[i];
   P[i] = -c[i] / denominator;
   Q[i] =  d[i] / denominator;
   for ( i = 1; i < n; i++ )
   {
      denominator = b[i] + a[i] * P[i-1];
      if ( abs( denominator ) < SMALL ) return false;
      P[i] =  -c[i]                   / denominator;
      Q[i] = ( d[i] - a[i] * Q[i-1] ) / denominator;
   }

   // Backward pass
   X[n-1] = Q[n-1];
   for ( i = n - 2; i >= 0; i-- ) X[i] = P[i] * X[i+1] + Q[i];

   return true;
}

//======================================================================

int main()              // Solve Ax = b, where A is tri-diagonal
{
   vector<double> lower    = {  0.0, -1.0, -2.0, -3.0, -4.0 };    // lower diagonal
   vector<double> diagonal = {  9.0,  8.0,  7.0,  6.0,  5.0 };    // main diagonal
   vector<double> upper    = { -1.0, -2.0, -3.0, -4.0,  0.0 };    // upper diagonal
   vector<double> rhs      = { 41.0, 21.0,  7.0, -1.0, -3.0 };    // RHS
   vector<double> X;

   int n = rhs.size();

   if ( tdma( lower, diagonal, upper, rhs, X ) )
   {
      cout << "X\tLHS\tRHS\n";
      for ( int i = 0; i < n; i++ )
      {
         double lhs = diagonal[i] * X[i];
         if ( i > 0     ) lhs += lower[i] * X[i-1];
         if ( i < n - 1 ) lhs += upper[i] * X[i+1];
         cout << X[i] << '\t' << lhs << '\t' << rhs[i] << '\n';
      }
   }
   else
   {
      cout << "Unable to solve\n";
   }
}

X	LHS	RHS
5	41	41
4	21	21
3	7	7
2	-1	-1
1	-3	-3

[nurulhudaismail]

Thank you. But when I run it using C++ Shell, there is an error which I don't understand as below:

In function 'bool tdma(const std::vector<double>&, const std::vector<double>&, const std::vector<double>&, const std::vector<double>&, std::vector<double>&)':
33:29: error: call of overloaded 'abs(double&)' is ambiguous
33:29: note: candidates are:
In file included from /usr/include/c++/4.9/cstdlib:72:0,
from /usr/include/c++/4.9/ext/string_conversions.h:41,
from /usr/include/c++/4.9/bits/basic_string.h:2850,
from /usr/include/c++/4.9/string:52,
from /usr/include/c++/4.9/bits/locale_classes.h:40,
from /usr/include/c++/4.9/bits/ios_base.h:41,
from /usr/include/c++/4.9/ios:42,
from /usr/include/c++/4.9/ostream:38,
from /usr/include/c++/4.9/iostream:39,
from 1:
/usr/include/stdlib.h:775:12: note: int abs(int)
extern int abs (int __x) __THROW __attribute__ ((__const__)) __wur;
^
In file included from /usr/include/c++/4.9/ext/string_conversions.h:41:0,
from /usr/include/c++/4.9/bits/basic_string.h:2850,
from /usr/include/c++/4.9/string:52,
from /usr/include/c++/4.9/bits/locale_classes.h:40,
from /usr/include/c++/4.9/bits/ios_base.h:41,
from /usr/include/c++/4.9/ios:42,
from /usr/include/c++/4.9/ostream:38,
from /usr/include/c++/4.9/iostream:39,
from 1:
/usr/include/c++/4.9/cstdlib:174:3: note: long long int std::abs(long long int)
abs(long long __x) { return __builtin_llabs (__x); }
^
/usr/include/c++/4.9/cstdlib:166:3: note: long int std::abs(long int)
abs(long __i) { return __builtin_labs(__i); }
^

[nurulhudaismail]

One more thing, if we look into the code, the matrix is written manually. What if the matrix is generated matrix, then, I extract the tridiagonal components to be a tridiagonal matrix where off-tridiagonal components are zero. Then, I compute the inversion of this tridiagonal matrix?

Any response is most appreciated. Thank you.

[lastchance]

The problem with "abs" should be resolved if you include the header cmath; i.e.
#include <cmath>
at the top. My compiler must have pulled in this header, or at least the requisite abs() bit, with one of the other headers. EDIT: I have amended the code to fix this.


It is unusual to want the actual inverse of a tridiagonal matrix; you usually just want to solve a single tridiagonal system. If you want the actual inverse then simply apply the TDMA to each of the columns of the identity matrix in turn; they will form the columns of the inverse matrix. Alternative, just do Gauss-Jordan elimination simultaneously on the original matrix and identity matrix.

[nurulhudaismail]

thank you lastchance. Basically, in my system:
1. I generate a full-rank complex matrix, A.

2. Then, I need to extract the tridiagonal components of A and put it in a full-rank matrix which consists of only tridiagonal components while the rest are zero value, called as matrix B.

3. Finally, I need to compute the inversion of B.

Is it possible to do this process in C++?

[lastchance]

Is it possible to do this process in C++?

Sure.

Store the non-zero diagonals of the matrix as rank-1 vectors, as in my code - there is no point in storing a lot of zeroes.

Amend the tridiagonal routine so that the "Backward pass" is done for each basis vector, i.e. the columns of the identity matrix. The resulting columns can be reassembled as the inverse matrix and returned instead of the rank-1 vector X.

Over to you.

[nurulhudaismail]

Thank you again. I will try to do it first.

[nurulhudaismail]

Dear lastchance(2887),

I realize that this technique will give me vector x where x = inv(A)b. However, what I am interested most is to compute the inversion of A.

Is anyone here knows how to compute the inversion of a tridiagonal matrix without using a linear equation, Ax=b?

Thanks in advance.