Hi,

I'm trying to write the equivalent of this matlab code in F#, with ILP64 MKL:

S = temp + (I - diag(diag(temp)))

where S, temp are matrices in Rn, I is the identity matrix in Rn, and diag(diag(temp)) is the matrix containing only the main diagonal of temp

In other words, given the matrix

[ a b c

d e f

g h i ]

I want the result

[1 b c

d 1 f

g h 1 ]

for any square sparse matrix.

I think I can perform part of this operation using mkl_dcsrmv, which computes

y := alpha*A*x + beta*y

where x and y are the values array of the identity matrix, alpha is -1.0 and beta is 1.0. According to the manual, I can operate on the diagonal by setting the first element of matdescra to 'D'.

Given the matrix

[1 2 3

4 5 6

7 8 9 ]

my calculation correctly returns the vector [-5 -14 -3] for a General matrix, but when I set matdescra(1) to D, I get [1 1 1] where I expect [0 -4 -8]

It's as though setting matdescra(1) to 'D' returns, not the main diagonal of my matrix, but [0 0 0].

I include some F# code that reproduces the issue, error handling and memory cleanup elided.

let test_dcsrmv (structure:char) =

let p x = PinnedArray.of_array x

let mutable transOption = 'N'

let mutable m = 3

let mutable k = 3

let mutable alpha = -1.0

let mutable beta = 1.0

let madescra = [| structure; '_'; 'N'; 'F'|] |> PinnedArray.of_array

let indx = p [|1; 2; 3;

1; 2; 3;

1; 2; 3; |]

let vals = p [|

1.0; 2.0; 3.0;

4.0; 5.0; 6.0;

7.0; 8.0; 9.0; |]

let pntrb = p [|1; 4; 7;|]

let pntre = p [|4; 7; 10|]

let x = p [| 1.0; 1.0; 1.0 |]

let y = [| 1.0; 1.0; 1.0 |]

let y_handle = p y

mkl_dcsrmv(&&transOption, &&m, &&k, &&alpha, madescra.Ptr, vals.Ptr, indx.Ptr, pntrb.Ptr, pntre.Ptr, x.Ptr, &&beta, y_handle.Ptr)

y

let general = test_dcsrmv 'G' // returns [-5 -14 -23]

let diagonal = test_dcsrmv 'D' // returns [1 1 1]