Multi-threading with arpack

[ajos6183]

Hi @wo80,

I have started using the arpack solver and it is very good. Thank you for the work you have done. The only minor issue is that i am not able to run it in multiple threads at the same time. This is my code below. Am i doing something wrong or is it not thread safe?

'EigenProblem
Friend Function EigenProblem(Ke As CSparse.Storage.CompressedColumnStorage(Of Double), Kg As CSparse.Storage.CompressedColumnStorage(Of Double),
                             PositiveKg As Boolean, BucklingSign As BucklingSign, LC As String) As Tuple(Of Double, Double())
    'Arpack solves Ke.x = λ.Kg.x, so -Kg is required as input
    If PositiveKg = True Then
        For i = 0 To Kg.Values.Count - 1
            Kg.Values(i) *= -1
        Next
    End If
    Dim Eig As CSparse.Double.Solver.Arpack = New CSparse.Double.Solver.Arpack(Kg.ToSparseMatrix, Ke.ToSparseMatrix) With {.Tolerance = 10 ^ -5, .ComputeEigenVectors = True, .Iterations = 1000}
    Dim Solution = Eig.SolveGeneralized(Math.Min(Ke.ColumnCount - 2, 10), CSparse.Interop.ARPACK.Job.LargestMagnitude)
    Dim EigenValueIndex As Integer = -1
    For i = 0 To Solution.EigenValues.Count - 1
        Select Case BucklingSign
            Case BucklingSign.Positive
                If Solution.EigenValues(i).Real > 0 Then
                    EigenValueIndex = i
                    Exit For
                End If
            Case BucklingSign.Negative
                If Solution.EigenValues(i).Real < 0 Then
                    EigenValueIndex = i
                    Exit For
                End If
            Case BucklingSign.Any
                EigenValueIndex = 0
                Exit For
            Case Else : Throw New NotImplementedException
        End Select
    Next
    If EigenValueIndex = -1 Then
        Return Tuple.Create(Of Double, Double())(1.0 / Solution.EigenValues.Last.Real, Solution.EigenVectors.Column(Solution.EigenVectors.ColumnCount - 1))
    Else
        Dim EigenValue As Double = 1.0 / Solution.EigenValues(EigenValueIndex).Real
        Dim EigenVector As Double() = Solution.EigenVectors.Column(EigenValueIndex)
        Return Tuple.Create(Of Double, Double())(EigenValue, EigenVector)
    End If
End Function

Friend Enum BucklingSign
    Positive
    Negative
    Any
End Enum

[wo80]

Thanks for reminder. The native code is not thread-safe and this should be mentioned.

Unfortunately, f2c introduces a lot of local static variables, even with -a switch enabled. Open any file in https://github.com/wo80/vs-arpack/tree/master/src/ARPACK/arpack-ng/SRC and you'll see the problem.

I think, changing those variables to non-static should be ok, but it has to be tested. I'll look into it, if I find the time.

[ajos6183]

Thank you for the clarification. In this case, i will consider using your feast solver. I didn't know MKL was free for use and redistribution but it seems so!

[wo80]

So, after having a first look, changing all local static variables to be non-static will break ARPACK. The global state is used for reverse communication (though this doesn't seem to be consistent - some local static variables are really just local).

It might still be possible to make ARPACK thread-safe, see https://stackoverflow.com/questions/3889222/is-arpack-thread-safe

[ajos6183]

The explanation on stackoverflow is beyond me because i don't understand the particular details of FORTRAN/C/C++. I am a structural engineer by training. Thank you for looking into it further! But it is not fair for me to ask you to invest more time in it. Please only look into it further if it is of interest to you or needed by you.

Also i found it strange how the FEAST solver needs a range for the eigenvalues. It doesn't really work for my usage in buckling problems (if i tested it correctly). So, i will stick to using the ARPACK solver making sure that i call it once at a time.

[wo80]

The FEAST generalized eigenproblem solver had a bug (see https://github.com/wo80/csparse-interop/commit/b825d54ba032558bdfa3bd0f5680c5882693df90), so this might have messed up your tests. But I agree, the interface is a bit inconvenient.

The latest version of MKL (2019) has added a new extended eigensolver, which seems to be pretty fast: TestExtendedEigensolver.cs

Regarding ARPACK: I don't think that I'll look into it any time soon. I opened an issue for the vs-arpack project, so please close this one, if you think we're done here.

[ajos6183]

OK thanks again!