> with(linalg):
> A := matrix(5,5,[0,1,1,0,0,
> 0,0,0,1,1, 0,0,0,1,1, 0,0,0,0,1,
> 0,0,0,0,0]);

                          [0    1    1    0    0]
                          [                     ]
                          [0    0    0    1    1]
                          [                     ]
                     A := [0    0    0    1    1]
                          [                     ]
                          [0    0    0    0    1]
                          [                     ]
                          [0    0    0    0    0]

> evalm(A + A^2 + A^3 + A^4);

                       [0    1    1    2    4]
                       [                     ]
                       [0    0    0    1    2]
                       [                     ]
                       [0    0    0    1    2]
                       [                     ]
                       [0    0    0    0    1]
                       [                     ]
                       [0    0    0    0    0]

> A := matrix(2,3);

                    A := array(1 .. 2, 1 .. 3, [])

> B := matrix(3,4);

                    B := array(1 .. 3, 1 .. 4, [])

> LS := evalm(transpose(A &* B));

  LS :=

        [A[1, 1] B[1, 1] + A[1, 2] B[2, 1] + A[1, 3] B[3, 1] ,

        A[2, 1] B[1, 1] + A[2, 2] B[2, 1] + A[2, 3] B[3, 1]]

        [A[1, 1] B[1, 2] + A[1, 2] B[2, 2] + A[1, 3] B[3, 2] ,

        A[2, 1] B[1, 2] + A[2, 2] B[2, 2] + A[2, 3] B[3, 2]]

        [A[1, 1] B[1, 3] + A[1, 2] B[2, 3] + A[1, 3] B[3, 3] ,

        A[2, 1] B[1, 3] + A[2, 2] B[2, 3] + A[2, 3] B[3, 3]]

        [A[1, 1] B[1, 4] + A[1, 2] B[2, 4] + A[1, 3] B[3, 4] ,

        A[2, 1] B[1, 4] + A[2, 2] B[2, 4] + A[2, 3] B[3, 4]]

> RS := evalm(transpose(B) &* transpose(A));

  RS :=

        [A[1, 1] B[1, 1] + A[1, 2] B[2, 1] + A[1, 3] B[3, 1] ,

        A[2, 1] B[1, 1] + A[2, 2] B[2, 1] + A[2, 3] B[3, 1]]

        [A[1, 1] B[1, 2] + A[1, 2] B[2, 2] + A[1, 3] B[3, 2] ,

        A[2, 1] B[1, 2] + A[2, 2] B[2, 2] + A[2, 3] B[3, 2]]

        [A[1, 1] B[1, 3] + A[1, 2] B[2, 3] + A[1, 3] B[3, 3] ,

        A[2, 1] B[1, 3] + A[2, 2] B[2, 3] + A[2, 3] B[3, 3]]

        [A[1, 1] B[1, 4] + A[1, 2] B[2, 4] + A[1, 3] B[3, 4] ,

        A[2, 1] B[1, 4] + A[2, 2] B[2, 4] + A[2, 3] B[3, 4]]

> equal(map(expand,LS), map(expand,RS));
> 

                                 true

> I2 := array(identity, 1..2, 1..2);

              I2 := array(identity, 1 .. 2, 1 .. 2, [])

> print(I2);

                               [1    0]
                               [      ]
                               [0    1]

> evalm(I2 &* A);

                   [A[1, 1]    A[1, 2]    A[1, 3]]
                   [                             ]
                   [A[2, 1]    A[2, 2]    A[2, 3]]

> e3 := vector([0,0,1]);

                           e3 := [0, 0, 1]

> c3 := evalm(A &* e3);

                       c3 := [A[1, 3], A[2, 3]]

> evalm(c3 &* A);

  [
  [A[1, 3] A[1, 1] + A[2, 3] A[2, 1],

                                                  2          2]
        A[1, 3] A[1, 2] + A[2, 3] A[2, 2], A[1, 3]  + A[2, 3] ]

> C := matrix(3,1,[a,b,c]);

                                    [a]
                                    [ ]
                               C := [b]
                                    [ ]
                                    [c]

> B := transpose(A);

                           [A[1, 1]    A[2, 1]]
                           [                  ]
                      B := [A[1, 2]    A[2, 2]]
                           [                  ]
                           [A[1, 3]    A[2, 3]]

> blockmatrix(1,2,[B,C]);

                      [A[1, 1]    A[2, 1]    a]
                      [                       ]
                      [A[1, 2]    A[2, 2]    b]
                      [                       ]
                      [A[1, 3]    A[2, 3]    c]

> ?diagonal
> D3 := diag(alpha,beta,gamma);

                          [alpha     0        0  ]
                          [                      ]
                    D3 := [  0      beta      0  ]
                          [                      ]
                          [  0       0      gamma]

> evalm(A &* D3);

           [A[1, 1] alpha    A[1, 2] beta    A[1, 3] gamma]
           [                                              ]
           [A[2, 1] alpha    A[2, 2] beta    A[2, 3] gamma]

> evalm(diag(2,3) &* A);

                [2 A[1, 1]    2 A[1, 2]    2 A[1, 3]]
                [                                   ]
                [3 A[2, 1]    3 A[2, 2]    3 A[2, 3]]

> T := matrix(2,2,[0,1,1,1]);

                                 [0    1]
                            T := [      ]
                                 [1    1]

> Tinv := evalm(T^(-1));

                                  [-1    1]
                          Tinv := [       ]
                                  [ 1    0]

> evalm(T &* Tinv);

                               [1    0]
                               [      ]
                               [0    1]

> evalm(Tinv &* T);

                               [1    0]
                               [      ]
                               [0    1]

> T2 := matrix(2,2,[0,0,0,1]);

                                  [0    0]
                            T2 := [      ]
                                  [0    1]

> inverse(T2);
Error, (in inverse) singular matrix
>