//  array2d.edp
//
//  Discussion:
//
//    Demonstrate how 2D matrices can be defined and manipulated.
//
//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/freefem++/array2d/array2d.edp
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    02 July 2015
//
//  Author:
//
//    John Burkardt
//
cout << "\n";
cout << "array2d:\n";
cout << "  FreeFem++ version\n";
cout << "  Demonstrate how 2D matrices can be defined and manipulated.\n";

load "lapack"
//
//  Define a 3x3 matrix A 
//
//  [ 0 1 0 ]
//  [ 0 0 1 ]
//  [ 1 1 1 ]
//
real[int,int] A(3,3);

A = 1.0; 
A(0,0) = 0.0;
A(0,2) = 0.0;
A(1,0) = 0.0;
A(1,1) = 0.0;
cout << "\n";
cout << "  Matrix A = " << A << endl;
//
//  Define a 3x3 matrix B:
//
//  [ 3 1 1 ]
//  [ 1 1 3 ]
//  [ 3 3 3 ]
//
real[int,int] B(3,3);
//B = 3.0; 
//B(0,1) = 1.0;
//B(0,2) = 1.0;
//B(1,0) = 1.0;
//B(1,0) = 1.0;
B = [ [ 3.0, 1.0, 1.0 ], [ 1.0, 1.0, 3.0 ], [ 3.0, 3.0, 3.0 ] ];
cout << "\n";
cout << "  Matrix B = " << B << endl;
//
//  Compute A*B.
//  This operation requires the "load lapack" statement above.
//
real[int,int] AxB(3,3);
AxB = A * B;
cout << "\n";
cout << "  Matrix AxB = product A*B = " << AxB << endl;
//
//  Compute A transpose.
//
real[int,int] AT(3,3);
AT = A';
cout << "\n";
cout << "  Matrix AT = A' = transpose A'  " << AT << endl;
//
//  Define a matrix C.

//
//  C = 1 1 0 1
//      1 1 1 1
//      1 0 1 1
//
real[int,int] C(3,4);
C = 1.0;
C(2,1) = 0.0;
C(0,2) = 0.0;
cout << "\n";
cout << "  Matrix C = " << C << endl;
cout << "  First row of C = C(0,:) = " << C(0,:) << endl;
//
//  Define a vector d, and copy row 0 of C into it.
//
real[int] d(4);
d = C(0,:);
cout << "\n";
cout << "  d = C(0,:)\n";
cout << "  Size of d = d.n = " << d.n << "\n";
cout << "  Entries of d = " << d << "\n";
//
//  Define vectors e and f, and create a matrix G as their outer product.
//
real[int] e(3);
e=[1,2,3];
real[int] f(4);
f=[4,5,6,7];
real[int,int] G(3,4);
G = e * f';
cout << "\n";
cout << "  Vector e = " << e << "\n";
cout << "  Vector f = " << f << "\n";
cout << "  Matrix G = e * f' = " << G << "\n";
//
//  Define vectors h and i, and use them to pre- and post-multiply
//  the 3x4 matrix C.
//
real[int] h(3);
h = [ 1.0, 1.0, 1.0 ];
real[int] i(4);
i = [ 1.0, 1.0, 1.0, 1.0 ];
real[int] j(4);
j = h' * C;
real[int] k(3);
k = C * i;
cout << "\n";
cout << "  Premultiply j = h' * C = " << j << "\n";
cout << "  Postmultiply k = C * i = " << k << "\n";
//
//  (Rectangular) Matrix multiplication.
//
real[int,int] L(3,4);
L = A * G;
cout << "\n";
cout << "  L = A * G = " << L << endl;
//
//  Matrix addition.
//
real[int,int] M(3,3);
M = A + B;
cout << "\n";
cout << "  M = A + B = " << M << endl;
//
//  Scalar multiplication.
//
M = 2 * M;
cout << "  M = 2 * M = " << M << endl;
//
//  Scalar addition.
//
//M = M .+ 10;
//cout << "  M = M + 10 = " << M << endl;
//
//  Elementwise multiplication.
//
M = A .* B;
cout << "  M = A .* B = " << M << endl;
//
//  Elementwise division
//
M = A ./ B;
cout << "  M = A ./ B = " << M << endl;
//
//  Inverse.
//
M = A ^ -1;
cout << "  M = A ^ -1 = " << M << endl;
//
//  Dimensions.
//
cout << "  Dimensions of A = (A.m,A.n) = " << A.m << "," << A.n << "\n";
//
//  Solution operator is multiplication by inverse:
//
real[int] x1(3);
real[int] x2(3);
real[int] y2(3);
real[int,int] Ainv(3,3);
x1 = [ 1, 2, 3];
cout << "\n";
cout << "  Solution X1 = " << x1 << "\n";
y2 = A * x1;
cout << "  Y2 = A * X1 = " << y2 << "\n";
//x2 = ( A ^ -1 ) * y2;
Ainv = ( A ^ - 1 );
x2 = Ainv * y2;
cout << "  X2 = A^-1 * Y2 = " << x2 << "\n";
//
//  Maximum, doesn't work!
//
//cout << "  max ( M ) = " << max ( M ) << "\n";
//
//  Norms?  Doesn't work.
//
//cout << || M ||_1 << "\n";
//
//  Terminate.
//
cout << "\n";
cout << "array2d:\n";
cout << "  Normal end of execution.\n";

exit ( 0 );