Sparse Matrix Type Operations

Morse storage is used by default. First line: "n m s nbcoef" where s tells if the matrix is symmetric or not and a line of the form " i j a_ij " where (i,j) \in {1,...,n}x{1,...,m} appears for each nonzero coefficient.

Observe that for a matrix A the number of rows is A.n and the number of columns is A.m

download example: sparsematrix.edp or **return to Matrices and Arrays**

load "lapack" real[int,int] tab(1,7); tab=1; tab(3:5,0)=0; cout << " the single array in double array format tab = " << tab << endl; matrix tab1; tab1=tab; cout << " the sparse vector tab = " << tab1 << endl; int N=3,M=4; real[int,int] A(N,M),B(M,N),C(N,M),AxB(N,N),AmC(N,M),BxA(M,M),AxBt(N,N); A=1; A(0,0)=0; A(0,2)=0; A(1,0)=0; A(1,1)=0; B=3; B(0,1)=0; B(0,2)=0; B(1,0)=0; B(1,0)=0; C=-1; C(0:1,2)=0; // // Matrix operation product // AxB=A*B;// works if lapack is loaded BxA=B*A; cout << " double array A (dim 3x4) = " << A << endl; cout << " double array B (dim 4x3) = " << B << endl; cout << " double array product A*B (dim 3x3)= " << AxB << endl; cout << " double array product B*A (dim 4x4)= " << BxA << endl; cout << " double array A " <<A<< endl; //B=A'; //cout << " double array transpose A' only works the square case " <<B<< endl; //cout <<"non square case the size of B is changed"<<endl; matrix SA,SB,SAxSB,SAt,SC,SAmSC; SA=A; SB=B; SC=C; cout << " SA = " << SA << endl; cout << " SB = " << SB << endl; cout << " C = " << C << endl; cout << " SC = " << SC << endl; SAmSC=SA+(-SC); AmC=A-C; SAxSB=SA*SB; SAt=SA'; cout << " sparse matrix diff SA-SC = " << SAmSC << endl; cout << " the same as matrix difference A-C = " << AmC << endl; cout << " sparse matrix prod SA*SB = " << SAxSB << endl; cout << " the same as matrix product A*B = " << AxB << endl; cout << " sparse matrix transpose SA' = " << SAt << endl; real[int,int] T(4,7),Tt(7,4); T=10; T(0,:)=1;T(:,3)=0; cout << " the matrix T = " << T << endl; matrix ST,STpSA,STt; ST=T; STt=ST'; cout << " the sparse matrix ST = " << ST << endl; cout << " the sparse matrix STt = " << STt << endl; cout << " +++++++++++++++++++++++++++++++++++++++++ " << endl; cout << " adding sparse matrices of different sizes " << endl; cout << " takes largest of both sizes and add entries with common indices "<<endl; cout << " +++++++++++++++++++++++++++++++++++++++++ " << endl; cout << " the sparse matrix SA = " << SA << endl; cout << " the sparse matrix ST = " << ST << endl; STpSA=ST+SA; cout << " the sparse matrix ST+SA = " << STpSA << endl; STt=ST'; cout << " the sparse matrix ST = " << ST << endl; cout << " the sparse matrix ST transpose IS WORKING in sparse matrix format= " << STt << endl; real[int] v(7); v=1; real[int] test=ST*v; cout << " ST = " << ST<< endl; cout << " v = " << v<< endl; cout << "4x7 sparse matrix times a 7x1 vector ST*v = " << test << endl;