Page **5** of **7**

### Choice of matrix library and interface

Posted: **November 22nd, 2011, 5:47 am**

by **Cuchulainn**

I'll clean up the code and present it soon. Maube we can use it as a baseline benchmark case? btw we could also use Cholesky and LU as tests for dense matrices?

### Choice of matrix library and interface

Posted: **November 22nd, 2011, 2:52 pm**

by **Polter**

Cuch, I can try comparing with Eigen -- regarding LU decompositions, it hasPartialPivLU (LU decomposition of a matrix with partial pivoting) and FullPivLU (LU decomposition of a matrix with complete pivoting) which offer different performance/accuracy/matrix-requirements trade-offs/characteristics (see below).As for Cholesky decompositions, those two are available:Standard Cholesky decomposition (LL^T)Robust Cholesky decomposition of a matrix with pivoting (LDL^T).

http://eigen.tuxfamily.org/dox/Tutorial ... tions.html

### Choice of matrix library and interface

Posted: **November 22nd, 2011, 3:30 pm**

by **Cuchulainn**

QuoteOriginally posted by: PolterCuch, I can try comparing with Eigen -- regarding LU decompositions, it hasPartialPivLU (LU decomposition of a matrix with partial pivoting) and FullPivLU (LU decomposition of a matrix with complete pivoting) which offer different performance/accuracy/matrix-requirements trade-offs/characteristics (see below).As for Cholesky decompositions, those two are available:Standard Cholesky decomposition (LL^T)Robust Cholesky decomposition of a matrix with pivoting (LDL^T).

http://eigen.tuxfamily.org/dox/Tutorial ... s.htmlThat would be very good. I would look at my LU code again to fix size_t stuff and so on. I have not used pivoting (matices are positive definite) yet and I assume symmetric matrices for CGM.Having a sparse CGM solver is useful because I can use Rothe's method to discretise BS PDE to an elliptic PDE to solve at each time level. The matrix is typically O(1600X1600). Maybe we should use this triad for Thijs as well?

### Choice of matrix library and interface

Posted: **November 22nd, 2011, 4:28 pm**

by **Cuchulainn**

QuoteOriginally posted by: outrunI would be interested is raw storage access performance measures. Writing reading to dense or sparse matrices of various density percentages and sizes. Also for some structured sparse matrices like banded + a full column.Would that make sense? A bit like old harddisk performance tests?Yesp. Try Boost GraphAnd you can reduce the bandwidth using Cuthill-McKee (also BGL).

http://www.boost.org/doc/libs/1_44_0/li ... ering.html

### Choice of matrix library and interface

Posted: **November 22nd, 2011, 7:50 pm**

by **Polter**

QuoteOriginally posted by: CuchulainnQuoteOriginally posted by: PolterCuch, I can try comparing with Eigen -- regarding LU decompositions, it hasPartialPivLU (LU decomposition of a matrix with partial pivoting) and FullPivLU (LU decomposition of a matrix with complete pivoting) which offer different performance/accuracy/matrix-requirements trade-offs/characteristics (see below).As for Cholesky decompositions, those two are available:Standard Cholesky decomposition (LL^T)Robust Cholesky decomposition of a matrix with pivoting (LDL^T).

http://eigen.tuxfamily.org/dox/Tutorial ... s.htmlThat would be very good. I would look at my LU code again to fix size_t stuff and so on. I have not used pivoting (matices are positive definite) yet and I assume symmetric matrices for CGM.Having a sparse CGM solver is useful because I can use Rothe's method to discretise BS PDE to an elliptic PDE to solve at each time level. The matrix is typically O(1600X1600). Maybe we should use this triad for Thijs as well?Great! Yeah, we can try solving some kind of simple linear system, perhaps indeed something along the lines of your CGM example (will wait for your presentation of the cleaned up code, of course). As I've mentioned, I can look into porting into Eigen and benchmarking against its linear algebra decompositions' implementations. Just a simple benchmark, V1, solely dense algebra. Sparse matrices and/or PDEs are "interesting future research directions" / V2, no time to look at it right now. If the code-relevant-for-benchmarking doesn't fit on a single screen of a traditional 80x25 terminal, I'd say let's consider that a V3 and leave to future generations!

### Choice of matrix library and interface

Posted: **November 22nd, 2011, 10:19 pm**

by **Cuchulainn**

Ideally, this time-sharing but it is not very interoperable these days

### Choice of matrix library and interface

Posted: **November 22nd, 2011, 11:24 pm**

by **Polter**

QuoteOriginally posted by: CuchulainnIdeally, this time-sharing but it is not very interoperable these daysCuch, LOL! Somehow, you made me think of this:yes, it's all your fault!

### Choice of matrix library and interface

Posted: **November 26th, 2011, 8:59 am**

by **Cuchulainn**

Here is the code for CGM as we discussed. It is in uBLAS.Now to compare with the other libraries. Feedback, please.

### Choice of matrix library and interface

Posted: **November 26th, 2011, 1:51 pm**

by **Cuchulainn**

Agree on all counts. At some stage we will have enough code/momentum to do all this.I suggest envionments;1 Dev (e.g. CGM and double)2 beta and feedback (CGM with <T>)3. Acceptance testingThen the choice is more flexibility (e.g. 95% of cases is double anyways).This environment set works in practice (especially with real OS with access control s/w) but we have the s/w tools here for this?

### Choice of matrix library and interface

Posted: **December 1st, 2011, 8:24 am**

by **quartz**

I was playing a bit with Cuch's nice CGM and making it work on both uBlas&eigen, in a very crude way.But I was wondering, does anyone remember the design rationale in uBlas for prod() in place of operator*? There might be something about the type of intermediate temporaries in the case of mixed types, but surely that cannot be all of the story! Google is not my friend today.There are even libraries retrofitting operator* onto uBlas for convenience anyway, and I'm not so sure why should one be against it, readability&mantainability is important indeed.(*this* is not a note in favour of eigen/against uBlas as the deciding factors are primarily others like e.g. accuracy&speed, but it'd be good being able to stay library independent anyway... and hardly prod(other_matrix_expr A, other_vec_expr x){return A*x;} is really the way to go..).

### Choice of matrix library and interface

Posted: **December 1st, 2011, 9:09 am**

by **Cuchulainn**

QuoteOriginally posted by: quartzI was playing a bit with Cuch's nice CGM and making it work on both uBlas&eigen, in a very crude way.But I was wondering, does anyone remember the design rationale in uBlas for prod() in place of operator*? There might be something about the type of intermediate temporaries in the case of mixed types, but surely that cannot be all of the story! Google is not my friend today.There are even libraries retrofitting operator* onto uBlas for convenience anyway, and I'm not so sure why should one be against it, readability&mantainability is important indeed.(*this* is not a note in favour of eigen/against uBlas as the deciding factors are primarily others like e.g. accuracy&speed, but it'd be good being able to stay library independent anyway... and hardly prod(other_matrix_expr A, other_vec_expr x){return A*x;} is really the way to go..).The prod(..) is a bit awkward, indeed and probably not portable. My hunch is that * may introduce temporary variables if the user is not careful (d = a*b + c). And I don't know if expression templates would charm numericists.What I like is to use a 'push model' in which the desired result M is modified when we call prod(const input& A, output& M). It means I create M once at a higher level, saves memory stuff.