Solve a CppAD Sparse Lower Triangular System

solve_lower_cppad

@(@\newcommand{\B}[1]{{\bf #1}} \newcommand{\R}[1]{{\rm #1}}@)@Solve a CppAD Sparse Lower Triangular System
Syntax

#include "utility.hpp"

msg = solve_lower_cppad(A, x, b)

Prototype


template <class Scalar> std::string
solve_lower_cppad(
     const CPPAD_SPARSE(Scalar)&  A ,
     VECTOR(Scalar)&              x ,
     const VECTOR(Scalar)&        b )

A
Is the a square lower triangular matrix in row-major order. We use n to denote the number for rows and columns in A ; i.e., n is equal to both A.nr() and A.nc() .

x
This vector has size n and satisfies the equation



     A * x = b

b
This is a vector with size n .

msg
If the return value msg is empty, no error has occurred. Otherwise, it is one of the following



     "A is not in row-major order"

     "A is not lower triangular"

     "A is not invertible"

.

Example


bool test_solve_lower_cppad(void)
{    bool ok = true;
     double eps99 = 99.0 * std::numeric_limits<double>::epsilon();

     //     [ 1 0 0 ]
     // A = [ 2 3 0 ]
     //     [ 4 5 6 ]
     size_t nr  = 3;
     size_t nc  = 3;
     size_t nnz = 6;
     CppAD::sparse_rc< VECTOR(size_t) > pattern(nr, nc, nnz);
     size_t k = 0;
     for(size_t i = 0; i < nr; i++)
     {    for(size_t j = 0; j <= i; j++)
               pattern.set(k++, i, j);
     }
     assert( k == nnz );
     CPPAD_SPARSE(double) A(pattern);
     for(k = 0; k < nnz; k++)
          A.set(k, double(k+1));
     //
     // b = [1, 8, 32]^T
     VECTOR(double) b(nr);
     b[0] = 1.0;
     b[1] = 8.0;
     b[2] = 32.0;
     //
     // solve A * x = b
     VECTOR(double) x(nr);
     std::string msg = solve_lower_cppad<double>(A, x, b);
     //
     // check result
     ok &= msg == "";
     ok &= CppAD::NearEqual(x[0], 1.0, eps99, eps99);
     ok &= CppAD::NearEqual(x[1], 2.0, eps99, eps99);
     ok &= CppAD::NearEqual(x[2], 3.0, eps99, eps99);
     //
     return ok;
}

Input File: src/utility.hpp