// file: $isip/doc/examples/class/math/matrix/math_matrix_example_01/example.cc
// version: $Id: example.cc 5955 2000-12-17 18:40:27Z hamaker $
//

// isip include files
//
#include <VectorFloat.h>
#include <MatrixFloat.h>

// main program starts here:
//  this example implements a sequence of matrix operations known
//  as a quadratic form: scalar = x_t * u_t * e * u * x
//
int main() {
  
  // declare the vector x
  //
  VectorFloat x;
  x.assign(L"1.0, 2.0, 3.0");
  
  // declare the required matrices. matrix u is declared to be a FULL
  // matrix (the default type), which means all possible matrix
  // locations have storage available. matrix e is DIAGONAL, which
  // means only elements along the diagonal are allowed to be
  // non-zero. an attempt to set an off-diagonal element will result
  // in an error, so use the non-FULL types with care.
  //
  MatrixFloat u;
  u.assign(3, 3, L"1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0",
	   Integral::FULL);
  
  MatrixFloat e;
  e.assign(3, 3, L"1.0, 2.0, 1.0", Integral::DIAGONAL);
  
  // compute the transpose of u.
  //
  MatrixFloat u_t;
  u_t.transpose(u);
  
  // compute the inner product: u_t * e * u
  //
  MatrixFloat w;
  w.mult(u_t, e);
  w.mult(u);
  
  // compute the outer product: x_t * w * x
  //  the vmult function multiples a vector to a matrix and stores the
  //  result in a matrix. mathematically speaking, the input vector is
  //  really a 1xn matrix and the output vector is really a nx1
  //  matrix.
  //
  VectorFloat y;
  w.vmult(y, x);
  Float result = y.dotProduct(x);
  
  // output the result
  //
  x.debug(L"x = ");
  Console::put(L"");
  u.debug(L"u = ");
  Console::put(L"");
  e.debug(L"e = ");
  Console::put(L"");
  result.debug(L"x_t * u_t * e * u * x = ");

  // exit gracefully
  //
  Integral::exit();
}