1% Test file for Sparse Matrices and the Linear Algebra Package for
2% Sparse Matrices.
3
4% Author: Stephen Scowcroft.                    Date: June 1995.
5
6% Firstly, the matrices need to be created.
7
8% This is the standard way to create a sparse matrix.
9
10% Create a sparse matrix.
11
12sparse mat1(5,5);
13
14%Fill the sparse matrix with data
15
16mat1(1,1):=2;
17mat1(2,2):=4;
18mat1(3,3):=6;
19mat1(4,4):=8;
20mat1(5,5):=10;
21
22sparse mat4(5,5);
23
24mat4(1,1):=x;
25mat4(2,2):=x;
26mat4(3,3):=x;
27mat4(4,4):=x;
28mat4(5,5):=x;
29
30% A small function to automatically fill a sparse matrix with data.
31
32procedure makematsp(nam,row);
33 begin;
34  sparse nam(row,row);
35    for i := 1:row do <<nam(i,i):=i>>
36 end;
37
38clear mat2;
39makematsp(mat2,100);
40
41% Matrices created in the standard Matrix way.
42
43zz1:=mat((1,2),(3,4));
44zz2:=mat((x,x),(x,x));
45zz3:=mat((i+1,i+2,i+3),(4,5,2),(1,i,0));
46
47% I have taken advantage of the Linear Algebra Package (Matt Rebbeck)
48% in order to create some Sparse Matrices.
49
50mat3:=diagonal(zz1,zz1,zz1);
51mat5:=band_matrix({1,3,1},100)$
52mat6:=diagonal(zz3,zz3);
53mat7:=band_matrix({a,b,c},4);
54
55% These are then "translated" into the Sparse operator using the
56% function transmat.
57% This is a destructive function in the sense that the matrices are no
58% longer of type 'matrix but are now 'sparse.
59
60transmat mat3;
61transmat mat5;
62transmat mat6;
63transmat mat7;
64
65poly  := x^7+x^5+4*x^4+5*x^3+12;
66poly1 := x^2+x*y^3+x*y*z^3+y*x+2+y*3;
67
68% Firstly some basic matrix operations.
69% These are the same as the present matrix package
70
71mat1^-1;
72mat4^-1;
73mat2 + mat5$
74mat2 - mat5$
75mat1-mat1;
76mat4 + mat1;
77mat4 * mat1;
78
792*mat1 + (3*mat4 + mat1);
80% It is also possible to combine both 'matrix and 'sparse type matrices
81% in these operations.
82
83pp:=band_matrix({1,3,1},100)$
84mat5*pp;
85
86mat5^2$
87
88det(mat1);
89det(mat4);
90trace(mat1);
91trace(mat4);
92
93rank(mat1);
94rank mat5;
95
96tp(mat3);
97
98spmateigen(mat3,eta);
99
100% Next, tests for the Linear Algebra Package for Sparse Matrices.
101
102%Basic matrix manipulations.
103
104spadd_columns(mat1,1,2,5*y);
105spadd_rows(mat1,1,2,x);
106
107spadd_to_columns(mat1,3,1000);
108spadd_to_columns(mat5,{1,2,3},y)$
109spadd_to_rows(mat1,2,1000);
110spadd_to_rows(mat5,{1,2,3},x)$
111
112spaugment_columns(mat3,2);
113spaugment_columns(mat1,{1,2,5});
114spstack_rows(mat1,3);
115spstack_rows(mat1,{1,3,5});
116
117spchar_poly(mat1,x);
118
119spcol_dim(mat2);
120sprow_dim(mat1);
121
122spcopy_into(mat7,mat1,2,2);
123spcopy_into(mat7,mat1,5,5);
124spcopy_into(zz1,mat1,1,1);
125
126spdiagonal(3);
127% spdiagonal can take both a list of arguments or just the arguments.
128spdiagonal({mat2,mat5})$
129spdiagonal(mat2,mat5)$
130% spdiagonal can also take a mixture of 'sparse and 'matrix types.
131spdiagonal(zz1,mat4,zz1);
132
133spextend(mat1,3,2,x);
134
135spfind_companion(mat5,x);
136
137spget_columns(mat1,1);
138spget_columns(mat1,{1,2});
139spget_rows(mat1,3);
140spget_rows(mat1,{1,3});
141
142sphermitian_tp(mat6);
143
144% matrix_augment and matrix_stack can take both a list of arguments
145% or just the arguments.
146
147spmatrix_augment({mat1,mat1});
148spmatrix_augment(mat5,mat2,mat5)$
149spmatrix_stack(mat2,mat2)$
150
151spminor(mat1,2,3);
152
153spmult_columns(mat1,3,y);
154spmult_columns(mat2,{2,3,4},100)$
155spmult_rows(mat2,2,x);
156spmult_rows(mat1,{1,3,5},10);
157
158sppivot(mat3,3,3);
159sprows_pivot(mat3,1,1,{2,4});
160
161spremove_columns(mat1,3);
162spremove_columns(mat3,{2,3,4});
163spremove_rows(mat1,2);
164spremove_rows(mat2,{1,3})$
165spremove_rows(mat1,{1,2,3,4,5});
166
167spswap_cols(mat1,2,4);
168spswap_rows(mat5,1,2)$
169spswap_entries(mat1,{1,1},{5,5});
170
171
172
173% Constructors - functions that create matrices.
174
175spband_matrix(x,500)$
176spband_matrix({x,y,z},6000)$
177
178spblock_matrix(1,2,{mat1,mat1});
179spblock_matrix(2,3,{mat3,mat6,mat3,mat6,mat3,mat6});
180
181spchar_matrix(mat3,x);
182
183cfmat := spcoeff_matrix({y+4*+-5*w=10,y-z=20,y+4+3*z,w+x+50});
184first cfmat * second cfmat;
185third cfmat;
186
187spcompanion(poly,x);
188
189sphessian(poly1,{w,x,y,z});
190
191spjacobian({x^4,x*y^2,x*y*z^3},{w,x,y,z});
192
193spjordan_block(x,500)$
194
195spmake_identity(1000)$
196
197on rounded; % makes output easier to read.
198ch := spcholesky(mat1);
199tp first ch - second ch;
200tmp := first ch * second ch;
201tmp - mat1;
202off rounded;
203
204% The gram schmidt functions takes a list of vectors.
205% These vectors are matrices of type 'sparse with column dimension 1.
206
207%Create the "vectors".
208sparse a(4,1);
209sparse b(4,1);
210sparse c(4,1);
211sparse d(4,1);
212
213%Fill the "vectors" with data.
214a(1,1):=1;
215b(1,1):=1;
216b(2,1):=1;
217c(1,1):=1;
218c(2,1):=1;
219c(3,1):=1;
220d(1,1):=1;
221d(2,1):=1;
222d(3,1):=1;
223d(4,1):=1;
224
225spgram_schmidt({{a},{b},{c},{d}});
226
227on rounded; % again, makes large quotients a bit more readable.
228% The algorithm used for splu_decom sometimes swaps the rows of the
229% input matrix so that (given matrix A, splu_decom(A) = {L,U,vec}),
230% we find L*U does not equal A but a row equivalent of it. The call
231% spconvert(A,vec) will return this row equivalent
232% (ie: L*U = convert(A,vec)).
233lu := splu_decom(mat5)$
234tmp := first lu * second lu$
235tmp1 := spconvert(mat5,third lu);
236tmp - tmp1;
237% and the complex case..
238on complex;
239lu1 := splu_decom(mat6);
240mat6;
241tmp := first lu1 * second lu1;
242tmp1 := spconvert(mat6,third lu1);
243tmp - tmp1;
244off complex;
245
246mat3inv := sppseudo_inverse(mat3);
247mat3 * mat3inv;
248
249
250% Predicates.
251
252matrixp(mat1);
253matrixp(poly);
254
255squarep(mat2);
256squarep(mat3);
257
258symmetricp(mat1);
259symmetricp(mat3);
260
261sparsematp(mat1);
262sparsematp(poly);
263
264off rounded;
265
266end;
267