2

我有一个方阵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]
      );
}

我想构造对角矩阵D = 1/sqrt(a_ii)。似乎我必须提取对角线，然后对每个元素进行操作。我希望这个矩阵非常大且稀疏，如果这改变了答案。

4

3 回答 3

2

这是使用LinearAlgebra模块的解决方案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);
于 2017-08-19T00:11:34.917 回答
1

你试过这种方法吗？

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
}

<TiO>-IDE（到目前为止， ATM还没有完全发挥该LinearAlgebra软件包的功能，因此无法实时显示结果，但希望您会喜欢前进的道路）

于 2017-08-18T20:57:25.617 回答
1

下面是一些在 1.15 版中适用于稀疏对角数组的代码，没有线性代数库支持：

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);
于 2017-08-19T17:07:25.333 回答