logo.jpg™

Sparse Matrix Library for C++

www.rejonesconsulting.com

sparse.h contains:

·         a complete Sparse matrix package comparable to the Matrix package above

·         a single excellent solver using the Generalized Minimum Residual method. 

·         a few preconditioning tools. (this area is under development).  See About Preconditioning and About Uncouple

·         some text displays of sparse matrices

·         a eigenvalue solvers, and optionally eigenvectors, for symmetric, very sparse systems.

 

You may download this library here: sparse.h.

Sparse.h includes this library: matrix.h, which you must download separately.

 

To go directly to solving sparse systems, scroll down or click here.

To go directly to sparse eigenvalues/eigenvectors, scroll further down, or click here.

 

A note about the implementation: A Sparse matrix is implemented as a list of “Node”s.  A Node represents one row of a Sparse matrix.  Each element of a row is represented by a column index and a coefficient value. •When normalized, the diagonal element is placed first in the list.  Certain operations require that the matrix be normalized.  Most users will never have to be aware of the Node class, or the need to normalize the Sparse matrix.  Developers who want to add algorithms can look at the Node class description in sparse.h for details.

Quick Example

 

#include “matrix.h”

//declare a 1000 by 1000 sparse matrix

  Sparse S(1000);

 

//set up whatever values are needed in S…

  for (int i=0; i<1000; i++)

   for (j=…whatever…)

    S(i,j) = …whatever…

 

//set up right hand side vector, r…

  Vector r(1000);

  for (int i=0; i<1000; i++)

   r[i] = …whatever…

 

//solve…

  Vector x=S.solve_gmres(r);

  //or x=S.solve_gmres(r,50,0.00001);

Creating Sparse Matrices

 

Sparse matrices have a single dimension, which is both the number of rows and the number of columns

 

  Sparse S; //no rows, no columns

  Sparse S(1000); //1000 rows, 1000 columns, diagonal elements=0.0

  Sparse S(1000,1.0); //diagonal elements=1.0

  Sparse S(D); //diagonal of S assigned from Diagonal matrix D.

  Sparse S(S2);    //S becomes a copy of S2(**)

  Sparse S = S2;   //S becomes a copy of S2(**)

  S.clear();       //all of S’s memory is deleted

                   //(and S now has no rows or columns.)

  S1.resize(500); //resizes to smaller (deletes data)

                  //or larger (zero fills)

 

(**) Note: any excessive storage space in S2 will be eliminated during the copy or assignment.

Index Operations

 

The dimensions of the sparse matrix can be accessed as follows:

  int i=S.dim(); //or size() or dim1() or dim2()

    //all are equivalent

 

The element at row i, column j, can be set or accessed like this:

  S(i,j) = 3.0;  //sets the sparse matrix element (i,j) to 3.0

  S.set(i,j,3.0);  //optional alternate call

  //Note: the element will be added to the list if needed

 

  double x = S(i,j);  //gets the element at (2,3).

  //Note this form inserts an element at (i,j) with the value 0.0

  //if the element is not already present.

  double x = S.get(i,j); 

  //Note that the “get” form does not insert an element at (i,j).

  //It just returns 0.0 if the element is not present in the list.

 

To see if a particular element is actually present, without causing any change to the matrix:

  bool pres=S.present(i,j);

Element-Wise Operations

 

  S.multiply_elements_by(2.0); //multiplies ALL elements by 2.0

  S.multiply_diagonal_by(2.0); //multiplies diagonal elements ONLY by 2.0

  S.add_to_diagonal(2.0); // //adds 2.0 ONLY to diagonal elements

Whole Matrix Operations

 

  Sparse S1,S2,S3;  //define as needed

  Vector V1,V2;     //define as needed

 

  S1.transpose(S2); //puts transpose of S1 into S2

 

  Vector V2=S1*V1;  //Sparse matrix times vector

  Sparse S3=S1*S2;  //Sparse matrix times sparse matrix

 

  S1.add(S2);       //Adds S2 into S1

  S1.subtract(S2);  //Subtracts S2 from S1

Comparison Operations

 

  if (S1==S2) …  //Corresponding elements of S1 and S2

           //must be exactly the same, but incidental

           //details of ordering within element lists are ignored

 

  if (S1.approximate(S2,0.00001)) … //Corresponding elements of S1 and S2

           //must agree with an absolute (not relative) difference

           //of the given amount

 

  If (S1.equivalent(S2,0.00001)) … //similar to approximate(),

           //but this method effectively sorts each matrix' rows

           //into increasing row moment, and then compares

           //the apparently corresponding rows.

             //Note: this method is NOT fool-proof

Miscellaneous Operations

 

  double x=S.maxabs(); //returns absolute value of element with the

                       //largest absolute value

 

  int i=S.bandedness(); //returns maximum off-diagonal distance

                        //for any element

 

  int i=S.band_average(); //returns average off-diagonal distance

 

There are many other operations in sparse.h.

Textual Displays of a Sparse Matrix

 

  S.print(); //prints a subset of S

 

  S.print_details(); //prints all structure information for up to first 21 rows

  S.print_star_magnitudes(); //prints “star magnitudes” for up to 60 rows or columns

  S.print_star_map(); //compactly displays entire matrix in a

      //40 row by 60 column display by grouping elements and 

      //displaying the largest magnitude of any element in the group.

Solving Sparse Matrix Problems

 

See Quick Demo above for a quick introduction to setting up and solving a sparse matrix problem.

  Vector B ... //set up corresponding right hand

  int m=30; //number of internal iterations before restart...default 50

  double accuracy=1.0E06; //default is 1.0E-7

  zero=1.0E-10; //default is 1.0E-10

  Vector X=S.solve_gmres(B,m,accuracy,zero);

  //the process iterates until the norm of the residual is smaller

  //than “accuracy” times the norm of B,

  //or the norm of the change in the solution is less than “zero”

  //times the norm of the initial solution to the system.

 

To use the default parameters, just call:

  Vector X=S.solve_gmres(B);

About Preconditioning

 

Most realistic sparse systems will require some kind of preconditioning for even the GMRES method to converge.   The point of preconditioning is usually to “partially solve” the system, preferably using methods appropriate to the particular situation.  Preconditioning methods are mostly beyond the present scope of this package.  However, we do offer some rescaling and reordering utilities, and one generic method, called “uncouple()”.

 

Rescaling. Every row of the sparse system is scaled to a unit norm.  Just this small step will transform some example sparse problems (such as certain ones  offered by the Matrix Market site) from hard to easy.  If a row contains only near-zero values, then it is  cleared completely.

  Sparse S ... //set up sparse matrix

  Vector B ... //set up corresponding right hand

  S.scale_rows(B);

 

Row renumbering.  For well set-up problem this process is usually counter-productive.  So we only proceed with the result of renumbering if the total off-diagonal mass for the entire matrix is reduced by it.  The process is straightforward: we compute the first moment of the coefficients about column zero. Then we sort the rows into increasing order of that moment.  In many cases this will restore the natural order to a matrix which is not presenting in a reasonable ordering.  You do not need to call S.scale_rows first.)

 

  Sparse S ... //set up sparse matrix

  Vector B ... //set up corresponding right hand

  S.reorder_rows(B);

 

Generic “uncoupling”.  By “partially solve” above, we mean to drive the Sparse matrix, S, toward a more nearly diagonal form, and of course, to carry along the right hand side vector, B, as appropriate.  The difficulty is that almost any row or column operation required to move to a more diagonalized form is likely to cause the creation of new non-zero elements.  This is called infill.  We have to limit any predconditioning process to the amount of infill which is tolerable.  For example, in Sparse::uncouple() currently we allow the amount of storage space for a given row to approximately double before we stop the process for that row.  (Note: uncouple can be called repeatedly, but the storage used may double again each time.)

 

A more detailed description of the process in uncouple() is available: About_Uncouple.  Here we give a shorter desctiption.  In  uncouple(), we look locally around each row of the matrix for opportunities to reduce the off-diagonal norm (of the unit-scaled version of the row).  Suppose we are processing row i.  Then we will look above and below row  i by an arbitrary heuristically-determined number of rows, and compute the projection of the off-diagonal elements of row i onto each such nearby row j.  We then try subtracting that projection from row i, and see if the off-diagonal mass norm (of the unit-scaled version of the modified row) is reduced.  If the reduction is not trivial, we accept this change.  Then we check to see if we have exceeded our infill budget for this row.  If so, we are done with row i and proceed to the next row.  This method is usually not very helpful.  But sometimes it transforms a hard problem into an easy one.  So it is worth a try unless any infill is intolerable.

 

  Sparse S ... //set up sparse matrix

  Vector B ... //set up corresponding right hand

  S.uncouple(B); //precondition.. changes S and B!

  Vector x=S.solve_gmres(B);  //solve

Computing Eigenvalues/Eigenvectors of a Sparse Symmetric Matrix

 

Sparse.h can compute Eigenvalues, and optionally Eigenvectors, of a sparse symmetirc matrix.  However, the matrix needs to be very sparse for our method to work well.  (Perhaps two elements per row.)  Note that we assume that your matrix is properly set up as symmetric.  Symmetry is essential, but is not checked.  Once your matrix, S, is set up, just call:

 

  SparseEig E(S);

 

to create an eigensolver object, E.  Then to compute the eigenvalues, call

 

  Vector eigensolve(T ctol, T etol, bool dovectors);

 

where

·         ctol is the converge criterion,

·         etol is an absolute tolerance below which any element can be set to zero. 

·         “dovectors” is a Boolean that tells whether or not to compute eigenvectors.  Computation of eigenvectors is fairly expensive in time and space, and should not be done if not necessary.

 

For example,

 

  Vector x = E.eigensolve(0.00001, 0.000001, true);

 

If you ask for eigenvectors to be computed, you can retrieve them one at a time with the call

 

  Vector x=E.getEVector(k);

 

where k is the integer index (zero to n-1) of the eigenvalue corresponding to the desired eigenvector.

 

If you do not ask for eigenvectors to be computed, you can still compute one by calling

 

  Vector x = E.computeEVector(S, T eigenvalue, T etol);

 

For example,

 

  Vector x = E.computeEVector(S, 1.0, 0.000001);

 

This would be a very slow way to compute all of the the eigenvalues, but can be fast to compute just one.

Go back to our  main page.