Efficient Band Matrix Computations

Listing 4: Solving an SPD linear system

#include <iostream.h>
#include <vector.h>
#include "spdband.h"
#include "spdbandf.h"

int main()
{
    const int n = 9;
    vector<float> diag(n, 4.0f), superDiag(n-1, -1.0f);
    vector<float> zeros(n-2, 0.0f), offdiag(n-3, -1.0f);

    SPDBandMatrix<float> B(n);
    B.upperBandWidth() = 3;
    B.lowerBandWidth() = -3;

    B.diagonal(0) = diag;
    B.diagonal(1) = superDiag;
    B.diagonal(2) = zeros;
    B.diagonal(3) = offdiag;

    cout << "Calling factor ..." << endl;
    SPDBandMatrixFactor<float> choleskyFactor(n);
    choleskyFactor.factor(B);

    cout << "Solving ..." << endl;
    vector<float> x = choleskyFactor.solve(diag);
    for (int i=0; i<n; i++)
        cout << x[i] << endl;

    return 0;
}
//End of File

OUTPUT

Calling factor ...
Solving ...
3.81928
4.96386
5.39759
6.31325
6.63855
6.31325
5.39759
4.96386
3.81928