• Matrix::value_type and VecX::value_type must be the same type
• the multiplication operator must be defined for Matrix::value_type
• the division operator must be defined for Matrix::value_type
• the addition operator must be defined for Matrix::value_type

  // Select a lower triangular matrix with packed storage
  typedef matrix< double, 
                  triangle<lower>,
                   packed<>, 
                   row_major >::type Matrix;
  // The external vector uses memory from an "outside" source,
  // in this case vector b below will use memory from array db
  typedef external_vec<double> Vector;
  const int N = 3;

  Matrix A(N);
  set_diagonal(A, 1);

  //Fill the matrix...
  A(0,0) = 1;
  A(1,0) = 2;  A(1,1) = 4;
  A(2,0) = 3;  A(2,1) = 5; A(2,2) = 7;

  double db[] = { 7, 46, 124};
  Vector b(db, N);

  cout << "A in packed form" << endl;
  print_row(A);// print the matrix in row-wise fashion

  cout << "b:" << endl;
  print_vector(b);

  // This is also known as backward or forward substitution
  // In this case, since A is a lower triangular matrix
  // so this call performs a forward substitution
  // computes: A^-1 * b -> b
  tri_solve(A, b);

  cout << "A^-1 * b:" << endl;

  print_vector(b);