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