diff --git a/MathNet.Numerics.NativeProviders.sln b/MathNet.Numerics.NativeProviders.sln index 141b7074..a83444e8 100644 --- a/MathNet.Numerics.NativeProviders.sln +++ b/MathNet.Numerics.NativeProviders.sln @@ -20,6 +20,10 @@ Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "Numerics", "src\Numerics\Nu EndProject Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "UnitTests-MKL", "src\UnitTests\UnitTests-MKL.csproj", "{3515A344-AB5F-41C7-A14C-04A79B3FFAB1}" EndProject +Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "CUDA", "src\NativeProviders\Windows\CUDA\CUDAWrapper.vcxproj", "{5A52B796-7F41-4C90-8DE2-F3F391C4482C}" +EndProject +Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "UnitTests-CUDA", "src\UnitTests\UnitTests-CUDA.csproj", "{E79C0395-01DC-4BC9-B86C-ED45790892C5}" +EndProject Global GlobalSection(SolutionConfigurationPlatforms) = preSolution Debug|Any CPU = Debug|Any CPU @@ -30,6 +34,14 @@ Global Release|Mixed Platforms = Release|Mixed Platforms Release|Win32 = Release|Win32 Release|x64 = Release|x64 + Release-CUDA|Any CPU = Release-CUDA|Any CPU + Release-CUDA|Mixed Platforms = Release-CUDA|Mixed Platforms + Release-CUDA|Win32 = Release-CUDA|Win32 + Release-CUDA|x64 = Release-CUDA|x64 + Release-MKL|Any CPU = Release-MKL|Any CPU + Release-MKL|Mixed Platforms = Release-MKL|Mixed Platforms + Release-MKL|Win32 = Release-MKL|Win32 + Release-MKL|x64 = Release-MKL|x64 Release-Signed|Any CPU = Release-Signed|Any CPU Release-Signed|Mixed Platforms = Release-Signed|Mixed Platforms Release-Signed|Win32 = Release-Signed|Win32 @@ -50,6 +62,18 @@ Global {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release|Win32.Build.0 = Release|Win32 {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release|x64.ActiveCfg = Release|x64 {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release|x64.Build.0 = Release|x64 + {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release-CUDA|Any CPU.ActiveCfg = Release|Win32 + {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release-CUDA|Mixed Platforms.ActiveCfg = Release|Win32 + {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release-CUDA|Mixed Platforms.Build.0 = Release|Win32 + {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release-CUDA|Win32.ActiveCfg = Release|Win32 + {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release-CUDA|x64.ActiveCfg = Release|x64 + {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release-MKL|Any CPU.ActiveCfg = Release|Win32 + {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release-MKL|Mixed Platforms.ActiveCfg = Release|Win32 + {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release-MKL|Mixed Platforms.Build.0 = Release|Win32 + {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release-MKL|Win32.ActiveCfg = Release|Win32 + {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release-MKL|Win32.Build.0 = Release|Win32 + {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release-MKL|x64.ActiveCfg = Release|x64 + {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release-MKL|x64.Build.0 = Release|x64 {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release-Signed|Any CPU.ActiveCfg = Release|Win32 {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release-Signed|Mixed Platforms.ActiveCfg = Release|Win32 {C0B0DBA9-7FB0-4C87-BDB1-3EED19DC2B8F}.Release-Signed|Mixed Platforms.Build.0 = Release|Win32 @@ -65,6 +89,14 @@ Global {2362B8AC-C52B-45E4-A1BF-C682A4DB4220}.Release|Mixed Platforms.ActiveCfg = Release|Win32 {2362B8AC-C52B-45E4-A1BF-C682A4DB4220}.Release|Win32.ActiveCfg = Release|Win32 {2362B8AC-C52B-45E4-A1BF-C682A4DB4220}.Release|x64.ActiveCfg = Release|x64 + {2362B8AC-C52B-45E4-A1BF-C682A4DB4220}.Release-CUDA|Any CPU.ActiveCfg = Release|Win32 + {2362B8AC-C52B-45E4-A1BF-C682A4DB4220}.Release-CUDA|Mixed Platforms.ActiveCfg = Release|Win32 + {2362B8AC-C52B-45E4-A1BF-C682A4DB4220}.Release-CUDA|Win32.ActiveCfg = Release|Win32 + {2362B8AC-C52B-45E4-A1BF-C682A4DB4220}.Release-CUDA|x64.ActiveCfg = Release|x64 + {2362B8AC-C52B-45E4-A1BF-C682A4DB4220}.Release-MKL|Any CPU.ActiveCfg = Release|Win32 + {2362B8AC-C52B-45E4-A1BF-C682A4DB4220}.Release-MKL|Mixed Platforms.ActiveCfg = Release|Win32 + {2362B8AC-C52B-45E4-A1BF-C682A4DB4220}.Release-MKL|Win32.ActiveCfg = Release|Win32 + {2362B8AC-C52B-45E4-A1BF-C682A4DB4220}.Release-MKL|x64.ActiveCfg = Release|x64 {2362B8AC-C52B-45E4-A1BF-C682A4DB4220}.Release-Signed|Any CPU.ActiveCfg = Release|Win32 {2362B8AC-C52B-45E4-A1BF-C682A4DB4220}.Release-Signed|Mixed Platforms.ActiveCfg = Release|Win32 {2362B8AC-C52B-45E4-A1BF-C682A4DB4220}.Release-Signed|Win32.ActiveCfg = Release|Win32 @@ -87,6 +119,22 @@ Global {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release|Win32.Build.0 = Release|Any CPU {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release|x64.ActiveCfg = Release|Any CPU {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release|x64.Build.0 = Release|Any CPU + {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-CUDA|Any CPU.ActiveCfg = Release|Any CPU + {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-CUDA|Any CPU.Build.0 = Release|Any CPU + {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-CUDA|Mixed Platforms.ActiveCfg = Release|Any CPU + {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-CUDA|Mixed Platforms.Build.0 = Release|Any CPU + {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-CUDA|Win32.ActiveCfg = Release|Any CPU + {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-CUDA|Win32.Build.0 = Release|Any CPU + {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-CUDA|x64.ActiveCfg = Release|Any CPU + {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-CUDA|x64.Build.0 = Release|Any CPU + {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-MKL|Any CPU.ActiveCfg = Release|Any CPU + {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-MKL|Any CPU.Build.0 = Release|Any CPU + {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-MKL|Mixed Platforms.ActiveCfg = Release|Any CPU + {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-MKL|Mixed Platforms.Build.0 = Release|Any CPU + {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-MKL|Win32.ActiveCfg = Release|Any CPU + {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-MKL|Win32.Build.0 = Release|Any CPU + {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-MKL|x64.ActiveCfg = Release|Any CPU + {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-MKL|x64.Build.0 = Release|Any CPU {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-Signed|Any CPU.ActiveCfg = Release-Signed|Any CPU {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-Signed|Any CPU.Build.0 = Release-Signed|Any CPU {B7CAE5F4-A23F-4438-B5BE-41226618B695}.Release-Signed|Mixed Platforms.ActiveCfg = Release-Signed|Any CPU @@ -105,12 +153,90 @@ Global {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release|Mixed Platforms.Build.0 = Release|Any CPU {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release|Win32.ActiveCfg = Release|Any CPU {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release|x64.ActiveCfg = Release|Any CPU + {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release-CUDA|Any CPU.ActiveCfg = Release|Any CPU + {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release-CUDA|Any CPU.Build.0 = Release|Any CPU + {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release-CUDA|Mixed Platforms.ActiveCfg = Release|Any CPU + {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release-CUDA|Mixed Platforms.Build.0 = Release|Any CPU + {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release-CUDA|Win32.ActiveCfg = Release|Any CPU + {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release-CUDA|x64.ActiveCfg = Release|Any CPU + {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release-MKL|Any CPU.ActiveCfg = Release|Any CPU + {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release-MKL|Any CPU.Build.0 = Release|Any CPU + {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release-MKL|Mixed Platforms.ActiveCfg = Release|Any CPU + {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release-MKL|Mixed Platforms.Build.0 = Release|Any CPU + {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release-MKL|Win32.ActiveCfg = Release|Any CPU + {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release-MKL|Win32.Build.0 = Release|Any CPU + {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release-MKL|x64.ActiveCfg = Release|Any CPU + {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release-MKL|x64.Build.0 = Release|Any CPU {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release-Signed|Any CPU.ActiveCfg = Release|Any CPU {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release-Signed|Any CPU.Build.0 = Release|Any CPU {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release-Signed|Mixed Platforms.ActiveCfg = Release|Any CPU {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release-Signed|Mixed Platforms.Build.0 = Release|Any CPU {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release-Signed|Win32.ActiveCfg = Release|Any CPU {3515A344-AB5F-41C7-A14C-04A79B3FFAB1}.Release-Signed|x64.ActiveCfg = Release|Any CPU + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Debug|Any CPU.ActiveCfg = Debug|Win32 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Debug|Mixed Platforms.ActiveCfg = Debug|Win32 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Debug|Mixed Platforms.Build.0 = Debug|Win32 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Debug|Win32.ActiveCfg = Debug|Win32 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Debug|Win32.Build.0 = Debug|Win32 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Debug|x64.ActiveCfg = Debug|x64 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Debug|x64.Build.0 = Debug|x64 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release|Any CPU.ActiveCfg = Release|Win32 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release|Mixed Platforms.ActiveCfg = Release|Win32 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release|Mixed Platforms.Build.0 = Release|Win32 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release|Win32.ActiveCfg = Release|Win32 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release|Win32.Build.0 = Release|Win32 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release|x64.ActiveCfg = Release|x64 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release|x64.Build.0 = Release|x64 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release-CUDA|Any CPU.ActiveCfg = Release|Win32 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release-CUDA|Mixed Platforms.ActiveCfg = Release|Win32 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release-CUDA|Mixed Platforms.Build.0 = Release|Win32 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release-CUDA|Win32.ActiveCfg = Release|Win32 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release-CUDA|Win32.Build.0 = Release|Win32 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release-CUDA|x64.ActiveCfg = Release|x64 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release-CUDA|x64.Build.0 = Release|x64 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release-MKL|Any CPU.ActiveCfg = Release|Win32 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release-MKL|Mixed Platforms.ActiveCfg = Release|Win32 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release-MKL|Mixed Platforms.Build.0 = Release|Win32 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release-MKL|Win32.ActiveCfg = Release|Win32 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release-MKL|x64.ActiveCfg = Release|x64 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release-Signed|Any CPU.ActiveCfg = Release|Win32 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release-Signed|Mixed Platforms.ActiveCfg = Release|Win32 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release-Signed|Mixed Platforms.Build.0 = Release|Win32 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release-Signed|Win32.ActiveCfg = Release|Win32 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release-Signed|Win32.Build.0 = Release|Win32 + {5A52B796-7F41-4C90-8DE2-F3F391C4482C}.Release-Signed|x64.ActiveCfg = Release|Win32 + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Debug|Any CPU.ActiveCfg = Debug|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Debug|Any CPU.Build.0 = Debug|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Debug|Mixed Platforms.ActiveCfg = Debug|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Debug|Mixed Platforms.Build.0 = Debug|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Debug|Win32.ActiveCfg = Debug|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Debug|x64.ActiveCfg = Debug|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release|Any CPU.ActiveCfg = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release|Any CPU.Build.0 = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release|Mixed Platforms.ActiveCfg = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release|Mixed Platforms.Build.0 = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release|Win32.ActiveCfg = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release|x64.ActiveCfg = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-CUDA|Any CPU.ActiveCfg = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-CUDA|Any CPU.Build.0 = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-CUDA|Mixed Platforms.ActiveCfg = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-CUDA|Mixed Platforms.Build.0 = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-CUDA|Win32.ActiveCfg = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-CUDA|Win32.Build.0 = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-CUDA|x64.ActiveCfg = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-CUDA|x64.Build.0 = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-MKL|Any CPU.ActiveCfg = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-MKL|Any CPU.Build.0 = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-MKL|Mixed Platforms.ActiveCfg = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-MKL|Mixed Platforms.Build.0 = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-MKL|Win32.ActiveCfg = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-MKL|x64.ActiveCfg = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-Signed|Any CPU.ActiveCfg = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-Signed|Any CPU.Build.0 = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-Signed|Mixed Platforms.ActiveCfg = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-Signed|Mixed Platforms.Build.0 = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-Signed|Win32.ActiveCfg = Release|Any CPU + {E79C0395-01DC-4BC9-B86C-ED45790892C5}.Release-Signed|x64.ActiveCfg = Release|Any CPU EndGlobalSection GlobalSection(SolutionProperties) = preSolution HideSolutionNode = FALSE diff --git a/MathNet.Numerics.sln b/MathNet.Numerics.sln index 58787f08..4eb3bebe 100644 --- a/MathNet.Numerics.sln +++ b/MathNet.Numerics.sln @@ -9,8 +9,10 @@ Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "Readme", "Readme", "{C2F374 CONTRIBUTORS.md = CONTRIBUTORS.md LICENSE.md = LICENSE.md README.md = README.md + RELEASENOTES-CUDA.md = RELEASENOTES-CUDA.md RELEASENOTES-Data.md = RELEASENOTES-Data.md RELEASENOTES-MKL.md = RELEASENOTES-MKL.md + RELEASENOTES-OpenBLAS.md = RELEASENOTES-OpenBLAS.md RELEASENOTES.md = RELEASENOTES.md EndProjectSection EndProject @@ -27,6 +29,7 @@ Project("{2150E333-8FDC-42A3-9474-1A3956D46DE8}") = "Build", "Build", "{A4A66FA9 docs\tools\build-docs.fsx = docs\tools\build-docs.fsx build.fsx = build.fsx build\MathNet.Numerics.Extension.nuspec = build\MathNet.Numerics.Extension.nuspec + build\MathNet.Numerics.MKL.Win.targets = build\MathNet.Numerics.MKL.Win.targets build\MathNet.Numerics.nuspec = build\MathNet.Numerics.nuspec paket.dependencies = paket.dependencies paket.lock = paket.lock diff --git a/RELEASENOTES-CUDA.md b/RELEASENOTES-CUDA.md new file mode 100644 index 00000000..0ab166d4 --- /dev/null +++ b/RELEASENOTES-CUDA.md @@ -0,0 +1,3 @@ +### 0.1.0-alpha - TBA +* With Nvidia CUDA 7.0.28 +* Initial version diff --git a/RELEASENOTES-OpenBLAS.md b/RELEASENOTES-OpenBLAS.md new file mode 100644 index 00000000..2951cc2c --- /dev/null +++ b/RELEASENOTES-OpenBLAS.md @@ -0,0 +1,2 @@ +### 0.1.0-alpha - TBA +* Initial version diff --git a/build.fsx b/build.fsx index 505d73fd..f0a7eade 100644 --- a/build.fsx +++ b/build.fsx @@ -51,6 +51,13 @@ let mklPackageVersion = mklRelease.NugetVersion let mklReleaseNotes = mklRelease.Notes |> List.map (fun l -> l.Replace("*","").Replace("`","")) |> toLines trace (sprintf " Math.NET Numerics MKL Provider v%s" mklPackageVersion) +let cudaRelease = LoadReleaseNotes "RELEASENOTES-CUDA.md" +let cudaBuildPart = "0" +let cudaAssemblyVersion = cudaRelease.AssemblyVersion + "." + cudaBuildPart +let cudaPackageVersion = cudaRelease.NugetVersion +let cudaReleaseNotes = cudaRelease.Notes |> List.map (fun l -> l.Replace("*","").Replace("`","")) |> toLines +trace (sprintf " Math.NET Numerics CUDA Provider v%s" cudaPackageVersion) + let dataRelease = LoadReleaseNotes "RELEASENOTES-Data.md" let dataBuildPart = "0" let dataAssemblyVersion = dataRelease.AssemblyVersion + "." + dataBuildPart @@ -187,10 +194,10 @@ let coreSignedBundle = // NATIVE PROVIDER PACKAGES -let mklWin32Pack = - { Id = "MathNet.Numerics.MKL.Win-x86" +let mklWinPack = + { Id = "MathNet.Numerics.MKL.Win" Version = mklPackageVersion - Title = "Math.NET Numerics - MKL Native Provider for Windows (x86)" + Title = "Math.NET Numerics - MKL Native Provider for Windows (x64 and x86)" Summary = "" Description = "Intel MKL native libraries for Math.NET Numerics. Requires an Intel MKL license if redistributed." ReleaseNotes = mklReleaseNotes @@ -198,15 +205,32 @@ let mklWin32Pack = Authors = [ "Christoph Ruegg"; "Marcus Cuda"; "Jurgen Van Gael" ] Dependencies = [ { FrameworkVersion="" - Dependencies=[ "MathNet.Numerics", "2.4.0" ] } ] + Dependencies=[ "MathNet.Numerics", "3.6.0" ] } ] Files = - [ @"..\..\out\MKL\Windows\x86\libiomp5md.dll", Some "content", None; - @"..\..\out\MKL\Windows\x86\MathNet.Numerics.MKL.dll", Some "content", None ] } + [ @"MathNet.Numerics.MKL.Win.targets", Some "build", None; + @"..\..\out\MKL\Windows\x64\libiomp5md.dll", Some "build\x64", None; + @"..\..\out\MKL\Windows\x64\MathNet.Numerics.MKL.dll", Some "build\x64", None; + @"..\..\out\MKL\Windows\x86\libiomp5md.dll", Some "build\x86", None; + @"..\..\out\MKL\Windows\x86\MathNet.Numerics.MKL.dll", Some "build\x86", None ] } + +let mklWin32Pack = + { mklWinPack with + Id = "MathNet.Numerics.MKL.Win-x86" + Title = "Math.NET Numerics - MKL Native Provider for Windows (x86)" + Dependencies = + [ { FrameworkVersion="" + Dependencies=[ "MathNet.Numerics", "2.4.0" ] } ] + Files = + [ @"..\..\out\MKL\Windows\x86\libiomp5md.dll", Some "content", None; + @"..\..\out\MKL\Windows\x86\MathNet.Numerics.MKL.dll", Some "content", None ] } let mklWin64Pack = - { mklWin32Pack with + { mklWinPack with Id = "MathNet.Numerics.MKL.Win-x64" Title = "Math.NET Numerics - MKL Native Provider for Windows (x64)" + Dependencies = + [ { FrameworkVersion="" + Dependencies=[ "MathNet.Numerics", "2.4.0" ] } ] Files = [ @"..\..\out\MKL\Windows\x64\libiomp5md.dll", Some "content", None; @"..\..\out\MKL\Windows\x64\MathNet.Numerics.MKL.dll", Some "content", None ] } @@ -241,7 +265,7 @@ let mklWinBundle = Title = "Math.NET Numerics MKL Native Provider for Windows" ReleaseNotesFile = "RELEASENOTES-MKL.md" FsLoader = false - Packages = [ mklWin32Pack; mklWin64Pack ] } + Packages = [ mklWinPack; mklWin32Pack; mklWin64Pack ] } let mklLinuxBundle = { Id = "MathNet.Numerics.MKL.Linux" @@ -251,6 +275,37 @@ let mklLinuxBundle = FsLoader = false Packages = [ mklLinux32Pack; mklLinux64Pack ] } +let cudaWinPack = + { Id = "MathNet.Numerics.CUDA.Win" + Version = cudaPackageVersion + Title = "Math.NET Numerics - CUDA Native Provider for Windows (x64 and x86)" + Summary = "" + Description = "Nvidia CUDA native libraries for Math.NET Numerics." + ReleaseNotes = cudaReleaseNotes + Tags = "math numeric statistics probability integration interpolation linear algebra matrix fft native cuda gpu" + Authors = [ "Matthew A Johnson"; "Christoph Ruegg" ] + Dependencies = + [ { FrameworkVersion="" + Dependencies=[ "MathNet.Numerics", "3.7.0" ] } ] + Files = + [ @"MathNet.Numerics.CUDA.Win.targets", Some "build", None; + @"..\..\out\CUDA\Windows\x64\cublas64_70.dll", Some "content", None; + @"..\..\out\CUDA\Windows\x64\cudart64_70.dll", Some "content", None; + @"..\..\out\CUDA\Windows\x64\cusolver64_70.dll", Some "content", None; + @"..\..\out\CUDA\Windows\x64\MathNet.Numerics.CUDA.dll", Some "content", None + @"..\..\out\CUDA\Windows\x86\cublas32_70.dll", Some "content", None; + @"..\..\out\CUDA\Windows\x86\cudart32_70.dll", Some "content", None; + @"..\..\out\CUDA\Windows\x86\cusolver32_70.dll", Some "content", None; + @"..\..\out\CUDA\Windows\x86\MathNet.Numerics.CUDA.dll", Some "content", None ] } + +let cudaWinBundle = + { Id = "MathNet.Numerics.CUDA.Win" + Version = mklPackageVersion + Title = "Math.NET Numerics CUDA Native Provider for Windows" + ReleaseNotesFile = "RELEASENOTES-CUDA.md" + FsLoader = false + Packages = [ cudaWinPack ] } + // DATA EXTENSION PACKAGES @@ -309,7 +364,7 @@ Target "Clean" (fun _ -> CleanDirs [ "out/lib-debug/Net35"; "out/lib-debug/Net40"; "out/lib-debug/Profile7"; "out/lib-debug/Profile47"; "out/lib-debug/Profile78"; "out/lib-debug/Profile259"; "out/lib-debug/Profile328" ] CleanDirs [ "out/test-debug/Net35"; "out/test-debug/Net40"; "out/test-debug/Profile7"; "out/test-debug/Profile47"; "out/test-debug/Profile78"; "out/test-debug/Profile259"; "out/test-debug/Profile328" ] CleanDirs [ "out/lib-signed/Net40"; "out/test-signed/Net40" ] // Signed Build - CleanDirs [ "out/MKL"; "out/ATLAS" ] // Native Providers + CleanDirs [ "out/MKL"; "out/ATLAS"; "out/CUDA"; "out/OpenBLAS" ] // Native Providers CleanDirs [ "out/Data" ]) // Data Extensions Target "ApplyVersion" (fun _ -> @@ -328,7 +383,11 @@ Target "ApplyVersion" (fun _ -> ReplaceInFile (regex_replace @"\d+\.\d+\.\d+\.\d+" mklAssemblyVersion >> regex_replace @"\d+,\d+,\d+,\d+" (replace "." "," mklAssemblyVersion)) - "src/NativeProviders/Common/resource.rc") + "src/NativeProviders/MKL/resource.rc" + ReplaceInFile + (regex_replace @"\d+\.\d+\.\d+\.\d+" cudaAssemblyVersion + >> regex_replace @"\d+,\d+,\d+,\d+" (replace "." "," cudaAssemblyVersion)) + "src/NativeProviders/CUDA/resource.rc") Target "Prepare" DoNothing "Start" @@ -344,8 +403,8 @@ Target "Prepare" DoNothing let buildConfig config subject = MSBuild "" (if hasBuildParam "incremental" then "Build" else "Rebuild") [ "Configuration", config ] subject |> ignore let build subject = buildConfig "Release" subject let buildSigned subject = buildConfig "Release-Signed" subject -let nativeWin32Build subject = MSBuild "" (if hasBuildParam "incremental" then "Build" else "Rebuild") [("Configuration","Release"); ("Platform","Win32")] subject |> ignore -let nativeWin64Build subject = MSBuild "" (if hasBuildParam "incremental" then "Build" else "Rebuild") [("Configuration","Release"); ("Platform","x64")] subject |> ignore +let buildConfig32 config subject = MSBuild "" (if hasBuildParam "incremental" then "Build" else "Rebuild") [("Configuration", config); ("Platform","Win32")] subject |> ignore +let buildConfig64 config subject = MSBuild "" (if hasBuildParam "incremental" then "Build" else "Rebuild") [("Configuration", config); ("Platform","x64")] subject |> ignore Target "BuildMain" (fun _ -> build !! "MathNet.Numerics.sln") Target "BuildNet35" (fun _ -> build !! "MathNet.Numerics.Net35Only.sln") @@ -360,13 +419,18 @@ Target "Build" DoNothing =?> ("BuildMain", not (hasBuildParam "all" || hasBuildParam "release" || hasBuildParam "net35" || hasBuildParam "signed")) ==> "Build" -Target "MklWin32Build" (fun _ -> nativeWin32Build !! "MathNet.Numerics.NativeProviders.sln") -Target "MklWin64Build" (fun _ -> nativeWin64Build !! "MathNet.Numerics.NativeProviders.sln") - +Target "MklWin32Build" (fun _ -> buildConfig32 "Release-MKL" !! "MathNet.Numerics.NativeProviders.sln") +Target "MklWin64Build" (fun _ -> buildConfig64 "Release-MKL" !! "MathNet.Numerics.NativeProviders.sln") Target "MklWinBuild" DoNothing "Prepare" ==> "MklWin32Build" ==> "MklWinBuild" "Prepare" ==> "MklWin64Build" ==> "MklWinBuild" +Target "CudaWin32Build" (fun _ -> buildConfig32 "Release-CUDA" !! "MathNet.Numerics.NativeProviders.sln") +Target "CudaWin64Build" (fun _ -> buildConfig64 "Release-CUDA" !! "MathNet.Numerics.NativeProviders.sln") +Target "CudaWinBuild" DoNothing +"Prepare" ==> "CudaWin32Build" ==> "CudaWinBuild" +"Prepare" ==> "CudaWin64Build" ==> "CudaWinBuild" + Target "DataBuild" (fun _ -> build !! "MathNet.Numerics.Data.sln") "Prepare" ==> "DataBuild" diff --git a/build/MathNet.Numerics.CUDA.Win.targets b/build/MathNet.Numerics.CUDA.Win.targets new file mode 100644 index 00000000..74be39a2 --- /dev/null +++ b/build/MathNet.Numerics.CUDA.Win.targets @@ -0,0 +1,74 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + $(BuildDependsOn); + CopyMathNetInteropFiles; + + + $(CleanDependsOn); + CleanMathNetInteropFiles; + + + diff --git a/build/MathNet.Numerics.MKL.Win.targets b/build/MathNet.Numerics.MKL.Win.targets new file mode 100644 index 00000000..870de2b4 --- /dev/null +++ b/build/MathNet.Numerics.MKL.Win.targets @@ -0,0 +1,74 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + $(BuildDependsOn); + CopyMathNetInteropFiles; + + + $(CleanDependsOn); + CleanMathNetInteropFiles; + + + diff --git a/src/NativeProviders/ATLAS/resource.h b/src/NativeProviders/ATLAS/resource.h new file mode 100644 index 00000000..27e2900c --- /dev/null +++ b/src/NativeProviders/ATLAS/resource.h @@ -0,0 +1,14 @@ +//{{NO_DEPENDENCIES}} +// Microsoft Visual C++ generated include file. +// Used by resource.rc + +// Next default values for new objects +// +#ifdef APSTUDIO_INVOKED +#ifndef APSTUDIO_READONLY_SYMBOLS +#define _APS_NEXT_RESOURCE_VALUE 101 +#define _APS_NEXT_COMMAND_VALUE 40001 +#define _APS_NEXT_CONTROL_VALUE 1001 +#define _APS_NEXT_SYMED_VALUE 101 +#endif +#endif diff --git a/src/NativeProviders/ATLAS/resource.rc b/src/NativeProviders/ATLAS/resource.rc new file mode 100644 index 00000000..b49bd076 --- /dev/null +++ b/src/NativeProviders/ATLAS/resource.rc @@ -0,0 +1,101 @@ +// Microsoft Visual C++ generated resource script. +// +#include "resource.h" + +#define APSTUDIO_READONLY_SYMBOLS +///////////////////////////////////////////////////////////////////////////// +// +// Generated from the TEXTINCLUDE 2 resource. +// +#include "windows.h" + +///////////////////////////////////////////////////////////////////////////// +#undef APSTUDIO_READONLY_SYMBOLS + +///////////////////////////////////////////////////////////////////////////// +// English (United States) resources + +#if !defined(AFX_RESOURCE_DLL) || defined(AFX_TARG_ENU) +LANGUAGE LANG_ENGLISH, SUBLANG_ENGLISH_US +#pragma code_page(1252) + +#ifdef APSTUDIO_INVOKED +///////////////////////////////////////////////////////////////////////////// +// +// TEXTINCLUDE +// + +1 TEXTINCLUDE +BEGIN + "resource.h\0" +END + +2 TEXTINCLUDE +BEGIN + "#include ""windows.h""\r\n" + "\0" +END + +3 TEXTINCLUDE +BEGIN + "\r\n" + "\0" +END + +#endif // APSTUDIO_INVOKED + + +///////////////////////////////////////////////////////////////////////////// +// +// Version +// + +VS_VERSION_INFO VERSIONINFO + FILEVERSION 0,1,0,0 + PRODUCTVERSION 0,1,0,0 + FILEFLAGSMASK 0x17L +#ifdef _DEBUG + FILEFLAGS 0x1L +#else + FILEFLAGS 0x0L +#endif + FILEOS 0x4L + FILETYPE 0x2L + FILESUBTYPE 0x0L +BEGIN + BLOCK "StringFileInfo" + BEGIN + BLOCK "040904b0" + BEGIN + VALUE "Comments", "http://numerics.mathdotnet.com/" + VALUE "CompanyName", "Math.NET" + VALUE "FileDescription", "MathNET Numerics ATLAS Native Provider" + VALUE "FileVersion", "0.1.0.0" + VALUE "InternalName", "Math.NET" + VALUE "LegalCopyright", "Copyright (C) Math.NET 2009-2015" + VALUE "OriginalFilename", "MathNet.Numerics.ATLAS" + VALUE "ProductName", "Math.NET Numerics" + VALUE "ProductVersion", "0.1.0.0" + END + END + BLOCK "VarFileInfo" + BEGIN + VALUE "Translation", 0x409, 1200 + END +END + +#endif // English (United States) resources +///////////////////////////////////////////////////////////////////////////// + + + +#ifndef APSTUDIO_INVOKED +///////////////////////////////////////////////////////////////////////////// +// +// Generated from the TEXTINCLUDE 3 resource. +// + + +///////////////////////////////////////////////////////////////////////////// +#endif // not APSTUDIO_INVOKED + diff --git a/src/NativeProviders/CUDA/blas.cpp b/src/NativeProviders/CUDA/blas.cpp new file mode 100644 index 00000000..b99559b6 --- /dev/null +++ b/src/NativeProviders/CUDA/blas.cpp @@ -0,0 +1,168 @@ +#include +#include "cublas_v2.h" +#include "cuda_runtime.h" +#include "wrapper_common.h" + +template +void cuda_axpy(const cublasHandle_t blasHandle, const int n, const T alpha, const T x[], int incX, T y[], int incY, AXPY axpy) +{ + T *d_X = NULL; + T *d_Y = NULL; + cudaMalloc((void**)&d_X, n*sizeof(T)); + cudaMalloc((void**)&d_Y, n*sizeof(T)); + + cublasSetVector(n, sizeof(T), x, incX, d_X, incX); + cublasSetVector(n, sizeof(T), y, incY, d_Y, incY); + + axpy(blasHandle, n, &alpha, d_X, incX, d_Y, incX); + + cublasGetVector(n, sizeof(T), d_Y, incY, y, incY); + + cudaFree(d_X); + cudaFree(d_Y); +} + +template +void cuda_scal(const cublasHandle_t blasHandle, const int n, const T alpha, T x[], int incX, SCAL scal) +{ + T *d_X = NULL; + cudaMalloc((void**)&d_X, n*sizeof(T)); + + cublasSetVector(n, sizeof(T), x, incX, d_X, incX); + + scal(blasHandle, n, &alpha, d_X, incX); + + cublasGetVector(n, sizeof(T), d_X, incX, x, incX); + + cudaFree(d_X); +} + +template +void cuda_dot(const cublasHandle_t blasHandle, const int n, const T x[], int incX, const T y[], int incY, T* result, DOT dot) +{ + T *d_X = NULL; + T *d_Y = NULL; + cudaMalloc((void**)&d_X, n*sizeof(T)); + cudaMalloc((void**)&d_Y, n*sizeof(T)); + + cublasSetVector(n, sizeof(T), x, incX, d_X, incX); + cublasSetVector(n, sizeof(T), y, incY, d_Y, incY); + + dot(blasHandle, n, d_X, incX, d_Y, incY, result); + + cudaFree(d_X); + cudaFree(d_Y); +} + +template +void cuda_gemm(const cublasHandle_t handle, const cublasOperation_t transa, const cublasOperation_t transb, int m, int n, int k, const T alpha, const T A[], int lda, const T B[], int ldb, const T beta, T C[], int ldc, GEMM gemm) +{ + T *d_A = NULL; + cudaMalloc((void**)&d_A, m*k*sizeof(T)); + cublasSetMatrix(m, k, sizeof(T), A, m, d_A, m); + + T *d_B = NULL; + cudaMalloc((void**)&d_B, k*n*sizeof(T)); + cublasSetMatrix(k, n, sizeof(T), B, k, d_B, k); + + T *d_C = NULL; + cudaMalloc((void**)&d_C, m*n*sizeof(T)); + cublasSetMatrix(m, n, sizeof(T), C, m, d_C, m); + + gemm(handle, transa, transb, m, n, k, &alpha, d_A, lda, d_B, ldb, &beta, d_C, ldc); + + cublasGetMatrix(m, n, sizeof(T), d_C, m, C, m); + + cudaFree(d_A); + cudaFree(d_B); + cudaFree(d_C); +} + +extern "C" { + + DLLEXPORT void s_axpy(const cublasHandle_t blasHandle, const int n, const float alpha, const float x[], float y[]){ + cuda_axpy(blasHandle, n, alpha, x, 1, y, 1, cublasSaxpy); + } + + DLLEXPORT void d_axpy(const cublasHandle_t blasHandle, const int n, const double alpha, const double x[], double y[]){ + cuda_axpy(blasHandle, n, alpha, x, 1, y, 1, cublasDaxpy); + } + + DLLEXPORT void c_axpy(const cublasHandle_t blasHandle, const int n, const cuComplex alpha, const cuComplex x[], cuComplex y[]){ + cuda_axpy(blasHandle, n, alpha, x, 1, y, 1, cublasCaxpy); + } + + DLLEXPORT void z_axpy(const cublasHandle_t blasHandle, const int n, const cuDoubleComplex alpha, const cuDoubleComplex x[], cuDoubleComplex y[]){ + cuda_axpy(blasHandle, n, alpha, x, 1, y, 1, cublasZaxpy); + } + + DLLEXPORT void s_scale(const cublasHandle_t blasHandle, const int n, const float alpha, float x[]){ + cuda_scal(blasHandle, n, alpha, x, 1, cublasSscal); + } + + DLLEXPORT void d_scale(const cublasHandle_t blasHandle, const int n, const double alpha, double x[]){ + cuda_scal(blasHandle, n, alpha, x, 1, cublasDscal); + } + + DLLEXPORT void c_scale(const cublasHandle_t blasHandle, const int n, const cuComplex alpha, cuComplex x[]){ + cuda_scal(blasHandle, n, alpha, x, 1, cublasCscal); + } + + DLLEXPORT void z_scale(const cublasHandle_t blasHandle, const int n, const cuDoubleComplex alpha, cuDoubleComplex x[]){ + cuda_scal(blasHandle, n, alpha, x, 1, cublasZscal); + } + + DLLEXPORT float s_dot_product(const cublasHandle_t blasHandle, const int n, const float x[], const float y[]){ + float ret; + cuda_dot(blasHandle, n, x, 1, y, 1, &ret, cublasSdot); + return ret; + } + + DLLEXPORT double d_dot_product(const cublasHandle_t blasHandle, const int n, const double x[], const double y[]){ + double ret; + cuda_dot(blasHandle, n, x, 1, y, 1, &ret, cublasDdot); + return ret; + } + + DLLEXPORT cuComplex c_dot_product(const cublasHandle_t blasHandle, const int n, const cuComplex x[], const cuComplex y[]){ + cuComplex ret; + cuda_dot(blasHandle, n, x, 1, y, 1, &ret, cublasCdotu); + return ret; + } + + DLLEXPORT cuDoubleComplex z_dot_product(const cublasHandle_t blasHandle, const int n, const cuDoubleComplex x[], const cuDoubleComplex y[]){ + cuDoubleComplex ret; + cuda_dot(blasHandle, n, x, 1, y, 1, &ret, cublasZdotu); + return ret; + } + + DLLEXPORT void s_matrix_multiply(const cublasHandle_t blasHandle, cublasOperation_t transA, cublasOperation_t transB, const int m, const int n, const int k, const float alpha, const float x[], const float y[], const float beta, float c[]){ + int lda = transA == CUBLAS_OP_N ? m : k; + int ldb = transB == CUBLAS_OP_N ? k : n; + + cuda_gemm(blasHandle, transA, transB, m, n, k, alpha, x, lda, y, ldb, beta, c, m, cublasSgemm); + } + + DLLEXPORT void d_matrix_multiply(const cublasHandle_t blasHandle, cublasOperation_t transA, cublasOperation_t transB, const int m, const int n, const int k, const double alpha, const double x[], const double y[], const double beta, double c[]){ + int lda = transA == CUBLAS_OP_N ? m : k; + int ldb = transB == CUBLAS_OP_N ? k : n; + + cuda_gemm(blasHandle, transA, transB, m, n, k, alpha, x, lda, y, ldb, beta, c, m, cublasDgemm); + } + + DLLEXPORT void c_matrix_multiply(const cublasHandle_t blasHandle, cublasOperation_t transA, cublasOperation_t transB, const int m, const int n, const int k, const cuComplex alpha, const cuComplex x[], const cuComplex y[], const cuComplex beta, cuComplex c[]){ + int lda = transA == CUBLAS_OP_N ? m : k; + int ldb = transB == CUBLAS_OP_N ? k : n; + + cuda_gemm(blasHandle, transA, transB, m, n, k, alpha, x, lda, y, ldb, beta, c, m, cublasCgemm); + } + + DLLEXPORT void z_matrix_multiply(const cublasHandle_t blasHandle, cublasOperation_t transA, cublasOperation_t transB, const int m, const int n, const int k, const cuDoubleComplex alpha, const cuDoubleComplex x[], const cuDoubleComplex y[], const cuDoubleComplex beta, cuDoubleComplex c[]){ + int lda = transA == CUBLAS_OP_N ? m : k; + int ldb = transB == CUBLAS_OP_N ? k : n; + + cuda_gemm(blasHandle, transA, transB, m, n, k, alpha, x, lda, y, ldb, beta, c, m, cublasZgemm); + } + +} + diff --git a/src/NativeProviders/CUDA/capabilities.cpp b/src/NativeProviders/CUDA/capabilities.cpp new file mode 100644 index 00000000..75fb1d59 --- /dev/null +++ b/src/NativeProviders/CUDA/capabilities.cpp @@ -0,0 +1,79 @@ +#include "wrapper_common.h" +#include "cublas_v2.h" +#include "cusolverDn.h" + +#ifdef __cplusplus +extern "C" { +#endif /* __cplusplus */ + + /* + Capability is supported if >0 + + Actual number can be increased over time to indicate + extensions/revisions (that do not break compatibility) + */ + DLLEXPORT int query_capability(const int capability) + { + switch (capability) + { + + // SANITY CHECKS + case 0: return 0; + case 1: return -1; + + // PLATFORM + case 8: +#ifdef _M_IX86 + return 1; +#else + return 0; +#endif + case 9: +#ifdef _M_X64 + return 1; +#else + return 0; +#endif + case 10: +#ifdef _M_IA64 + return 1; +#else + return 0; +#endif + + // COMMON/SHARED + case 64: return 1; // revision + + // LINEAR ALGEBRA + case 128: return 1; // basic dense linear algebra + + // OPTIMIZATION + case 256: return 0; // basic optimization + + // FFT + case 384: return 0; // basic FFT + + default: return 0; // unknown or not supported + + } + } + + DLLEXPORT cublasStatus_t createBLASHandle(cublasHandle_t *blasHandle){ + return cublasCreate(blasHandle); + } + + DLLEXPORT cublasStatus_t destroyBLASHandle(cublasHandle_t blasHandle){ + return cublasDestroy(blasHandle); + } + + DLLEXPORT cusolverStatus_t createSolverHandle(cusolverDnHandle_t *solverHandle){ + return cusolverDnCreate(solverHandle); + } + + DLLEXPORT cusolverStatus_t destroySolverHandle(cusolverDnHandle_t solverHandle){ + return cusolverDnDestroy(solverHandle); + } + +#ifdef __cplusplus +} +#endif /* __cplusplus */ diff --git a/src/NativeProviders/CUDA/lapack.cpp b/src/NativeProviders/CUDA/lapack.cpp new file mode 100644 index 00000000..a7bbbfbd --- /dev/null +++ b/src/NativeProviders/CUDA/lapack.cpp @@ -0,0 +1,1041 @@ +#include + +#include "lapack_common.h" +#include "wrapper_common.h" +#include "cublas_v2.h" +#include "cusolverDn.h" +#include "cuda_runtime.h" + +template +inline int lu_factor(cusolverDnHandle_t solverHandle, int m, T a[], int ipiv[], GETRF getrf, GETRFBSIZE getrfbsize) +{ + int info = 0; + + T* d_A = NULL; + cudaMalloc((void**)&d_A, m*m*sizeof(T)); + cublasSetMatrix(m, m, sizeof(T), a, m, d_A, m); + + int* d_I = NULL; + cudaMalloc((void**)&d_I, m*sizeof(int)); + + T* work = NULL; + int lwork = 0; + getrfbsize(solverHandle, m, m, a, m, &lwork); + cudaMalloc((void**)&work, sizeof(T)*lwork); + + int* d_info = NULL; + cudaMalloc((void**)&d_info, sizeof(int)); + + getrf(solverHandle, m, m, d_A, m, work, d_I, d_info); + + cudaMemcpy(&info, d_info, sizeof(int), cudaMemcpyDeviceToHost); + + cublasGetMatrix(m, m, sizeof(T), d_A, m, a, m); + cublasGetVector(m, sizeof(int), d_I, 1, ipiv, 1); + + shift_ipiv_down(m, ipiv); + + cudaFree(d_A); + cudaFree(d_I); + cudaFree(d_info); + cudaFree(work); + + return info; +}; + +template +inline int lu_inverse(cusolverDnHandle_t solverHandle, cublasHandle_t blasHandle, int n, T a[], GETRF getrf, GETRIBATCHED getribatched, GETRFBSIZE getrfbsize) +{ + int info = 0; + + int* d_I = NULL; + cudaMalloc((void**)&d_I, n*sizeof(int)); + + T* d_A = NULL; + cudaMalloc((void**)&d_A, n*n*sizeof(T)); + cublasSetMatrix(n, n, sizeof(T), a, n, d_A, n); + + T* work = NULL; + int lwork = 0; + getrfbsize(solverHandle, n, n, d_A, n, &lwork); + cudaMalloc((void**)&work, sizeof(T)*lwork); + + int* d_info = NULL; + cudaMalloc((void**)&d_info, sizeof(int)); + + getrf(solverHandle, n, n, d_A, n, work, d_I, d_info); + cudaMemcpy(&info, d_info, sizeof(int), cudaMemcpyDeviceToHost); + + cudaFree(work); + + if (info != 0) + { + cudaFree(d_A); + cudaFree(d_I); + cudaFree(d_info); + return info; + } + + T* d_C = NULL; + cudaMalloc((void**)&d_C, n*n*sizeof(T)); + + const T **d_Aarray = NULL; + cudaMalloc((void**)&d_Aarray, sizeof(T*)); + cudaMemcpy(d_Aarray, &d_A, sizeof(T*), cudaMemcpyHostToDevice); + + T **d_Carray = NULL; + cudaMalloc((void**)&d_Carray, sizeof(T*)); + cudaMemcpy(d_Carray, &d_C, sizeof(T*), cudaMemcpyHostToDevice); + + getribatched(blasHandle, n, d_Aarray, n, d_I, d_Carray, n, d_info, 1); + cudaMemcpy(&info, d_info, sizeof(int), cudaMemcpyDeviceToHost); + + cublasGetMatrix(n, n, sizeof(T), d_C, n, a, n); + + cudaFree(d_A); + cudaFree(d_I); + cudaFree(d_C); + cudaFree(d_info); + cudaFree(d_Aarray); + cudaFree(d_Carray); + + return info; +}; + +template +inline int lu_inverse_factored(cublasHandle_t blasHandle, int n, T a[], int ipiv[], GETRI getri) +{ + int info = 0; + + shift_ipiv_up(n, ipiv); + + T* d_A = NULL; + cudaMalloc((void**)&d_A, n*n*sizeof(T)); + cublasSetMatrix(n, n, sizeof(T), a, n, d_A, n); + + T* d_C = NULL; + cudaMalloc((void**)&d_C, n*n*sizeof(T)); + + int* d_I = NULL; + cudaMalloc((void**)&d_I, n*sizeof(int)); + cublasSetVector(n, sizeof(int), ipiv, 1, d_I, 1); + + int* d_info = NULL; + cudaMalloc((void**)&d_info, sizeof(int)); + + const T **d_Aarray = NULL; + cudaMalloc((void**)&d_Aarray, sizeof(T*)); + cudaMemcpy(d_Aarray, &d_A, sizeof(T*), cudaMemcpyHostToDevice); + + T **d_Carray = NULL; + cudaMalloc((void**)&d_Carray, sizeof(T*)); + cudaMemcpy(d_Carray, &d_C, sizeof(T*), cudaMemcpyHostToDevice); + + getri(blasHandle, n, d_Aarray, n, d_I, d_Carray, n, d_info, 1); + cudaMemcpy(&info, d_info, sizeof(int), cudaMemcpyDeviceToHost); + + cublasGetMatrix(n, n, sizeof(T), d_C, n, a, n); + cublasGetVector(n, sizeof(int), d_I, 1, ipiv, 1); + + shift_ipiv_down(n, ipiv); + + cudaFree(d_A); + cudaFree(d_I); + cudaFree(d_C); + cudaFree(d_info); + cudaFree(d_Aarray); + cudaFree(d_Carray); + + return info; +} + +template +inline int lu_solve_factored(cusolverDnHandle_t solverHandle, int n, int nrhs, T a[], int ipiv[], T b[], GETRS getrs) +{ + int info = 0; + + shift_ipiv_up(n, ipiv); + + T* d_A = NULL; + cudaMalloc((void**)&d_A, n*n*sizeof(T)); + cublasSetMatrix(n, n, sizeof(T), a, n, d_A, n); + + T* d_B = NULL; + cudaMalloc((void**)&d_B, n*nrhs*sizeof(T)); + cublasSetMatrix(n, nrhs, sizeof(T), b, n, d_B, n); + + int* d_I = NULL; + cudaMalloc((void**)&d_I, n*sizeof(int)); + cublasSetVector(n, sizeof(int), ipiv, 1, d_I, 1); + + int* d_info = NULL; + cudaMalloc((void**)&d_info, sizeof(int)); + + getrs(solverHandle, CUBLAS_OP_N, n, nrhs, d_A, n, d_I, d_B, n, d_info); + cudaMemcpy(&info, d_info, sizeof(int), cudaMemcpyDeviceToHost); + + cublasGetMatrix(n, nrhs, sizeof(T), d_B, n, b, n); + + shift_ipiv_down(n, ipiv); + + cudaFree(d_A); + cudaFree(d_B); + cudaFree(d_I); + cudaFree(d_info); + + return info; +} + +template +inline int lu_solve(cusolverDnHandle_t solverHandle, int n, int nrhs, T a[], T b[], GETRF getrf, GETRS getrs, GETRFBSIZE getrfbsize) +{ + int info = 0; + + int* d_I = NULL; + cudaMalloc((void**)&d_I, n*sizeof(int)); + + T* d_A = NULL; + cudaMalloc((void**)&d_A, n*n*sizeof(T)); + cublasSetMatrix(n, n, sizeof(T), a, n, d_A, n); + + T* work = NULL; + int lwork = 0; + getrfbsize(solverHandle, n, n, a, n, &lwork); + cudaMalloc((void**)&work, sizeof(T)*lwork); + + int* d_info = NULL; + cudaMalloc((void**)&d_info, sizeof(int)); + + getrf(solverHandle, n, n, d_A, n, work, d_I, d_info); + cudaMemcpy(&info, d_info, sizeof(int), cudaMemcpyDeviceToHost); + + cudaFree(work); + + if (info != 0) + { + cudaFree(d_I); + cudaFree(d_A); + cudaFree(d_info); + return info; + } + + T* d_B = NULL; + cudaMalloc((void**)&d_B, n*nrhs*sizeof(T)); + cublasSetMatrix(n, nrhs, sizeof(T), b, n, d_B, n); + + getrs(solverHandle, CUBLAS_OP_N, n, nrhs, d_A, n, d_I, d_B, n, d_info); + cudaMemcpy(&info, d_info, 1, cudaMemcpyDeviceToHost); + + cublasGetMatrix(n, nrhs, sizeof(T), d_B, n, b, n); + + cudaFree(d_A); + cudaFree(d_B); + cudaFree(d_I); + cudaFree(d_info); + + return info; +} + + +template +inline int cholesky_factor(cusolverDnHandle_t solverHandle, int n, T a[], POTRF potrf, POTRFBSIZE potrfbsize) +{ + int info = 0; + + T* d_A = NULL; + cudaMalloc((void**)&d_A, n*n*sizeof(T)); + cublasSetMatrix(n, n, sizeof(T), a, n, d_A, n); + + T* work = NULL; + int lWork = 0; + potrfbsize(solverHandle, CUBLAS_FILL_MODE_LOWER, n, d_A, n, &lWork); + cudaMalloc((void**)&work, sizeof(T)*lWork); + + int* d_info = NULL; + cudaMalloc((void**)&d_info, sizeof(int)); + + potrf(solverHandle, CUBLAS_FILL_MODE_LOWER, n, d_A, n, work, lWork, d_info); + cudaMemcpy(&info, d_info, sizeof(int), cudaMemcpyDeviceToHost); + + cublasGetMatrix(n, n, sizeof(T), d_A, n, a, n); + + T zero = T(); + + for (int i = 0; i < n; ++i) + { + int index = i * n; + + for (int j = 0; j < n && i > j; ++j) + { + a[index + j] = zero; + } + } + + cudaFree(d_A); + cudaFree(d_info); + cudaFree(work); + + return info; +} + +template +inline int cholesky_solve(cusolverDnHandle_t solverHandle, int n, int nrhs, T a[], T b[], POTRF potrf, POTRS potrs, POTRFBSIZE potrfbsize) +{ + int info = 0; + + T* d_A = NULL; + cudaMalloc((void**)&d_A, n*n*sizeof(T)); + cublasSetMatrix(n, n, sizeof(T), a, n, d_A, n); + + T* work = NULL; + int lWork = 0; + potrfbsize(solverHandle, CUBLAS_FILL_MODE_LOWER, n, d_A, n, &lWork); + cudaMalloc((void**)&work, sizeof(T)*lWork); + + int* d_info = NULL; + cudaMalloc((void**)&d_info, sizeof(int)); + + potrf(solverHandle, CUBLAS_FILL_MODE_LOWER, n, d_A, n, work, lWork, d_info); + cudaMemcpy(&info, d_info, sizeof(int), cudaMemcpyDeviceToHost); + + cudaFree(work); + + if (info != 0) + { + cudaFree(d_A); + cudaFree(d_info); + return info; + } + + T* d_B = NULL; + cudaMalloc((void**)&d_B, n*nrhs*sizeof(T)); + cublasSetMatrix(n, nrhs, sizeof(T), b, n, d_B, n); + + potrs(solverHandle, CUBLAS_FILL_MODE_LOWER, n, nrhs, d_A, n, d_B, n, d_info); + cudaMemcpy(&info, d_info, sizeof(int), cudaMemcpyDeviceToHost); + + cublasGetMatrix(n, nrhs, sizeof(T), d_B, n, b, n); + + cudaFree(d_A); + cudaFree(d_B); + cudaFree(d_info); + + return info; +} + +template +inline int cholesky_solve_factored(cusolverDnHandle_t solverHandle, int n, int nrhs, T a[], T b[], POTRS potrs) +{ + int info = 0; + + T* d_A = NULL; + cudaMalloc((void**)&d_A, n*n*sizeof(T)); + cublasSetMatrix(n, n, sizeof(T), a, n, d_A, n); + + T* d_B = NULL; + cudaMalloc((void**)&d_B, n*nrhs*sizeof(T)); + cublasSetMatrix(n, nrhs, sizeof(T), b, n, d_B, n); + + int* d_info = NULL; + cudaMalloc((void**)&d_info, sizeof(int)); + + potrs(solverHandle, CUBLAS_FILL_MODE_LOWER, n, nrhs, d_A, n, d_B, n, d_info); + cudaMemcpy(&info, d_info, sizeof(int), cudaMemcpyDeviceToHost); + + cublasGetMatrix(n, nrhs, sizeof(T), d_B, n, b, n); + + cudaFree(d_A); + cudaFree(d_B); + cudaFree(d_info); + + return info; +} + +//template +//inline int qr_factor(int m, int n, T r[], T tau[], T q[], T work[], int len, GEQRF geqrf, ORGQR orgqr) +//{ +// int info = 0; +// geqrf(&m, &n, r, &m, tau, work, &len, &info); +// +// for (int i = 0; i < m; ++i) +// { +// for (int j = 0; j < m && j < n; ++j) +// { +// if (i > j) +// { +// q[j * m + i] = r[j * m + i]; +// } +// } +// } +// +// //compute the q elements explicitly +// if (m <= n) +// { +// orgqr(&m, &m, &m, q, &m, tau, work, &len, &info); +// } +// else +// { +// orgqr(&m, &m, &n, q, &m, tau, work, &len, &info); +// } +// +// return info; +//} +// +//template +//inline int qr_thin_factor(int m, int n, T q[], T tau[], T r[], T work[], int len, GEQRF geqrf, ORGQR orgqr) +//{ +// int info = 0; +// geqrf(&m, &n, q, &m, tau, work, &len, &info); +// +// for (int i = 0; i < n; ++i) +// { +// for (int j = 0; j < n; ++j) +// { +// if (i <= j) +// { +// r[j * n + i] = q[j * m + i]; +// } +// } +// } +// +// orgqr(&m, &n, &n, q, &m, tau, work, &len, &info); +// return info; +//} +// +//template +//inline int qr_solve(int m, int n, int bn, T a[], T b[], T x[], T work[], int len, GELS gels) +//{ +// T* clone_a = Clone(m, n, a); +// T* clone_b = Clone(m, bn, b); +// char N = 'N'; +// int info = 0; +// gels(&N, &m, &n, &bn, clone_a, &m, clone_b, &m, work, &len, &info); +// copyBtoX(m, n, bn, clone_b, x); +// delete[] clone_a; +// delete[] clone_b; +// return info; +//} + +//template +//inline int qr_solve_factored(cusolverDnHandle_t solverHandle, cublasHandle_t blasHandle, int m, int n, int bn, T r[], T b[], T tau[], T x[], T work[], int len, ORMQR ormqr, TRSM trsm) +//{ +// T* clone_b = Clone(m, bn, b); +// char side = 'L'; +// char tran = 'T'; +// int info = 0; +// ormqr(solverHandle, &side, &tran, &m, &bn, &n, r, &m, tau, clone_b, &m, work, &len, &info); +// trsm(blasHandle, CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, n, bn, 1.0, r, m, clone_b, m); +// +// copyBtoX(m, n, bn, clone_b, x); +// delete[] clone_b; +// return info; +//} + +//template +//inline int complex_qr_solve_factored(int m, int n, int bn, T r[], T b[], T tau[], T x[], T work[], int len, UNMQR unmqr, TRSM trsm) +//{ +// T* clone_b = Clone(m, bn, b); +// char side = 'L'; +// char tran = 'C'; +// int info = 0; +// unmqr(&side, &tran, &m, &bn, &n, r, &m, tau, clone_b, &m, work, &len, &info); +// T one = 1.0f; +// trsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, n, bn, &one, r, m, clone_b, m); +// copyBtoX(m, n, bn, clone_b, x); +// delete[] clone_b; +// return info; +//} + +template +inline int svd_factor(cusolverDnHandle_t solverHandle, bool compute_vectors, int m, int n, T a[], T s[], T u[], T v[], GESVD gesvd, GESVDBSIZE gesvdbsize) +{ + int info = 0; + int dim_s = std::min(m, n); + + T* d_A = NULL; + cudaMalloc((void**)&d_A, m*n*sizeof(T)); + cublasSetMatrix(m, n, sizeof(T), a, m, d_A, m); + + T* d_S = NULL; + cudaMalloc((void**)&d_S, dim_s*sizeof(T)); + + T* d_U = NULL; + cudaMalloc((void**)&d_U, m*m*sizeof(T)); + + T* d_V = NULL; + cudaMalloc((void**)&d_V, n*n*sizeof(T)); + + T* work = NULL; + int lWork = 0; + gesvdbsize(solverHandle, m, n, &lWork); + cudaMalloc((void**)&work, lWork*sizeof(T)); + + T* rwork = NULL; + cudaMalloc((void**)&rwork, 5 * dim_s * sizeof(T)); + + int* d_info = NULL; + cudaMalloc((void**)&d_info, sizeof(int)); + + char job = compute_vectors ? 'A' : 'N'; + gesvd(solverHandle, job, job, m, n, d_A, m, d_S, d_U, m, d_V, n, work, lWork, rwork, d_info); + cudaMemcpy(&info, d_info, sizeof(int), cudaMemcpyDeviceToHost); + + cublasGetVector(dim_s, sizeof(T), d_S, 1, s, 1); + cublasGetMatrix(m, m, sizeof(T), d_U, m, u, m); + cublasGetMatrix(n, n, sizeof(T), d_V, n, v, n); + + cudaFree(d_A); + cudaFree(d_S); + cudaFree(d_U); + cudaFree(d_V); + cudaFree(work); + cudaFree(rwork); + cudaFree(d_info); + + return info; +} + +template +inline int complex_svd_factor(cusolverDnHandle_t solverHandle, bool compute_vectors, int m, int n, T a[], T s[], T u[], T v[], GESVD gesvd, GESVDBSIZE gesvdbsize) +{ + int info = 0; + int dim_s = std::min(m, n); + + T* d_A = NULL; + cudaMalloc((void**)&d_A, m*n*sizeof(T)); + cublasSetMatrix(m, n, sizeof(T), a, m, d_A, m); + + R* s_local = new R[dim_s]; + R* d_S = NULL; + cudaMalloc((void**)&d_S, dim_s*sizeof(R)); + + T* d_U = NULL; + cudaMalloc((void**)&d_U, m*m*sizeof(T)); + + T* d_V = NULL; + cudaMalloc((void**)&d_V, n*m*sizeof(T)); + + T* work = NULL; + int lWork = 0; + gesvdbsize(solverHandle, m, n, &lWork); + cudaMalloc((void**)&work, lWork*sizeof(T)); + + R* rwork = NULL; + cudaMalloc((void**)&rwork, 5 * dim_s * sizeof(R)); + + int* d_info = NULL; + cudaMalloc((void**)&d_info, sizeof(int)); + + char job = compute_vectors ? 'A' : 'N'; + gesvd(solverHandle, job, job, m, n, d_A, m, d_S, d_U, m, d_V, n, work, lWork, rwork, d_info); + cudaMemcpy(&info, d_info, sizeof(int), cudaMemcpyDeviceToHost); + + cublasGetVector(dim_s, sizeof(R), d_S, 1, s_local, 1); + cublasGetMatrix(m, m, sizeof(T), d_U, m, u, m); + cublasGetMatrix(n, n, sizeof(T), d_V, n, v, n); + + for (int index = 0; index < dim_s; ++index) + { + s[index].x = s_local[index]; + } + + delete[] s_local; + cudaFree(d_A); + cudaFree(d_S); + cudaFree(d_U); + cudaFree(d_V); + cudaFree(work); + cudaFree(rwork); + cudaFree(d_info); + + return info; +} + +//template +//inline int eigen_factor(int n, T a[], T vectors[], R values[], T d[], GEES gees, TREVC trevc) +//{ +// T* clone_a = Clone(n, n, a); +// T* wr = new T[n]; +// T* wi = new T[n]; +// +// int sdim; +// int info = gees(LAPACK_COL_MAJOR, 'V', 'N', nullptr, n, clone_a, n, &sdim, wr, wi, vectors, n); +// if (info != 0) +// { +// delete[] clone_a; +// delete[] wr; +// delete[] wi; +// return info; +// } +// +// int m; +// info = trevc(LAPACK_COL_MAJOR, 'R', 'B', nullptr, n, clone_a, n, nullptr, n, vectors, n, n, &m); +// if (info != 0) +// { +// delete[] clone_a; +// delete[] wr; +// delete[] wi; +// return info; +// } +// +// for (int index = 0; index < n; ++index) +// { +// values[index] = R(wr[index], wi[index]); +// } +// +// for (int i = 0; i < n; ++i) +// { +// int in = i * n; +// d[in + i] = wr[i]; +// +// if (wi[i] > 0) +// { +// d[in + n + i] = wi[i]; +// } +// else if (wi[i] < 0) +// { +// d[in - n + i] = wi[i]; +// } +// } +// +// delete[] clone_a; +// delete[] wr; +// delete[] wi; +// return info; +//} +// +//template +//inline int eigen_complex_factor(int n, T a[], T vectors[], cuDoubleComplex values[], T d[], GEES gees, TREVC trevc) +//{ +// T* clone_a = Clone(n, n, a); +// T* w = new T[n]; +// +// int sdim; +// int info = gees(LAPACK_COL_MAJOR, 'V', 'N', nullptr, n, clone_a, n, &sdim, w, vectors, n); +// if (info != 0) +// { +// delete[] clone_a; +// delete[] w; +// return info; +// } +// +// int m; +// info = trevc(LAPACK_COL_MAJOR, 'R', 'B', nullptr, n, clone_a, n, nullptr, n, vectors, n, n, &m); +// if (info != 0) +// { +// delete[] clone_a; +// delete[] w; +// return info; +// } +// +// for (int i = 0; i < n; ++i) +// { +// values[i] = w[i]; +// d[i * n + i] = w[i]; +// } +// +// delete[] clone_a; +// delete[] w; +// return info; +//} +// +//template +//inline int sym_eigen_factor(int n, T a[], T vectors[], cuDoubleComplex values[], T d[], SYEV syev) +//{ +// T* clone_a = Clone(n, n, a); +// R* w = new R[n]; +// +// int info = syev(LAPACK_COL_MAJOR, 'V', 'U', n, clone_a, n, w); +// if (info != 0) +// { +// delete[] clone_a; +// delete[] w; +// return info; +// } +// +// memcpy(vectors, clone_a, n*n*sizeof(T)); +// +// for (int index = 0; index < n; ++index) +// { +// values[index] = cuDoubleComplex(w[index]); +// } +// +// for (int j = 0; j < n; ++j) +// { +// int jn = j*n; +// +// for (int i = 0; i < n; ++i) +// { +// if (i == j) +// { +// d[jn + i] = w[i]; +// } +// } +// } +// +// delete[] clone_a; +// delete[] w; +// return info; +//} + +#define sgetrf cusolverDnSgetrf +#define dgetrf cusolverDnDgetrf +#define cgetrf cusolverDnCgetrf +#define zgetrf cusolverDnZgetrf +#define sgetrfbsize cusolverDnSgetrf_bufferSize +#define dgetrfbsize cusolverDnDgetrf_bufferSize +#define cgetrfbsize cusolverDnCgetrf_bufferSize +#define zgetrfbsize cusolverDnZgetrf_bufferSize + +#define sgetrs cusolverDnSgetrs +#define dgetrs cusolverDnDgetrs +#define cgetrs cusolverDnCgetrs +#define zgetrs cusolverDnZgetrs + +#define spotrf cusolverDnSpotrf +#define dpotrf cusolverDnDpotrf +#define cpotrf cusolverDnCpotrf +#define zpotrf cusolverDnZpotrf +#define spotrfbsize cusolverDnSpotrf_bufferSize +#define dpotrfbsize cusolverDnDpotrf_bufferSize +#define cpotrfbsize cusolverDnCpotrf_bufferSize +#define zpotrfbsize cusolverDnZpotrf_bufferSize + +#define spotrs cusolverDnSpotrs +#define dpotrs cusolverDnDpotrs +#define cpotrs cusolverDnCpotrs +#define zpotrs cusolverDnZpotrs + +#define sgeqrf cusolverDnSgeqrf +#define dgeqrf cusolverDnDgeqrf +#define cgeqrf cusolverDnCgeqrf +#define zgeqrf cusolverDnZgeqrf + +#define sormqr cusolverDnSormqr +#define dormqr cusolverDnDormqr + +#define sgesvd cusolverDnSgesvd +#define dgesvd cusolverDnDgesvd +#define cgesvd cusolverDnCgesvd +#define zgesvd cusolverDnZgesvd +#define sgesvdbsize cusolverDnSgesvd_bufferSize +#define dgesvdbsize cusolverDnDgesvd_bufferSize +#define cgesvdbsize cusolverDnCgesvd_bufferSize +#define zgesvdbsize cusolverDnZgesvd_bufferSize + +#define sgetribatched cublasSgetriBatched +#define dgetribatched cublasDgetriBatched +#define cgetribatched cublasCgetriBatched +#define zgetribatched cublasZgetriBatched + +extern "C" { + + DLLEXPORT int s_lu_factor(cusolverDnHandle_t solverHandle, int m, float a[], int ipiv[]) + { + return lu_factor(solverHandle, m, a, ipiv, sgetrf, sgetrfbsize); + } + + DLLEXPORT int d_lu_factor(cusolverDnHandle_t solverHandle, int m, double a[], int ipiv[]) + { + return lu_factor(solverHandle, m, a, ipiv, dgetrf, dgetrfbsize); + } + + DLLEXPORT int c_lu_factor(cusolverDnHandle_t solverHandle, int m, cuComplex a[], int ipiv[]) + { + return lu_factor(solverHandle, m, a, ipiv, cgetrf, cgetrfbsize); + } + + DLLEXPORT int z_lu_factor(cusolverDnHandle_t solverHandle, int m, cuDoubleComplex a[], int ipiv[]) + { + return lu_factor(solverHandle, m, a, ipiv, zgetrf, zgetrfbsize); + } + + DLLEXPORT int s_lu_inverse(cusolverDnHandle_t solverHandle, cublasHandle_t blasHandle, int n, float a[]) + { + return lu_inverse(solverHandle, blasHandle, n, a, sgetrf, sgetribatched, sgetrfbsize); + } + + DLLEXPORT int d_lu_inverse(cusolverDnHandle_t solverHandle, cublasHandle_t blasHandle, int n, double a[]) + { + return lu_inverse(solverHandle, blasHandle, n, a, dgetrf, dgetribatched, dgetrfbsize); + } + + DLLEXPORT int c_lu_inverse(cusolverDnHandle_t solverHandle, cublasHandle_t blasHandle, int n, cuComplex a[]) + { + return lu_inverse(solverHandle, blasHandle, n, a, cgetrf, cgetribatched, cgetrfbsize); + } + + DLLEXPORT int z_lu_inverse(cusolverDnHandle_t solverHandle, cublasHandle_t blasHandle, int n, cuDoubleComplex a[]) + { + return lu_inverse(solverHandle, blasHandle, n, a, zgetrf, zgetribatched, zgetrfbsize); + } + + DLLEXPORT int s_lu_inverse_factored(cublasHandle_t blasHandle, int n, float a[], int ipiv[]) + { + return lu_inverse_factored(blasHandle, n, a, ipiv, sgetribatched); + } + + DLLEXPORT int d_lu_inverse_factored(cublasHandle_t blasHandle, int n, double a[], int ipiv[]) + { + return lu_inverse_factored(blasHandle, n, a, ipiv, dgetribatched); + } + + DLLEXPORT int c_lu_inverse_factored(cublasHandle_t blasHandle, int n, cuComplex a[], int ipiv[]) + { + return lu_inverse_factored(blasHandle, n, a, ipiv, cgetribatched); + } + + DLLEXPORT int z_lu_inverse_factored(cublasHandle_t blasHandle, int n, cuDoubleComplex a[], int ipiv[]) + { + return lu_inverse_factored(blasHandle, n, a, ipiv, zgetribatched); + } + + DLLEXPORT int s_lu_solve_factored(cusolverDnHandle_t solverHandle, int n, int nrhs, float a[], int ipiv[], float b[]) + { + return lu_solve_factored(solverHandle, n, nrhs, a, ipiv, b, sgetrs); + } + + DLLEXPORT int d_lu_solve_factored(cusolverDnHandle_t solverHandle, int n, int nrhs, double a[], int ipiv[], double b[]) + { + return lu_solve_factored(solverHandle, n, nrhs, a, ipiv, b, dgetrs); + } + + DLLEXPORT int c_lu_solve_factored(cusolverDnHandle_t solverHandle, int n, int nrhs, cuComplex a[], int ipiv[], cuComplex b[]) + { + return lu_solve_factored(solverHandle, n, nrhs, a, ipiv, b, cgetrs); + } + + DLLEXPORT int z_lu_solve_factored(cusolverDnHandle_t solverHandle, int n, int nrhs, cuDoubleComplex a[], int ipiv[], cuDoubleComplex b[]) + { + return lu_solve_factored(solverHandle, n, nrhs, a, ipiv, b, zgetrs); + } + + DLLEXPORT int s_lu_solve(cusolverDnHandle_t solverHandle, int n, int nrhs, float a[], float b[]) + { + return lu_solve(solverHandle, n, nrhs, a, b, sgetrf, sgetrs, sgetrfbsize); + } + + DLLEXPORT int d_lu_solve(cusolverDnHandle_t solverHandle, int n, int nrhs, double a[], double b[]) + { + return lu_solve(solverHandle, n, nrhs, a, b, dgetrf, dgetrs, dgetrfbsize); + } + + DLLEXPORT int c_lu_solve(cusolverDnHandle_t solverHandle, int n, int nrhs, cuComplex a[], cuComplex b[]) + { + return lu_solve(solverHandle, n, nrhs, a, b, cgetrf, cgetrs, cgetrfbsize); + } + + DLLEXPORT int z_lu_solve(cusolverDnHandle_t solverHandle, int n, int nrhs, cuDoubleComplex a[], cuDoubleComplex b[]) + { + return lu_solve(solverHandle, n, nrhs, a, b, zgetrf, zgetrs, zgetrfbsize); + } + + DLLEXPORT int s_cholesky_factor(cusolverDnHandle_t solverHandle, int n, float a[]) + { + return cholesky_factor(solverHandle, n, a, spotrf, spotrfbsize); + } + + DLLEXPORT int d_cholesky_factor(cusolverDnHandle_t solverHandle, int n, double a[]) + { + return cholesky_factor(solverHandle, n, a, dpotrf, dpotrfbsize); + } + + DLLEXPORT int c_cholesky_factor(cusolverDnHandle_t solverHandle, int n, cuComplex a[]) + { + return cholesky_factor(solverHandle, n, a, cpotrf, cpotrfbsize); + } + + DLLEXPORT int z_cholesky_factor(cusolverDnHandle_t solverHandle, int n, cuDoubleComplex a[]) + { + return cholesky_factor(solverHandle, n, a, zpotrf, zpotrfbsize); + } + + DLLEXPORT int s_cholesky_solve(cusolverDnHandle_t solverHandle, int n, int nrhs, float a[], float b[]) + { + return cholesky_solve(solverHandle, n, nrhs, a, b, spotrf, spotrs, spotrfbsize); + } + + DLLEXPORT int d_cholesky_solve(cusolverDnHandle_t solverHandle, int n, int nrhs, double a[], double b[]) + { + return cholesky_solve(solverHandle, n, nrhs, a, b, dpotrf, dpotrs, dpotrfbsize); + } + + DLLEXPORT int c_cholesky_solve(cusolverDnHandle_t solverHandle, int n, int nrhs, cuComplex a[], cuComplex b[]) + { + return cholesky_solve(solverHandle, n, nrhs, a, b, cpotrf, cpotrs, cpotrfbsize); + } + + DLLEXPORT int z_cholesky_solve(cusolverDnHandle_t solverHandle, int n, int nrhs, cuDoubleComplex a[], cuDoubleComplex b[]) + { + return cholesky_solve(solverHandle, n, nrhs, a, b, zpotrf, zpotrs, zpotrfbsize); + } + + DLLEXPORT int s_cholesky_solve_factored(cusolverDnHandle_t solverHandle, int n, int nrhs, float a[], float b[]) + { + return cholesky_solve_factored(solverHandle, n, nrhs, a, b, spotrs); + } + + DLLEXPORT int d_cholesky_solve_factored(cusolverDnHandle_t solverHandle, int n, int nrhs, double a[], double b[]) + { + return cholesky_solve_factored(solverHandle, n, nrhs, a, b, dpotrs); + } + + DLLEXPORT int c_cholesky_solve_factored(cusolverDnHandle_t solverHandle, int n, int nrhs, cuComplex a[], cuComplex b[]) + { + return cholesky_solve_factored(solverHandle, n, nrhs, a, b, cpotrs); + } + + DLLEXPORT int z_cholesky_solve_factored(cusolverDnHandle_t solverHandle, int n, int nrhs, cuDoubleComplex a[], cuDoubleComplex b[]) + { + return cholesky_solve_factored(solverHandle, n, nrhs, a, b, zpotrs); + } + + // MJ: I am fairly certain that it would be straightforward to implement ?orgqr and ?gels but I'm focusing on getting the low-hanging fruit working first + /*DLLEXPORT int s_qr_factor(int m, int n, float r[], float tau[], float q[], float work[], int len) + { + return qr_factor(m, n, r, tau, q, work, len, sgeqrf, sorgqr); + } + + DLLEXPORT int s_qr_thin_factor(int m, int n, float q[], float tau[], float r[], float work[], int len) + { + return qr_thin_factor(m, n, q, tau, r, work, len, sgeqrf, sorgqr); + } + + DLLEXPORT int d_qr_factor(int m, int n, double r[], double tau[], double q[], double work[], int len) + { + return qr_factor(m, n, r, tau, q, work, len, dgeqrf, dorgqr); + } + + DLLEXPORT int d_qr_thin_factor(int m, int n, double q[], double tau[], double r[], double work[], int len) + { + return qr_thin_factor(m, n, q, tau, r, work, len, dgeqrf, dorgqr); + } + + DLLEXPORT int c_qr_factor(int m, int n, cuComplex r[], cuComplex tau[], cuComplex q[], cuComplex work[], int len) + { + return qr_factor(m, n, r, tau, q, work, len, cgeqrf, cungqr); + } + + DLLEXPORT int c_qr_thin_factor(int m, int n, cuComplex q[], cuComplex tau[], cuComplex r[], cuComplex work[], int len) + { + return qr_thin_factor(m, n, q, tau, r, work, len, cgeqrf, cungqr); + } + + DLLEXPORT int z_qr_factor(int m, int n, cuDoubleComplex r[], cuDoubleComplex tau[], cuDoubleComplex q[], cuDoubleComplex work[], int len) + { + return qr_factor(m, n, r, tau, q, work, len, zgeqrf, zungqr); + } + + DLLEXPORT int z_qr_thin_factor(int m, int n, cuDoubleComplex q[], cuDoubleComplex tau[], cuDoubleComplex r[], cuDoubleComplex work[], int len) + { + return qr_thin_factor(m, n, q, tau, r, work, len, zgeqrf, zungqr); + } + + DLLEXPORT int s_qr_solve(int m, int n, int bn, float a[], float b[], float x[], float work[], int len) + { + return qr_solve(m, n, bn, a, b, x, work, len, sgels); + } + + DLLEXPORT int d_qr_solve(int m, int n, int bn, double a[], double b[], double x[], double work[], int len) + { + return qr_solve(m, n, bn, a, b, x, work, len, dgels); + } + + DLLEXPORT int c_qr_solve(int m, int n, int bn, cuComplex a[], cuComplex b[], cuComplex x[], cuComplex work[], int len) + { + return qr_solve(m, n, bn, a, b, x, work, len, cgels); + } + + DLLEXPORT int z_qr_solve(int m, int n, int bn, cuDoubleComplex a[], cuDoubleComplex b[], cuDoubleComplex x[], cuDoubleComplex work[], int len) + { + return qr_solve(m, n, bn, a, b, x, work, len, zgels); + }*/ + + //DLLEXPORT int s_qr_solve_factored(cusolverDnHandle_t solverHandle, cublasHandle_t blasHandle, int m, int n, int bn, float r[], float b[], float tau[], float x[], float work[], int len) + //{ + // return qr_solve_factored(solverHandle, blasHandle, m, n, bn, r, b, tau, x, work, len, sormqr, cublasStrsm); + //} + + //DLLEXPORT int d_qr_solve_factored(cusolverDnHandle_t solverHandle, cublasHandle_t blasHandle, int m, int n, int bn, double r[], double b[], double tau[], double x[], double work[], int len) + //{ + // return qr_solve_factored(solverHandle, blasHandle, m, n, bn, r, b, tau, x, work, len, dormqr, cublasDtrsm); + //} + + //DLLEXPORT int c_qr_solve_factored(int m, int n, int bn, cuComplex r[], cuComplex b[], cuComplex tau[], cuComplex x[], cuComplex work[], int len) + //{ + // return complex_qr_solve_factored(m, n, bn, r, b, tau, x, work, len, cunmqr, cublasCtrsm); + //} + + //DLLEXPORT int z_qr_solve_factored(int m, int n, int bn, cuDoubleComplex r[], cuDoubleComplex b[], cuDoubleComplex tau[], cuDoubleComplex x[], cuDoubleComplex work[], int len) + //{ + // return complex_qr_solve_factored(m, n, bn, r, b, tau, x, work, len, zunmqr, cublasZtrsm); + //} + + DLLEXPORT int s_svd_factor(cusolverDnHandle_t solverHandle, bool compute_vectors, int m, int n, float a[], float s[], float u[], float v[]) + { + return svd_factor(solverHandle, compute_vectors, m, n, a, s, u, v, sgesvd, sgesvdbsize); + } + + DLLEXPORT int d_svd_factor(cusolverDnHandle_t solverHandle, bool compute_vectors, int m, int n, double a[], double s[], double u[], double v[]) + { + return svd_factor(solverHandle, compute_vectors, m, n, a, s, u, v,dgesvd, dgesvdbsize); + } + + DLLEXPORT int c_svd_factor(cusolverDnHandle_t solverHandle, bool compute_vectors, int m, int n, cuComplex a[], cuComplex s[], cuComplex u[], cuComplex v[]) + { + return complex_svd_factor(solverHandle, compute_vectors, m, n, a, s, u, v, cgesvd, cgesvdbsize); + } + + DLLEXPORT int z_svd_factor(cusolverDnHandle_t solverHandle, bool compute_vectors, int m, int n, cuDoubleComplex a[], cuDoubleComplex s[], cuDoubleComplex u[], cuDoubleComplex v[]) + { + return complex_svd_factor(solverHandle, compute_vectors, m, n, a, s, u, v, zgesvd, zgesvdbsize); + } + + /*DLLEXPORT int s_eigen(bool isSymmetric, int n, float a[], float vectors[], cuDoubleComplex values[], float d[]) + { + if (isSymmetric) + { + return sym_eigen_factor(n, a, vectors, values, d, LAPACKE_ssyev); + } + else + { + return eigen_factor(n, a, vectors, values, d, LAPACKE_sgees, LAPACKE_strevc); + } + } + + DLLEXPORT int d_eigen(bool isSymmetric, int n, double a[], double vectors[], cuDoubleComplex values[], double d[]) + { + if (isSymmetric) + { + return sym_eigen_factor(n, a, vectors, values, d, LAPACKE_dsyev); + } + else + { + return eigen_factor(n, a, vectors, values, d, LAPACKE_dgees, LAPACKE_dtrevc); + } + } + + DLLEXPORT int c_eigen(bool isSymmetric, int n, cuComplex a[], cuComplex vectors[], cuDoubleComplex values[], cuComplex d[]) + { + if (isSymmetric) + { + return sym_eigen_factor(n, a, vectors, values, d, LAPACKE_cheev); + } + else + { + return eigen_complex_factor(n, a, vectors, values, d, LAPACKE_cgees, LAPACKE_ctrevc); + } + } + + DLLEXPORT int z_eigen(bool isSymmetric, int n, cuDoubleComplex a[], cuDoubleComplex vectors[], cuDoubleComplex values[], cuDoubleComplex d[]) + { + if (isSymmetric) + { + return sym_eigen_factor(n, a, vectors, values, d, LAPACKE_zheev); + } + else + { + return eigen_complex_factor(n, a, vectors, values, d, LAPACKE_zgees, LAPACKE_ztrevc); + } + }*/ +} diff --git a/src/NativeProviders/CUDA/memory.c b/src/NativeProviders/CUDA/memory.c new file mode 100644 index 00000000..e69de29b diff --git a/src/NativeProviders/CUDA/resource.h b/src/NativeProviders/CUDA/resource.h new file mode 100644 index 00000000..27e2900c --- /dev/null +++ b/src/NativeProviders/CUDA/resource.h @@ -0,0 +1,14 @@ +//{{NO_DEPENDENCIES}} +// Microsoft Visual C++ generated include file. +// Used by resource.rc + +// Next default values for new objects +// +#ifdef APSTUDIO_INVOKED +#ifndef APSTUDIO_READONLY_SYMBOLS +#define _APS_NEXT_RESOURCE_VALUE 101 +#define _APS_NEXT_COMMAND_VALUE 40001 +#define _APS_NEXT_CONTROL_VALUE 1001 +#define _APS_NEXT_SYMED_VALUE 101 +#endif +#endif diff --git a/src/NativeProviders/CUDA/resource.rc b/src/NativeProviders/CUDA/resource.rc new file mode 100644 index 00000000..30aa52c5 --- /dev/null +++ b/src/NativeProviders/CUDA/resource.rc @@ -0,0 +1,101 @@ +// Microsoft Visual C++ generated resource script. +// +#include "resource.h" + +#define APSTUDIO_READONLY_SYMBOLS +///////////////////////////////////////////////////////////////////////////// +// +// Generated from the TEXTINCLUDE 2 resource. +// +#include "windows.h" + +///////////////////////////////////////////////////////////////////////////// +#undef APSTUDIO_READONLY_SYMBOLS + +///////////////////////////////////////////////////////////////////////////// +// English (United States) resources + +#if !defined(AFX_RESOURCE_DLL) || defined(AFX_TARG_ENU) +LANGUAGE LANG_ENGLISH, SUBLANG_ENGLISH_US +#pragma code_page(1252) + +#ifdef APSTUDIO_INVOKED +///////////////////////////////////////////////////////////////////////////// +// +// TEXTINCLUDE +// + +1 TEXTINCLUDE +BEGIN + "resource.h\0" +END + +2 TEXTINCLUDE +BEGIN + "#include ""windows.h""\r\n" + "\0" +END + +3 TEXTINCLUDE +BEGIN + "\r\n" + "\0" +END + +#endif // APSTUDIO_INVOKED + + +///////////////////////////////////////////////////////////////////////////// +// +// Version +// + +VS_VERSION_INFO VERSIONINFO + FILEVERSION 0,1,0,0 + PRODUCTVERSION 0,1,0,0 + FILEFLAGSMASK 0x17L +#ifdef _DEBUG + FILEFLAGS 0x1L +#else + FILEFLAGS 0x0L +#endif + FILEOS 0x4L + FILETYPE 0x2L + FILESUBTYPE 0x0L +BEGIN + BLOCK "StringFileInfo" + BEGIN + BLOCK "040904b0" + BEGIN + VALUE "Comments", "http://numerics.mathdotnet.com/" + VALUE "CompanyName", "Math.NET" + VALUE "FileDescription", "MathNET Numerics CUDA Native Provider" + VALUE "FileVersion", "0.1.0.0" + VALUE "InternalName", "Math.NET" + VALUE "LegalCopyright", "Copyright (C) Math.NET 2009-2015" + VALUE "OriginalFilename", "MathNet.Numerics.CUDA" + VALUE "ProductName", "Math.NET Numerics" + VALUE "ProductVersion", "0.1.0.0" + END + END + BLOCK "VarFileInfo" + BEGIN + VALUE "Translation", 0x409, 1200 + END +END + +#endif // English (United States) resources +///////////////////////////////////////////////////////////////////////////// + + + +#ifndef APSTUDIO_INVOKED +///////////////////////////////////////////////////////////////////////////// +// +// Generated from the TEXTINCLUDE 3 resource. +// + + +///////////////////////////////////////////////////////////////////////////// +#endif // not APSTUDIO_INVOKED + diff --git a/src/NativeProviders/MKL/resource.h b/src/NativeProviders/MKL/resource.h new file mode 100644 index 00000000..27e2900c --- /dev/null +++ b/src/NativeProviders/MKL/resource.h @@ -0,0 +1,14 @@ +//{{NO_DEPENDENCIES}} +// Microsoft Visual C++ generated include file. +// Used by resource.rc + +// Next default values for new objects +// +#ifdef APSTUDIO_INVOKED +#ifndef APSTUDIO_READONLY_SYMBOLS +#define _APS_NEXT_RESOURCE_VALUE 101 +#define _APS_NEXT_COMMAND_VALUE 40001 +#define _APS_NEXT_CONTROL_VALUE 1001 +#define _APS_NEXT_SYMED_VALUE 101 +#endif +#endif diff --git a/src/NativeProviders/MKL/resource.rc b/src/NativeProviders/MKL/resource.rc new file mode 100644 index 00000000..6604b1a7 --- /dev/null +++ b/src/NativeProviders/MKL/resource.rc @@ -0,0 +1,101 @@ +// Microsoft Visual C++ generated resource script. +// +#include "resource.h" + +#define APSTUDIO_READONLY_SYMBOLS +///////////////////////////////////////////////////////////////////////////// +// +// Generated from the TEXTINCLUDE 2 resource. +// +#include "windows.h" + +///////////////////////////////////////////////////////////////////////////// +#undef APSTUDIO_READONLY_SYMBOLS + +///////////////////////////////////////////////////////////////////////////// +// English (United States) resources + +#if !defined(AFX_RESOURCE_DLL) || defined(AFX_TARG_ENU) +LANGUAGE LANG_ENGLISH, SUBLANG_ENGLISH_US +#pragma code_page(1252) + +#ifdef APSTUDIO_INVOKED +///////////////////////////////////////////////////////////////////////////// +// +// TEXTINCLUDE +// + +1 TEXTINCLUDE +BEGIN + "resource.h\0" +END + +2 TEXTINCLUDE +BEGIN + "#include ""windows.h""\r\n" + "\0" +END + +3 TEXTINCLUDE +BEGIN + "\r\n" + "\0" +END + +#endif // APSTUDIO_INVOKED + + +///////////////////////////////////////////////////////////////////////////// +// +// Version +// + +VS_VERSION_INFO VERSIONINFO + FILEVERSION 1,7,0,0 + PRODUCTVERSION 1,7,0,0 + FILEFLAGSMASK 0x17L +#ifdef _DEBUG + FILEFLAGS 0x1L +#else + FILEFLAGS 0x0L +#endif + FILEOS 0x4L + FILETYPE 0x2L + FILESUBTYPE 0x0L +BEGIN + BLOCK "StringFileInfo" + BEGIN + BLOCK "040904b0" + BEGIN + VALUE "Comments", "http://numerics.mathdotnet.com/" + VALUE "CompanyName", "Math.NET" + VALUE "FileDescription", "MathNET Numerics MKL Native Provider" + VALUE "FileVersion", "1.7.0.0" + VALUE "InternalName", "Math.NET" + VALUE "LegalCopyright", "Copyright (C) Math.NET 2009-2015" + VALUE "OriginalFilename", "MathNet.Numerics.MKL" + VALUE "ProductName", "Math.NET Numerics" + VALUE "ProductVersion", "1.7.0.0" + END + END + BLOCK "VarFileInfo" + BEGIN + VALUE "Translation", 0x409, 1200 + END +END + +#endif // English (United States) resources +///////////////////////////////////////////////////////////////////////////// + + + +#ifndef APSTUDIO_INVOKED +///////////////////////////////////////////////////////////////////////////// +// +// Generated from the TEXTINCLUDE 3 resource. +// + + +///////////////////////////////////////////////////////////////////////////// +#endif // not APSTUDIO_INVOKED + diff --git a/src/NativeProviders/Windows/ATLAS/ATLASWrapper.vcxproj b/src/NativeProviders/Windows/ATLAS/ATLASWrapper.vcxproj index afdea811..55eeaefa 100644 --- a/src/NativeProviders/Windows/ATLAS/ATLASWrapper.vcxproj +++ b/src/NativeProviders/Windows/ATLAS/ATLASWrapper.vcxproj @@ -148,7 +148,7 @@ - + diff --git a/src/NativeProviders/Windows/ATLAS/ATLASWrapper.vcxproj.filters b/src/NativeProviders/Windows/ATLAS/ATLASWrapper.vcxproj.filters index d6464ae4..518bf6df 100644 --- a/src/NativeProviders/Windows/ATLAS/ATLASWrapper.vcxproj.filters +++ b/src/NativeProviders/Windows/ATLAS/ATLASWrapper.vcxproj.filters @@ -15,7 +15,7 @@ - + Resource Files diff --git a/src/NativeProviders/Windows/ATLASEx/ATLASWrapper.vcproj b/src/NativeProviders/Windows/ATLASEx/ATLASWrapper.vcproj index d38c82fd..64823410 100644 --- a/src/NativeProviders/Windows/ATLASEx/ATLASWrapper.vcproj +++ b/src/NativeProviders/Windows/ATLASEx/ATLASWrapper.vcproj @@ -206,7 +206,7 @@ UniqueIdentifier="{67DA6AB6-F800-4c08-8B7A-83BB121AAD01}" > diff --git a/src/NativeProviders/Windows/ATLASEx/ATLASWrapper.vcxproj b/src/NativeProviders/Windows/ATLASEx/ATLASWrapper.vcxproj index 11109832..5c811855 100644 --- a/src/NativeProviders/Windows/ATLASEx/ATLASWrapper.vcxproj +++ b/src/NativeProviders/Windows/ATLASEx/ATLASWrapper.vcxproj @@ -15,7 +15,7 @@ - + diff --git a/src/NativeProviders/Windows/ATLASEx/ATLASWrapper.vcxproj.filters b/src/NativeProviders/Windows/ATLASEx/ATLASWrapper.vcxproj.filters index 58436d58..1cfd00b1 100644 --- a/src/NativeProviders/Windows/ATLASEx/ATLASWrapper.vcxproj.filters +++ b/src/NativeProviders/Windows/ATLASEx/ATLASWrapper.vcxproj.filters @@ -23,7 +23,7 @@ - + Resource Files diff --git a/src/NativeProviders/Windows/CUDA/CUDAWrapper.vcxproj b/src/NativeProviders/Windows/CUDA/CUDAWrapper.vcxproj new file mode 100644 index 00000000..476349ff --- /dev/null +++ b/src/NativeProviders/Windows/CUDA/CUDAWrapper.vcxproj @@ -0,0 +1,174 @@ + + + + + Debug + Win32 + + + Debug + x64 + + + Release + Win32 + + + Release + x64 + + + + + + + + + + + + + {5A52B796-7F41-4C90-8DE2-F3F391C4482C} + CUDA + CUDA + + + + DynamicLibrary + true + v120 + MultiByte + + + DynamicLibrary + true + v120 + MultiByte + + + DynamicLibrary + false + v120 + true + MultiByte + + + DynamicLibrary + false + v120 + true + MultiByte + + + + + + + + + + + + + + + + + + + $(ProjectDir)..\..\..\..\out\CUDA\Windows\x86\ + $(Platform)\$(Configuration)\ + MathNet.Numerics.CUDA + + + $(Platform)\$(Configuration)\ + MathNet.Numerics.CUDA + $(ProjectDir)..\..\..\..\out\CUDA\Windows\x64\ + + + $(ProjectDir)..\..\..\..\out\CUDA\Windows\x86\ + $(Platform)\$(Configuration)\ + MathNet.Numerics.CUDA + + + $(Platform)\$(Configuration)\ + MathNet.Numerics.CUDA + $(ProjectDir)..\..\..\..\out\CUDA\Windows\x64\ + + + + Level3 + Disabled + true + $(CUDA_PATH)\include;$(ProjectDir)..\..\Common;$(ProjectDir)..\..\CUDA;%(AdditionalIncludeDirectories) + + + true + cublas.lib;cublas_device.lib;%(AdditionalDependencies) + $(CUDA_PATH)\lib\Win32;%(AdditionalLibraryDirectories) + + + + + Level3 + Disabled + $(CUDA_PATH)\include;$(ProjectDir)..\..\Common;$(ProjectDir)..\..\CUDA;%(AdditionalIncludeDirectories) + _WINDOWS;%(PreprocessorDefinitions) + MultiThreadedDebug + + + true + cudart.lib;cusolver.lib;cublas.lib;cublas_device.lib;%(AdditionalDependencies) + $(CUDA_PATH)\lib\x64;%(AdditionalLibraryDirectories) + + + copy "$(CUDA_PATH)\bin\cublas64_70.dll" $(OutputPath) +copy "$(CUDA_PATH)\bin\cusolver64_70.dll" $(OutputPath) +copy "$(CUDA_PATH)\bin\cudart64_70.dll" $(OutputPath) + + + + + Level3 + MaxSpeed + true + true + true + $(CUDA_PATH)\include;$(ProjectDir)..\..\Common;$(ProjectDir)..\..\CUDA;%(AdditionalIncludeDirectories) + + + true + true + true + cublas.lib;%(AdditionalDependencies) + $(CUDA_PATH)\lib\Win32;%(AdditionalLibraryDirectories) + + + + + Level3 + MaxSpeed + true + true + $(CUDA_PATH)\include;$(ProjectDir)..\..\Common;$(ProjectDir)..\..\CUDA;%(AdditionalIncludeDirectories) + _WINDOWS;%(PreprocessorDefinitions) + /Qvec-report:1 + MultiThreaded + + + true + true + true + cudart.lib;cusolver.lib;cublas.lib;%(AdditionalDependencies) + $(CUDA_PATH)\lib\x64;%(AdditionalLibraryDirectories) + + + copy "$(CUDA_PATH)\bin\cublas64_70.dll" $(OutputPath) +copy "$(CUDA_PATH)\bin\cusolver64_70.dll" $(OutputPath) +copy "$(CUDA_PATH)\bin\cudart64_70.dll" $(OutputPath) + + + + + + \ No newline at end of file diff --git a/src/NativeProviders/Windows/CUDA/CUDAWrapper.vcxproj.filters b/src/NativeProviders/Windows/CUDA/CUDAWrapper.vcxproj.filters new file mode 100644 index 00000000..c67ec930 --- /dev/null +++ b/src/NativeProviders/Windows/CUDA/CUDAWrapper.vcxproj.filters @@ -0,0 +1,36 @@ + + + + + {4FC737F1-C7A5-4376-A066-2A32D752A2FF} + cpp;c;cc;cxx;def;odl;idl;hpj;bat;asm;asmx + + + {93995380-89BD-4b04-88EB-625FBE52EBFB} + h;hh;hpp;hxx;hm;inl;inc;xsd + + + {67DA6AB6-F800-4c08-8B7A-83BB121AAD01} + rc;ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe;resx;tiff;tif;png;wav;mfcribbon-ms + + + + + Resource Files + + + + + Source Files + + + Source Files + + + Source Files + + + Source Files + + + \ No newline at end of file diff --git a/src/NativeProviders/Windows/MKL/MKLWrapper.vcproj b/src/NativeProviders/Windows/MKL/MKLWrapper.vcproj index 0ce61789..7b46f402 100644 --- a/src/NativeProviders/Windows/MKL/MKLWrapper.vcproj +++ b/src/NativeProviders/Windows/MKL/MKLWrapper.vcproj @@ -368,7 +368,7 @@ UniqueIdentifier="{67DA6AB6-F800-4c08-8B7A-83BB121AAD01}" > diff --git a/src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj b/src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj index 26d4874d..b2d46580 100644 --- a/src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj +++ b/src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj @@ -297,7 +297,7 @@ - + diff --git a/src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj.filters b/src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj.filters index 9ccf129e..2a00620e 100644 --- a/src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj.filters +++ b/src/NativeProviders/Windows/MKL/MKLWrapper.vcxproj.filters @@ -35,7 +35,7 @@ - + Resource Files diff --git a/src/Numerics/Control.cs b/src/Numerics/Control.cs index fa1028fc..0e9570fd 100644 --- a/src/Numerics/Control.cs +++ b/src/Numerics/Control.cs @@ -78,6 +78,10 @@ namespace MathNet.Numerics case "MKL": LinearAlgebraProvider = new Providers.LinearAlgebra.Mkl.MklLinearAlgebraProvider(); break; + + case "CUDA": + LinearAlgebraProvider = new Providers.LinearAlgebra.Cuda.CudaLinearAlgebraProvider(); + break; #endif default: LinearAlgebraProvider = new ManagedLinearAlgebraProvider(); @@ -127,6 +131,11 @@ namespace MathNet.Numerics { LinearAlgebraProvider = new Providers.LinearAlgebra.Mkl.MklLinearAlgebraProvider(consistency, precision, accuracy); } + + public static void UseNativeCUDA() + { + LinearAlgebraProvider = new Providers.LinearAlgebra.Cuda.CudaLinearAlgebraProvider(); + } #endif /// diff --git a/src/Numerics/Numerics.csproj b/src/Numerics/Numerics.csproj index b783234e..b451cf78 100644 --- a/src/Numerics/Numerics.csproj +++ b/src/Numerics/Numerics.csproj @@ -157,6 +157,12 @@ + + + + + + diff --git a/src/Numerics/Properties/AssemblyInfo.cs b/src/Numerics/Properties/AssemblyInfo.cs index 2281025c..930f5179 100644 --- a/src/Numerics/Properties/AssemblyInfo.cs +++ b/src/Numerics/Properties/AssemblyInfo.cs @@ -76,6 +76,7 @@ using System.Runtime.InteropServices; #else [assembly: InternalsVisibleTo("MathNet.Numerics.UnitTests")] [assembly: InternalsVisibleTo("MathNet.Numerics.UnitTestsMKL")] +[assembly: InternalsVisibleTo("MathNet.Numerics.UnitTestsCUDA")] [assembly: InternalsVisibleTo("Performance")] #endif diff --git a/src/Numerics/Properties/Resources.Designer.cs b/src/Numerics/Properties/Resources.Designer.cs index 0c7fae4d..443ba150 100644 --- a/src/Numerics/Properties/Resources.Designer.cs +++ b/src/Numerics/Properties/Resources.Designer.cs @@ -1,15 +1,13 @@ //------------------------------------------------------------------------------ // // This code was generated by a tool. -// Runtime Version:4.0.30319.34209 +// Runtime Version:4.0.30319.34014 // // Changes to this file may cause incorrect behavior and will be lost if // the code is regenerated. // //------------------------------------------------------------------------------ -using System.Reflection; - namespace MathNet.Numerics.Properties { using System; @@ -40,18 +38,10 @@ namespace MathNet.Numerics.Properties { [global::System.ComponentModel.EditorBrowsableAttribute(global::System.ComponentModel.EditorBrowsableState.Advanced)] public static global::System.Resources.ResourceManager ResourceManager { get { -#if NET45REFLECTION - if (object.ReferenceEquals(resourceMan, null)) - { - global::System.Resources.ResourceManager temp = new global::System.Resources.ResourceManager("MathNet.Numerics.Properties.Resources", typeof(Resources).GetTypeInfo().Assembly); - resourceMan = temp; - } -#else if (object.ReferenceEquals(resourceMan, null)) { global::System.Resources.ResourceManager temp = new global::System.Resources.ResourceManager("MathNet.Numerics.Properties.Resources", typeof(Resources).Assembly); resourceMan = temp; } -#endif return resourceMan; } } @@ -935,6 +925,15 @@ namespace MathNet.Numerics.Properties { } } + /// + /// Looks up a localized string similar to User work buffers are not supported by this provider.. + /// + public static string UserWorkBufferNotSupported { + get { + return ResourceManager.GetString("UserWorkBufferNotSupported", resourceCulture); + } + } + /// /// Looks up a localized string similar to Vectors can not be empty and must have at least one element.. /// diff --git a/src/Numerics/Properties/Resources.resx b/src/Numerics/Properties/Resources.resx index d46d44c8..e1360503 100644 --- a/src/Numerics/Properties/Resources.resx +++ b/src/Numerics/Properties/Resources.resx @@ -412,4 +412,7 @@ Vectors can not be empty and must have at least one element. + + User work buffers are not supported by this provider. + \ No newline at end of file diff --git a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex.cs b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex.cs new file mode 100644 index 00000000..31abfbc5 --- /dev/null +++ b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex.cs @@ -0,0 +1,705 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2013 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +#if NATIVE + +using System; +using System.Numerics; +using System.Security; +using MathNet.Numerics.LinearAlgebra.Factorization; +using MathNet.Numerics.Properties; + +namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda +{ + /// + /// NVidia's CUDA Toolkit linear algebra provider. + /// + public partial class CudaLinearAlgebraProvider + { + /// + /// Computes the dot product of x and y. + /// + /// The vector x. + /// The vector y. + /// The dot product of x and y. + /// This is equivalent to the DOT BLAS routine. + [SecuritySafeCritical] + public override Complex DotProduct(Complex[] x, Complex[] y) + { + if (y == null) + { + throw new ArgumentNullException("y"); + } + + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (x.Length != y.Length) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength); + } + + return SafeNativeMethods.z_dot_product(_blasHandle, x.Length, x, y); + } + + /// + /// Adds a scaled vector to another: result = y + alpha*x. + /// + /// The vector to update. + /// The value to scale by. + /// The vector to add to . + /// The result of the addition. + /// This is similar to the AXPY BLAS routine. + [SecuritySafeCritical] + public override void AddVectorToScaledVector(Complex[] y, Complex alpha, Complex[] x, Complex[] result) + { + if (y == null) + { + throw new ArgumentNullException("y"); + } + + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (y.Length != x.Length) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength); + } + + if (!ReferenceEquals(y, result)) + { + Array.Copy(y, 0, result, 0, y.Length); + } + + if (alpha == Complex.Zero) + { + return; + } + + SafeNativeMethods.z_axpy(_blasHandle, y.Length, alpha, x, result); + } + + /// + /// Scales an array. Can be used to scale a vector and a matrix. + /// + /// The scalar. + /// The values to scale. + /// This result of the scaling. + /// This is similar to the SCAL BLAS routine. + [SecuritySafeCritical] + public override void ScaleArray(Complex alpha, Complex[] x, Complex[] result) + { + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (!ReferenceEquals(x, result)) + { + Array.Copy(x, 0, result, 0, x.Length); + } + + if (alpha == Complex.One) + { + return; + } + + SafeNativeMethods.z_scale(_blasHandle, x.Length, alpha, result); + } + + /// + /// Multiples two matrices. result = x * y + /// + /// The x matrix. + /// The number of rows in the x matrix. + /// The number of columns in the x matrix. + /// The y matrix. + /// The number of rows in the y matrix. + /// The number of columns in the y matrix. + /// Where to store the result of the multiplication. + /// This is a simplified version of the BLAS GEMM routine with alpha + /// set to Complex.One and beta set to Complex.Zero, and x and y are not transposed. + public override void MatrixMultiply(Complex[] x, int rowsX, int columnsX, Complex[] y, int rowsY, int columnsY, Complex[] result) + { + MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, Complex.One, x, rowsX, columnsX, y, rowsY, columnsY, Complex.Zero, result); + } + + /// + /// Multiplies two matrices and updates another with the result. c = alpha*op(a)*op(b) + beta*c + /// + /// How to transpose the matrix. + /// How to transpose the matrix. + /// The value to scale matrix. + /// The a matrix. + /// The number of rows in the matrix. + /// The number of columns in the matrix. + /// The b matrix + /// The number of rows in the matrix. + /// The number of columns in the matrix. + /// The value to scale the matrix. + /// The c matrix. + [SecuritySafeCritical] + public override void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex alpha, Complex[] a, int rowsA, int columnsA, Complex[] b, int rowsB, int columnsB, Complex beta, Complex[] c) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (c == null) + { + throw new ArgumentNullException("c"); + } + + var m = transposeA == Transpose.DontTranspose ? rowsA : columnsA; + var n = transposeB == Transpose.DontTranspose ? columnsB : rowsB; + var k = transposeA == Transpose.DontTranspose ? columnsA : rowsA; + var l = transposeB == Transpose.DontTranspose ? rowsB : columnsB; + + if (c.Length != m*n) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + + if (k != l) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + + SafeNativeMethods.z_matrix_multiply(_blasHandle, transposeA.ToCUDA(), transposeB.ToCUDA(), m, n, k, alpha, a, b, beta, c); + } + + /// + /// Computes the LUP factorization of A. P*A = L*U. + /// + /// An by matrix. The matrix is overwritten with the + /// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of (the diagonal is always Complex.One + /// for the L factor). The upper triangular factor U is stored on and above the diagonal of . + /// The order of the square matrix . + /// On exit, it contains the pivot indices. The size of the array must be . + /// This is equivalent to the GETRF LAPACK routine. + [SecuritySafeCritical] + public override void LUFactor(Complex[] data, int order, int[] ipiv) + { + if (data == null) + { + throw new ArgumentNullException("data"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (data.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "data"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + Solver(SafeNativeMethods.z_lu_factor(_solverHandle, order, data, ipiv)); + } + + /// + /// Computes the inverse of matrix using LU factorization. + /// + /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . + /// This is equivalent to the GETRF and GETRI LAPACK routines. + [SecuritySafeCritical] + public override void LUInverse(Complex[] a, int order) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + Solver(SafeNativeMethods.z_lu_inverse(_solverHandle, _blasHandle, order, a)); + } + + /// + /// Computes the inverse of a previously factored matrix. + /// + /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . + /// The pivot indices of . + /// This is equivalent to the GETRI LAPACK routine. + [SecuritySafeCritical] + public override void LUInverseFactored(Complex[] a, int order, int[] ipiv) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + BLAS(SafeNativeMethods.z_lu_inverse_factored(_blasHandle, order, a, ipiv)); + } + + /// + /// Computes the inverse of matrix using LU factorization. + /// + /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . + /// Not supported. Should be left null. + /// This is equivalent to the GETRF and GETRI LAPACK routines. + [SecuritySafeCritical] + public override void LUInverse(Complex[] a, int order, Complex[] work) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (work != null) + { + throw new ArgumentException(Resources.UserWorkBufferNotSupported); + } + + Solver(SafeNativeMethods.z_lu_inverse(_solverHandle, _blasHandle, order, a)); + } + + /// + /// Computes the inverse of a previously factored matrix. + /// + /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . + /// The pivot indices of . + /// Not supported. Should be left null. + /// This is equivalent to the GETRI LAPACK routine. + [SecuritySafeCritical] + public override void LUInverseFactored(Complex[] a, int order, int[] ipiv, Complex[] work) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + if (work != null) + { + throw new ArgumentException(Resources.UserWorkBufferNotSupported); + } + + BLAS(SafeNativeMethods.z_lu_inverse_factored(_blasHandle, order, a, ipiv)); + } + + /// + /// Solves A*X=B for X using LU factorization. + /// + /// The number of columns of B. + /// The square matrix A. + /// The order of the square matrix . + /// On entry the B matrix; on exit the X matrix. + /// This is equivalent to the GETRF and GETRS LAPACK routines. + [SecuritySafeCritical] + public override void LUSolve(int columnsOfB, Complex[] a, int order, Complex[] b) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (b.Length != columnsOfB*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + Solver(SafeNativeMethods.z_lu_solve(_solverHandle, order, columnsOfB, a, b)); + } + + /// + /// Solves A*X=B for X using a previously factored A matrix. + /// + /// The number of columns of B. + /// The factored A matrix. + /// The order of the square matrix . + /// The pivot indices of . + /// On entry the B matrix; on exit the X matrix. + /// This is equivalent to the GETRS LAPACK routine. + [SecuritySafeCritical] + public override void LUSolveFactored(int columnsOfB, Complex[] a, int order, int[] ipiv, Complex[] b) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + if (b.Length != columnsOfB*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + Solver(SafeNativeMethods.z_lu_solve_factored(_solverHandle, order, columnsOfB, a, ipiv, b)); + } + + /// + /// Computes the Cholesky factorization of A. + /// + /// On entry, a square, positive definite matrix. On exit, the matrix is overwritten with the + /// the Cholesky factorization. + /// The number of rows or columns in the matrix. + /// This is equivalent to the POTRF LAPACK routine. + [SecuritySafeCritical] + public override void CholeskyFactor(Complex[] a, int order) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (order < 1) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "order"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + Solver(SafeNativeMethods.z_cholesky_factor(_solverHandle, order, a)); + } + + /// + /// Solves A*X=B for X using Cholesky factorization. + /// + /// The square, positive definite matrix A. + /// The number of rows and columns in A. + /// On entry the B matrix; on exit the X matrix. + /// The number of columns in the B matrix. + /// This is equivalent to the POTRF add POTRS LAPACK routines. + /// + [SecuritySafeCritical] + public override void CholeskySolve(Complex[] a, int orderA, Complex[] b, int columnsB) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (b.Length != orderA*columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + Solver(SafeNativeMethods.z_cholesky_solve(_solverHandle, orderA, columnsB, a, b)); + } + + /// + /// Solves A*X=B for X using a previously factored A matrix. + /// + /// The square, positive definite matrix A. + /// The number of rows and columns in A. + /// On entry the B matrix; on exit the X matrix. + /// The number of columns in the B matrix. + /// This is equivalent to the POTRS LAPACK routine. + [SecuritySafeCritical] + public override void CholeskySolveFactored(Complex[] a, int orderA, Complex[] b, int columnsB) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (b.Length != orderA*columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + Solver(SafeNativeMethods.z_cholesky_solve_factored(_solverHandle, orderA, columnsB, a, b)); + } + + /// + /// Computes the singular value decomposition of A. + /// + /// Compute the singular U and VT vectors or not. + /// On entry, the M by N matrix to decompose. On exit, A may be overwritten. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The singular values of A in ascending value. + /// If is true, on exit U contains the left + /// singular vectors. + /// If is true, on exit VT contains the transposed + /// right singular vectors. + /// This is equivalent to the GESVD LAPACK routine. + [SecuritySafeCritical] + public override void SingularValueDecomposition(bool computeVectors, Complex[] a, int rowsA, int columnsA, Complex[] s, Complex[] u, Complex[] vt) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (s == null) + { + throw new ArgumentNullException("s"); + } + + if (u == null) + { + throw new ArgumentNullException("u"); + } + + if (vt == null) + { + throw new ArgumentNullException("vt"); + } + + if (u.Length != rowsA*rowsA) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); + } + + if (vt.Length != columnsA*columnsA) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); + } + + if (s.Length != Math.Min(rowsA, columnsA)) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); + } + + SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, null); + } + + /// + /// Solves A*X=B for X using the singular value decomposition of A. + /// + /// On entry, the M by N matrix to decompose. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + public override void SvdSolve(Complex[] a, int rowsA, int columnsA, Complex[] b, int columnsB, Complex[] x) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (b.Length != rowsA*columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (x.Length != columnsA*columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + var s = new Complex[Math.Min(rowsA, columnsA)]; + var u = new Complex[rowsA*rowsA]; + var vt = new Complex[columnsA*columnsA]; + + var clone = new Complex[a.Length]; + a.Copy(clone); + SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt, null); + SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); + } + + /// + /// Computes the singular value decomposition of A. + /// + /// Compute the singular U and VT vectors or not. + /// On entry, the M by N matrix to decompose. On exit, A may be overwritten. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The singular values of A in ascending value. + /// If is true, on exit U contains the left + /// singular vectors. + /// If is true, on exit VT contains the transposed + /// right singular vectors. + /// User work buffers are not supported. Should be null. + /// This is equivalent to the GESVD LAPACK routine. + [SecuritySafeCritical] + public override void SingularValueDecomposition(bool computeVectors, Complex[] a, int rowsA, int columnsA, Complex[] s, Complex[] u, Complex[] vt, Complex[] work) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (s == null) + { + throw new ArgumentNullException("s"); + } + + if (u == null) + { + throw new ArgumentNullException("u"); + } + + if (vt == null) + { + throw new ArgumentNullException("vt"); + } + + if (work != null) + { + throw new ArgumentException(Resources.UserWorkBufferNotSupported); + } + + if (u.Length != rowsA*rowsA) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); + } + + if (vt.Length != columnsA*columnsA) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); + } + + if (s.Length != Math.Min(rowsA, columnsA)) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); + } + + if (columnsA > rowsA || !computeVectors) // see remarks http://docs.nvidia.com/cuda/cusolver/index.html#cuds-lt-t-gt-gesvd + base.SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, new Complex[rowsA]); + else Solver(SafeNativeMethods.z_svd_factor(_solverHandle, computeVectors, rowsA, columnsA, a, s, u, vt)); + } + } +} + +#endif diff --git a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex32.cs b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex32.cs new file mode 100644 index 00000000..81679c15 --- /dev/null +++ b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Complex32.cs @@ -0,0 +1,705 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2013 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +#if NATIVE + +using System; +using System.Numerics; +using System.Security; +using MathNet.Numerics.LinearAlgebra.Factorization; +using MathNet.Numerics.Properties; + +namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda +{ + /// + /// NVidia's CUDA Toolkit linear algebra provider. + /// + public partial class CudaLinearAlgebraProvider + { + /// + /// Computes the dot product of x and y. + /// + /// The vector x. + /// The vector y. + /// The dot product of x and y. + /// This is equivalent to the DOT BLAS routine. + [SecuritySafeCritical] + public override Complex32 DotProduct(Complex32[] x, Complex32[] y) + { + if (y == null) + { + throw new ArgumentNullException("y"); + } + + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (x.Length != y.Length) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength); + } + + return SafeNativeMethods.c_dot_product(_blasHandle, x.Length, x, y); + } + + /// + /// Adds a scaled vector to another: result = y + alpha*x. + /// + /// The vector to update. + /// The value to scale by. + /// The vector to add to . + /// The result of the addition. + /// This is similar to the AXPY BLAS routine. + [SecuritySafeCritical] + public override void AddVectorToScaledVector(Complex32[] y, Complex32 alpha, Complex32[] x, Complex32[] result) + { + if (y == null) + { + throw new ArgumentNullException("y"); + } + + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (y.Length != x.Length) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength); + } + + if (!ReferenceEquals(y, result)) + { + Array.Copy(y, 0, result, 0, y.Length); + } + + if (alpha == Complex32.Zero) + { + return; + } + + SafeNativeMethods.c_axpy(_blasHandle, y.Length, alpha, x, result); + } + + /// + /// Scales an array. Can be used to scale a vector and a matrix. + /// + /// The scalar. + /// The values to scale. + /// This result of the scaling. + /// This is similar to the SCAL BLAS routine. + [SecuritySafeCritical] + public override void ScaleArray(Complex32 alpha, Complex32[] x, Complex32[] result) + { + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (!ReferenceEquals(x, result)) + { + Array.Copy(x, 0, result, 0, x.Length); + } + + if (alpha == Complex32.One) + { + return; + } + + SafeNativeMethods.c_scale(_blasHandle, x.Length, alpha, result); + } + + /// + /// Multiples two matrices. result = x * y + /// + /// The x matrix. + /// The number of rows in the x matrix. + /// The number of columns in the x matrix. + /// The y matrix. + /// The number of rows in the y matrix. + /// The number of columns in the y matrix. + /// Where to store the result of the multiplication. + /// This is a simplified version of the BLAS GEMM routine with alpha + /// set to Complex32.One and beta set to Complex32.Zero, and x and y are not transposed. + public override void MatrixMultiply(Complex32[] x, int rowsX, int columnsX, Complex32[] y, int rowsY, int columnsY, Complex32[] result) + { + MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, Complex32.One, x, rowsX, columnsX, y, rowsY, columnsY, Complex32.Zero, result); + } + + /// + /// Multiplies two matrices and updates another with the result. c = alpha*op(a)*op(b) + beta*c + /// + /// How to transpose the matrix. + /// How to transpose the matrix. + /// The value to scale matrix. + /// The a matrix. + /// The number of rows in the matrix. + /// The number of columns in the matrix. + /// The b matrix + /// The number of rows in the matrix. + /// The number of columns in the matrix. + /// The value to scale the matrix. + /// The c matrix. + [SecuritySafeCritical] + public override void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, Complex32 alpha, Complex32[] a, int rowsA, int columnsA, Complex32[] b, int rowsB, int columnsB, Complex32 beta, Complex32[] c) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (c == null) + { + throw new ArgumentNullException("c"); + } + + var m = transposeA == Transpose.DontTranspose ? rowsA : columnsA; + var n = transposeB == Transpose.DontTranspose ? columnsB : rowsB; + var k = transposeA == Transpose.DontTranspose ? columnsA : rowsA; + var l = transposeB == Transpose.DontTranspose ? rowsB : columnsB; + + if (c.Length != m*n) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + + if (k != l) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + + SafeNativeMethods.c_matrix_multiply(_blasHandle, transposeA.ToCUDA(), transposeB.ToCUDA(), m, n, k, alpha, a, b, beta, c); + } + + /// + /// Computes the LUP factorization of A. P*A = L*U. + /// + /// An by matrix. The matrix is overwritten with the + /// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of (the diagonal is always Complex32.One + /// for the L factor). The upper triangular factor U is stored on and above the diagonal of . + /// The order of the square matrix . + /// On exit, it contains the pivot indices. The size of the array must be . + /// This is equivalent to the GETRF LAPACK routine. + [SecuritySafeCritical] + public override void LUFactor(Complex32[] data, int order, int[] ipiv) + { + if (data == null) + { + throw new ArgumentNullException("data"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (data.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "data"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + Solver(SafeNativeMethods.c_lu_factor(_solverHandle, order, data, ipiv)); + } + + /// + /// Computes the inverse of matrix using LU factorization. + /// + /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . + /// This is equivalent to the GETRF and GETRI LAPACK routines. + [SecuritySafeCritical] + public override void LUInverse(Complex32[] a, int order) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + Solver(SafeNativeMethods.c_lu_inverse(_solverHandle, _blasHandle, order, a)); + } + + /// + /// Computes the inverse of a previously factored matrix. + /// + /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . + /// The pivot indices of . + /// This is equivalent to the GETRI LAPACK routine. + [SecuritySafeCritical] + public override void LUInverseFactored(Complex32[] a, int order, int[] ipiv) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + BLAS(SafeNativeMethods.c_lu_inverse_factored(_blasHandle, order, a, ipiv)); + } + + /// + /// Computes the inverse of matrix using LU factorization. + /// + /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . + /// Not supported. Should be left null. + /// This is equivalent to the GETRF and GETRI LAPACK routines. + [SecuritySafeCritical] + public override void LUInverse(Complex32[] a, int order, Complex32[] work) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (work != null) + { + throw new ArgumentException(Resources.UserWorkBufferNotSupported); + } + + Solver(SafeNativeMethods.c_lu_inverse(_solverHandle, _blasHandle, order, a)); + } + + /// + /// Computes the inverse of a previously factored matrix. + /// + /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . + /// The pivot indices of . + /// Not supported. Should be left null. + /// This is equivalent to the GETRI LAPACK routine. + [SecuritySafeCritical] + public override void LUInverseFactored(Complex32[] a, int order, int[] ipiv, Complex32[] work) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + if (work != null) + { + throw new ArgumentException(Resources.UserWorkBufferNotSupported); + } + + BLAS(SafeNativeMethods.c_lu_inverse_factored(_blasHandle, order, a, ipiv)); + } + + /// + /// Solves A*X=B for X using LU factorization. + /// + /// The number of columns of B. + /// The square matrix A. + /// The order of the square matrix . + /// On entry the B matrix; on exit the X matrix. + /// This is equivalent to the GETRF and GETRS LAPACK routines. + [SecuritySafeCritical] + public override void LUSolve(int columnsOfB, Complex32[] a, int order, Complex32[] b) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (b.Length != columnsOfB*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + Solver(SafeNativeMethods.c_lu_solve(_solverHandle, order, columnsOfB, a, b)); + } + + /// + /// Solves A*X=B for X using a previously factored A matrix. + /// + /// The number of columns of B. + /// The factored A matrix. + /// The order of the square matrix . + /// The pivot indices of . + /// On entry the B matrix; on exit the X matrix. + /// This is equivalent to the GETRS LAPACK routine. + [SecuritySafeCritical] + public override void LUSolveFactored(int columnsOfB, Complex32[] a, int order, int[] ipiv, Complex32[] b) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + if (b.Length != columnsOfB*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + Solver(SafeNativeMethods.c_lu_solve_factored(_solverHandle, order, columnsOfB, a, ipiv, b)); + } + + /// + /// Computes the Cholesky factorization of A. + /// + /// On entry, a square, positive definite matrix. On exit, the matrix is overwritten with the + /// the Cholesky factorization. + /// The number of rows or columns in the matrix. + /// This is equivalent to the POTRF LAPACK routine. + [SecuritySafeCritical] + public override void CholeskyFactor(Complex32[] a, int order) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (order < 1) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "order"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + Solver(SafeNativeMethods.c_cholesky_factor(_solverHandle, order, a)); + } + + /// + /// Solves A*X=B for X using Cholesky factorization. + /// + /// The square, positive definite matrix A. + /// The number of rows and columns in A. + /// On entry the B matrix; on exit the X matrix. + /// The number of columns in the B matrix. + /// This is equivalent to the POTRF add POTRS LAPACK routines. + /// + [SecuritySafeCritical] + public override void CholeskySolve(Complex32[] a, int orderA, Complex32[] b, int columnsB) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (b.Length != orderA*columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + Solver(SafeNativeMethods.c_cholesky_solve(_solverHandle, orderA, columnsB, a, b)); + } + + /// + /// Solves A*X=B for X using a previously factored A matrix. + /// + /// The square, positive definite matrix A. + /// The number of rows and columns in A. + /// On entry the B matrix; on exit the X matrix. + /// The number of columns in the B matrix. + /// This is equivalent to the POTRS LAPACK routine. + [SecuritySafeCritical] + public override void CholeskySolveFactored(Complex32[] a, int orderA, Complex32[] b, int columnsB) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (b.Length != orderA*columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + Solver(SafeNativeMethods.c_cholesky_solve_factored(_solverHandle, orderA, columnsB, a, b)); + } + + /// + /// Computes the singular value decomposition of A. + /// + /// Compute the singular U and VT vectors or not. + /// On entry, the M by N matrix to decompose. On exit, A may be overwritten. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The singular values of A in ascending value. + /// If is true, on exit U contains the left + /// singular vectors. + /// If is true, on exit VT contains the transposed + /// right singular vectors. + /// This is equivalent to the GESVD LAPACK routine. + [SecuritySafeCritical] + public override void SingularValueDecomposition(bool computeVectors, Complex32[] a, int rowsA, int columnsA, Complex32[] s, Complex32[] u, Complex32[] vt) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (s == null) + { + throw new ArgumentNullException("s"); + } + + if (u == null) + { + throw new ArgumentNullException("u"); + } + + if (vt == null) + { + throw new ArgumentNullException("vt"); + } + + if (u.Length != rowsA*rowsA) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); + } + + if (vt.Length != columnsA*columnsA) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); + } + + if (s.Length != Math.Min(rowsA, columnsA)) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); + } + + SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, null); + } + + /// + /// Solves A*X=B for X using the singular value decomposition of A. + /// + /// On entry, the M by N matrix to decompose. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + public override void SvdSolve(Complex32[] a, int rowsA, int columnsA, Complex32[] b, int columnsB, Complex32[] x) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (b.Length != rowsA*columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (x.Length != columnsA*columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + var s = new Complex32[Math.Min(rowsA, columnsA)]; + var u = new Complex32[rowsA*rowsA]; + var vt = new Complex32[columnsA*columnsA]; + + var clone = new Complex32[a.Length]; + a.Copy(clone); + SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt, null); + SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); + } + + /// + /// Computes the singular value decomposition of A. + /// + /// Compute the singular U and VT vectors or not. + /// On entry, the M by N matrix to decompose. On exit, A may be overwritten. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The singular values of A in ascending value. + /// If is true, on exit U contains the left + /// singular vectors. + /// If is true, on exit VT contains the transposed + /// right singular vectors. + /// Not supported. Should be left null. + /// This is equivalent to the GESVD LAPACK routine. + [SecuritySafeCritical] + public override void SingularValueDecomposition(bool computeVectors, Complex32[] a, int rowsA, int columnsA, Complex32[] s, Complex32[] u, Complex32[] vt, Complex32[] work) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (s == null) + { + throw new ArgumentNullException("s"); + } + + if (u == null) + { + throw new ArgumentNullException("u"); + } + + if (vt == null) + { + throw new ArgumentNullException("vt"); + } + + if (work != null) + { + throw new ArgumentException(Resources.UserWorkBufferNotSupported); + } + + if (u.Length != rowsA*rowsA) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); + } + + if (vt.Length != columnsA*columnsA) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); + } + + if (s.Length != Math.Min(rowsA, columnsA)) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); + } + + if (columnsA > rowsA || !computeVectors) // see remarks http://docs.nvidia.com/cuda/cusolver/index.html#cuds-lt-t-gt-gesvd + base.SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, new Complex32[rowsA]); + else Solver(SafeNativeMethods.c_svd_factor(_solverHandle, computeVectors, rowsA, columnsA, a, s, u, vt)); + } + } +} + +#endif diff --git a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Double.cs b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Double.cs new file mode 100644 index 00000000..6e73ed46 --- /dev/null +++ b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Double.cs @@ -0,0 +1,705 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2013 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +#if NATIVE + +using System; +using System.Numerics; +using System.Security; +using MathNet.Numerics.LinearAlgebra.Factorization; +using MathNet.Numerics.Properties; + +namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda +{ + /// + /// NVidia's CUDA Toolkit linear algebra provider. + /// + public partial class CudaLinearAlgebraProvider + { + /// + /// Computes the dot product of x and y. + /// + /// The vector x. + /// The vector y. + /// The dot product of x and y. + /// This is equivalent to the DOT BLAS routine. + [SecuritySafeCritical] + public override double DotProduct(double[] x, double[] y) + { + if (y == null) + { + throw new ArgumentNullException("y"); + } + + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (x.Length != y.Length) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength); + } + + return SafeNativeMethods.d_dot_product(_blasHandle, x.Length, x, y); + } + + /// + /// Adds a scaled vector to another: result = y + alpha*x. + /// + /// The vector to update. + /// The value to scale by. + /// The vector to add to . + /// The result of the addition. + /// This is similar to the AXPY BLAS routine. + [SecuritySafeCritical] + public override void AddVectorToScaledVector(double[] y, double alpha, double[] x, double[] result) + { + if (y == null) + { + throw new ArgumentNullException("y"); + } + + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (y.Length != x.Length) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength); + } + + if (!ReferenceEquals(y, result)) + { + Array.Copy(y, 0, result, 0, y.Length); + } + + if (alpha == 0.0) + { + return; + } + + SafeNativeMethods.d_axpy(_blasHandle, y.Length, alpha, x, result); + } + + /// + /// Scales an array. Can be used to scale a vector and a matrix. + /// + /// The scalar. + /// The values to scale. + /// This result of the scaling. + /// This is similar to the SCAL BLAS routine. + [SecuritySafeCritical] + public override void ScaleArray(double alpha, double[] x, double[] result) + { + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (!ReferenceEquals(x, result)) + { + Array.Copy(x, 0, result, 0, x.Length); + } + + if (alpha == 1.0) + { + return; + } + + SafeNativeMethods.d_scale(_blasHandle, x.Length, alpha, result); + } + + /// + /// Multiples two matrices. result = x * y + /// + /// The x matrix. + /// The number of rows in the x matrix. + /// The number of columns in the x matrix. + /// The y matrix. + /// The number of rows in the y matrix. + /// The number of columns in the y matrix. + /// Where to store the result of the multiplication. + /// This is a simplified version of the BLAS GEMM routine with alpha + /// set to 1.0 and beta set to 0.0, and x and y are not transposed. + public override void MatrixMultiply(double[] x, int rowsX, int columnsX, double[] y, int rowsY, int columnsY, double[] result) + { + MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, 1.0, x, rowsX, columnsX, y, rowsY, columnsY, 0.0, result); + } + + /// + /// Multiplies two matrices and updates another with the result. c = alpha*op(a)*op(b) + beta*c + /// + /// How to transpose the matrix. + /// How to transpose the matrix. + /// The value to scale matrix. + /// The a matrix. + /// The number of rows in the matrix. + /// The number of columns in the matrix. + /// The b matrix + /// The number of rows in the matrix. + /// The number of columns in the matrix. + /// The value to scale the matrix. + /// The c matrix. + [SecuritySafeCritical] + public override void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, double alpha, double[] a, int rowsA, int columnsA, double[] b, int rowsB, int columnsB, double beta, double[] c) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (c == null) + { + throw new ArgumentNullException("c"); + } + + var m = transposeA == Transpose.DontTranspose ? rowsA : columnsA; + var n = transposeB == Transpose.DontTranspose ? columnsB : rowsB; + var k = transposeA == Transpose.DontTranspose ? columnsA : rowsA; + var l = transposeB == Transpose.DontTranspose ? rowsB : columnsB; + + if (c.Length != m*n) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + + if (k != l) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + + SafeNativeMethods.d_matrix_multiply(_blasHandle, transposeA.ToCUDA(), transposeB.ToCUDA(), m, n, k, alpha, a, b, beta, c); + } + + /// + /// Computes the LUP factorization of A. P*A = L*U. + /// + /// An by matrix. The matrix is overwritten with the + /// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of (the diagonal is always 1.0 + /// for the L factor). The upper triangular factor U is stored on and above the diagonal of . + /// The order of the square matrix . + /// On exit, it contains the pivot indices. The size of the array must be . + /// This is equivalent to the GETRF LAPACK routine. + [SecuritySafeCritical] + public override void LUFactor(double[] data, int order, int[] ipiv) + { + if (data == null) + { + throw new ArgumentNullException("data"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (data.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "data"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + Solver(SafeNativeMethods.d_lu_factor(_solverHandle, order, data, ipiv)); + } + + /// + /// Computes the inverse of matrix using LU factorization. + /// + /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . + /// This is equivalent to the GETRF and GETRI LAPACK routines. + [SecuritySafeCritical] + public override void LUInverse(double[] a, int order) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + Solver(SafeNativeMethods.d_lu_inverse(_solverHandle, _blasHandle, order, a)); + } + + /// + /// Computes the inverse of a previously factored matrix. + /// + /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . + /// The pivot indices of . + /// This is equivalent to the GETRI LAPACK routine. + [SecuritySafeCritical] + public override void LUInverseFactored(double[] a, int order, int[] ipiv) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + BLAS(SafeNativeMethods.d_lu_inverse_factored(_blasHandle, order, a, ipiv)); + } + + /// + /// Computes the inverse of matrix using LU factorization. + /// + /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . + /// Not supported. Should be left null. + /// This is equivalent to the GETRF and GETRI LAPACK routines. + [SecuritySafeCritical] + public override void LUInverse(double[] a, int order, double[] work) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (work != null) + { + throw new ArgumentException(Resources.UserWorkBufferNotSupported); + } + + Solver(SafeNativeMethods.d_lu_inverse(_solverHandle, _blasHandle, order, a)); + } + + /// + /// Computes the inverse of a previously factored matrix. + /// + /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . + /// The pivot indices of . + /// Not supported. Should be left null. + /// This is equivalent to the GETRI LAPACK routine. + [SecuritySafeCritical] + public override void LUInverseFactored(double[] a, int order, int[] ipiv, double[] work) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + if (work != null) + { + throw new ArgumentException(Resources.UserWorkBufferNotSupported); + } + + BLAS(SafeNativeMethods.d_lu_inverse_factored(_blasHandle, order, a, ipiv)); + } + + /// + /// Solves A*X=B for X using LU factorization. + /// + /// The number of columns of B. + /// The square matrix A. + /// The order of the square matrix . + /// On entry the B matrix; on exit the X matrix. + /// This is equivalent to the GETRF and GETRS LAPACK routines. + [SecuritySafeCritical] + public override void LUSolve(int columnsOfB, double[] a, int order, double[] b) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (b.Length != columnsOfB*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + Solver(SafeNativeMethods.d_lu_solve(_solverHandle, order, columnsOfB, a, b)); + } + + /// + /// Solves A*X=B for X using a previously factored A matrix. + /// + /// The number of columns of B. + /// The factored A matrix. + /// The order of the square matrix . + /// The pivot indices of . + /// On entry the B matrix; on exit the X matrix. + /// This is equivalent to the GETRS LAPACK routine. + [SecuritySafeCritical] + public override void LUSolveFactored(int columnsOfB, double[] a, int order, int[] ipiv, double[] b) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + if (b.Length != columnsOfB*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + Solver(SafeNativeMethods.d_lu_solve_factored(_solverHandle, order, columnsOfB, a, ipiv, b)); + } + + /// + /// Computes the Cholesky factorization of A. + /// + /// On entry, a square, positive definite matrix. On exit, the matrix is overwritten with the + /// the Cholesky factorization. + /// The number of rows or columns in the matrix. + /// This is equivalent to the POTRF LAPACK routine. + [SecuritySafeCritical] + public override void CholeskyFactor(double[] a, int order) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (order < 1) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "order"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + Solver(SafeNativeMethods.d_cholesky_factor(_solverHandle, order, a)); + } + + /// + /// Solves A*X=B for X using Cholesky factorization. + /// + /// The square, positive definite matrix A. + /// The number of rows and columns in A. + /// On entry the B matrix; on exit the X matrix. + /// The number of columns in the B matrix. + /// This is equivalent to the POTRF add POTRS LAPACK routines. + /// + [SecuritySafeCritical] + public override void CholeskySolve(double[] a, int orderA, double[] b, int columnsB) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (b.Length != orderA*columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + Solver(SafeNativeMethods.d_cholesky_solve(_solverHandle, orderA, columnsB, a, b)); + } + + /// + /// Solves A*X=B for X using a previously factored A matrix. + /// + /// The square, positive definite matrix A. + /// The number of rows and columns in A. + /// On entry the B matrix; on exit the X matrix. + /// The number of columns in the B matrix. + /// This is equivalent to the POTRS LAPACK routine. + [SecuritySafeCritical] + public override void CholeskySolveFactored(double[] a, int orderA, double[] b, int columnsB) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (b.Length != orderA*columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + Solver(SafeNativeMethods.d_cholesky_solve_factored(_solverHandle, orderA, columnsB, a, b)); + } + + /// + /// Computes the singular value decomposition of A. + /// + /// Compute the singular U and VT vectors or not. + /// On entry, the M by N matrix to decompose. On exit, A may be overwritten. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The singular values of A in ascending value. + /// If is true, on exit U contains the left + /// singular vectors. + /// If is true, on exit VT contains the transposed + /// right singular vectors. + /// This is equivalent to the GESVD LAPACK routine. + [SecuritySafeCritical] + public override void SingularValueDecomposition(bool computeVectors, double[] a, int rowsA, int columnsA, double[] s, double[] u, double[] vt) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (s == null) + { + throw new ArgumentNullException("s"); + } + + if (u == null) + { + throw new ArgumentNullException("u"); + } + + if (vt == null) + { + throw new ArgumentNullException("vt"); + } + + if (u.Length != rowsA*rowsA) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); + } + + if (vt.Length != columnsA*columnsA) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); + } + + if (s.Length != Math.Min(rowsA, columnsA)) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); + } + + SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, null); + } + + /// + /// Solves A*X=B for X using the singular value decomposition of A. + /// + /// On entry, the M by N matrix to decompose. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + public override void SvdSolve(double[] a, int rowsA, int columnsA, double[] b, int columnsB, double[] x) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (b.Length != rowsA*columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (x.Length != columnsA*columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + var s = new double[Math.Min(rowsA, columnsA)]; + var u = new double[rowsA*rowsA]; + var vt = new double[columnsA*columnsA]; + + var clone = new double[a.Length]; + a.Copy(clone); + SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt); + SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); + } + + /// + /// Computes the singular value decomposition of A. + /// + /// Compute the singular U and VT vectors or not. + /// On entry, the M by N matrix to decompose. On exit, A may be overwritten. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The singular values of A in ascending value. + /// If is true, on exit U contains the left + /// singular vectors. + /// If is true, on exit VT contains the transposed + /// right singular vectors. + /// Not supported. Should be left null. + /// This is equivalent to the GESVD LAPACK routine. + [SecuritySafeCritical] + public override void SingularValueDecomposition(bool computeVectors, double[] a, int rowsA, int columnsA, double[] s, double[] u, double[] vt, double[] work) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (s == null) + { + throw new ArgumentNullException("s"); + } + + if (u == null) + { + throw new ArgumentNullException("u"); + } + + if (vt == null) + { + throw new ArgumentNullException("vt"); + } + + if (work != null) + { + throw new ArgumentException(Resources.UserWorkBufferNotSupported); + } + + if (u.Length != rowsA*rowsA) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); + } + + if (vt.Length != columnsA*columnsA) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); + } + + if (s.Length != Math.Min(rowsA, columnsA)) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); + } + + if (columnsA > rowsA || !computeVectors) // see remarks http://docs.nvidia.com/cuda/cusolver/index.html#cuds-lt-t-gt-gesvd + base.SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, new double[rowsA]); + else Solver (SafeNativeMethods.d_svd_factor(_solverHandle, computeVectors, rowsA, columnsA, a, s, u, vt)); + } + } +} + +#endif diff --git a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Single.cs b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Single.cs new file mode 100644 index 00000000..0ac4e34f --- /dev/null +++ b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.Single.cs @@ -0,0 +1,705 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2013 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +#if NATIVE + +using System; +using System.Numerics; +using System.Security; +using MathNet.Numerics.LinearAlgebra.Factorization; +using MathNet.Numerics.Properties; + +namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda +{ + /// + /// NVidia's CUDA Toolkit linear algebra provider. + /// + public partial class CudaLinearAlgebraProvider + { + /// + /// Computes the dot product of x and y. + /// + /// The vector x. + /// The vector y. + /// The dot product of x and y. + /// This is equivalent to the DOT BLAS routine. + [SecuritySafeCritical] + public override float DotProduct(float[] x, float[] y) + { + if (y == null) + { + throw new ArgumentNullException("y"); + } + + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (x.Length != y.Length) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength); + } + + return SafeNativeMethods.s_dot_product(_blasHandle, x.Length, x, y); + } + + /// + /// Adds a scaled vector to another: result = y + alpha*x. + /// + /// The vector to update. + /// The value to scale by. + /// The vector to add to . + /// The result of the addition. + /// This is similar to the AXPY BLAS routine. + [SecuritySafeCritical] + public override void AddVectorToScaledVector(float[] y, float alpha, float[] x, float[] result) + { + if (y == null) + { + throw new ArgumentNullException("y"); + } + + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (y.Length != x.Length) + { + throw new ArgumentException(Resources.ArgumentVectorsSameLength); + } + + if (!ReferenceEquals(y, result)) + { + Array.Copy(y, 0, result, 0, y.Length); + } + + if (alpha == 0.0f) + { + return; + } + + SafeNativeMethods.s_axpy(_blasHandle, y.Length, alpha, x, result); + } + + /// + /// Scales an array. Can be used to scale a vector and a matrix. + /// + /// The scalar. + /// The values to scale. + /// This result of the scaling. + /// This is similar to the SCAL BLAS routine. + [SecuritySafeCritical] + public override void ScaleArray(float alpha, float[] x, float[] result) + { + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (!ReferenceEquals(x, result)) + { + Array.Copy(x, 0, result, 0, x.Length); + } + + if (alpha == 1.0f) + { + return; + } + + SafeNativeMethods.s_scale(_blasHandle, x.Length, alpha, result); + } + + /// + /// Multiples two matrices. result = x * y + /// + /// The x matrix. + /// The number of rows in the x matrix. + /// The number of columns in the x matrix. + /// The y matrix. + /// The number of rows in the y matrix. + /// The number of columns in the y matrix. + /// Where to store the result of the multiplication. + /// This is a simplified version of the BLAS GEMM routine with alpha + /// set to 1.0f and beta set to 0.0f, and x and y are not transposed. + public override void MatrixMultiply(float[] x, int rowsX, int columnsX, float[] y, int rowsY, int columnsY, float[] result) + { + MatrixMultiplyWithUpdate(Transpose.DontTranspose, Transpose.DontTranspose, 1.0f, x, rowsX, columnsX, y, rowsY, columnsY, 0.0f, result); + } + + /// + /// Multiplies two matrices and updates another with the result. c = alpha*op(a)*op(b) + beta*c + /// + /// How to transpose the matrix. + /// How to transpose the matrix. + /// The value to scale matrix. + /// The a matrix. + /// The number of rows in the matrix. + /// The number of columns in the matrix. + /// The b matrix + /// The number of rows in the matrix. + /// The number of columns in the matrix. + /// The value to scale the matrix. + /// The c matrix. + [SecuritySafeCritical] + public override void MatrixMultiplyWithUpdate(Transpose transposeA, Transpose transposeB, float alpha, float[] a, int rowsA, int columnsA, float[] b, int rowsB, int columnsB, float beta, float[] c) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (c == null) + { + throw new ArgumentNullException("c"); + } + + var m = transposeA == Transpose.DontTranspose ? rowsA : columnsA; + var n = transposeB == Transpose.DontTranspose ? columnsB : rowsB; + var k = transposeA == Transpose.DontTranspose ? columnsA : rowsA; + var l = transposeB == Transpose.DontTranspose ? rowsB : columnsB; + + if (c.Length != m*n) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + + if (k != l) + { + throw new ArgumentException(Resources.ArgumentMatrixDimensions); + } + + SafeNativeMethods.s_matrix_multiply(_blasHandle, transposeA.ToCUDA(), transposeB.ToCUDA(), m, n, k, alpha, a, b, beta, c); + } + + /// + /// Computes the LUP factorization of A. P*A = L*U. + /// + /// An by matrix. The matrix is overwritten with the + /// the LU factorization on exit. The lower triangular factor L is stored in under the diagonal of (the diagonal is always 1.0f + /// for the L factor). The upper triangular factor U is stored on and above the diagonal of . + /// The order of the square matrix . + /// On exit, it contains the pivot indices. The size of the array must be . + /// This is equivalent to the GETRF LAPACK routine. + [SecuritySafeCritical] + public override void LUFactor(float[] data, int order, int[] ipiv) + { + if (data == null) + { + throw new ArgumentNullException("data"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (data.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "data"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + Solver(SafeNativeMethods.s_lu_factor(_solverHandle, order, data, ipiv)); + } + + /// + /// Computes the inverse of matrix using LU factorization. + /// + /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . + /// This is equivalent to the GETRF and GETRI LAPACK routines. + [SecuritySafeCritical] + public override void LUInverse(float[] a, int order) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + Solver(SafeNativeMethods.s_lu_inverse(_solverHandle, _blasHandle, order, a)); + } + + /// + /// Computes the inverse of a previously factored matrix. + /// + /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . + /// The pivot indices of . + /// This is equivalent to the GETRI LAPACK routine. + [SecuritySafeCritical] + public override void LUInverseFactored(float[] a, int order, int[] ipiv) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + BLAS(SafeNativeMethods.s_lu_inverse_factored(_blasHandle, order, a, ipiv)); + } + + /// + /// Computes the inverse of matrix using LU factorization. + /// + /// The N by N matrix to invert. Contains the inverse On exit. + /// The order of the square matrix . + /// Not supported. Should be left null. + /// This is equivalent to the GETRF and GETRI LAPACK routines. + [SecuritySafeCritical] + public override void LUInverse(float[] a, int order, float[] work) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (work != null) + { + throw new ArgumentException(Resources.UserWorkBufferNotSupported); + } + + Solver(SafeNativeMethods.s_lu_inverse(_solverHandle, _blasHandle, order, a)); + } + + /// + /// Computes the inverse of a previously factored matrix. + /// + /// The LU factored N by N matrix. Contains the inverse On exit. + /// The order of the square matrix . + /// The pivot indices of . + /// Not supported. This should be left null. + /// This is equivalent to the GETRI LAPACK routine. + [SecuritySafeCritical] + public override void LUInverseFactored(float[] a, int order, int[] ipiv, float[] work) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + if (work != null) + { + throw new ArgumentException(Resources.UserWorkBufferNotSupported); + } + + BLAS(SafeNativeMethods.s_lu_inverse_factored(_blasHandle, order, a, ipiv)); + } + + /// + /// Solves A*X=B for X using LU factorization. + /// + /// The number of columns of B. + /// The square matrix A. + /// The order of the square matrix . + /// On entry the B matrix; on exit the X matrix. + /// This is equivalent to the GETRF and GETRS LAPACK routines. + [SecuritySafeCritical] + public override void LUSolve(int columnsOfB, float[] a, int order, float[] b) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (b.Length != columnsOfB*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + Solver(SafeNativeMethods.s_lu_solve(_solverHandle, order, columnsOfB, a, b)); + } + + /// + /// Solves A*X=B for X using a previously factored A matrix. + /// + /// The number of columns of B. + /// The factored A matrix. + /// The order of the square matrix . + /// The pivot indices of . + /// On entry the B matrix; on exit the X matrix. + /// This is equivalent to the GETRS LAPACK routine. + [SecuritySafeCritical] + public override void LUSolveFactored(int columnsOfB, float[] a, int order, int[] ipiv, float[] b) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (ipiv == null) + { + throw new ArgumentNullException("ipiv"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + if (ipiv.Length != order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "ipiv"); + } + + if (b.Length != columnsOfB*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + Solver(SafeNativeMethods.s_lu_solve_factored(_solverHandle, order, columnsOfB, a, ipiv, b)); + } + + /// + /// Computes the Cholesky factorization of A. + /// + /// On entry, a square, positive definite matrix. On exit, the matrix is overwritten with the + /// the Cholesky factorization. + /// The number of rows or columns in the matrix. + /// This is equivalent to the POTRF LAPACK routine. + [SecuritySafeCritical] + public override void CholeskyFactor(float[] a, int order) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (order < 1) + { + throw new ArgumentException(Resources.ArgumentMustBePositive, "order"); + } + + if (a.Length != order*order) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "a"); + } + + Solver(SafeNativeMethods.s_cholesky_factor(_solverHandle, order, a)); + } + + /// + /// Solves A*X=B for X using Cholesky factorization. + /// + /// The square, positive definite matrix A. + /// The number of rows and columns in A. + /// On entry the B matrix; on exit the X matrix. + /// The number of columns in the B matrix. + /// This is equivalent to the POTRF add POTRS LAPACK routines. + /// + [SecuritySafeCritical] + public override void CholeskySolve(float[] a, int orderA, float[] b, int columnsB) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (b.Length != orderA*columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + Solver(SafeNativeMethods.s_cholesky_solve(_solverHandle, orderA, columnsB, a, b)); + } + + /// + /// Solves A*X=B for X using a previously factored A matrix. + /// + /// The square, positive definite matrix A. + /// The number of rows and columns in A. + /// On entry the B matrix; on exit the X matrix. + /// The number of columns in the B matrix. + /// This is equivalent to the POTRS LAPACK routine. + [SecuritySafeCritical] + public override void CholeskySolveFactored(float[] a, int orderA, float[] b, int columnsB) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (b.Length != orderA*columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (ReferenceEquals(a, b)) + { + throw new ArgumentException(Resources.ArgumentReferenceDifferent); + } + + Solver(SafeNativeMethods.s_cholesky_solve_factored(_solverHandle, orderA, columnsB, a, b)); + } + + /// + /// Computes the singular value decomposition of A. + /// + /// Compute the singular U and VT vectors or not. + /// On entry, the M by N matrix to decompose. On exit, A may be overwritten. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The singular values of A in ascending value. + /// If is true, on exit U contains the left + /// singular vectors. + /// If is true, on exit VT contains the transposed + /// right singular vectors. + /// This is equivalent to the GESVD LAPACK routine. + [SecuritySafeCritical] + public override void SingularValueDecomposition(bool computeVectors, float[] a, int rowsA, int columnsA, float[] s, float[] u, float[] vt) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (s == null) + { + throw new ArgumentNullException("s"); + } + + if (u == null) + { + throw new ArgumentNullException("u"); + } + + if (vt == null) + { + throw new ArgumentNullException("vt"); + } + + if (u.Length != rowsA*rowsA) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); + } + + if (vt.Length != columnsA*columnsA) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); + } + + if (s.Length != Math.Min(rowsA, columnsA)) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); + } + + SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, null); + } + + /// + /// Solves A*X=B for X using the singular value decomposition of A. + /// + /// On entry, the M by N matrix to decompose. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The B matrix. + /// The number of columns of B. + /// On exit, the solution matrix. + public override void SvdSolve(float[] a, int rowsA, int columnsA, float[] b, int columnsB, float[] x) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (b == null) + { + throw new ArgumentNullException("b"); + } + + if (x == null) + { + throw new ArgumentNullException("x"); + } + + if (b.Length != rowsA*columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + if (x.Length != columnsA*columnsB) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "b"); + } + + var s = new float[Math.Min(rowsA, columnsA)]; + var u = new float[rowsA*rowsA]; + var vt = new float[columnsA*columnsA]; + + var clone = new float[a.Length]; + a.Copy(clone); + SingularValueDecomposition(true, clone, rowsA, columnsA, s, u, vt, null); + SvdSolveFactored(rowsA, columnsA, s, u, vt, b, columnsB, x); + } + + /// + /// Computes the singular value decomposition of A. + /// + /// Compute the singular U and VT vectors or not. + /// On entry, the M by N matrix to decompose. On exit, A may be overwritten. + /// The number of rows in the A matrix. + /// The number of columns in the A matrix. + /// The singular values of A in ascending value. + /// If is true, on exit U contains the left + /// singular vectors. + /// If is true, on exit VT contains the transposed + /// right singular vectors. + /// Not supported. Should be left null. + /// This is equivalent to the GESVD LAPACK routine. + [SecuritySafeCritical] + public override void SingularValueDecomposition(bool computeVectors, float[] a, int rowsA, int columnsA, float[] s, float[] u, float[] vt, float[] work) + { + if (a == null) + { + throw new ArgumentNullException("a"); + } + + if (s == null) + { + throw new ArgumentNullException("s"); + } + + if (u == null) + { + throw new ArgumentNullException("u"); + } + + if (vt == null) + { + throw new ArgumentNullException("vt"); + } + + if (work != null) + { + throw new ArgumentException(Resources.UserWorkBufferNotSupported); + } + + if (u.Length != rowsA*rowsA) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "u"); + } + + if (vt.Length != columnsA*columnsA) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "vt"); + } + + if (s.Length != Math.Min(rowsA, columnsA)) + { + throw new ArgumentException(Resources.ArgumentArraysSameLength, "s"); + } + + if (columnsA > rowsA || !computeVectors) // see remarks http://docs.nvidia.com/cuda/cusolver/index.html#cuds-lt-t-gt-gesvd + base.SingularValueDecomposition(computeVectors, a, rowsA, columnsA, s, u, vt, new float[rowsA]); + else Solver(SafeNativeMethods.s_svd_factor(_solverHandle, computeVectors, rowsA, columnsA, a, s, u, vt)); + } + } +} + +#endif diff --git a/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.cs b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.cs new file mode 100644 index 00000000..fa11f02e --- /dev/null +++ b/src/Numerics/Providers/LinearAlgebra/Cuda/CudaLinearAlgebraProvider.cs @@ -0,0 +1,202 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://numerics.mathdotnet.com +// http://github.com/mathnet/mathnet-numerics +// http://mathnetnumerics.codeplex.com +// +// Copyright (c) 2009-2015 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +using System; + +#if NATIVE + +namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda +{ + /// + /// NVidia's CUDA Toolkit linear algebra provider. + /// + public partial class CudaLinearAlgebraProvider : ManagedLinearAlgebraProvider, IDisposable + { + private int _nativeRevision; + private bool _nativeIX86; + private bool _nativeX64; + private bool _nativeIA64; + private IntPtr _blasHandle; + private IntPtr _solverHandle; + + + /// + /// Constructor. + /// + [CLSCompliant(false)] + public CudaLinearAlgebraProvider() + { + } + + /// + /// Initialize and verify that the provided is indeed available. + /// If calling this method fails, consider to fall back to alternatives like the managed provider. + /// + public override void InitializeVerify() + { + int a, b, linearAlgebra; + try + { + // Load the native library + NativeProviderLoader.TryLoad(SafeNativeMethods.DllName); + + a = SafeNativeMethods.query_capability(0); + b = SafeNativeMethods.query_capability(1); + + _nativeIX86 = SafeNativeMethods.query_capability(8) > 0; + _nativeX64 = SafeNativeMethods.query_capability(9) > 0; + _nativeIA64 = SafeNativeMethods.query_capability(10) > 0; + + _nativeRevision = SafeNativeMethods.query_capability(64); + linearAlgebra = SafeNativeMethods.query_capability(128); + } + catch (DllNotFoundException e) + { + throw new NotSupportedException("Cuda Native Provider not found.", e); + } + catch (BadImageFormatException e) + { + throw new NotSupportedException("Cuda Native Provider found but failed to load. Please verify that the platform matches (x64 vs x32, Windows vs Linux).", e); + } + catch (EntryPointNotFoundException e) + { + throw new NotSupportedException("Cuda Native Provider does not support capability querying and is therefore not compatible. Consider upgrading to a newer version.", e); + } + + if (a != 0 || b != -1 || linearAlgebra <=0 || _nativeRevision < 1) + { + throw new NotSupportedException("Cuda Native Provider too old or not compatible. Consider upgrading to a newer version."); + } + + BLAS(SafeNativeMethods.createBLASHandle(ref _blasHandle)); + Solver(SafeNativeMethods.createSolverHandle(ref _solverHandle)); + } + + private void BLAS(int status) + { + switch (status) + { + case 0: // CUBLAS_STATUS_SUCCESS + return; + + case 1: // CUBLAS_STATUS_NOT_INITIALIZED + throw new Exception("The CUDA Runtime initialization failed"); + + case 2: // CUSOLVER_STATUS_ALLOC_FAILED + throw new OutOfMemoryException("The resources could not be allocated"); + + case 7: // CUBLAS_STATUS_INVALID_VALUE + throw new ArgumentException("Invalid value"); + + case 8: // CUBLAS_STATUS_ARCH_MISMATCH + throw new NotSupportedException("The device does not support this opeation."); + + case 11: // CUBLAS_STATUS_MAPPING_ERROR + throw new Exception("Mapping error."); + + case 13: // CUBLAS_STATUS_EXECUTION_FAILED + throw new Exception("Execution failed"); + + case 14: // CUBLAS_STATUS_INTERNAL_ERROR + throw new Exception("Internal error"); + + case 15: // CUBLAS_STATUS_NOT_SUPPORTED + throw new NotSupportedException(); + + case 16: // CUBLAS_STATUS_LICENSE_ERROR + throw new Exception("License error"); + + default: + throw new Exception("Unrecognized cuBLAS status code: " + status); + } + } + + private void Solver(int status) + { + switch (status) + { + case 0: // CUSOLVER_STATUS_SUCCESS + return; + + case 1: // CUSOLVER_STATUS_NOT_INITIALIZED + throw new Exception("The library was not initialized"); + + case 2: // CUSOLVER_STATUS_ALLOC_FAILED + throw new OutOfMemoryException("The resources could not be allocated"); + + case 3: // CUSOLVER_STATUS_INVALID_VALUE + throw new ArgumentException("Invalid value"); + + case 4: // CUSOLVER_STATUS_ARCH_MISMATCH + throw new NotSupportedException("The device does not support compute capability 2.0 and above"); + + case 5: // CUSOLVER_STATUS_MAPPING_ERROR + throw new Exception("Mapping error"); + + case 6: // CUSOLVER_STATUS_EXECUTION_FAILED + throw new NonConvergenceException("Execution failed"); + + case 7: //CUSOLVER_STATUS_INTERNAL_ERROR + throw new Exception("Internal error"); + + case 8: // CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED + throw new ArgumentException("Matrix type not supported"); + + case 9: // CUSOLVER_STATUS_NOT_SUPPORTED + throw new NotSupportedException(); + + case 10: // CUSOLVER_STATUS_ZERO_PIVOT + throw new Exception("Zero pivot"); + + case 11: //CUSOLVER_STATUS_INVALID_LICENSE + throw new Exception("Invalid license"); + + default: + throw new Exception("Unrecognized cuSolverDn status code: " + status); + + + } + } + + public override string ToString() + { + return string.Format("Nvidia CUDA ({1}; revision {0})", _nativeRevision, _nativeIX86 ? "x86" : _nativeX64 ? "x64" : _nativeIA64 ? "IA64" : "unknown"); + } + + + public void Dispose() + { + BLAS(SafeNativeMethods.destroyBLASHandle(_blasHandle)); + Solver(SafeNativeMethods.destroySolverHandle(_solverHandle)); + } + } +} + +#endif diff --git a/src/Numerics/Providers/LinearAlgebra/Cuda/SafeNativeMethods.cs b/src/Numerics/Providers/LinearAlgebra/Cuda/SafeNativeMethods.cs new file mode 100644 index 00000000..ba58a483 --- /dev/null +++ b/src/Numerics/Providers/LinearAlgebra/Cuda/SafeNativeMethods.cs @@ -0,0 +1,378 @@ +// +// Math.NET Numerics, part of the Math.NET Project +// http://mathnet.opensourcedotnet.info +// +// Copyright (c) 2009-2014 Math.NET +// +// Permission is hereby granted, free of charge, to any person +// obtaining a copy of this software and associated documentation +// files (the "Software"), to deal in the Software without +// restriction, including without limitation the rights to use, +// copy, modify, merge, publish, distribute, sublicense, and/or sell +// copies of the Software, and to permit persons to whom the +// Software is furnished to do so, subject to the following +// conditions: +// +// The above copyright notice and this permission notice shall be +// included in all copies or substantial portions of the Software. +// +// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +// OTHER DEALINGS IN THE SOFTWARE. +// + +#if NATIVE + +using System; +using System.Numerics; +using System.Runtime.InteropServices; +using System.Security; + +namespace MathNet.Numerics.Providers.LinearAlgebra.Cuda +{ + /// + /// P/Invoke methods to the native math libraries. + /// + [SuppressUnmanagedCodeSecurity] + [SecurityCritical] + internal static class SafeNativeMethods + { + // ReSharper disable InconsistentNaming + + /// + /// Name of the native DLL. + /// + const string _DllName = "MathNet.Numerics.CUDA.dll"; + internal static string DllName { get { return _DllName; } } + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int query_capability(int capability); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int createBLASHandle(ref IntPtr blasHandle); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int destroyBLASHandle(IntPtr blasHandle); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int createSolverHandle(ref IntPtr solverHandle); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int destroySolverHandle(IntPtr solverHandle); + + #region BLAS + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void s_axpy(IntPtr blasHandle, int n, float alpha, float[] x, [In, Out] float[] y); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void d_axpy(IntPtr blasHandle, int n, double alpha, double[] x, [In, Out] double[] y); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void c_axpy(IntPtr blasHandle, int n, Complex32 alpha, Complex32[] x, [In, Out] Complex32[] y); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void z_axpy(IntPtr blasHandle, int n, Complex alpha, Complex[] x, [In, Out] Complex[] y); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void s_scale(IntPtr blasHandle, int n, float alpha, [Out] float[] x); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void d_scale(IntPtr blasHandle, int n, double alpha, [Out] double[] x); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void c_scale(IntPtr blasHandle, int n, Complex32 alpha, [In, Out] Complex32[] x); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void z_scale(IntPtr blasHandle, int n, Complex alpha, [In, Out] Complex[] x); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern float s_dot_product(IntPtr blasHandle, int n, float[] x, float[] y); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern double d_dot_product(IntPtr blasHandle, int n, double[] x, double[] y); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern Complex32 c_dot_product(IntPtr blasHandle, int n, Complex32[] x, Complex32[] y); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern Complex z_dot_product(IntPtr blasHandle, int n, Complex[] x, Complex[] y); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void s_matrix_multiply(IntPtr blasHandle, int transA, int transB, int m, int n, int k, float alpha, float[] x, float[] y, float beta, [In, Out] float[] c); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void d_matrix_multiply(IntPtr blasHandle, int transA, int transB, int m, int n, int k, double alpha, double[] x, double[] y, double beta, [In, Out] double[] c); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void c_matrix_multiply(IntPtr blasHandle, int transA, int transB, int m, int n, int k, Complex32 alpha, Complex32[] x, Complex32[] y, Complex32 beta, [In, Out] Complex32[] c); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern void z_matrix_multiply(IntPtr blasHandle, int transA, int transB, int m, int n, int k, Complex alpha, Complex[] x, Complex[] y, Complex beta, [In, Out] Complex[] c); + + internal static int ToCUDA(this Transpose transpose) + { + switch (transpose) + { + case Transpose.DontTranspose: + return 0; + + case Transpose.Transpose: + return 1; + + case Transpose.ConjugateTranspose: + return 2; + + default: + throw new ArgumentException("Unsupported transpose: " + transpose); + } + } + + #endregion BLAS + + #region LAPACK + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern float s_matrix_norm(byte norm, int rows, int columns, [In] float[] a, [In, Out] float[] work); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern double d_matrix_norm(byte norm, int rows, int columns, [In] double[] a, [In, Out] double[] work); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern float c_matrix_norm(byte norm, int rows, int columns, [In] Complex32[] a, [In, Out] float[] work); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern double z_matrix_norm(byte norm, int rows, int columns, [In] Complex[] a, [In, Out] double[] work); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_cholesky_factor(IntPtr solverHandle, int n, [In, Out] float[] a); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_cholesky_factor(IntPtr solverHandle, int n, [In, Out] double[] a); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_cholesky_factor(IntPtr solverHandle, int n, [In, Out] Complex32[] a); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_cholesky_factor(IntPtr solverHandle, int n, [In, Out] Complex[] a); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_lu_factor(IntPtr solverHandle, int n, [In, Out] float[] a, [In, Out] int[] ipiv); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_lu_factor(IntPtr solverHandle, int n, [In, Out] double[] a, [In, Out] int[] ipiv); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_lu_factor(IntPtr solverHandle, int n, [In, Out] Complex32[] a, [In, Out] int[] ipiv); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_lu_factor(IntPtr solverHandle, int n, [In, Out] Complex[] a, [In, Out] int[] ipiv); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_lu_inverse(IntPtr solverHandle, IntPtr blasHandle, int n, [In, Out] float[] a); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_lu_inverse(IntPtr solverHandle, IntPtr blasHandle, int n, [In, Out] double[] a); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_lu_inverse(IntPtr solverHandle, IntPtr blasHandle, int n, [In, Out] Complex32[] a); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_lu_inverse(IntPtr solverHandle, IntPtr blasHandle, int n, [In, Out] Complex[] a); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_lu_inverse_factored(IntPtr blasHandle, int n, [In, Out] float[] a, [In, Out] int[] ipiv); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_lu_inverse_factored(IntPtr blasHandle, int n, [In, Out] double[] a, [In, Out] int[] ipiv); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_lu_inverse_factored(IntPtr blasHandle, int n, [In, Out] Complex32[] a, [In, Out] int[] ipiv); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_lu_inverse_factored(IntPtr blasHandle, int n, [In, Out] Complex[] a, [In, Out] int[] ipiv); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_lu_solve_factored(IntPtr solverHandle, int n, int nrhs, float[] a, [In, Out] int[] ipiv, [In, Out] float[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_lu_solve_factored(IntPtr solverHandle, int n, int nrhs, double[] a, [In, Out] int[] ipiv, [In, Out] double[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_lu_solve_factored(IntPtr solverHandle, int n, int nrhs, Complex32[] a, [In, Out] int[] ipiv, [In, Out] Complex32[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_lu_solve_factored(IntPtr solverHandle, int n, int nrhs, Complex[] a, [In, Out] int[] ipiv, [In, Out] Complex[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_lu_solve(IntPtr solverHandle, int n, int nrhs, float[] a, [In, Out] float[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_lu_solve(IntPtr solverHandle, int n, int nrhs, double[] a, [In, Out] double[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_lu_solve(IntPtr solverHandle, int n, int nrhs, Complex32[] a, [In, Out] Complex32[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_lu_solve(IntPtr solverHandle, int n, int nrhs, Complex[] a, [In, Out] Complex[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_cholesky_solve(IntPtr solverHandle, int n, int nrhs, float[] a, [In, Out] float[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_cholesky_solve(IntPtr solverHandle, int n, int nrhs, double[] a, [In, Out] double[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_cholesky_solve(IntPtr solverHandle, int n, int nrhs, Complex32[] a, [In, Out] Complex32[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_cholesky_solve(IntPtr solverHandle, int n, int nrhs, Complex[] a, [In, Out] Complex[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_cholesky_solve_factored(IntPtr solverHandle, int n, int nrhs, float[] a, [In, Out] float[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_cholesky_solve_factored(IntPtr solverHandle, int n, int nrhs, double[] a, [In, Out] double[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_cholesky_solve_factored(IntPtr solverHandle, int n, int nrhs, Complex32[] a, [In, Out] Complex32[] b); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_cholesky_solve_factored(IntPtr solverHandle, int n, int nrhs, Complex[] a, [In, Out] Complex[] b); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern int s_qr_factor(int m, int n, [In, Out] float[] r, [In, Out] float[] tau, [In, Out] float[] q, [In, Out] float[] work, int len); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern int d_qr_factor(int m, int n, [In, Out] double[] r, [In, Out] double[] tau, [In, Out] double[] q, [In, Out] double[] work, int len); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern int c_qr_factor(int m, int n, [In, Out] Complex32[] r, [In, Out] Complex32[] tau, [In, Out] Complex32[] q, [In, Out] Complex32[] work, int len); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern int z_qr_factor(int m, int n, [In, Out] Complex[] r, [In, Out] Complex[] tau, [In, Out] Complex[] q, [In, Out] Complex[] work, int len); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern int s_qr_thin_factor(int m, int n, [In, Out] float[] q, [In, Out] float[] tau, [In, Out] float[] r, [In, Out] float[] work, int len); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern int d_qr_thin_factor(int m, int n, [In, Out] double[] q, [In, Out] double[] tau, [In, Out] double[] r, [In, Out] double[] work, int len); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern int c_qr_thin_factor(int m, int n, [In, Out] Complex32[] q, [In, Out] Complex32[] tau, [In, Out] Complex32[] r, [In, Out] Complex32[] work, int len); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern int z_qr_thin_factor(int m, int n, [In, Out] Complex[] q, [In, Out] Complex[] tau, [In, Out] Complex[] r, [In, Out] Complex[] work, int len); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern int s_qr_solve(int m, int n, int bn, float[] r, float[] b, [In, Out] float[] x, [In, Out] float[] work, int len); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern int d_qr_solve(int m, int n, int bn, double[] r, double[] b, [In, Out] double[] x, [In, Out] double[] work, int len); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern int c_qr_solve(int m, int n, int bn, Complex32[] r, Complex32[] b, [In, Out] Complex32[] x, [In, Out] Complex32[] work, int len); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern int z_qr_solve(int m, int n, int bn, Complex[] r, Complex[] b, [In, Out] Complex[] x, [In, Out] Complex[] work, int len); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern int s_qr_solve_factored(int m, int n, int bn, float[] r, float[] b, float[] tau, [In, Out] float[] x, [In, Out] float[] work, int len); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern int d_qr_solve_factored(int m, int n, int bn, double[] r, double[] b, double[] tau, [In, Out] double[] x, [In, Out] double[] work, int len); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern int c_qr_solve_factored(int m, int n, int bn, Complex32[] r, Complex32[] b, Complex32[] tau, [In, Out] Complex32[] x, [In, Out] Complex32[] work, int len); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern int z_qr_solve_factored(int m, int n, int bn, Complex[] r, Complex[] b, Complex[] tau, [In, Out] Complex[] x, [In, Out] Complex[] work, int len); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int s_svd_factor(IntPtr solverHandle, [MarshalAs(UnmanagedType.U1)] bool computeVectors, int m, int n, [In, Out] float[] a, [In, Out] float[] s, [In, Out] float[] u, [In, Out] float[] v); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int d_svd_factor(IntPtr solverHandle, [MarshalAs(UnmanagedType.U1)] bool computeVectors, int m, int n, [In, Out] double[] a, [In, Out] double[] s, [In, Out] double[] u, [In, Out] double[] v); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int c_svd_factor(IntPtr solverHandle, [MarshalAs(UnmanagedType.U1)] bool computeVectors, int m, int n, [In, Out] Complex32[] a, [In, Out] Complex32[] s, [In, Out] Complex32[] u, [In, Out] Complex32[] v); + + [DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + internal static extern int z_svd_factor(IntPtr solverHandle, [MarshalAs(UnmanagedType.U1)] bool computeVectors, int m, int n, [In, Out] Complex[] a, [In, Out] Complex[] s, [In, Out] Complex[] u, [In, Out] Complex[] v); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern int s_eigen([MarshalAs(UnmanagedType.U1)] bool isSymmetric, int n, [In] float[] a, [In, Out] float[] vectors, [In, Out] Complex[] values, [In, Out] float[] d); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern int d_eigen([MarshalAs(UnmanagedType.U1)] bool isSymmetric, int n, [In] double[] a, [In, Out] double[] vectors, [In, Out] Complex[] values, [In, Out] double[] d); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern int c_eigen([MarshalAs(UnmanagedType.U1)] bool isSymmetric, int n, [In] Complex32[] a, [In, Out] Complex32[] vectors, [In, Out] Complex[] values, [In, Out] Complex32[] d); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern int z_eigen([MarshalAs(UnmanagedType.U1)] bool isSymmetric, int n, [In] Complex[] a, [In, Out] Complex[] vectors, [In, Out] Complex[] values, [In, Out] Complex[] d); + + #endregion LAPACK + + #region Vector Functions + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern void s_vector_add(int n, float[] x, float[] y, [In, Out] float[] result); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern void s_vector_subtract(int n, float[] x, float[] y, [In, Out] float[] result); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern void s_vector_multiply(int n, float[] x, float[] y, [In, Out] float[] result); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern void s_vector_divide(int n, float[] x, float[] y, [In, Out] float[] result); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern void d_vector_add(int n, double[] x, double[] y, [In, Out] double[] result); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern void d_vector_subtract(int n, double[] x, double[] y, [In, Out] double[] result); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern void d_vector_multiply(int n, double[] x, double[] y, [In, Out] double[] result); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern void d_vector_divide(int n, double[] x, double[] y, [In, Out] double[] result); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern void c_vector_add(int n, Complex32[] x, Complex32[] y, [In, Out] Complex32[] result); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern void c_vector_subtract(int n, Complex32[] x, Complex32[] y, [In, Out] Complex32[] result); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern void c_vector_multiply(int n, Complex32[] x, Complex32[] y, [In, Out] Complex32[] result); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern void c_vector_divide(int n, Complex32[] x, Complex32[] y, [In, Out] Complex32[] result); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern void z_vector_add(int n, Complex[] x, Complex[] y, [In, Out] Complex[] result); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern void z_vector_subtract(int n, Complex[] x, Complex[] y, [In, Out] Complex[] result); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern void z_vector_multiply(int n, Complex[] x, Complex[] y, [In, Out] Complex[] result); + + //[DllImport(_DllName, ExactSpelling = true, SetLastError = false, CallingConvention = CallingConvention.Cdecl)] + //internal static extern void z_vector_divide(int n, Complex[] x, Complex[] y, [In, Out] Complex[] result); + + #endregion Vector Functions + + // ReSharper restore InconsistentNaming + } +} + +#endif diff --git a/src/UnitTests/LinearAlgebraProviderTests/Complex/LinearAlgebraProviderTests.cs b/src/UnitTests/LinearAlgebraProviderTests/Complex/LinearAlgebraProviderTests.cs index 7508a520..d17a75cd 100644 --- a/src/UnitTests/LinearAlgebraProviderTests/Complex/LinearAlgebraProviderTests.cs +++ b/src/UnitTests/LinearAlgebraProviderTests/Complex/LinearAlgebraProviderTests.cs @@ -444,7 +444,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex var a = new Complex[matrix.RowCount*matrix.RowCount]; Array.Copy(matrix.Values, a, a.Length); +#if CUDA + Complex[] work = null; +#else var work = new Complex[matrix.RowCount]; +#endif Control.LinearAlgebraProvider.LUInverse(a, matrix.RowCount, work); AssertHelpers.AlmostEqualRelative(a[0], -0.454545454545454, 13); @@ -473,7 +477,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex Control.LinearAlgebraProvider.LUFactor(a, matrix.RowCount, ipiv); +#if CUDA + Complex[] work = null; +#else var work = new Complex[matrix.RowCount]; +#endif Control.LinearAlgebraProvider.LUInverseFactored(a, matrix.RowCount, ipiv, work); AssertHelpers.AlmostEqualRelative(a[0], -0.454545454545454, 13); @@ -1447,7 +1455,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex var s = new Complex[matrix.RowCount]; var u = new Complex[matrix.RowCount*matrix.RowCount]; var vt = new Complex[matrix.ColumnCount*matrix.ColumnCount]; +#if CUDA + Complex[] work = null; +#else var work = new Complex[100]; +#endif Control.LinearAlgebraProvider.SingularValueDecomposition(true, a, matrix.RowCount, matrix.ColumnCount, s, u, vt, work); @@ -1486,7 +1498,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex var s = new Complex[matrix.ColumnCount]; var u = new Complex[matrix.RowCount*matrix.RowCount]; var vt = new Complex[matrix.ColumnCount*matrix.ColumnCount]; +#if CUDA + Complex[] work = null; +#else var work = new Complex[100]; +#endif Control.LinearAlgebraProvider.SingularValueDecomposition(true, a, matrix.RowCount, matrix.ColumnCount, s, u, vt, work); @@ -1522,7 +1538,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex var s = new Complex[matrix.RowCount]; var u = new Complex[matrix.RowCount*matrix.RowCount]; var vt = new Complex[matrix.ColumnCount*matrix.ColumnCount]; +#if CUDA + Complex[] work = null; +#else var work = new Complex[100]; +#endif Control.LinearAlgebraProvider.SingularValueDecomposition(true, a, matrix.RowCount, matrix.ColumnCount, s, u, vt, work); diff --git a/src/UnitTests/LinearAlgebraProviderTests/Complex32/LinearAlgebraProviderTests.cs b/src/UnitTests/LinearAlgebraProviderTests/Complex32/LinearAlgebraProviderTests.cs index d13ea874..d0b8a4e3 100644 --- a/src/UnitTests/LinearAlgebraProviderTests/Complex32/LinearAlgebraProviderTests.cs +++ b/src/UnitTests/LinearAlgebraProviderTests/Complex32/LinearAlgebraProviderTests.cs @@ -448,7 +448,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex32 var a = new Complex32[matrix.RowCount*matrix.RowCount]; Array.Copy(matrix.Values, a, a.Length); +#if CUDA + Complex32[] work = null; +#else var work = new Complex32[matrix.RowCount]; +#endif Control.LinearAlgebraProvider.LUInverse(a, matrix.RowCount, work); AssertHelpers.AlmostEqualRelative(a[0], -0.454545454545454f, 5); @@ -477,7 +481,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex32 Control.LinearAlgebraProvider.LUFactor(a, matrix.RowCount, ipiv); +#if CUDA + Complex32[] work = null; +#else var work = new Complex32[matrix.RowCount]; +#endif Control.LinearAlgebraProvider.LUInverseFactored(a, matrix.RowCount, ipiv, work); AssertHelpers.AlmostEqualRelative(a[0], -0.454545454545454f, 5); @@ -1452,7 +1460,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex32 var s = new Complex32[matrix.RowCount]; var u = new Complex32[matrix.RowCount*matrix.RowCount]; var vt = new Complex32[matrix.ColumnCount*matrix.ColumnCount]; +#if CUDA + Complex32[] work = null; +#else var work = new Complex32[100]; +#endif Control.LinearAlgebraProvider.SingularValueDecomposition(true, a, matrix.RowCount, matrix.ColumnCount, s, u, vt, work); @@ -1491,7 +1503,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex32 var s = new Complex32[matrix.ColumnCount]; var u = new Complex32[matrix.RowCount*matrix.RowCount]; var vt = new Complex32[matrix.ColumnCount*matrix.ColumnCount]; +#if CUDA + Complex32[] work = null; +#else var work = new Complex32[100]; +#endif Control.LinearAlgebraProvider.SingularValueDecomposition(true, a, matrix.RowCount, matrix.ColumnCount, s, u, vt, work); @@ -1527,7 +1543,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Complex32 var s = new Complex32[matrix.RowCount]; var u = new Complex32[matrix.RowCount*matrix.RowCount]; var vt = new Complex32[matrix.ColumnCount*matrix.ColumnCount]; +#if CUDA + Complex32[] work = null; +#else var work = new Complex32[100]; +#endif Control.LinearAlgebraProvider.SingularValueDecomposition(true, a, matrix.RowCount, matrix.ColumnCount, s, u, vt, work); diff --git a/src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs b/src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs index 2e50b594..6546bc55 100644 --- a/src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs +++ b/src/UnitTests/LinearAlgebraProviderTests/Double/LinearAlgebraProviderTests.cs @@ -438,7 +438,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double var a = new double[matrix.RowCount*matrix.RowCount]; Array.Copy(matrix.Values, a, a.Length); +#if CUDA + double[] work = null; +#else var work = new double[matrix.RowCount]; +#endif Control.LinearAlgebraProvider.LUInverse(a, matrix.RowCount, work); AssertHelpers.AlmostEqualRelative(a[0], -0.454545454545454, 13); @@ -467,7 +471,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double Control.LinearAlgebraProvider.LUFactor(a, matrix.RowCount, ipiv); +#if CUDA + double[] work = null; +#else var work = new double[matrix.RowCount]; +#endif Control.LinearAlgebraProvider.LUInverseFactored(a, matrix.RowCount, ipiv, work); AssertHelpers.AlmostEqualRelative(a[0], -0.454545454545454, 13); @@ -1442,7 +1450,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double var s = new double[matrix.RowCount]; var u = new double[matrix.RowCount*matrix.RowCount]; var vt = new double[matrix.ColumnCount*matrix.ColumnCount]; +#if CUDA + double[] work = null; +#else var work = new double[100]; +#endif Control.LinearAlgebraProvider.SingularValueDecomposition(true, a, matrix.RowCount, matrix.ColumnCount, s, u, vt, work); @@ -1481,7 +1493,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double var s = new double[matrix.ColumnCount]; var u = new double[matrix.RowCount*matrix.RowCount]; var vt = new double[matrix.ColumnCount*matrix.ColumnCount]; +#if CUDA + double[] work = null; +#else var work = new double[100]; +#endif Control.LinearAlgebraProvider.SingularValueDecomposition(true, a, matrix.RowCount, matrix.ColumnCount, s, u, vt, work); @@ -1517,7 +1533,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Double var s = new double[matrix.RowCount]; var u = new double[matrix.RowCount*matrix.RowCount]; var vt = new double[matrix.ColumnCount*matrix.ColumnCount]; +#if CUDA + double[] work = null; +#else var work = new double[100]; +#endif Control.LinearAlgebraProvider.SingularValueDecomposition(true, a, matrix.RowCount, matrix.ColumnCount, s, u, vt, work); diff --git a/src/UnitTests/LinearAlgebraProviderTests/Single/LinearAlgebraProviderTests.cs b/src/UnitTests/LinearAlgebraProviderTests/Single/LinearAlgebraProviderTests.cs index 1ed7b3f0..fab970db 100644 --- a/src/UnitTests/LinearAlgebraProviderTests/Single/LinearAlgebraProviderTests.cs +++ b/src/UnitTests/LinearAlgebraProviderTests/Single/LinearAlgebraProviderTests.cs @@ -446,7 +446,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Single var a = new float[matrix.RowCount*matrix.RowCount]; Array.Copy(matrix.Values, a, a.Length); +#if CUDA + float[] work = null; +#else var work = new float[matrix.RowCount]; +#endif Control.LinearAlgebraProvider.LUInverse(a, matrix.RowCount, work); AssertHelpers.AlmostEqual(a[0], -0.454545454545454, 5); @@ -475,7 +479,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Single Control.LinearAlgebraProvider.LUFactor(a, matrix.RowCount, ipiv); +#if CUDA + float[] work = null; +#else var work = new float[matrix.RowCount]; +#endif Control.LinearAlgebraProvider.LUInverseFactored(a, matrix.RowCount, ipiv, work); AssertHelpers.AlmostEqual(a[0], -0.454545454545454, 5); @@ -1450,7 +1458,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Single var s = new float[matrix.RowCount]; var u = new float[matrix.RowCount*matrix.RowCount]; var vt = new float[matrix.ColumnCount*matrix.ColumnCount]; +#if CUDA + float[] work = null; +#else var work = new float[100]; +#endif Control.LinearAlgebraProvider.SingularValueDecomposition(true, a, matrix.RowCount, matrix.ColumnCount, s, u, vt, work); @@ -1489,7 +1501,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Single var s = new float[matrix.ColumnCount]; var u = new float[matrix.RowCount*matrix.RowCount]; var vt = new float[matrix.ColumnCount*matrix.ColumnCount]; +#if CUDA + float[] work = null; +#else var work = new float[100]; +#endif Control.LinearAlgebraProvider.SingularValueDecomposition(true, a, matrix.RowCount, matrix.ColumnCount, s, u, vt, work); @@ -1525,7 +1541,11 @@ namespace MathNet.Numerics.UnitTests.LinearAlgebraProviderTests.Single var s = new float[matrix.RowCount]; var u = new float[matrix.RowCount*matrix.RowCount]; var vt = new float[matrix.ColumnCount*matrix.ColumnCount]; +#if CUDA + float[] work = null; +#else var work = new float[100]; +#endif Control.LinearAlgebraProvider.SingularValueDecomposition(true, a, matrix.RowCount, matrix.ColumnCount, s, u, vt, work); diff --git a/src/UnitTests/UnitTests-CUDA.csproj b/src/UnitTests/UnitTests-CUDA.csproj new file mode 100644 index 00000000..074df81c --- /dev/null +++ b/src/UnitTests/UnitTests-CUDA.csproj @@ -0,0 +1,348 @@ + + + + 10.0 + Debug + AnyCPU + 8.0.30703 + 2.0 + {E79C0395-01DC-4BC9-B86C-ED45790892C5} + Library + Properties + MathNet.Numerics.UnitTests + MathNet.Numerics.UnitTestsCUDA + v4.5 + 512 + ..\..\ + + + TRACE;NATIVE;CUDA + ..\..\out\CUDA\Windows\ + ..\..\obj\CUDA\Windows\x86\ + ..\..\obj\CUDA\Windows\x86\ + true + pdbonly + prompt + MinimumRecommendedRules.ruleset + 1591 + AnyCPU + + + TRACE;DEBUG;NATIVE;CUDA + ..\..\out\CUDA\Windows\ + ..\..\obj\CUDA\Windows\x86\ + ..\..\obj\CUDA\Windows\x86\ + false + full + true + prompt + 4 + 1591 + AnyCPU + + + + + + + + + + + + + + + + + + + data\Codeplex-5667.csv + Always + + + data\Github-Cureos-1.csv + Always + + + data\Matlab\A.mat + Always + + + data\Matlab\collection-nocompress.mat + Always + + + data\Matlab\collection.mat + Always + + + data\Matlab\complex.mat + Always + + + data\Matlab\sparse-large.mat + Always + + + data\Matlab\sparse-small.mat + Always + + + data\Matlab\sparse_complex.mat + Always + + + data\Matlab\v.mat + Always + + + data\NIST\AtmWtAgt.dat + Always + + + data\NIST\Bennett5.dat + Always + + + data\NIST\BoxBOD.dat + Always + + + data\NIST\Chwirut1.dat + Always + + + data\NIST\Chwirut2.dat + Always + + + data\NIST\DanWood.dat + Always + + + data\NIST\Eckerle4.dat + Always + + + data\NIST\ENSO.dat + Always + + + data\NIST\Filip.dat + Always + + + data\NIST\Gauss1.dat + Always + + + data\NIST\Gauss2.dat + Always + + + data\NIST\Gauss3.dat + Always + + + data\NIST\Hahn1.dat + Always + + + data\NIST\Kirby2.dat + Always + + + data\NIST\Lanczos1.dat + Always + + + data\NIST\Lanczos2.dat + Always + + + data\NIST\Lanczos3.dat + Always + + + data\NIST\Lew.dat + Always + + + data\NIST\Longley.dat + Always + + + data\NIST\Lottery.dat + Always + + + data\NIST\Mavro.dat + Always + + + data\NIST\MGH09.dat + Always + + + data\NIST\MGH10.dat + Always + + + data\NIST\MGH17.dat + Always + + + data\NIST\Michelso.dat + Always + + + data\NIST\Misra1a.dat + Always + + + data\NIST\Misra1b.dat + Always + + + data\NIST\Misra1c.dat + Always + + + data\NIST\Misra1d.dat + Always + + + data\NIST\Nelson.dat + Always + + + data\NIST\NoInt1.dat + Always + + + data\NIST\NoInt2.dat + Always + + + data\NIST\Norris.dat + Always + + + data\NIST\NumAcc1.dat + Always + + + data\NIST\NumAcc2.dat + Always + + + data\NIST\NumAcc3.dat + Always + + + data\NIST\NumAcc4.dat + Always + + + data\NIST\Pontius.dat + Always + + + data\NIST\Rat42.dat + Always + + + data\NIST\Rat43.dat + Always + + + data\NIST\Roszman1.dat + Always + + + data\NIST\SiRstvt.dat + Always + + + data\NIST\SmLs01t.dat + Always + + + data\NIST\SmLs02t.dat + Always + + + data\NIST\SmLs03t.dat + Always + + + data\NIST\SmLs04t.dat + Always + + + data\NIST\SmLs05t.dat + Always + + + data\NIST\SmLs06t.dat + Always + + + data\NIST\SmLs07t.dat + Always + + + data\NIST\SmLs08t.dat + Always + + + data\NIST\SmLs09t.dat + Always + + + data\NIST\Thurber.dat + Always + + + data\NIST\Wampler1.dat + Always + + + data\NIST\Wampler2.dat + Always + + + data\NIST\Wampler3.dat + Always + + + data\NIST\Wampler4.dat + Always + + + data\NIST\Wampler5.dat + Always + + + + data\NIST\Meixner.dat + Always + + + + + + {b7cae5f4-a23f-4438-b5be-41226618b695} + Numerics + + + + + + ..\..\packages\NUnit\lib\nunit.framework.dll + True + True + + + \ No newline at end of file diff --git a/src/UnitTests/UnitTests-MKL.csproj b/src/UnitTests/UnitTests-MKL.csproj index bbd73bbf..07bd43ec 100644 --- a/src/UnitTests/UnitTests-MKL.csproj +++ b/src/UnitTests/UnitTests-MKL.csproj @@ -16,7 +16,7 @@ ..\..\ - TRACE;NATIVE + TRACE;NATIVE;MKL ..\..\out\MKL\Windows\ ..\..\obj\MKL\Windows\x86\ ..\..\obj\MKL\Windows\x86\ @@ -28,7 +28,7 @@ AnyCPU - TRACE;DEBUG;NATIVE + TRACE;DEBUG;NATIVE;MKL ..\..\out\MKL\Windows\ ..\..\obj\MKL\Windows\x86\ ..\..\obj\MKL\Windows\x86\ diff --git a/src/UnitTests/UseLinearAlgebraProvider.cs b/src/UnitTests/UseLinearAlgebraProvider.cs index fc4a42d6..00266687 100644 --- a/src/UnitTests/UseLinearAlgebraProvider.cs +++ b/src/UnitTests/UseLinearAlgebraProvider.cs @@ -39,7 +39,11 @@ namespace MathNet.Numerics.UnitTests public void BeforeTest(TestDetails testDetails) { #if !NET35 && NATIVE +#if MKL Control.UseNativeMKL(); +#elif CUDA + Control.UseNativeCUDA(); +#endif #endif }