Add preventive and corrective SCOPF (#67)

- migrate to MadNLP 0.6
- migrate to CUDA 4.*
- migrate to ExaPF 0.9
- add preventive and corrective SCOPF
- add script for real-time SCOPF

---------

Co-authored-by: Michel Schanen <mschanen@anl.gov>

.github/workflows/action.yml

@@ -39,7 +39,6 @@ jobs:
     env:
       CUDA_VISIBLE_DEVICES: 1
       JULIA_DEPOT_PATH: /scratch/github-actions/julia_depot_argos
-      JULIA_CUDA_USE_BINARYBUILDER: false
     runs-on: self-hosted
     strategy:
       matrix:

Project.toml

@@ -14,20 +14,18 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
-CUDA = "^3.8"
-CUDAKernels = "0.4"
-ExaPF = "~0.8"
+CUDA = "4.1"
+ExaPF = "~0.9.3"
 FiniteDiff = "2.7"
 Ipopt = "1"
-KernelAbstractions = "0.8"
-MadNLP = "~0.5.2"
+KernelAbstractions = "0.9"
+MadNLP = "0.7"
 MathOptInterface = "1"
-NLPModels = "0.18"
+NLPModels = ">= 0.18"
 julia = "1.6"
 
 [extras]
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
-CUDAKernels = "72cfdca4-0801-4ab0-bf6a-d52aa10adc57"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
 Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
@@ -36,4 +34,4 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "CUDA", "CUDAKernels", "DelimitedFiles", "FiniteDiff", "Ipopt", "LazyArtifacts", "Random"]
+test = ["Test", "CUDA", "DelimitedFiles", "FiniteDiff", "Ipopt", "LazyArtifacts", "Random"]

docs/make.jl

@@ -19,6 +19,7 @@ makedocs(
     ),
     modules = [Argos],
     repo = "https://github.com/exanauts/Argos.jl/blob/{commit}{path}#{line}",
+    doctest=false,
     checkdocs = :exports,
     pages = [
         "Home" => "index.md",

lib/ArgosCUDA.jl/Project.toml

@@ -6,7 +6,6 @@ version = "0.1.0"
 [deps]
 Argos = "ef244971-cf80-42b0-9762-2c2c832df5d5"
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
-CUDAKernels = "72cfdca4-0801-4ab0-bf6a-d52aa10adc57"
 CUSOLVERRF = "a8cc9031-bad2-4722-94f5-40deabb4245c"
 ExaPF = "0cf0e50c-a82e-488f-ac7e-41ffdff1b8aa"
 KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"

lib/ArgosCUDA.jl/src/ArgosCUDA.jl

@@ -7,8 +7,8 @@ using CUDA
 using CUDA.CUSPARSE
 using CUSOLVERRF
 
-using KernelAbstractions
-using CUDAKernels
+import KernelAbstractions as KA
+import KernelAbstractions: @kernel, @index
 
 using ExaPF
 const LS = ExaPF.LinearSolvers

lib/ArgosCUDA.jl/src/kernels.jl

@@ -5,8 +5,8 @@ end
 function Argos.transfer2tril!(hessvals::AbstractVector, H::CuSparseMatrixCSR, csc2tril)
     Hz = nonzeros(H)
     ndrange = (length(hessvals),)
-    ev = _map2vec_kernel!(CUDADevice())(hessvals, Hz, csc2tril; ndrange=ndrange)
-    wait(ev)
+    _map2vec_kernel!(CUDABackend())(hessvals, Hz, csc2tril; ndrange=ndrange)
+    KA.synchronize(CUDABackend())
 end
 
 
@@ -17,8 +17,8 @@ end
 function Argos.fixed!(dest::CuVector, ind_fixed, val::Number)
     length(ind_fixed) == 0 && return
     g_ind_fixed = ind_fixed |> CuArray
-    ev = _fixed_kernel!(CUDADevice())(dest, g_ind_fixed, val; ndrange=length(ind_fixed))
-    wait(ev)
+    _fixed_kernel!(CUDABackend())(dest, g_ind_fixed, val; ndrange=length(ind_fixed))
+    KA.synchronize(CUDABackend())
 end
 
 
@@ -30,8 +30,8 @@ function Argos.copy_index!(dest::CuVector{T}, src::CuVector{T}, idx) where T
     @assert length(dest) == length(idx)
     ndrange = (length(dest),)
     idx_d = CuVector(idx)
-    ev = _copy_index_kernel!(CUDADevice())(dest, src, idx_d; ndrange=ndrange)
-    wait(ev)
+    _copy_index_kernel!(CUDABackend())(dest, src, idx_d; ndrange=ndrange)
+    KA.synchronize(CUDABackend())
 end
 
 
@@ -43,8 +43,8 @@ end
 function Argos.fixed_diag!(dest::CuMatrix, ind_fixed, val::Number)
     length(ind_fixed) == 0 && return
     g_ind_fixed = ind_fixed |> CuArray
-    ev = _fixed_diag_kernel!(CUDADevice())(dest, g_ind_fixed, val; ndrange=length(ind_fixed))
-    wait(ev)
+    _fixed_diag_kernel!(CUDABackend())(dest, g_ind_fixed, val; ndrange=length(ind_fixed))
+    KA.synchronize(CUDABackend())
 end
 
 @kernel function _transfer_auglag_hessian!(dest, H, J, ρ, n, m)
@@ -70,8 +70,8 @@ function Argos.transfer_auglag_hessian!(
     @assert size(ρ, 1) == m
 
     ndrange = (n+m, n)
-    ev = _transfer_auglag_hessian!(CUDADevice())(dest, H, J, ρ, n, m, ndrange=ndrange, dependencies=Event(CUDADevice()))
-    wait(ev)
+    _transfer_auglag_hessian!(CUDABackend())(dest, H, J, ρ, n, m, ndrange=ndrange)
+    KA.synchronize(CUDABackend())
     return
 end
 
@@ -84,11 +84,11 @@ function Argos.set_batch_tangents!(seeds::CuMatrix, offset, n, n_batches)
     @assert offset + n_batches <= n
     ndrange = (n_batches)
     fill!(seeds, 0.0)
-    ev = _batch_tangents_kernel!(CUDADevice())(
+    _batch_tangents_kernel!(CUDABackend())(
         seeds, offset, n_batches;
-        ndrange=ndrange, dependencies=Event(CUDADevice()),
+        ndrange=ndrange,
     )
-    wait(ev)
+    KA.synchronize(CUDABackend())
     return
 end
 
@@ -115,11 +115,11 @@ function Argos.tgtmul!(
     k = size(z, 2)
     ndrange = (n, k)
     y .*= beta
-    ev = _tgtmul_1_kernel!(CUDADevice())(
+    _tgtmul_1_kernel!(CUDABackend())(
         y, A.rowPtr, A.colVal, A.nzVal, z, w, alpha, nz, nw;
-        ndrange=ndrange, dependencies=Event(CUDADevice()),
+        ndrange=ndrange,
     )
-    wait(ev)
+    KA.synchronize(CUDABackend())
 end
 
 
@@ -150,11 +150,11 @@ function Argos.tgtmul!(
     ndrange = (n, k)
     yx .*= beta
     yu .*= beta
-    ev = _tgtmul_2_kernel!(CUDADevice())(
+    _tgtmul_2_kernel!(CUDABackend())(
         yx, yu, A.rowPtr, A.colVal, A.nzVal, z, w, alpha, nz, nw;
-        ndrange=ndrange, dependencies=Event(CUDADevice()),
+        ndrange=ndrange,
     )
-    wait(ev)
+    KA.synchronize(CUDABackend())
 end
 
 
@@ -171,11 +171,11 @@ end
 
 function Argos.update!(K::Argos.HJDJ, A, D, Σ)
     m = size(A, 1)
-    ev = _scale_transpose_kernel!(CUDADevice())(
+    ev = _scale_transpose_kernel!(CUDABackend())(
         K.Jt.nzVal, A.rowPtr, A.colVal, A.nzVal, D, K.transperm,
         ndrange=(m, 1),
     )
-    wait(ev)
+    KA.synchronize(ev)
     spgemm!('N', 'N', 1.0, K.Jt, A, 0.0, K.JtJ, 'O')
     K.Σ .= Σ
 end

lib/ArgosCUDA.jl/src/reduction.jl

@@ -24,21 +24,21 @@ end
 
 function Argos.update!(K::Argos.HJDJ, A, D, Σ)
     m = size(A, 1)
-    ev = _scale_transpose_kernel!(CUDADevice())(
+    _scale_transpose_kernel!(CUDABackend())(
         K.Jt.nzVal, A.rowPtr, A.colVal, A.nzVal, D, K.transperm,
         ndrange=(m, 1),
     )
-    wait(ev)
+    KA.synchronize(CUDABackend())
     spgemm!('N', 'N', 1.0, K.Jt, A, 0.0, K.JtJ, 'O')
     K.Σ .= Σ
 end
 function Argos.update!(K::Argos.HJDJ, A, D)
     m = size(A, 1)
-    ev = _scale_transpose_kernel!(CUDADevice())(
+    _scale_transpose_kernel!(CUDABackend())(
         K.Jt.nzVal, A.rowPtr, A.colVal, A.nzVal, D, K.transperm,
         ndrange=(m, 1),
     )
-    wait(ev)
+    KA.synchronize(CUDABackend())
     spgemm!('N', 'N', 1.0, K.Jt, A, 0.0, K.JtJ, 'O')
 end
 
@@ -47,7 +47,13 @@ function MadNLP.set_aug_diagonal!(kkt::Argos.BieglerKKTSystem{T, VI, VT, MT}, ip
     pr_diag_h = kkt.etc[:pr_diag_host]::Vector{T}
     # Broadcast is not working as MadNLP array are allocated on the CPU,
     # whereas pr_diag is allocated on the GPU
-    pr_diag_h .= ips.zl./(ips.x.-ips.xl) .+ ips.zu./(ips.xu.-ips.x)
+    x = MadNLP.full(ips.x)
+    xl = MadNLP.full(ips.xl)
+    xu = MadNLP.full(ips.xu)
+    zl = MadNLP.full(ips.zl)
+    zu = MadNLP.full(ips.zu)
+
+    pr_diag_h .= zl ./ (x .- xl) .+ zu ./ (xu .- x)
     copyto!(kkt.pr_diag, pr_diag_h)
     fill!(kkt.du_diag, 0.0)
 end