sparsematrix.rlg

on errcont;



sparse_matrix sd(5,5);



length sd;


{5,5}


sd;


[0  0  0  0  0]
[             ]
[0  0  0  0  0]
[             ]
[0  0  0  0  0]
[             ]
[0  0  0  0  0]
[             ]
[0  0  0  0  0]



matrix_density sd;


0


for i := 1:5 do sd(i,i) := i;



sd;


[1  0  0  0  0]
[             ]
[0  2  0  0  0]
[             ]
[0  0  3  0  0]
[             ]
[0  0  0  4  0]
[             ]
[0  0  0  0  5]



matrix_density sd;


20


for i := 1:5 do write sd(i,i) := sd(i,i);


sd(1,1) := 1

sd(2,2) := 2

sd(3,3) := 3

sd(4,4) := 4

sd(5,5) := 5


densify sd;


[1  0  0  0  0]
[             ]
[0  2  0  0  0]
[             ]
[0  0  3  0  0]
[             ]
[0  0  0  4  0]
[             ]
[0  0  0  0  5]



matrix_density ws;


20


m := mat((a,b),(c,d));


     [a  b]
m := [    ]
     [c  d]



sparsify m;


[a  b]
[    ]
[c  d]



s := ws;


     [a  b]
s := [    ]
     [c  d]



% Implicit mapping:
abs s;


[abs(a)  abs(b)]
[              ]
[abs(c)  abs(d)]


depend {a,b,c,d}, x;


df(s,x);


[df(a,x)  df(b,x)]
[                ]
[df(c,x)  df(d,x)]



% Substitution:
sub(a=aa,s);


[aa  b]
[     ]
[c   d]



sparse_random_matrix sr(5,5);


matrix_density sr;


12


% Trace:
sparse_trace sr;


452

trace sr;


452


densify sr;


[422   0   0  0   0]
[                  ]
[ 0    0   0  0   0]
[                  ]
[ 0    0   0  0   0]
[                  ]
[ 0   428  0  30  0]
[                  ]
[ 0    0   0  0   0]



% Transpose:
densify sparse_tp sr;


[422  0  0   0   0]
[                 ]
[ 0   0  0  428  0]
[                 ]
[ 0   0  0   0   0]
[                 ]
[ 0   0  0  30   0]
[                 ]
[ 0   0  0   0   0]


tp sr;


[422  0  0   0   0]
[                 ]
[ 0   0  0  428  0]
[                 ]
[ 0   0  0   0   0]
[                 ]
[ 0   0  0  30   0]
[                 ]
[ 0   0  0   0   0]


densify tp sr;


[422  0  0   0   0]
[                 ]
[ 0   0  0  428  0]
[                 ]
[ 0   0  0   0   0]
[                 ]
[ 0   0  0  30   0]
[                 ]
[ 0   0  0   0   0]



% Determinant:
m := mat((1,2,3),(4,5,6),(7,8,9));


     [1  2  3]
     [       ]
m := [4  5  6]
     [       ]
     [7  8  9]


s := sparsify m;


     [1  2  3]
     [       ]
s := [4  5  6]
     [       ]
     [7  8  9]


det m;


0

sparse_det s;


0

det s;


0

m := mat((a,b,c),(d,e,f),(g,h,i));


     [a  b  c]
     [       ]
m := [d  e  f]
     [       ]
     [g  h  i]


s := sparsify m;


     [a  b  c]
     [       ]
s := [d  e  f]
     [       ]
     [g  h  i]


det m;


a*e*i - a*f*h - b*d*i + b*f*g + c*d*h - c*e*g

sparse_det s;


a*e*i - a*f*h - b*d*i + b*f*g + c*d*h - c*e*g

det s;


a*e*i - a*f*h - b*d*i + b*f*g + c*d*h - c*e*g

if det m = det s then "OK" else "***** ERROR *****";


OK


% Addition and scalar multiplication:
m1 := mat((a,b),(c,d));


      [a  b]
m1 := [    ]
      [c  d]


m2 := mat((e,f),(g,h));


      [e  f]
m2 := [    ]
      [g  h]


s1 := sparsify m1;


      [a  b]
s1 := [    ]
      [c  d]


s2 := sparsify m2;


      [e  f]
s2 := [    ]
      [g  h]


% symbolic;
% global '(ss1 ss2 ss);
% ss1 := sparse!-matsm 's1;
% ss2 := sparse!-matsm 's2;
% ss := sparse!-addm(ss1,ss2);
% sparse!-matpri sparse!-matsm!*1 ss;
m1 + m2;


[a + e  b + f]
[            ]
[c + g  d + h]


s1 + s2;


[a + e  b + f]
[            ]
[c + g  d + h]


2m1 + 3m2;


[2*a + 3*e  2*b + 3*f]
[                    ]
[2*c + 3*g  2*d + 3*h]


2s1 + 3s2;


[2*a + 3*e  2*b + 3*f]
[                    ]
[2*c + 3*g  2*d + 3*h]



% Matrix multiplication:
m3 := mat((a,b,c),(d,e,f));


      [a  b  c]
m3 := [       ]
      [d  e  f]


m4 := mat((g,h),(i,j),(k,l));


      [g  h]
      [    ]
m4 := [i  j]
      [    ]
      [k  l]


s3 := sparsify m3;


      [a  b  c]
s3 := [       ]
      [d  e  f]


s4 := sparsify m4;


      [g  h]
      [    ]
s4 := [i  j]
      [    ]
      [k  l]


m3*m4;


[a*g + b*i + c*k  a*h + b*j + c*l]
[                                ]
[d*g + e*i + f*k  d*h + e*j + f*l]


s3*s4;


[a*g + b*i + c*k  a*h + b*j + c*l]
[                                ]
[d*g + e*i + f*k  d*h + e*j + f*l]



m4*m3;


[a*g + d*h  b*g + e*h  c*g + f*h]
[                               ]
[a*i + d*j  b*i + e*j  c*i + f*j]
[                               ]
[a*k + d*l  b*k + e*l  c*k + f*l]


s4*s3;


[a*g + d*h  b*g + e*h  c*g + f*h]
[                               ]
[a*i + d*j  b*i + e*j  c*i + f*j]
[                               ]
[a*k + d*l  b*k + e*l  c*k + f*l]



2m1*m3*m4;


         2
mat((2*(a *g + a*b*i + a*c*k + b*d*g + b*e*i + b*f*k),

         2
     2*(a *h + a*b*j + a*c*l + b*d*h + b*e*j + b*f*l)),

                         2      2
    (2*(a*c*g + b*c*i + c *k + d *g + d*e*i + d*f*k),

                         2      2
     2*(a*c*h + b*c*j + c *l + d *h + d*e*j + d*f*l)))


2s1*s3*s4;


         2
mat((2*(a *g + a*b*i + a*c*k + b*d*g + b*e*i + b*f*k),

         2
     2*(a *h + a*b*j + a*c*l + b*d*h + b*e*j + b*f*l)),

                         2      2
    (2*(a*c*g + b*c*i + c *k + d *g + d*e*i + d*f*k),

                         2      2
     2*(a*c*h + b*c*j + c *l + d *h + d*e*j + d*f*l)))



m3*m4*m1*3;


         2                                      2
mat((3*(a *g + a*b*i + a*c*h + a*c*k + b*c*j + c *l),

                         2
     3*(a*b*g + a*d*h + b *i + b*c*k + b*d*j + c*d*l)),

    (3*(a*d*g + a*e*i + a*f*k + c*d*h + c*e*j + c*f*l),

                                 2
     3*(b*d*g + b*e*i + b*f*k + d *h + d*e*j + d*f*l)))


s3*s4*s1*3;


         2                                      2
mat((3*(a *g + a*b*i + a*c*h + a*c*k + b*c*j + c *l),

                         2
     3*(a*b*g + a*d*h + b *i + b*c*k + b*d*j + c*d*l)),

    (3*(a*d*g + a*e*i + a*f*k + c*d*h + c*e*j + c*f*l),

                                 2
     3*(b*d*g + b*e*i + b*f*k + d *h + d*e*j + d*f*l)))



% Positive integer powers:
m5 := mat((1,2),(3,4));


      [1  2]
m5 := [    ]
      [3  4]


s5 := sparsify m5;


      [1  2]
s5 := [    ]
      [3  4]


m5^10;


[4783807   6972050 ]
[                  ]
[10458075  15241882]


s5^10;


[4783807   6972050 ]
[                  ]
[10458075  15241882]



% Rank:
% REDUCE manual example
rank m3;


2

sparse_rank sparsify m3;


2

rank sparsify m3;


2


% Wikipedia example
m6 := mat((1,0,1),(0,1,1),(0,1,1));


      [1  0  1]
      [       ]
m6 := [0  1  1]
      [       ]
      [0  1  1]


rank m6;


2

sparse_rank sparsify m6;


2

rank sparsify m6;


2


% Wikipedia example
m6 := mat((1,1,0,2),(-1,-1,0,-2));


      [1   1   0  2 ]
m6 := [             ]
      [-1  -1  0  -2]


rank m6;


1

sparse_rank sparsify m6;


1

rank sparsify m6;


1


% Submatrices:
% m := mat((a,b,c),(d,e,f),(g,h,i));
% s := sparsify m;
m;


[a  b  c]
[       ]
[d  e  f]
[       ]
[g  h  i]


densify sparse_submatrix(s,2,2);


[a  c]
[    ]
[g  i]

        % not advertised!

% Cofactors:
sparse_cofactor(s,1,1);


e*i - f*h

cofactor(s,1,1);


e*i - f*h

% Determinant via cofactors of first row:
for j := 1 : 3 sum s(1,j)*cofactor(s,1,j);


a*e*i - a*f*h - b*d*i + b*f*g + c*d*h - c*e*g

if ws = det s then "OK" else "***** ERROR *****";


OK


% Inverses and non-positive integer powers:
m5^0;


[1  0]
[    ]
[0  1]


s5^0;


[1  0]
[    ]
[0  1]


m5^(-1);


[-2    1  ]
[         ]
[ 3    -1 ]
[---  ----]
[ 2    2  ]


s5^(-1);


[-2    1  ]
[         ]
[ 3    -1 ]
[---  ----]
[ 2    2  ]


m5^(-1)*m5;


[1  0]
[    ]
[0  1]


s5^(-1)*s5;


[1  0]
[    ]
[0  1]


if ws = s5^0 then "OK" else "***** ERROR *****";


OK

m5*m5^(-1);


[1  0]
[    ]
[0  1]


s5*s5^(-1);


[1  0]
[    ]
[0  1]


if ws = s5^0 then "OK" else "***** ERROR *****";


OK

m5^(-2);


[ 11     -5 ]
[----   ----]
[ 2      2  ]
[           ]
[ -15    7  ]
[-----  --- ]
[  4     4  ]


s5^(-2);


[ 11     -5 ]
[----   ----]
[ 2      2  ]
[           ]
[ -15    7  ]
[-----  --- ]
[  4     4  ]



% 1*1 zero matrix:
matrix m0(1,1);


m0^0;


[1]


sparsify m0;


[0]


sparse_matrix s0(1,1);


s0;


[0]


densify s0;


[0]


s0^0;


[1]


m0^(-1);


***** Zero divisor 

s0^(-1);


***** Zero divisor 

m5*m0;


[0  0]
[    ]
[0  0]


s5*s0;


[0  0]
[    ]
[0  0]


m5/m0;


***** Zero divisor 

s5/s0;


***** Zero divisor 


% 1*1 nonzero matrix:
m0(1,1) := 42;


m0(1,1) := 42

s0(1,1) := 42;


s0(1,1) := 42

m5*m0;


[42   84 ]
[        ]
[126  168]


s5*s0;


[42   84 ]
[        ]
[126  168]


m5/m0;


[ 1     1  ]
[----  ----]
[ 42    21 ]
[          ]
[ 1     2  ]
[----  ----]
[ 14    21 ]


s5/s0;


[ 1     1  ]
[----  ----]
[ 42    21 ]
[          ]
[ 1     2  ]
[----  ----]
[ 14    21 ]



% Explicit mapping:
m := mat((x^2,x^5),(x^4,x^5));


     [ 2   5]
     [x   x ]
m := [      ]
     [ 4   5]
     [x   x ]


s := sparsify m;


     [ 2   5]
     [x   x ]
s := [      ]
     [ 4   5]
     [x   x ]


map(int(~w,x), m);


[  3     6 ]
[ x     x  ]
[----  ----]
[ 3     6  ]
[          ]
[  5     6 ]
[ x     x  ]
[----  ----]
[ 5     6  ]


map(int(~w,x), s);


[  3     6 ]
[ x     x  ]
[----  ----]
[ 3     6  ]
[          ]
[  5     6 ]
[ x     x  ]
[----  ----]
[ 5     6  ]


map(1 + ~w, m);


[ 2       5    ]
[x  + 1  x  + 1]
[              ]
[ 4       5    ]
[x  + 1  x  + 1]


map(1 + ~w, s);


[ 2       5    ]
[x  + 1  x  + 1]
[              ]
[ 4       5    ]
[x  + 1  x  + 1]



% Nullspace:
m := mat((1,2,3,4),(5,6,7,8));


     [1  2  3  4]
m := [          ]
     [5  6  7  8]


s := sparsify m;


     [1  2  3  4]
s := [          ]
     [5  6  7  8]


nullspace mat2list m;


{{1,0,-3,2},{0,1,-2,1}}

sparse_nullspace s;


{{1,0,-3,2},{0,1,-2,1}}

nullspace s;


{{1,0,-3,2},{0,1,-2,1}}


sparse_random_matrix s(5,10);


matrix_density s;


14

m := densify s;


     [0  835  718  0   0    0    0   655  0  0]
     [                                        ]
     [0   0    0   0   0    0   100   0   0  0]
     [                                        ]
m := [0   0    0   0   0    0    0    0   0  0]
     [                                        ]
     [0   0    0   0   0    0    0   563  0  0]
     [                                        ]
     [0   0    0   0  867  552   0    0   0  0]


matrix_density m;


14

nullspace mat2list m;


{{1,0,0,0,0,0,0,0,0,0},

 {0,718,-835,0,0,0,0,0,0,0},

 {0,0,0,1,0,0,0,0,0,0},

 {0,0,0,0,184,-289,0,0,0,0},

 {0,0,0,0,0,0,0,0,1,0},

 {0,0,0,0,0,0,0,0,0,1}}

if second length m - length ws = rank m then "OK" else "***** ERROR *****";


OK

nullspace s;


{{1,0,0,0,0,0,0,0,0,0},

 {0,718,-835,0,0,0,0,0,0,0},

 {0,0,0,1,0,0,0,0,0,0},

 {0,0,0,0,184,-289,0,0,0,0},

 {0,0,0,0,0,0,0,0,1,0},

 {0,0,0,0,0,0,0,0,0,1}}

if second length s - length ws = rank s then "OK" else "***** ERROR *****";


OK


% Eigenvalues and eigenvectors:

% REDUCE manual example:
m := mat((2,-1,1),(0,1,1),(-1,1,1));


     [2   -1  1]
     [         ]
m := [0   1   1]
     [         ]
     [-1  1   1]


s := sparsify m;


     [2   -1  1]
     [         ]
s := [0   1   1]
     [         ]
     [-1  1   1]


mateigen(m, eta);


{{eta - 1,2,

  [arbcomplex(1)]
  [             ]
  [arbcomplex(1)]
  [             ]
  [      0      ]

  },

 {eta - 2,1,

  [      0      ]
  [             ]
  [arbcomplex(2)]
  [             ]
  [arbcomplex(2)]

  }}

sparse_mateigen(s, eta);


{{eta - 1,2,

  [arbcomplex(3)]
  [             ]
  [arbcomplex(3)]
  [             ]
  [      0      ]

  },{eta - 2,1,

     [      0      ]
     [             ]
     [arbcomplex(4)]
     [             ]
     [arbcomplex(4)]

     }}

off sparse_matrix_dense_print;

          % to show matrix type
sparse_mateigen(s, eta);



*** 3 × 1 sparse matrix — 2 nonzero elements: 
{{eta - 1,2,

  ?(1,1) := arbcomplex(5)

  ?(2,1) := arbcomplex(5)

*** 3 × 1 sparse matrix — 2 nonzero elements: 
  },{eta - 2,1,

     ?(2,1) := arbcomplex(6)

     ?(3,1) := arbcomplex(6)
     }}


% Lipschutz examples using generic mateigen:
s := sparsify mat((4,1,-1),(2,5,-2),(1,1,2));


*** 3 × 3 sparse matrix — 9 nonzero elements: 

s(1,1) := 4

s(1,2) := 1

s(1,3) := -1

s(2,1) := 2

s(2,2) := 5

s(2,3) := -2

s(3,1) := 1

s(3,2) := 1

s(3,3) := 2

mateigen(s, lam);



*** 3 × 1 sparse matrix — 3 nonzero elements: 
{{lam - 3,2,

  ?(1,1) := arbcomplex(8) - arbcomplex(7)

  ?(2,1) := arbcomplex(7)

  ?(3,1) := arbcomplex(8)

*** 3 × 1 sparse matrix — 3 nonzero elements: 
  },{lam - 5,1,

     ?(1,1) := arbcomplex(9)

     ?(2,1) := 2*arbcomplex(9)

     ?(3,1) := arbcomplex(9)
     }}

s := sparsify mat((3,-1,1),(7,-5,1),(6,-6,2));


*** 3 × 3 sparse matrix — 9 nonzero elements: 

s(1,1) := 3

s(1,2) := -1

s(1,3) := 1

s(2,1) := 7

s(2,2) := -5

s(2,3) := 1

s(3,1) := 6

s(3,2) := -6

s(3,3) := 2

mateigen(s, lam);



*** 3 × 1 sparse matrix — 2 nonzero elements: 
{{lam - 2,2,

  ?(1,1) := arbcomplex(10)

  ?(2,1) := arbcomplex(10)

*** 3 × 1 sparse matrix — 2 nonzero elements: 
  },{lam + 4,1,

     ?(2,1) := arbcomplex(11)

     ?(3,1) := arbcomplex(11)
     }}


% Mixed matrix types in algebraic expressions:


m := mat((1,2),(3,4));


     [1  2]
m := [    ]
     [3  4]


s := sparsify m;


*** 2 × 2 sparse matrix — 4 nonzero elements: 

s(1,1) := 1

s(1,2) := 2

s(2,1) := 3

s(2,2) := 4


sparse_matrix_auto_convert_type();


dense

(2s+m)*m;


[21  30]
[      ]
[45  66]


(2m+s)^2;


[63   90 ]
[        ]
[135  198]



sparse_matrix_auto_convert_type sparse;


dense

(2s+m)*m;


*** 2 × 2 sparse matrix — 4 nonzero elements: 

?(1,1) := 21

?(1,2) := 30

?(2,1) := 45

?(2,2) := 66

(2m+s)^2;


*** 2 × 2 sparse matrix — 4 nonzero elements: 

?(1,1) := 63

?(1,2) := 90

?(2,1) := 135

?(2,2) := 198


sparse_matrix_auto_convert_type none;


sparse

(2s+m)*m;


***** 

[1  2]
[    ]
[3  4]

 invalid as scalar

(2m+s)^2;



*** 2 × 2 sparse matrix — 4 nonzero elements: 
***** 

?(1,1) := 1

?(1,2) := 2

?(2,1) := 3

?(2,2) := 4
 invalid as scalar


;

end;