Create / manipulate diagonal matrices in Chapel

Question

I have a square matrix A

use LinearAlgebra;    
proc main() {
  var A = Matrix(
       [4.0, 0.8, 1.1, 0.0, 2.0]
      ,[0.8, 9.0, 1.3, 1.0, 0.0]
      ,[1.1, 1.3, 1.0, 0.5, 1.7]
      ,[0.0, 1.0, 0.5, 4.0, 1.5]
      ,[2.0, 0.0, 1.7, 1.5, 16.0]
      );
}

And I want to construct the diagonal matrix D = 1/sqrt(a_ii). It seems like I have have to extract the diagonal, then operate on each element. I expect this matrix is be very large and sparse, if that changes the answer.

ben-albrecht · Accepted Answer · 2017-08-21T19:47:43.930

2

Here's a solution using the LinearAlgebra module in 1.16 (pre-release):

use LinearAlgebra;

var A = Matrix(
               [4.0, 0.8, 1.1, 0.0, 2.0],
               [0.8, 9.0, 1.3, 1.0, 0.0],
               [1.1, 1.3, 1.0, 0.5, 1.7],
               [0.0, 1.0, 0.5, 4.0, 1.5],
               [2.0, 0.0, 1.7, 1.5, 16.0]
              );

var S = sqrt(1.0/diag(A));

// Type required because of promotion-flatting
// See Linear Algebra documentation for more details..
var B: A.type = diag(S);

writeln(B);

edited Aug 21 '17 at 19:47

answered Aug 19 '17 at 00:11

ben-albrecht

1,785
10
23

1

oooohhh, slick! – Brian Dolan Aug 19 '17 at 00:12
2

Will the conversion step from **`A( 2D, sparse ) -> S( 1D, dense )`**, and the trailing-end flattening-promotion guarantee the final representation of **`B( 2D, diagonal )`** yet remain a sparse representation? – user3666197 Aug 19 '17 at 06:28
1

I mistakenly missed that part of the question. The `LinearAlgebra` module does not yet officially sparse arrays, so there can be on guarantees that its routines will just *work* with sparse arrays. That said, once they do, I'd expect the above code to work as well (where `Matrix` might be something like `sparseMatrix`). – ben-albrecht Aug 19 '17 at 14:28
2

The type of `S` can be omitted once [#7057](https://github.com/chapel-lang/chapel/pull/7057) is merged. I'll update this answer once that is live on the `1.16 (pre-release)`. – ben-albrecht Aug 19 '17 at 14:56
1

So, if I follow your note above, Ben, the declaration below **`var D: [A.domain];`** ought precisely work right-enough for this very reason -- replicating the sparse-**`A`**'s-domain ( so taking here care of having an identical sparse-domain allocation-map for **`D`** ( no values set into `D` so far ) and next, the **`forall{...}`** computes and sets those and only those sparse- **`D`** diagonal elements, that were originally non-zero elements ( there were by-definition none others than non-zero, right? ) on the sparse-**`A`**'s diagonal, correct? – user3666197 Aug 19 '17 at 16:07

score 1 · Answer 2 · answered Aug 18 '17 at 20:57

1

Did you try this approach?

use Math;
var D: [A.domain];
forall      i in D.dim( 1 ) {
        D[i,i] = 1 / Math.sqrt( A[i,i] ); // ought get fused-DIV!0 protection
}

_{( A.T.M. <TiO>-IDE has not so far fully functional the LinearAlgebra package, so cannot show you the results live, yet hope you would enjoy the way forwards )}

answered Aug 18 '17 at 20:57

user3666197

1
6
50
92

Although I would prefer the Linear Algebra package because of the new sparse formats coming out. – Brian Dolan Aug 18 '17 at 21:04

score 1 · Answer 3 · answered Aug 19 '17 at 17:07

Here's some code that works with a sparse diagonal array in version 1.15 today without linear algebra library support:

config const n = 10;   // problem size; override with --n=1000 on command-line

const D = {1..n, 1..n},                        // dense/conceptual matrix size
      Diag: sparse subdomain(D) = genDiag(n);  // sparse diagonal matrix

// iterator that yields indices describing the diagonal
iter genDiag(n) {
  for i in 1..n do
    yield (i,i);
}

// sparse diagonal matrix
var DiagMat: [Diag] real;

// assign sparse matrix elements in parallel
forall ((r,c), elem) in zip(Diag, DiagMat) do
  elem = r + c/10.0;

// print sparse matrix elements serially
for (ind, elem) in zip(Diag, DiagMat) do
  writeln("A[", ind, "] is ", elem);

// dense array
var Dense: [D] real;

// assign from sparse to dense
forall ij in D do
  Dense[ij] = DiagMat[ij];

// print dense array
writeln(Dense);

Note that I think that, ultimately, one would want to create a specialized 'SparseDiagonal' domain map since it's a common case and one that can be represented much more efficiently than a general sparse array. — Brad, Aug 19 '17 at 17:08

Create / manipulate diagonal matrices in Chapel

3 Answers3

Did you try this approach?